|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 from taskinit import * 00002 from getazel import * 00003 import pylab as pl 00004 import os 00005 import datetime 00006 from matplotlib.dates import DateFormatter 00007 # plotting azimuth and elevation of source(s) 00008 00009 def plotsrcazel(planet='', srcfile='', date='', obs='ALMA', plotsun=False, plottype='', saveplot='', elimit='3.0deg'): 00010 ''' 00011 plot Azimuth/elevation track of a source (or list of sources) 00012 on given date. For solar system objects known to Measures(this excludes minor planets, Jovian moons, etc), 00013 name is sufficient but for other sources 00014 RAs and Decs are required. 00015 (requires getazel.py) 00016 00017 Input: 00018 planet - a planet (optional) 00019 srcfile - (optional) An ascii file with source_name coordinates (J2000) seperated by 00020 a white space. A line with '#' will be skipped. 00021 The solar system can also be include without coodinates. 00022 00023 0423-013 4h23m0s -1d20m0s 00024 #3c273 12h29m0s 2d3m0s 00025 3c279 12h56m0s -5d47m0s 00026 Mars 00027 00028 date - 'YYYY/MM/DD' if not specified it will be prompted to enter later 00029 obs - observatory name (ALMA,EVLA,etc) 00030 plotsun - True will plot the Sun also 00031 plottype - Az, El, or both (if left blank it will ask to enter later) 00032 saveplot - save to the plot to ps file with the name specified in this paramter 00033 elimit - draw horizontal line to indicate elevation limit. 00034 00035 Examples: 00036 plotsrcazel('Uranus') # enter date, plottype when prompted 00037 plotsrcazel(planet='Uranus', date='2012/06/01', plottype='both') 00038 plotsrcazel(planet='Uranus', srcfile='calandtarget.txt') #plot Uranus and sources defined 00039 # in calandtarget.txt 00040 00041 version 2012.04.20 TT 00042 ''' 00043 00044 00045 #def plotazel(date=''): 00046 # date = '2010/02/01' 00047 inpd ="" 00048 if date=='': 00049 inpd=raw_input('date? YYYY/MM/DD or hit return to use today\'s date:') 00050 if inpd == "": 00051 date=qa.time('today', form=['ymd','no_time']) 00052 else: 00053 date=inpd 00054 else: 00055 date=date 00056 print "Use date:", date 00057 00058 intz=raw_input('Show in UTC, CLT,or LST? ') 00059 tz='UTC' 00060 if intz!='': 00061 check_intz = [intz.upper()==tl for tl in ['UTC','CLT','LST']] 00062 if not any(check_intz): 00063 raise Exception, "Input error: should be 'UTC','CLT',or 'LST'" 00064 tz=intz.upper() 00065 00066 # read source list from a file 00067 if srcfile!="": 00068 if not os.path.exists(srcfile): 00069 raise IOError, "%s does not exist!" % srcfile 00070 readfromfile = True 00071 else: 00072 readfromfile = False 00073 insrc = True 00074 # plot planets 00075 knownsrcs=['SUN', 'MOON', 'MERCURY','VENUS','MARS','JUPITER','SATURN','URANUS','NEPTUNE','PLUTO'] 00076 plotplanets=False 00077 if planet!='': 00078 for ksrc in knownsrcs: 00079 if planet.upper() == ksrc: 00080 plotplanets= True 00081 #plotplanets = True 00082 00083 if plottype=='': 00084 # default plottype 00085 #plottype='az' 00086 plottype='both' 00087 00088 inptype=raw_input('plot type? (el, az, or both):') 00089 00090 if inptype=='el' or inptype=='az': 00091 plottype=inptype 00092 00093 00094 #observatory 00095 # should work any observatory name in 00096 # me.obslist() 00097 #obs='ALMA' 00098 #obs='VLA' 00099 #obs='NRO' 00100 if obs.upper()=='ALMA': 00101 elimit_hard=3.0 00102 else: 00103 elimit_hard=0.0 00104 #user specified EL limit (draws horizontal line) 00105 elimitv=qa.quantity(elimit)['value'] 00106 00107 # source list from a file or use hardcoded one in here 00108 # or can be read from source in casa database 00109 # (but maybe they are mostly northern hemisphere soruces???) 00110 00111 # read calibrator list as a file 00112 # format 00113 # name ra dec 00114 # 0522-364 5h23m0s -36d28m0s 00115 # 00116 # 00117 srclist={} 00118 srclist['name']=[] 00119 srclist['ra']=[] 00120 srclist['dec']=[] 00121 00122 if readfromfile: 00123 f = open('calibrators.list','r') 00124 while f: 00125 ln = f.readline() 00126 if len(ln) ==0: 00127 break 00128 if ln.rfind('#')==0: 00129 # skip comment 00130 pass 00131 else: 00132 #(name,ra,dec) = ln.split() 00133 src = ln.split() 00134 if len(src)==3: 00135 srclist['name'].append(src[0]) 00136 srclist['ra'].append(src[1]) 00137 srclist['dec'].append(src[2]) 00138 elif len(src)==1: 00139 # this part still does not work (need to modify getazel) 00140 knownEphemSrcs=me.listcodes(me.direction())['extra'].tolist() 00141 if any(i == src[0].upper() for i in knownEphemSrcs): 00142 srclist['name'].append(src[0]) 00143 srclist['ra'].append(src[0]) 00144 srclist['dec'].append(' ') 00145 if insrc: 00146 inpos = raw_input("Enter name xxhxxmxx.xs +/-xxdxxmxx.xs: (return to skip)") 00147 if inpos=='': 00148 if planet=='' and srcfile=='' and (not plotsun): 00149 raise Exception, "No source is specified!!!" 00150 pass 00151 else: 00152 inpos.replace("\'","") 00153 srcpos = inpos.split() 00154 srclist['name'].append(srcpos[0]) 00155 srclist['ra'].append(srcpos[1]) 00156 srclist['dec'].append(srcpos[2]) 00157 00158 if plotsun: 00159 # add sun 00160 srclist['name'].append('SUN') 00161 srclist['ra'].append('SUN') 00162 srclist['dec'].append('') 00163 00164 if plotplanets: 00165 srclist['name'].append(planet) 00166 srclist['ra'].append(planet) 00167 srclist['dec'].append('') 00168 00169 # 00170 last = False 00171 fig = pl.figure() 00172 fig.clf() 00173 tmfmt = DateFormatter('%H') 00174 seq = range(len(srclist['name'])) 00175 print "Number of sources to be plotted:", len(srclist['name']) 00176 for i in seq: 00177 if i == seq[-1] : 00178 last = True 00179 srcname=srclist['name'][i] 00180 srcradec=srclist['ra'][i]+' '+srclist['dec'][i] 00181 00182 srcCoord='J2000 '+srcradec 00183 casalog.post('Plotting for '+srcname) 00184 (az,el,tm,tutc) = getazel(obs,srcname,srcCoord,date,tz) 00185 00186 eld = el*180./pi 00187 azd = az*180./pi 00188 #elevation limit 00189 import numpy 00190 for iel in range(len(eld)): 00191 if eld[iel] < elimit_hard: 00192 eld[iel]=-1.0 00193 azd[iel]=360. 00194 00195 #print "max(eld)=",max(eld) 00196 #print "max(azd)=",max(azd) 00197 00198 tx = tm 00199 #for t in tx: 00200 # tq = qa.time(str(t)+'d',form='ymd') 00201 # tqymd=tq.split('/') 00202 # tqhms=tqymd[-1].split(':') 00203 # dt=datetime.datetime(int(tqymd[0]),int(tqymd[1]),int(tqymd[2]),int(tqhms[0]),int(tqhms[1]),int(tqhms[2])) 00204 # print "date", dt.year, dt.month, dt.day 00205 # print "time", dt.hour 00206 00207 #toff = int(min(tx)) 00208 #if tz=='LST': 00209 # tlab = tutc 00210 #tofflab = int(min(tlab)) 00211 #tx = tutc 00212 # toff = int(min(tx)) 00213 # tofflab = toff 00214 # print "tx=",tx 00215 # print "toff=",toff 00216 # print "tofflab=",tofflab 00217 #else: 00218 # tofflab = toff 00219 00220 #tx -= toff 00221 #tx = tx*24. 00222 #print "tx[0]=", tx[0], " tx[-1]=",tx[-1] 00223 #for itm in range(len(tx)): 00224 # if tx[itm] >24.0: 00225 # tx[itm]-=24.0 00226 #print "tx[0]=", tx[0], " tx[-1]=",tx[-1] 00227 #timlab = qa.time(qa.quantity(tofflab,'d'), form=['ymd', 'no_time']) 00228 #timlab += ' %s' % obs 00229 timlab = qa.time(qa.quantity(int(min(tx)),'d'), form=['ymd', 'no_time']) 00230 timlab += ' %s' % obs 00231 00232 defaultcolors=['b','g','r','c','m','y','k'] 00233 cindx = fmod(i,len(defaultcolors)) 00234 00235 if plottype=='both' or plottype=='el': 00236 #print "plotting El vs. time.." 00237 if i==0: 00238 if plottype=='both': 00239 ax=fig.add_subplot(2,1,1) 00240 if plottype=='el': 00241 #pl.subplot(1,1,1) 00242 ax=fig.add_subplot(1,1,1) 00243 ax.xaxis.set_major_formatter(tmfmt) 00244 if last: 00245 #ax.xlabel(tz) 00246 #ax.ylabel('EL(deg.)') 00247 ax.set_xlabel(tz) 00248 ax.set_ylabel('EL(deg.)') 00249 #pl.xlim(min(tx),max(tx)) 00250 ax.plot_date(tx,eld,marker='.',markersize=2.0,color=defaultcolors[cindx]) 00251 if last: 00252 if elimitv>elimit_hard: 00253 ax.axhline(elimitv, linestyle='dashed') 00254 # get max for label 00255 maxel = max(eld) 00256 maxidx = argmax(eld) 00257 if last: 00258 #if tz=='LST': 00259 # pl.xlim(tx[0],tx[-1]) 00260 #else: 00261 # pl.xlim(0,24.) 00262 ax.set_ylim(elimit_hard,90.) 00263 ax.set_title(timlab) 00264 ax.text(tx[maxidx],maxel*1.01,srcname,color=defaultcolors[cindx]) 00265 00266 if plottype=='both' or plottype=='az': 00267 #print "plotting Az vs. time.." 00268 if i==0: 00269 if plottype=='both': 00270 ax2=fig.add_subplot(2,1,2) 00271 if plottype=='az': 00272 ax2=fig.add_subplot(1,1,1) 00273 ax2.xaxis.set_major_formatter(tmfmt) 00274 ax2.plot_date(tx,azd,marker='.',markersize=2.0,color=defaultcolors[cindx]) 00275 maxaz = max(azd) 00276 maxazidx = argmax(azd) 00277 if last: 00278 ax2.set_xlabel(tz) 00279 ax2.set_ylabel('Az(deg.)') 00280 #if tz=='LST': 00281 # pl.xlim(tx[0],tx[-1]) 00282 #else: 00283 # pl.xlim(0,24.) 00284 #pl.xlim(0,24.) 00285 ax2.set_ylim(-180.,180.) 00286 #ax.text(tx[maxidx],*1.01,srcname,color=defaultcolors[cindx]) 00287 print srcname," max EL:", maxel 00288 00289 #pl.draw() 00290 if saveplot!='': 00291 fig.savefig(saveplot)
1.8.0
