|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 # setjy helper functions 00002 from casac import casac 00003 import os 00004 import sys 00005 import shutil 00006 import odict 00007 import numpy 00008 00009 class ss_setjy_helper: 00010 def __init__(self,imtool, vis, casalog=None): 00011 self.im = imtool 00012 self.vis = vis 00013 if not casalog: 00014 casalog = casac.logsink() 00015 self._casalog = casalog 00016 00017 def setSolarObjectJy(self,field,spw,scalebychan, timerange,observation, scan, useephemdir): 00018 """ 00019 Set flux density of a solar system object using Bryan Butler's new 00020 python model calculation code. 00021 A single time stamp (frist time stamp of MS after selections are applied) is 00022 currently used per execution. For flux observation done in a long time span 00023 may need to run multiple setjy with selections by time range (or scans). 00024 """ 00025 retval = True 00026 cleanupcomps = True # leave genenerated cl files 00027 00028 #from taskinit import * 00029 from taskinit import gentools 00030 qa = casac.quanta() 00031 00032 00033 (myms, mytb, mycl, myme) = gentools(['ms','tb','cl','me']) 00034 # prepare parameters need to pass to the Bryan's code 00035 # make ms selections 00036 # get source name 00037 # get time ranges 00038 # get spwused and get frequency ranges 00039 00040 sel={} 00041 sel['field']=field 00042 sel['spw']=spw 00043 sel['timerange']=timerange 00044 sel['observation']=str(observation) 00045 sel['scan']=scan 00046 00047 measframes=['REST','LSRK','LSRD','BARY','GEO','TOPO','GALACTO','LGROUP','CMB'] 00048 myms.open(self.vis) 00049 myms.msselect(sel,False) 00050 scansummary=myms.getscansummary() 00051 nscan=len(scansummary.keys()) 00052 fieldids=myms.msselectedindices()['field'] 00053 obsids=myms.msselectedindices()['observationid'] 00054 myms.close() 00055 00056 mytb.open(self.vis+'/OBSERVATION') 00057 if len(obsids)==0: 00058 getrow=0 00059 else: 00060 getrow=obsids[0] 00061 observatory=mytb.getcell('TELESCOPE_NAME',getrow) 00062 mytb.close() 00063 00064 mytb.open(self.vis+'/FIELD') 00065 if len(fieldids)==0: 00066 fieldids = range(mytb.nrows()) 00067 # frame reference for field position 00068 phdir_info=mytb.getcolkeyword("PHASE_DIR","MEASINFO") 00069 if phdir_info.has_key('Ref'): 00070 fieldref=phdir_info['Ref'] 00071 elif phdir_info.has_key('TabRefTypes'): 00072 colnames=mytb.colnames() 00073 for col in colnames: 00074 if col=='PhaseDir_Ref': 00075 fieldrefind=mytb.getcell(col,fieldids[0]) 00076 fieldref=phdir_info['TabRefTypes'][fieldrefind] 00077 else: 00078 fieldref='' 00079 srcnames={} 00080 fielddirs={} 00081 ftimes={} 00082 for fid in fieldids: 00083 #srcnames.append(mytb.getcell('NAME',int(fid))) 00084 srcnames[fid]=(mytb.getcell('NAME',int(fid))) 00085 fielddirs[fid]=(mytb.getcell('PHASE_DIR',int(fid))) 00086 ftimes[fid]=(mytb.getcell('TIME',int(fid))) 00087 mytb.close() 00088 # need to get a list of time 00089 # but for now for test just get a time centroid of the all scans 00090 # get time range 00091 # Also, this needs to be done per source 00092 00093 # gather freq info from Spw table 00094 mytb.open(self.vis+'/SPECTRAL_WINDOW') 00095 nspw = mytb.nrows() 00096 freqcol=mytb.getvarcol('CHAN_FREQ') 00097 freqwcol=mytb.getvarcol('CHAN_WIDTH') 00098 fmeasrefcol=mytb.getcol('MEAS_FREQ_REF') 00099 reffreqs=mytb.getcol('REF_FREQUENCY') 00100 mytb.close() 00101 00102 # store all parameters need to call solar_system_fd for all sources 00103 inparams={} 00104 validfids = [] # keep track of valid fid that has data (after selection) 00105 # if same source name ....need handle this... 00106 for fid in fieldids: 00107 sel['field']=str(fid) 00108 myms.open(self.vis) 00109 #reset the selection 00110 try: 00111 status=myms.msselect(sel) 00112 except: 00113 #skip this field 00114 continue 00115 if not status: 00116 continue 00117 00118 validfids.append(fid) 00119 trange=myms.range('time') 00120 if not inparams.has_key(srcnames[fid]): 00121 inparams[srcnames[fid]]={} 00122 00123 #tc = (trange['time'][0]+trange['time'][1])/2. #in sec. 00124 # use first timestamp to be consistent with the ALMA Control 00125 # old setjy (Butler-JPL-Horizons 2010) seems to be using 00126 # time in FIELD... but here is first selected time in main table 00127 tc = trange['time'][0] #in sec. 00128 if inparams[srcnames[fid]].has_key('mjd'): 00129 inparams[srcnames[fid]]['mjds'][0].append([myme.epoch('utc',qa.quantity(tc,'s'))['m0']['value']]) 00130 else: 00131 inparams[srcnames[fid]]['mjds']=[myme.epoch('utc',qa.quantity(tc,'s'))['m0']['value']] 00132 00133 # somehow it gives you duplicated ids .... so need to uniquify 00134 selspws= list(set(myms.msselectedindices()['spw'])) 00135 # make sure it is int rather than numpy.int32, etc. 00136 selspws = [int(ispw) for ispw in selspws] 00137 inparams[srcnames[fid]]['spwids']= selspws if len(selspws)!=0 else range(nspw) 00138 00139 #create a list of freq ranges with selected spws 00140 # should worry about freq order??? 00141 freqlist=[] 00142 framelist=[] 00143 for idx in inparams[srcnames[fid]]['spwids']: 00144 freqs=freqcol['r'+str(idx+1)] 00145 freqws=freqwcol['r'+str(idx+1)] 00146 fmeasref=fmeasrefcol[idx] 00147 00148 #if scalebychan=T, this has to be min, max determined from 00149 # chan_freq(channel center)+/- chan width. 00150 if scalebychan: 00151 # pack into list of list of list (freqlist[nf][nspw]) 00152 perspwfreqlist=[] 00153 for nf in range(len(freqs)): 00154 fl = freqs[nf][0]-freqws[nf][0]/2. 00155 fh = freqs[nf][0]+freqws[nf][0]/2. 00156 perspwfreqlist.append([fl,fh]) 00157 freqlist.append(perspwfreqlist) 00158 else: 00159 if (len(freqs)==1): 00160 fl = freqs[0][0]-freqws[0][0]/2. 00161 fh = freqs[0][0]+freqws[0][0]/2. 00162 freqlist.append([float(fl),float(fh)]) 00163 else: 00164 freqlist.append([float(min(freqs)),float(max(freqs))]) 00165 00166 framelist.append(measframes[fmeasref]) 00167 inparams[srcnames[fid]]['freqlist']=freqlist 00168 inparams[srcnames[fid]]['framelist']=framelist 00169 inparams[srcnames[fid]]['reffreqs']=reffreqs 00170 myms.close() 00171 00172 00173 # call Bryan's code 00174 # errcode: list of list - inner list - each values for range of freqs 00175 # flluxes: list of list 00176 # fluxerrs: 00177 # size: [majoraxis, minoraxis, pa] 00178 # direction: direction for each time stamp 00179 # 00180 import solar_system_setjy as ss_setjy 00181 retdict={} 00182 #for src in srcnames: 00183 for vfid in validfids: 00184 src=srcnames[vfid] 00185 mjds=inparams[src]['mjds'] 00186 fluxes=[] 00187 # call solar_system_fd() per spw (for scalebychan freqlist has an extra dimention) 00188 nspwused=len(inparams[src]['freqlist']) 00189 # warning for many channels but it is really depends on the source 00190 if scalebychan: 00191 maxnf=0 00192 for ispw in range(nspwused): 00193 nf = inparams[src]['freqlist'][ispw] 00194 maxnf = max(nf,maxnf) 00195 if maxnf >= 3840 and src.upper()!="MARS": # mars shoulde be ok 00196 self._casalog.post("Processing %s spw(s) with at least some of them are a few 1000 channels or more. This may takes \ 00197 many minutes (>3min per spw for 3840 channels) in some cases. Please be patient." % nspwused,"WARN") 00198 00199 for i in range(nspwused): # corresponds to n spw 00200 if type(freqlist[0][0])==list: 00201 infreqs=inparams[src]['freqlist'][i] 00202 else: 00203 infreqs=[inparams[src]['freqlist'][i]] 00204 self._casalog.post("Calling solar_system_fd: %s for spw%s freqs=%s" % (src, i,freqlist[i]),'DEBUG1') 00205 (errcodes, subfluxes, fluxerrs, sizes, dirs)=\ 00206 ss_setjy.solar_system_fd(source_name=src, MJDs=mjds, frequencies=infreqs, observatory=observatory, casalog=self._casalog) 00207 # for old code 00208 #(errcodes, subfluxes, fluxerrs, sizes, dirs)=\ 00209 # ss_setjy.solar_system_fd(source_name=src, MJDs=mjds, frequencies=infreqs, casalog=self._casalog) 00210 # for nf freq ranges, nt mjds 00211 # errcodes[nf][nt], subfluxes[nf][nt], fluxerrs[nf][nt], sizes[nt], dirs[nt] 00212 self._casalog.post("+++++ solar_system_fd() returned values +++++", 'DEBUG1') 00213 self._casalog.post(" fluxes(fds)=%s" % subfluxes, 'DEBUG1') 00214 self._casalog.post(" sizes=%s" % sizes, 'DEBUG1') 00215 self._casalog.post(" directions=%s\n" % dirs, 'DEBUG1') 00216 00217 # packed fluxes for all spws 00218 fluxes.append(subfluxes) 00219 # fluxes has fluxes[nf][nt][nspw] 00220 00221 # ------------------------------------------------------------------------ 00222 # For testing with hardcoded values without calling solar_system_fd()... 00223 #errcodes=[[0,0,0,0,0]] 00224 #fluxes=[[26.40653147,65.23839313,65.23382757,65.80638802,69.33396562]] 00225 #fluxerrs=[[0.0,0.0,0.0,0.0,0.0]] 00226 #sizes=[[3.6228991032674371,3.6228991032674371,0.0]] 00227 #dirs=[{'m0': {'unit': 'rad', 'value': 0.0}, 'm1': {'unit': 'rad', 'value': 0.0}, 'refer': 'J2000', 'type': 'direction'}] 00228 # ------------------------------------------------------------------------ 00229 # local params for selected src 00230 framelist=inparams[src]['framelist'] 00231 freqlist=inparams[src]['freqlist'] 00232 reffreqs=inparams[src]['reffreqs'] 00233 spwids=inparams[src]['spwids'] 00234 00235 clrecs=odict.odict() 00236 labels = [] 00237 # loop for over for multiple directions (=multiple MJDs) for a given src 00238 for i in range(len(dirs)): # this is currently only length of 1 since no mutliple timpstamps were used 00239 # check errcode - error code would be there per flux per time (dir) 00240 reterr=testerrs(errcodes[i],src) 00241 if reterr == 2: 00242 continue 00243 00244 #dirstring = [dirs[i]['refer'],qa.tos(dirs[i]['m0']),qa.tos(dirs[i]['m1'])] 00245 if useephemdir: 00246 dirstring = [dirs[i]['refer'],qa.tos(dirs[i]['m0']),qa.tos(dirs[i]['m1'])] 00247 else: 00248 #dirstring = [fieldref, str(fielddirs[fieldids[0]][0][0])+'rad', str(fielddirs[fieldids[0]][1][0])+'rad'] 00249 # extract field direction of first id of the selected field ids 00250 dirstring = [fieldref, "%.18frad" % (fielddirs[fieldids[0]][0][0]), "%.18frad" % (fielddirs[fieldids[0]][1][0])] 00251 #print "dirstring=",dirstring 00252 00253 # setup componentlists 00254 # need to set per dir 00255 # if scalebychan=F, len(freqs) corresponds to nspw selected 00256 # Need to put in for-loop to create cl for each time stamp? or scan? 00257 00258 #clpath='/tmp/' 00259 clpath='./' 00260 self.clnamelist=[] 00261 for j in range(len(freqlist)): # loop over nspw 00262 freqlabel = '%.3fGHz' % (reffreqs[int(spwids[j])]/1.e9) 00263 tmlabel = '%.1fd' % (tc/86400.) 00264 clabel = src+'_spw'+str(spwids[j])+'_'+freqlabel+'_'+tmlabel 00265 clname = clpath+clabel+'.cl' 00266 00267 if(os.path.exists(clname)): 00268 shutil.rmtree(clname) 00269 #print "addcomponent...for j=",j," flux=",fluxes[j][i] 00270 if scalebychan: 00271 index= 2.0 00272 sptype = 'spectral index' 00273 else: 00274 index= 0.0 00275 sptype = 'constant' 00276 # adjust to change in returned flux shape. An extra [] now seems to be gone. 2012-09-27 00277 iflux=fluxes[j][i][0] 00278 self._casalog.post("addcomponent with flux=%s at frequency=%s" %\ 00279 # (fluxes[j][i][0],str(reffreqs[int(spwids[j])]/1.e9)+'GHz'), 'INFO1') 00280 (iflux,str(reffreqs[int(spwids[j])]/1.e9)+'GHz'), 'INFO1') 00281 #print "addcomponent with flux=%s at frequency=%s" %\ 00282 # (iflux,str(reffreqs[int(spwids[j])]/1.e9)+'GHz') 00283 # i - time stamps = 0 for now, j = a freq range 00284 infreq=freqlist[j][0][0] if type(freqlist[j][0])==list else freqlist[j][0] 00285 00286 #mycl.addcomponent(flux=fluxes[j][i][0],fluxunit='Jy', polarization="Stokes", dir=dirstring, 00287 mycl.addcomponent(iflux,fluxunit='Jy', polarization="Stokes", dir=dirstring, 00288 shape='disk', majoraxis=str(sizes[i][0])+'arcsec', minoraxis=str(sizes[i][1])+'arcsec', 00289 # positionangle=str(sizes[i][2])+'deg', freq=[framelist[j],str(reffreqs[int(spwids[j])])+'Hz'], 00290 # positionangle=str(sizes[i][2])+'deg', freq=[framelist[j],str((freqlist[j][0]+freqlist[j][1])/2)+'Hz'], 00291 positionangle=str(sizes[i][2])+'deg', freq=[framelist[j],str(infreq)+'Hz'], 00292 spectrumtype=sptype, index=index, label=clabel) 00293 # if it's list of fluxes try to put in tabular form 00294 if type(fluxes[j][i]) ==list and len(fluxes[j][i])> 1: 00295 #print "framelist[j]=",framelist[j] 00296 if type(freqlist[j][0])==list and len(freqlist[j][0])>1: 00297 freqs=[] 00298 for fr in freqlist[j]: 00299 freqs.append((fr[1]+fr[0])/2) 00300 else: 00301 freqs=freqlist[j] 00302 clind = mycl.length() - 1 00303 mycl.setspectrum(which=clind, type='tabular', tabularfreqs=freqs, tabularflux=fluxes[j][0], 00304 tabularframe=framelist[j]) 00305 mycl.rename(clname) 00306 #put in a record for log output 00307 clrecs[clabel] = mycl.torecord() 00308 mycl.close(False) # False for not to send a warning message 00309 mycl.done() 00310 00311 # finally, put the componentlist as model 00312 self.im.selectvis(spw=spwids[j],field=field,observation=observation,time=timerange) 00313 self.im.ft(complist=clname) 00314 #debug: set locally saved 2010-version component list instead 00315 #cl2010='mod_setjy_spw0_Titan_230.543GHz55674.1d.cl' 00316 #print "Using complist=",cl2010 00317 #self.im.ft(complist=cl2010) 00318 00319 #if cleanupcomps: 00320 # shutil.rmtree(clname) 00321 self.clnamelist.append(clname) 00322 # end of for loop over fields 00323 msg="Using channel dependent " if scalebychan else "Using spw dependent " 00324 00325 self._casalog.post(msg+" flux densities") 00326 self._reportoLog(clrecs,self._casalog) 00327 self._updateHistory(clrecs,self.vis) 00328 00329 return retval 00330 00331 def getclnamelist(self): 00332 return self.clnamelist 00333 00334 def _reportoLog(self,clrecs,casalog): 00335 """ 00336 send model parameters to log 00337 """ 00338 #print "clrecs=", clrecs 00339 for ky in clrecs.keys(): 00340 # expect only one component for each 00341 comp = clrecs[ky]['component0'] 00342 srcn = ky.split('_')[0] 00343 ispw = ky.split('_')[1] 00344 msg=" direction set in the componentlist: RA=%s rad, Dec%s rad" %\ 00345 (float('%.6g' % comp['shape']['direction']['m0']['value']), 00346 float('%.6g' % comp['shape']['direction']['m1']['value'])) 00347 casalog.post(msg,'INFO2') 00348 00349 msg=" %s: %s Flux:[I=%s,Q=%s,U=%s,V=%s] +/- [I=%s,Q=%s,U=%s,V=%s] Jy" %\ 00350 (srcn, ispw, float('%.5g' % comp['flux']['value'][0]), 00351 float('%.5g' % comp['flux']['value'][1]), 00352 float('%.5g' % comp['flux']['value'][2]), 00353 float('%.5g' % comp['flux']['value'][3]), 00354 float('%.5g' % comp['flux']['error'][0]), 00355 float('%.5g' % comp['flux']['error'][1]), 00356 float('%.5g' % comp['flux']['error'][2]), 00357 float('%.5g' % comp['flux']['error'][3])) 00358 casalog.post(msg, 'INFO') 00359 00360 def _updateHistory(self,clrecs,vis): 00361 """ 00362 Update history table when setSolarObjectJy is run 00363 """ 00364 # 00365 from taskinit import gentools 00366 (mytb,) = gentools(['tb']) 00367 mytb.open(vis+'/HISTORY',nomodify=False) 00368 nrow = mytb.nrows() 00369 lasttime=mytb.getcol('TIME')[nrow-1] 00370 rown=nrow 00371 # 00372 for ky in clrecs.keys(): 00373 # expect only one component for each 00374 comp = clrecs[ky]['component0'] 00375 srcn = ky.split('_')[0] 00376 ispw = ky.split('_')[1] 00377 mytb.addrows(1) 00378 appl='setjy' 00379 msg = (" %s: spw %s Flux:[I=%s,Q=%s,U=%s,V=%s] +/- [I=%s,Q=%s,U=%s,V=%s] Jy" % 00380 (srcn, ispw, comp['flux']['value'][0], comp['flux']['value'][1], 00381 comp['flux']['value'][2], comp['flux']['value'][3], 00382 comp['flux']['error'][0], comp['flux']['error'][1],comp['flux']['error'][2], 00383 comp['flux']['error'][3])) 00384 00385 origin='setjy' 00386 priority='INFO' 00387 time=lasttime 00388 emptystrarr=numpy.array(['']) 00389 mytb.putcell('APP_PARAMS', rown, ['']) 00390 mytb.putcell('CLI_COMMAND',rown, ['']) 00391 mytb.putcell('APPLICATION',rown, appl) 00392 mytb.putcell('MESSAGE',rown, msg) 00393 mytb.putcell('OBSERVATION_ID',rown, -1) 00394 mytb.putcell('ORIGIN',rown,origin) 00395 mytb.putcell('PRIORITY',rown,priority) 00396 mytb.putcell('TIME',rown, time) 00397 rown += 1 00398 mytb.close() 00399 00400 def testerrs(errcode,srcname): 00401 """ 00402 Check error codes from solar_system_fd 00403 return = 0 all ok 00404 return = 1 partly ok 00405 return = 2 all bad - should not proceed to set component 00406 """ 00407 from taskinit import casalog 00408 errcount = 0 00409 if type(errcode)!=list: 00410 errcode=[errcode] 00411 for ec in errcode: 00412 if ec != 0: 00413 errcount += 1 00414 if ec == 1: 00415 casalog.post("The model for %s is not supported" % srcname, 'WARN') 00416 elif ec == 2: 00417 casalog.post("Unsupported freuquency range",'WARN') 00418 elif ec == 3: 00419 casalog.post("Tb model not found",'WARN') 00420 elif ec == 4: 00421 casalog.post("The ephemeris table is not found or the time is out of range",'WARN') 00422 if errcount == len(errcode): 00423 return 2 00424 if errcount != 0: 00425 return 1 00426 else: 00427 return 0
1.8.0
