#
# Example Jython script. 
# This small script retrieves some observations
# from the archive, extracts metadata information
# and executes some tasks and operations on the
# level 2 products. The results are stored in a formatted
# file (see also  asciiTableTool). Basic user interaction
# is required.
# MSP, February 2012
#
from herschel.share.fltdyn.time import *
from java.util import Date
from java.text import DateFormat
from herschel.share.fltdyn.time import FineTime
from java.text import SimpleDateFormat
from herschel.share.unit import Angle
# Set a few 'variables'
df=SimpleDateFormat("yyyy-MM-dd'T'HH:mm:ssZ")
# Put the 1st January 2000, midnight as a fine time
j2000="2000-01-01T00:00:00Z"
date2000=FineTime(df.parse(j2000.split("Z")[0]+"-0000"))
# Unit conversions
secTorad=Angle.SECONDS_ARC.to(Angle.RADIANS)
degTorad=Angle.DEGREES.to(Angle.RADIANS)
# Raw input allows you to do some nice interaction
y_n = ""
y_n = raw_input("Do you want graphics output? (y/n)")
# 
if ( y_n.rfind("y") != -1 ) | ( y_n.rfind("Y") != -1 ):
 print "Graphics output required"
 GRAPHICS=True
else:
 print "Graphics output not required"
 GRAPHICS=False
# 
file=open("/Users/msanchez/hcss_work/workshop_2011/output.csv",'w')
# 
# Write the header of the output text file
#
line="obsid,target,ra_nominal,dec_nominal,ra_fitted,dec_fitted,pm_ra,pm_dec,\
ra_offset,dec_offset,delta_y,delta_z\n"
#
file.write(line)
#
# List of observations to process
list = Int1d([\
1342181645,\
1342181646,\
1342181645,\
1342181651\
])
#
# Initialise the acces to the archive
archive = HsaReadPool()
store = ProductStorage()
store.register(archive)
#
# 'for' loop, until obsId has passed through all values within «list«
for obsId in list:
 # Explicitly call the 'garbage collector'
 System.gc()
 print "obsId =", obsId
 # Get the observation context from the archive
 obs = getObservation(obsId, useHsa = True)
 # Get some values from the products' metadata
 ranom=obs.meta["raNominal"].value
 decnom=obs.meta["decNominal"].value
 # getting proper motion data (can be done with obs.meta["pmRA"] as well)
 pmRa=obs.properMotion[0]
 pmDec=obs.properMotion[1]
 # Getting the start date 
 start=obs.startDate
 # Ellapsed time in years!
 ellaps=start.subtract(date2000)/(365*24*3600e6)
 # RA-DEC correction due to proper motion
 corr_ra=pmRa*ellaps/COS(degTorad*decnom)
 corr_dec=pmDec*ellaps
 # 
 print "target nominal RA = %lf, DEC = %lf" % (ranom, decnom)
 #
 print "PM correction RA = %lf, DEC = %lf" % (corr_ra, corr_dec)
 # Recover the level 2 "blue" map image ... use the copy/paste utility to learn!!
 # this is a "SimpleImage" object
 #
# if obsId == 1342193137:
#  myIma=obs.refs["level2"].product.refs["HPPPMAPB"].product
# else:
#  myIma=obs.refs["level2"].product.refs["HPPPMAPB"].product.refs[0].product
 myIma=obs.refs["level2"].product.refs["HPPPMAPB"].product

 #
 # Choosing graphics
 #
 if GRAPHICS:    
   Display(myIma)
 # 
 # Convert equatorial coordinates to Y, Z
 #  
 www = myIma.getWcs()
 cpos = www.getPixelCoordinates(ranom,decnom)
 # This line does a nice source fitting
 parameters = sourceFitting(elongated=False,slope=False,image=myIma,\
                            minX=cpos[1]-13.,minY=cpos[0]-13.,width=40.0,height=40.0)
 # The "parameters" object has two attributes, centerRA, centerDEC
 print "Fitted source centre RA = %lf, DEC = %lf" % (parameters.centerRA, parameters.centerDec)
 # Pass from Equinox 2000, Epoch 2011.x to Epoch 2000
 fittedRa_2000 = parameters.centerRA - corr_ra/3600
 fittedDec_2000 = parameters.centerDec - corr_dec/3600
 # Compute the ra, dec ofsets (in degrees)
 delta_ra = (ranom  - fittedRa_2000 ) * 3600
 delta_dec = (decnom - fittedDec_2000) * 3600
 #
 print "Astrometrical error  RA = %lf, DEC = %lf" % (delta_ra, delta_dec)
 # Convert this into delta_Y, delta_Z offsets
 posAng=obs.meta["posAngle"].value
 phi=degTorad*posAng
 print "Position angle = ",posAng
 delta_y=COS(phi)*delta_ra*COS(degTorad*decnom)-SIN(phi)*delta_dec
 delta_z=SIN(phi)*delta_ra*COS(degTorad*decnom)+COS(phi)*delta_dec
 #
 print "Astrometrical error  Y = %lf, Z = %lf" % (delta_y, delta_z)
 #
 line="%s,%s,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf\n" % (obsId,obs.meta["object"].value,ranom, decnom, \
                                                 parameters.centerRA, parameters.centerDec,\
                                                 corr_ra,corr_dec,delta_ra,delta_dec,delta_y,delta_z)
 #
 file.write(line)
# You can "close" the loop using a pass sentence - does nothing and can
# serve as a delimiter.
pass
# Don't forget to close the file, otherwise you'll get an empty file!!
file.close()