#!/usr/local/python-2.5.1/bin/python


'''
Set path above to your local Python path, please.

ColourCutter (2007 ACR) provides selected columns from high-resolution
optical and IR survey object catalogues based on position and colour.
(To make different selections from a single catalogue please use
direct ADQL queries).
See http://www2.astrogrid.org/science/documentation/workbench-data-
           analysis/science-examples/colourcutter
for a more detailed description of the scientific workings (details of
implimentation need to be updated as of 2007 Nov 21).

See http://www2.astrogrid.org/science/documentation/workbench-advanced
           /advanced-usage/scripting-user-s-guide/pydoc
for more information on scripting, in particular how to install the
initial module and create the configuration file which gives access to
MySpace.

Anita Richards 2007 Nov 21
'''
################################################################
# You must have the ACR running (e.g. the AstroGrid VO Desktop)

# USER INPUTS

# PLEASE USE EXACT NAMES/FORMAT AS IN EXAMPLES
# See end of this file for links to more information about surveys
#
# Chose name of desired survey_1 and survey_2 from list.
# For each catalogue, select two bands b1, b2 for colour selection
# INT-WFS: 'U' 'B' 'V' 'g' 'r' 'i' 'z' 'Harris_R' 'RGO_I'
# SDSSDR5: 'u' 'g' 'r' 'i' 'z'
# IPHAS:   'r' 'i' 'ha'
# TWOMASS: 'j_m' 'h_m' 'k_m'
# GLIMPSE: 'mag1' 'mag2' 'mag3' 'mag4' # The 4 IRAC bands
# UKIDSS_las: 'y' 'j_1' 'j_2' 'h' 'k' # merged point sources, AperMag3
# UKIDSS_gps: 'j' 'h2_1' 'h2_2' 'h2_3' 'h' 'k_1' 'k_2' 'k_3' # as above
#
# If you want to use UKIDSS, you must have an AstroGrid account in the
# ukidss.roe.ac.uk community (called 'ukidss' in ~/.python-acr).
#
# Insert limiting colour differences. The default example will
# give data with
# INT-WFS   i - z     > 1.1
# TWOMASS   j_m - k_m > 0.9
# e.g. replace 'INT-WFS' with 'SDSSDR5' and 'i' with 'u'

# EXAMPLES
survey_1  = 'INT-WFS'
band1_1   = 'i'
band2_1   = 'z'
magdiff_1 = 1.1

#survey_1 = 'GLIMPSE'
#band1_1 = 'mag1'
#band2_1 = 'mag2'
#magdiff_1 = 1.0

survey_2  = 'TWOMASS'
band1_2   = 'j_m'
band2_2   = 'k_m'
magdiff_2 = 0.9

# Enter search box in decimal degrees or sexagessimal (J2000)
# e.g. centre '56.75 +23.867' or '03 45 00.00 +23 52 01.2'
# or '03:45:00.00 23:52:01.2' 
# boxside is full width in degrees (no special treatment of curvature) 
# Please check coverage of selected surveys!

# INT-WFS example
centre   = '56.75 +23.867'
boxside  = 1.0

# GLIMPSE example
#centre = '18:10:00.0 -19:37:00'

# Enter cross-match radius. As many objects are extended on the scales
# of the survey resolution this may be greater than the position
# uncertainties, which are presently ignored, but could be included. 

radius = 2.0

# Output file names will be derived automatically from inputs e.g.
# survey_1:             INT_WFS_56.75+23.867_1.0_i-z1.1.vot
# survey_2:             TWOMASS_56.75+23.867_1.0_j_m-k_m0.9.vot
# crossmatched output: CCX_56.75+23.867_1.0_i-z1.1_j_m-k_m0.9.vot

#Enter the name of the MySpace directory for output files
ccoutdir = 'CCC'

# Send the final crossmatched colour table to TopCat? (Y or N)
# (if you have TopCat running).
toTopCat = 'Y'

# How to acknowledge your use of AstroGrid:
# 'This research has made use of data obtained using, or software provided 
# by, the UK's AstroGrid Virtual Observatory Project, which is funded by the 
# Science & Technology Facilities Council and through the EU's 
# Framework 6 programme'

# Use of any data discovered or accessed through AstroGrid should of
# course be mentioned as noted by the data providers.

# END OF USER INPUTS
####################################################################

# Define subroutines

# Make position and allow for RA wrapping

def makepos(centre, boxside):
    centre = centre.replace(':', ' ')
    global pos
    pos = centre.split()
    i=0
    while (i<>len(pos)):
        pos[i] = float(pos[i])
        i = i+1
    if (len(pos) > 2):
        rasex = pos[0:3]
        rasex.reverse()
        pos[0] = 15*reduce((lambda x, y: x/60.0 + y), rasex)
        if ('-' in str(pos[3])):
            pos[4] = -1.0*pos[4]
            pos[5] = -1.0*pos[5]
        decsex = pos[3:6]
        decsex.reverse()
        pos[1] = reduce((lambda x, y: x/60.0 + y), decsex)
    pos = pos[0:2]  
    pos.append(360.0)
    if (pos[0]-boxside/2.0 < 0.0):
        pos[2] = 360.0 + (pos[0] - boxside/2.0)
    pos.append(0.0)
    if (pos[0]+boxside/2.0 > 360.0):
        pos[3] = pos[0] + boxside/2.0 - 360.0
    return(pos)

# Make SQL

def makequery(survey, pos, boxside, band1, band2, magdiff):
    global query
    if (survey == 'INT-WFS'):
        query = 'SELECT o.ra, o."dec", o.objID, o.U, o.Err_U, o.g, o.Err_g, o.r, o.Err_r, o.i, o.Err_i, o.z, o.Err_z, o.B, o.Err_B, o.V, o.Err_V, o.Harris_R, o.Err_Harris_R, o.RGO_I, o.Err_RGO_I FROM obj as o where (((o.ra > %f) AND (o.ra < %f))  OR (o.ra > %f) OR (o.ra < %f)) AND (o."dec" > %f) AND (o."dec" < %f) AND (o.Err_%s > 0.0) AND (o.Err_%s < 0.5) AND ((o.%s - o.%s) > %f)' % (pos[0]-boxside/2.0, pos[0]+boxside/2.0, pos[2], pos[3], pos[1]-boxside/2.0, pos[1]+boxside/2.0, band2, band2, band1, band2, magdiff)
    if (survey == 'IPHAS'):
        query = 'SELECT o.ra, o."dec", o.objID, o.coremag_r, o.coremagerr_r, o.coremag_i,o.coremagerr_i, o.coremag_ha, o.coremagerr_ha from PhotoObjBest as o where (((o.ra > %f) AND (o.ra < %f))  OR (o.ra > %f) OR (o.ra < %f)) AND (o."dec" > %f) AND (o."dec" < %f) AND (o.coremagerr_%s > 0.0) AND (o.coremagerr_%s < 0.5) AND ((o.coremag_%s - o.coremag_%s) > %f)' % (pos[0]-boxside/2.0, pos[0]+boxside/2.0, pos[2], pos[3], pos[1]-boxside/2.0, pos[1]+boxside/2.0, band2, band2, band1, band2, magdiff)
    if (survey == 'TWOMASS'):
        query = 'SELECT o.ra, o."dec", o.designation, o.j_m, o.j_msigcom, o.h_m ,o.h_msigcom, o.k_m, o.k_msigcom FROM twomass_psc AS o where (((o.ra > %f) AND (o.ra < %f))  OR (o.ra > %f) OR (o.ra < %f)) AND (o."dec" > %f) AND (o."dec" < %f) AND (o.%ssigcom > 0.0) AND (o.%ssigcom < 0.5) AND ((o.%s - o.%s) > %f)' % (pos[0]-boxside/2.0, pos[0]+boxside/2.0, pos[2], pos[3], pos[1]-boxside/2.0, pos[1]+boxside/2.0, band2, band2, band1, band2, magdiff) 
    if (survey == 'SDSSDR5'):
        query = 'SELECT o.ra, o."dec", o.objID, o.dered_u, o.err_u, o.dered_g, o.err_g, o.dered_r, o.err_r, o.dered_i, o.err_i, o.dered_z, o.err_z FROM  PhotoObjAll AS o where (((o.ra > %f) AND (o.ra < %f))  OR (o.ra > %f) OR (o.ra < %f)) AND (o."dec" > %f) AND (o."dec" < %f) AND (o.err_%s > 0.0) AND (o.err_%s < 0.5) AND ((o.dered_%s - o.dered_%s) > %f)' % (pos[0]-boxside/2.0, pos[0]+boxside/2.0, pos[2], pos[3], pos[1]-boxside/2.0, pos[1]+boxside/2.0, band2, band2, band1, band2, magdiff)
    if (survey == 'GLIMPSE'):
        query = 'SELECT o.ra, o."dec", o.srcName, o.mag1, o.mag1_err, o.mag2, o.mag2_err, o.mag3, o.mag3_err, o.mag4, o.mag4_err FROM glimpse_hrc_inter AS o where (((o.ra > %f) AND (o.ra < %f))  OR (o.ra > %f) OR (o.ra < %f)) AND (o."dec" > %f) AND (o."dec" < %f) AND (o.%s_err > 0.0) AND (o.%s_err < 0.5) AND ((o.%s - o.%s) > %f)' % (pos[0]-boxside/2.0, pos[0]+boxside/2.0, pos[2], pos[3], pos[1]-boxside/2.0, pos[1]+boxside/2.0, band2, band2, band1, band2, magdiff)
    if (survey == 'UKIDSS_gps'):
        query = 'SELECT o.ra, o."dec", o.sourceID, o.jAperMag3, o.jAperMag3Err, o.h2_1AperMag3, o.h2_1AperMag3Err, o.h2_2AperMag3, o.h2_2AperMag3Err, o.h2_3AperMag3, o.h2_3AperMag3Err, o.hAperMag3, o.hAperMag3Err, o.k_1AperMag3, o.k_1AperMag3Err, o.k_2AperMag3, o.k_2AperMag3Err, o.k_3AperMag3, o.k_3AperMag3Err from gpsPointSource as o where (((o.ra > %f) AND (o.ra < %f))  OR (o.ra > %f) OR (o.ra < %f)) AND (o."dec" > %f) AND (o."dec" < %f) AND (o.%sAperMag3Err > 0.0) AND (o.%sAperMag3Err < 0.5) AND ((o.%sAperMag3 - o.%sAperMag3) > %f)' % (pos[0]-boxside/2.0, pos[0]+boxside/2.0, pos[2], pos[3], pos[1]-boxside/2.0, pos[1]+boxside/2.0, band2, band2, band1, band2, magdiff)
    if (survey == 'UKIDSS_las'):
        query = 'SELECT o.ra, o."dec", o.sourceID, o.yAperMag3, o.yAperMag3Err, o.j_1AperMag3, o.j_1AperMag3Err, o.j_2AperMag3, o.j_2AperMag3Err, o.hAperMag3, o.hAperMag3Err, o.kAperMag3, o.kAperMag3Err from lasPointSource as o where (((o.ra > %f) AND (o.ra < %f))  OR (o.ra > %f) OR (o.ra < %f)) AND (o."dec" > %f) AND (o."dec" < %f) AND (o.%sAperMag3Err > 0.0) AND (o.%sAperMag3Err < 0.5) AND ((o.%sAperMag3 - o.%sAperMag3) > %f)' % (pos[0]-boxside/2.0, pos[0]+boxside/2.0, pos[2], pos[3], pos[1]-boxside/2.0, pos[1]+boxside/2.0, band2, band2, band1, band2, magdiff)    
    return(query)

# Define database access routes
def setdb(survey):
    global dbivorn
    if(survey == 'INT-WFS'):
        dbivorn = 'ivo://uk.ac.cam.ast/INT-WFS/merged-object-catalogue/ceaApplication'
    if (survey == 'TWOMASS'):
        dbivorn = 'ivo://wfau.roe.ac.uk/twomass-dsa/ceaApplication'
    if (survey == 'IPHAS'):
        dbivorn = 'ivo://uk.ac.cam.ast/IPHAS/catalogue/ceaApplication'
    if (survey == 'SDSSDR5'):
        dbivorn = 'ivo://wfau.roe.ac.uk/sdssdr5-dsa/ceaApplication'
    if (survey == 'GLIMPSE'):
        dbivorn = 'ivo://wfau.roe.ac.uk/glimpse-dsa/ceaApplication'
    if (survey[0:6] == 'UKIDSS'):
        dbivorn = 'ivo://wfau.roe.ac.uk/ukidssDR2-dsa/ceaApplication'
    return(dbivorn)

# Define how to make the colours
def newcol(band1, band2):
    global newcoldef 
    if(survey == 'INT-WFS'):
        newcoldef = 'addcol %s-%s "(%s-%s)"' % (band1, band2, band1, band2) 
    if (survey == 'IPHAS'):
        newcoldef = 'addcol coremag_%s-coremag_%s "(coremag_%s-coremag_%s)"' % (band1, band2, band1, band2) 
    if (survey == 'TWOMASS'):    
        newcoldef = 'addcol %s-%s "(%s-%s)"' % (band1, band2, band1, band2) 
    if (survey == 'SDSSDR5'):
        newcoldef = 'addcol dered_%s-dered_%s "(dered_%s-dered_%s)"' % (band1, band2, band1, band2)
    if (survey == 'GLIMPSE'):
        newcoldef = 'addcol %s-%s "(%s-%s)"' % (band1, band2, band1, band2)
    if (survey[0:6] == 'UKIDSS'):
        newcoldef = 'addcol %sAperMag3-%sAperMag3 "(%sAperMag3-%sAperMag3)"' % (band1, band2, band1, band2)
    return (newcoldef)

####################################################################

# Translate position into decimal degrees if necessary and
# allow for RA wrapping
makepos(centre, boxside)

# Make position string for output filename
posstr = str(pos[1])[0:5]
if (pos[1] > 0.0):
    posstr = '+' + str(pos[1])[0:5]

posstr = str(pos[0])[0:5] + posstr + '_' + str(boxside) + '_' 

# Construct SQL and descriptions for output filename
# First survey
survey = survey_1
band1 = band1_1
band2 = band2_1
magdiff = magdiff_1

query_1 = makequery(survey, pos, boxside, band1, band2, magdiff)
dbivorn_1 = setdb(survey)
newcol_1 =  newcol(band1, band2)

ccoutdir = '#' + ccoutdir 
if (ccoutdir[-1] <> '/'):
    ccoutdir = ccoutdir + '/'
    
survey1vot =  ccoutdir + survey + '_' + posstr + band1 + '-' + band2 + str(magdiff) + '.vot'
ccoutputvot = ccoutdir + 'CCX_' + posstr + survey + band1 + '-' + band2 + str(magdiff) 

# Second survey
survey = survey_2
band1 = band1_2
band2 = band2_2
magdiff = magdiff_2

query_2 = makequery(survey, pos, boxside, band1, band2, magdiff)
dbivorn_2 = setdb(survey)
newcol_2 =  newcol(band1, band2)

survey2vot =  ccoutdir + survey + '_' + posstr + band1 + '-' + band2 + str(magdiff) + '.vot'
tempvot = ccoutdir + 'temp.vot'
ccoutputvot = ccoutputvot + survey + band1 + '-' + band2 + str(magdiff) + '.vot'

#########################################################
# Connect to ACR

from astrogrid import acr
from astrogrid import MySpace, DSA, Applications

if ((survey_1[0:6] == 'UKIDSS') or (survey_2[0:6] == 'UKIDSS')):
    acr.login('ukidss')
else:    
    acr.login()
import time

# Run first query
db = DSA(dbivorn_1)
job = db.query(query_1, ofmt='votable', saveAs=survey1vot)
print 'Querying ' + survey_1
while job.status()<>'COMPLETED':
    if (job.status() == 'ERROR'):
        print job.status()
        break
    else:
        time.sleep(10)

print survey_1 + ' query ' + job.status()        

# Run second query
db = DSA(dbivorn_2)
job = db.query(query_2, ofmt='votable', saveAs=survey2vot)
print 'Querying ' + survey_2
while job.status()<>'COMPLETED':
    if (job.status() == 'ERROR'):
        print job.status()
        break
    else:
        time.sleep(10)

print survey_2 + ' query ' + job.status()        

# Crossmatch

app = Applications('ivo://starlink.ac.uk/stilts', 'tmatch2')

app.inputs['tmatch2_in1']['value'] = survey1vot
app.inputs['tmatch2_in2']['value'] = survey2vot
app.inputs['tmatch2_values1']['value'] = 'ra dec'
app.inputs['tmatch2_values2']['value'] = 'ra dec'
app.inputs['tmatch2_params']['value'] = radius
app.inputs['tmatch2_ifmt1']['value'] = '(auto)'
app.inputs['tmatch2_ifmt2']['value'] = '(auto)'
app.inputs['tmatch2_ofmt']['value'] = 'votable-tabledata'
app.inputs.pop('tmatch2_icmd1')
app.inputs.pop('tmatch2_icmd2')
app.inputs.pop('tmatch2_ocmd')

app.outputs['tmatch2_out']['value'] = tempvot

job = app.submit()
print 'Crossmatching ' + survey_1 + ' and ' + survey_2 + ' selections'
while job.status()<>'COMPLETED':
    if (job.status() == 'ERROR'):
        print job.status()
        break
    else:
        time.sleep(10)

print 'Crossmatch ' + job.status()

# Make colours

app = Applications('ivo://starlink.ac.uk/stilts', 'tpipe')
app.inputs['tpipe_in']['value'] = tempvot
app.inputs['tpipe_cmd']['value'] = newcol_1 +'; ' + newcol_2
app.inputs['tpipe_ofmt']['value'] = 'votable-tabledata'
app.inputs['tpipe_ofmt']['value'] = 'vot'
app.outputs['tpipe_out']['value'] = ccoutputvot 
print 'Adding colour columns'
job = app.submit()
while job.status()<>'COMPLETED':
    if (job.status() == 'ERROR'):
        print job.status()
        break
    else:
        time.sleep(10)

print 'Addition of colour columns ' + job.status()

# Remove temp.vot
m = MySpace()
m.rm(tempvot)


print 'Jobs completed, results are in ' + ccoutputvot
# Send to Topcat?
if (toTopCat == 'Y'):
    print 'Preparing to send output to TopCat'
    acr.startplastic()
    acr.plastic.broadcast(ccoutputvot, 'Topcat')

# INT-WFS
# http://www.ast.cam.ac.uk/~wfcsur/science/overview/
# http://www.ast.cam.ac.uk/~wfcsur/science/fields/   # coverage

# IPHAS
# http://www2.astrogrid.org/science/science-examples-stars/iphas-catalogue-images/iphas-database-access # including coverage

# UKIDSS
# http://surveys.roe.ac.uk/wsa/index.html
# http://surveys.roe.ac.uk/wsa/dr2_gps.html  # coverage
# http://surveys.roe.ac.uk/wsa/dr2_las.html  # coverage

# GLIMPSE
# http://irsa.ipac.caltech.edu/data/SPITZER/GLIMPSE/

# 2MASS
# http://pegasus.phast.umass.edu/