# # 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()