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

# Interrogates TGAS for sources in the region of the Pleiades,
# (using either conesearch or TAP), and uses the results to
# estimate the distance of the Pleiades cluster,
# plotting the proper motions and parallax histogram.
#
# Works with JyStilts, but requires v3.0-9.
# (there's a bug in previous versions that prevents the subset plotting
# from working correctly).

from stilts import *

# Pleiades coordinates
ra=56.75
dec=24.12
radius=5

use_tap = False

# Here is how to do a TGAS cone search using the GACS TAP service
if use_tap:
   gacs_tap = 'http://gea.esac.esa.int/tap-server/tap'
   adql = """
      SELECT * from gaiadr1.tgas_source AS g
      WHERE 1=CONTAINS(POINT('ICRS', g.ra, g.dec),
                       CIRCLE('ICRS', %f, %f, %f))
   """ % (ra, dec, radius)
   plei = tapquery(tapurl=gacs_tap, adql=adql, sync=True)

# Here is an alternative way to do it using a Simple Cone Search service
else:
   tgas_cone = 'http://gaia.ari.uni-heidelberg.de/cone/tgas?'
   coneurl = ( tgas_cone + "RA=%f&DEC=%f&SR=%f" % ( ra, dec, radius ) )
   plei = tread(coneurl)

# Check how many rows we have
plei.mode_count()

# Define a subset tables within a certain region of proper motion space
comoving = plei.cmd_select('hypot(pmra-20,pmdec+45)<2')

# Refine it to a certain region of parallax space
cluster = comoving.cmd_select('parallax>6 && parallax<9')

# Check how many rows we have
print ("Comoving rows: %d; Cluster rows: %d"
       % (comoving.count_rows(), cluster.count_rows()))

# Plot the data in proper motion space.
# Note the plot is 'live' - you can interact using the mouse.
plot2plane(x='pmra', y='pmdec',
           layer1='mark', in1=plei, color1='red', leglabel1="all",
           layer2='mark', in2=comoving, color2='cyan', leglabel2="comoving",
           layer3='mark', in3=cluster, color3='blue', leglabel3="cluster",
           size=3,
           seq='1,2,3')

# Work out the mean and standard deviation of parallax for the cluster subset.
parallax_stats = cluster.cmd_keepcols("parallax").cmd_stats('Mean', 'StDev')[0]
p_mean = parallax_stats['Mean']
p_stdev = parallax_stats['StDev']
print "Pleaides parallax: %.2f +/- %.2f (+ systematics)" % (p_mean, p_stdev)

# Plot a histogram with the parallax data and a fitted Gaussian.
fitfunc = "%d*exp(-0.5 * pow((x-%f)/%f, 2))" % (25, p_mean, p_stdev)
fitlatex = '''%s e^{-\\frac{1}{2}
              \\left(\\frac{\\varpi-%.2f}{%.2f}
              \\right)^2}
           ''' % ('k', p_mean, p_stdev)
plot2plane(layer1='histogram', in1=cluster, x1='parallax', binsize1=0.2,
           layer2='function', fexpr2=fitfunc, thick2=2, color2='blue',
           antialias2=True,
           leglabel1='Cluster', leglabel2=fitlatex,
           legend=True, xlabel='\\varpi', legpos='.1,.9', texttype='latex',
           xcrowd=0.8, ycrowd=0.8,
           fontsize=18, grid=True )

# Plot the positions of the cluster members in 3d space.
plot2sphere(title='Pleiades 3d positions',
            layer1='mark', in1=cluster,
            lat1='ra', lon1='dec', r1="1000./parallax",
            size1=3, color1='blue', shading1='transparent', opaque1=2)