##### # qa - Quanta ##### # # Make a quantity var = qa.quantity(5.4,'km/s') #or var = qa.quantity('5.4km/s') print var #{'value': 5.4000000000000004, 'unit': 'km/s'} #make a vector of quantities qa.quantity([1.2,2.4,3.6],'deg') #{'unit': 'deg', 'value': array([ 1.2, 2.4, 3.6])} # Convert units print qa.convert(var,'pc/h') #{'value': 6.3000749403988587e-10, 'unit': 'pc/h'} qa.convertfreq('3mm','GHz') #{'unit': 'GHz', 'value': 99.930819333333332} # Define new units qa.define('ALMAunit','0.5 Jy') qa.convert('5 ALMAunit','Jy') #{'unit': 'Jy', 'value': 2.5} print qa.map('user') #== User ==== # ALMAunit (User defined) 5e-27 kg.s-2 # BEAM (dimensionless beam) 1 _ # DAYS (day) 86400 s # DEG (degree) 0.0174532925199 rad # DEGREES (degree) 0.0174532925199 rad # HZ (hertz) 1 s-1 # JY (jansky) 1e-26 kg.s-2 # KELVIN (kelvin) 1 K # KELVINS (kelvin) 1 K # KM (km) 1000 m # M (meter) 1 m # METERS (meter) 1 m # PASCAL (pascal) 1 m-1.kg.s-2 # PIXEL (dimensionless pixel) 1 _ # S (second) 1 s # SEC (second) 1 s # SECONDS (second) 1 s # VOLTS (volt) 1 m2.kg.s-3.A-1 # YEAR (year) 31557600 s # YEARS (year) 31557600 s # pix (pixel units) 1 # Do quantity math qa.add('5km','6000m') #{'unit': 'km', 'value': 11.0} qa.canonical(qa.add('5km','6000m')) # use canonical distance units #{'unit': 'm', 'value': 11000.0} qa.div('5m','3s') #{'unit': 'm/(s)', 'value': 1.6666666666666667} # Play with angles ang1=qa.quantity('5.7.12.345678') print ang1 #{'value': 5.1200960216666669, 'unit': 'deg'} qa.angle(ang1) # '+005.07.12' qa.angle(ang1,prec=7) #'+005.07.12.3' qa.angle(ang1,form='time') #'00:20:29' # Play with time now=qa.quantity('today') #{'unit': 'd', 'value': 54385.802151331016} qa.time(now,form='dmy') #'12-Oct-2007/19:15:06' qa.time(now,form='fits') #'2007-10-12T19:15:06' qa.time(now,form='ymd') #'2007/10/12/19:15:06' qa.time(now,form=['ymd','local']) #'2007/10/12/13:15:06' # Play with constants print qa.map('const') #== Constants ==== #pi 3.14.. 3.14159 #ee 2.71.. 2.71828 #c light vel. 2.99792e+08 m/s #G grav. const 6.67259e-11 N.m2/kg2 #h Planck const 6.62608e-34 J.s #HI HI line 1420.41 MHz #R gas const 8.31451 J/K/mol #NA Avogadro # 6.02214e+23 mol-1 #e electron charge 1.60218e-19 C #mp proton mass 1.67262e-27 kg #mp_me mp/me 1836.15 #mu0 permeability vac. 1.25664e-06 H/m #eps0 permittivity vac. 8.85419e-12 F/m #k Boltzmann const 1.38066e-23 J/K #F Faraday const 96485.3 C/mol #me electron mass 9.10939e-31 kg #re electron radius 2.8179e-15 m #a0 Bohr's radius 5.2918e-11 m #R0 solar radius 6.9599e+08 m #k2 IAU grav. const^2 0.000295912 AU3/d2/S0 boltzmann = qa.constants('k') print boltzmann {'value': 1.3806577987510647e-23, 'unit': 'J/K'} ##### # MEASURES ##### # Set frame information (when, where) #me.epoch('utc','now') # current time meetingend=me.epoch('utc','2007/10/18/17:00:00') # time at the end of this meeting print meetingend #{'type': 'epoch', 'm0': {'value': 54391.708333333336, 'unit': 'd'}, 'refer': 'UTC'} me.show(meetingend) #'54391.7 UTC' # show measure as a string qa.time(x,form='ymd') #'2007/10/18/17:00:00' #me.listcodes(me.epoch()) # to see list of available reference frames # for epoch #me.doframe - takes measure information (epoch, location) me.doframe(me.epoch('utc','now')) # set time to current for calculations #me.observatory('VLA') # specify location to a known observatory #me.obslist() # gives the current list of observatories known me.doframe(me.observatory('VLA')) # set location for calculations print me.showframe() #Frame: Epoch: 54385::19:24:15.0220 (TDB = 54385.8, UT1 = 54385.8, TT = 54385.8) # Position: [-1.60119e+06, -5.04198e+06, 3.55488e+06] # (Longitude = -1.87829 Latitude = 0.591675) # Direction: [0, 0, 1] # (J2000 = [0, 90] deg) # Get direction to a source # me.listcodes(me.direction()) # to see list of available reference frames # includes planets, SUN, MOON satdir=me.direction('Saturn') # the direction to Saturn; use direction for planets/comets caldir=me.source ('J181945.3-552120') # use .source for known calibrators (print me.sourcelist()) #convert Saturn's direction to Az-El me.measure(me.direction('Saturn'),'AZEL') # format it nicely me.dirshow(me.measure(me.direction('Saturn'),'AZEL')) # '[-79.6479, 4.51178] deg AZEL' #do it for J2000 me.dirshow(me.measure(me.direction('Saturn'),'J2000')) # '[156.889, 11.1405] deg J2000' qa.formxxx('156.889deg','hms') #'10:27:33.36' me.riseset('moon') # See in logger: #LAST of rise= 09:02:51.53 #LAST of set= 19:18:18.75 #UTC of rise= Epoch: 54384::14:54:03.3707 #UTC of set= Epoch: 54385::01:07:49.7606 #also me.riseset('J2000 19h30m00 -20d00m00') ##### # TABLES & MATPLOTLIB ##### #MSs are just sets of tables tb.open('ngc5921.ms') # this opens the Main table tb.colnames() tb.getcol('SCAN_NUMBER') # can do same thing with image files, calibration tables, etc #Ingest a tabular/CSV ASCII file # !more Vertex20070606_16_55.dat # !more Vertex_ascii.dat # show change to headers/can also use autoheader # Look at parameters to import tb.fromascii? tb.fromascii('Vertex_pointing.tb','Vertex_ascii.dat') tb.colnames() #grab the on-TE Elevation readings onteel=tb.getcol('ONTEEL') #grab the mid-point TE Elevation readings midteel=tb.getcol('MIDTEEL') #get difference and convert to arcseconds diffarcsec=(onteel-midteel)*3600 # turn on plotter pl.ion() # pl.plot(diffarcsec,label='diff') #get difference between the on-TE reading and he mid-point TE from the next reading diffarcsecb=(onteel[0:12950]-midteel[1:12951])*3600 # pl.plot(diffarcsecb,label='off by one') # Add plot information pl.title('ONTEEL-MIDTEEL') pl.legend() pl.xlabel('Time') pl.ylabel('Difference (arcseconds)') # !rm -rf outtable !more readcorrdump.py execfile 'readcorrdump.py' !more plotcorrdump.py pl.clf() execfile 'plotcorrdump.py' #See matplotlib site for full range ##### # COORDSYS ##### # Retrieve the coordinate system from an image ia.open('ngc5921_task.image') csys=ia.coordsys() #csys. to see all methods # utility csys.summary() # output goes to logger # shows header and axis information csys.names() csys.naxes() csys.axiscoordinatetypes() csys.referencepixel() csys.units() csys.projection('all') - see all projections # get/set # Note - these just change the header information csys.telescope() csys.observer() csys.setobserver('MCMULLIN') csys.summary() # see change in observer # set new coordinate system for the image ia.setcoordsys(csys.torecord()) ia.summary() # see if change is recorded csys.done() ia.close() # play with the increment !cp ngc5921_task.image original_task.image ia.open('ngc5921_task.image') csys=ia.coordsys() incr=csys.increment('n') incrn=incr['numeric'] for i in range(len(incrn)): incrn[i]=2*incrn[i] incr['numeric']=incrn csys.setincrement(value=incr) ia.setcoordsys(csys.torecord()) csys.done() ia.done() # Now look at ngc5921_task.image and original_task.image in the viewer # conversion #find the Right Ascension at a given pixel ra=csys.toworld(value=[128,128])['numeric'][0] qa.angle(str(ra)+'rad',form='time') dec=csys.toworld(value=[128,128])['numeric'][1] qa.angle(str(dec)+'rad',form='dms') # Change coordinate frame J2000->B1950 ia.open('ngc5921_task.image') # J2000 coordinate csys=ia.coordsys() csys.setreferencecode('B1950') csys.referencecode() ia.regrid(outfile='ngc5921_regridded',csys=csys.torecord()) #look at it in the viewer