# Run with jystilts (java -jar jystilts.jar)
# Jar file: http://www.starlink.ac.uk/stilts/jystilts.jar

# Crossmatches Gouliermis et al. 2006 (Vizier J/ApJS/166/549/table2)
# with full Gaia DR1 positions, to identify NGC 346 motion between
# J2014.6 and J2015.0, and supplement the Gouliermis table with
# Gaia G-band photometry.

from stilts import *

# VizieR dataset
goul_url = ('http://vizier.u-strasbg.fr/viz-bin/votable?' +
            '-source=J%2fApJS%2f166%2f549&-oc.form=dec&-out.meta=DhuL' +
            '&-out.add=_RAJ%2C_DEJ&-out.max=100000')
goul = tread(goul_url)
goul = (goul.cmd_colmeta('-name', 'ra_goul', '_raj2000')
            .cmd_colmeta('-name', 'dec_goul', '_dej2000'))

# Match with Gaia DR1
# Note the "in_" parameter is so named because of python syntax;
# it's equivalent to "in" in command-line stilts.
match = cdsskymatch(in_=goul, cdstable='gaia dr1', find='all', radius=1)
match = (match.cmd_colmeta('-name', 'ra_gaia', 'ra')
              .cmd_colmeta('-name', 'dec_gaia', 'dec'))

# Plot the results of the match.  We set the plot2sky arguments up
# as a dictionary so that we can do something similar but not quite
# the same later.
matchplot_kwargs = dict(
   layer_a='mark', in_a=goul, lon_a='ra_goul', lat_a='dec_goul',
                   size_a=2, shape_a='open_circle', color_a='red',
   layer_b='mark', in_b=match, lon_b='ra_gaia', lat_b='dec_gaia',
                   size_b=2, shape_b='filled_circle', color_b='cyan',
   layer_m='link2', in_m=match, lon1_m='ra_goul', lat1_m='dec_goul',
                                lon2_m='ra_gaia', lat2_m='dec_gaia',
                   color_m='blue',
   leglabel_a='Gouliermis', leglabel_b='Gaia', leglabel_m='match',
   legpos='.9,.9',
   seq='_a,_b,_m'
)
plot2sky(**matchplot_kwargs)

# Look at a histogram of angular distances.
plot2plane(layer1='histogram', in1=match, x1='angDist')

# Add deltaRa and deltaDec columns for the 2-d displacement of the
# matches.
match = (match
        .cmd_addcol('deltaRa', '-units', 'arcsec', '3600*(ra_gaia-ra_goul)')
        .cmd_addcol('deltaDec', '-units', 'arcsec', '3600*(dec_gaia-dec_goul)'))

# Plot this and see an overdensity.
plot2plane(layer1='mark', in1=match, x1='deltaRa', y1='deltaDec', grid=True)

# Identify the overdense region by eye to get a deltaRa/deltaDec
# interval for true matches.
true_match = match.cmd_select('deltaRa>0.5 && deltaRa<1.0 && ' +
                              'deltaDec>-0.35 && deltaDec<-0.25')

# Repeat the crossmatch result plot from above but this time including
# only the true matches.
matchplot_kwargs['in_m'] = true_match
plot2sky(**matchplot_kwargs)

# Plot a colour-magnitude diagram combining HST/ACS and Gaia photometry.
plot2plane(layer1='mark', in1=true_match,
           x1='Vmag-Imag', y1='phot_g_mean_mag', yflip=True)