Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/setjy_helper.py: 77%
444 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1# setjy helper functions
2import os
3import sys
4import shutil
5import numpy
7from casatasks import casalog as default_casalog
8from casatools import quanta, ms, table, componentlist, measures, calibrater, msmetadata
9from collections import OrderedDict as odict
11class ss_setjy_helper:
12 def __init__(self,imtool, vis, casalog=None):
13 self.im = imtool
14 self.vis = vis
15 if not casalog:
16 casalog = default_casalog
17 self._casalog = casalog
19 def setSolarObjectJy(self,field,spw,scalebychan, timerange,observation, scan, intent, useephemdir, usescratch=False):
20 """
21 Set flux density of a solar system object using Bryan Butler's new
22 python model calculation code.
23 A single time stamp (first time stamp of MS after selections are applied) is
24 currently used per execution. For flux observation done in a long time span
25 may need to run multiple setjy with selections by time range (or scans).
26 """
27 #retval = True
28 output = {}
29 cleanupcomps = True # leave generated cl files
31 qa = quanta()
32 myms = ms( )
33 mytb = table( )
34 mycl = componentlist( )
35 myme = measures( )
36 mycb = calibrater( )
37 mymsmd = msmetadata( )
39 # prepare parameters need to pass to the Bryan's code
40 # make ms selections
41 # get source name
42 # get time ranges
43 # get spwused and get frequency ranges
45 sel={}
46 sel['field']=field
47 sel['spw']=spw
48 sel['timerange']=timerange
49 sel['observation']=str(observation)
50 sel['scan']=scan
51 sel['scanintent']=intent
53 measframes=['REST','LSRK','LSRD','BARY','GEO','TOPO','GALACTO','LGROUP','CMB']
54 myms.open(self.vis)
55 myms.msselect(sel,False)
56 scansummary=myms.getscansummary()
57 nscan=len(scansummary.keys())
58 selectedids=myms.msselectedindices()
59 fieldids=selectedids['field']
60 obsids=selectedids['observationid']
61 myms.close()
63 mytb.open(os.path.join(self.vis,'OBSERVATION'))
64 if len(obsids)==0:
65 getrow=0
66 else:
67 getrow=obsids[0]
68 observatory=mytb.getcell('TELESCOPE_NAME',getrow)
69 mytb.close()
71 mytb.open(os.path.join(self.vis,'FIELD'))
72 colnames = mytb.colnames()
73 if len(fieldids)==0:
74 fieldids = range(mytb.nrows())
75 # frame reference for field position
76 phdir_info=mytb.getcolkeyword("PHASE_DIR","MEASINFO")
77 if 'Ref' in phdir_info:
78 fieldref=phdir_info['Ref']
79 elif 'TabRefTypes' in phdir_info:
80 colnames=mytb.colnames()
81 for col in colnames:
82 if col=='PhaseDir_Ref':
83 fieldrefind=mytb.getcell(col,fieldids[0])
84 fieldref=phdir_info['TabRefTypes'][fieldrefind]
85 else:
86 fieldref=''
87 srcnames={}
88 fielddirs={}
89 ftimes={}
90 ephemid={}
91 for fid in fieldids:
92 #srcnames.append(mytb.getcell('NAME',int(fid)))
93 srcnames[fid]=(mytb.getcell('NAME',int(fid)))
94 #fielddirs[fid]=(mytb.getcell('PHASE_DIR',int(fid)))
95 ftimes[fid]=(mytb.getcell('TIME',int(fid)))
96 if colnames.count('EPHEMERIS_ID'):
97 ephemid[fid]=(mytb.getcell('EPHEMERIS_ID',int(fid)))
98 else:
99 ephemid[fid]=-1
100 mytb.close() # close FIELD table
102 # need to get a list of time
103 # but for now for test just get a time centroid of the all scans
104 # get time range
105 # Also, this needs to be done per source
107 # gather freq info from Spw table
108 mytb.open(os.path.join(self.vis,'SPECTRAL_WINDOW'))
109 nspw = mytb.nrows()
110 freqcol=mytb.getvarcol('CHAN_FREQ')
111 freqwcol=mytb.getvarcol('CHAN_WIDTH')
112 fmeasrefcol=mytb.getcol('MEAS_FREQ_REF')
113 reffreqs=mytb.getcol('REF_FREQUENCY')
114 mytb.close()
116 # store all parameters need to call solar_system_fd for all sources
117 inparams={}
118 validfids = [] # keep track of valid fid that has data (after selection)
119 # if same source name ....need handle this...
120 for fid in fieldids:
121 sel['field']=str(fid)
122 myms.open(self.vis)
123 #reset the selection
124 try:
125 status=myms.msselect(sel)
126 except:
127 #skip this field
128 continue
129 if not status:
130 continue
133 validfids.append(fid)
134 trange=myms.range('time')
135 #if not srcnames[fid] in inparams:
136 # inparams[srcnames[fid]]={}
137 if not fid in inparams:
138 inparams[fid]={}
139 inparams[fid]['fieldname']=srcnames[fid]
141 #tc = (trange['time'][0]+trange['time'][1])/2. #in sec.
142 # use first timestamp to be consistent with the ALMA Control
143 # old setjy (Butler-JPL-Horizons 2010) seems to be using
144 # time in FIELD... but here is first selected time in main table
145 if ephemid[fid]==-1: # keep old behavior
146 tc = trange['time'][0] #in sec.
147 tmsg = "Time used for the model calculation (=first time stamp of the selected data) for field "
148 else: # must have ephem table attached...
149 tc = ftimes[fid]#in sec.
150 tmsg = "Time used for the model calculation for field "
152 #if 'mjd' in inparams[srcnames[fid]]:
153 # inparams[srcnames[fid]]['mjds'][0].append([myme.epoch('utc',qa.quantity(tc,'s'))['m0']['value']])
154 #else:
155 # inparams[srcnames[fid]]['mjds']=[myme.epoch('utc',qa.quantity(tc,'s'))['m0']['value']]
156 if 'mjd' in inparams[fid]:
157 inparams[fid]['mjds'][0].append([myme.epoch('utc',qa.quantity(tc,'s'))['m0']['value']])
158 else:
159 inparams[fid]['mjds']=[myme.epoch('utc',qa.quantity(tc,'s'))['m0']['value']]
160 usedtime = qa.time(qa.quantity(tc,'s'),form='ymd')[0]
161 self._casalog.post(tmsg+str(fid)+":"+usedtime)
163 # get a correct direction measure
164 mymsmd.open(self.vis)
165 dirdic = mymsmd.phasecenter(int(fid))
166 mymsmd.close()
168 # check if frame is J2000 and if not do the transform
169 if dirdic['refer'] != 'J2000':
170 myme.doframe(myme.observatory(observatory))
171 myme.doframe(myme.epoch('utc', qa.quantity(tc,'s')))
172 azeldir = myme.measure(dirdic,'AZELGEO')
173 dirdic=myme.measure(azeldir,'J2000')
175 fielddirs[fid]=(numpy.array([[dirdic['m0']['value']],[dirdic['m1']['value']]]))
177 # somehow it gives you duplicated ids .... so need to uniquify
178 selspws= list(set(myms.msselectedindices()['spw']))
179 # make sure it is int rather than numpy.int32, etc.
180 selspws = [int(ispw) for ispw in selspws]
181 #inparams[srcnames[fid]]['spwids']= selspws if len(selspws)!=0 else range(nspw)
182 inparams[fid]['spwids']= selspws if len(selspws)!=0 else range(nspw)
184 #create a list of freq ranges with selected spws
185 # should worry about freq order???
186 freqlist=[]
187 freqlistraw=[]
188 framelist=[]
189 #for idx in inparams[srcnames[fid]]['spwids']:
190 for idx in inparams[fid]['spwids']:
191 freqs=freqcol['r'+str(idx+1)]
192 freqws=freqwcol['r'+str(idx+1)]
193 fmeasref=fmeasrefcol[idx]
195 #if scalebychan=T, this has to be min, max determined from
196 # chan_freq(channel center)+/- chan width.
197 if scalebychan:
198 # pack into list of list of list (freqlist[nf][nspw])
199 perspwfreqlist=[]
200 for nf in range(len(freqs)):
201 fl = freqs[nf][0]-freqws[nf][0]/2.
202 fh = freqs[nf][0]+freqws[nf][0]/2.
203 perspwfreqlist.append([fl,fh])
204 freqlist.append(perspwfreqlist)
205 else:
206 if (len(freqs)==1):
207 fl = freqs[0][0]-freqws[0][0]/2.
208 fh = freqs[0][0]+freqws[0][0]/2.
209 freqlist.append([float(fl),float(fh)])
210 else:
211 freqlist.append([float(min(freqs)),float(max(freqs))])
213 framelist.append(measframes[fmeasref])
214 freqlistraw.append(freqs.transpose()[0].tolist()) # data chanfreq list
215 #inparams[srcnames[fid]]['freqlist']=freqlist
216 #inparams[srcnames[fid]]['freqlistraw']=freqlistraw
217 #inparams[srcnames[fid]]['framelist']=framelist
218 #inparams[srcnames[fid]]['reffreqs']=reffreqs
220 inparams[fid]['freqlist']=freqlist
221 inparams[fid]['freqlistraw']=freqlistraw
222 inparams[fid]['framelist']=framelist
223 inparams[fid]['reffreqs']=reffreqs
224 myms.close()
227 # call Bryan's code
228 # errcode: list of list - inner list - each values for range of freqs
229 # flluxes: list of list
230 # fluxerrs:
231 # size: [majoraxis, minoraxis, pa]
232 # direction: direction for each time stamp
233 #
234 from . import solar_system_setjy as SSsetjy
236 #import solar_system_setjy2 as SSsetjy
237 retdict={} # for returning flux densities?
238 ss_setjy=SSsetjy.solar_system_setjy()
239 #for src in srcnames:
240 self.clnamelist=[]
241 for vfid in validfids:
242 src=srcnames[vfid]
243 #mjds=inparams[src]['mjds']
244 mjds=inparams[vfid]['mjds']
245 fluxes=[]
246 # call solar_system_fd() per spw (for scalebychan freqlist has an extra dimention)
247 #nspwused=len(inparams[src]['freqlist'])
248 nspwused=len(inparams[vfid]['freqlist'])
250 # warning for many channels but it is really depends on the source
251 if scalebychan:
252 maxnf=0
253 for ispw in range(nspwused):
254 #nf = len(inparams[src]['freqlist'][ispw])
255 nf = len(inparams[vfid]['freqlist'][ispw])
256 maxnf = max(nf,maxnf)
257 if maxnf >= 3840 and src.upper()!="MARS": # mars shoulde be ok
258 self._casalog.post("Processing %s spw(s) and at least some of them are a few 1000 channels or more. This may take \
259 many minutes (>3min per spw for 3840 channels) in some cases. Please be patient." % nspwused,"WARN")
261 for i in range(nspwused): # corresponds to n spw
262 if type(freqlist[0][0])==list:
263 #infreqs=inparams[src]['freqlist'][i]
264 infreqs=inparams[vfid]['freqlist'][i]
265 else:
266 #infreqs=[inparams[src]['freqlist'][i]]
267 infreqs=[inparams[vfid]['freqlist'][i]]
268 self._casalog.post("Calling solar_system_fd: %s(%s) for spw%s freqs=%s" % (src, vfid,i,freqlist[i]),'DEBUG1')
269 # casalog.post( "calling solar system fd...")
270 (errcodes, subfluxes, fluxerrs, sizes, dirs)=\
271 ss_setjy.solar_system_fd(source_name=src, MJDs=mjds, frequencies=infreqs, observatory=observatory, casalog=self._casalog)
272 # for old code
273 #(errcodes, subfluxes, fluxerrs, sizes, dirs)=\
274 # ss_setjy.solar_system_fd(source_name=src, MJDs=mjds, frequencies=infreqs, casalog=self._casalog)
275 # for nf freq ranges, nt mjds
276 # errcodes[nf][nt], subfluxes[nf][nt], fluxerrs[nf][nt], sizes[nt], dirs[nt]
277 self._casalog.post("+++++ solar_system_fd() returned values +++++", 'DEBUG1')
278 self._casalog.post(" fluxes(fds)=%s" % subfluxes, 'DEBUG1')
279 self._casalog.post(" sizes=%s" % sizes, 'DEBUG1')
280 self._casalog.post(" directions=%s\n" % dirs, 'DEBUG1')
282 # packed fluxes for all spws
283 fluxes.append(subfluxes)
284 #casalog.post( "fluxes=" + fluxes)
285 # fluxes has fluxesi[nspw][nf][nt]
287 # ------------------------------------------------------------------------
288 # For testing with hardcoded values without calling solar_system_fd()...
289 #errcodes=[[0,0,0,0,0]]
290 #fluxes=[[26.40653147,65.23839313,65.23382757,65.80638802,69.33396562]]
291 #fluxerrs=[[0.0,0.0,0.0,0.0,0.0]]
292 #sizes=[[3.6228991032674371,3.6228991032674371,0.0]]
293 #dirs=[{'m0': {'unit': 'rad', 'value': 0.0}, 'm1': {'unit': 'rad', 'value': 0.0}, 'refer': 'J2000', 'type': 'direction'}]
294 # ------------------------------------------------------------------------
295 # local params for selected src
296 #framelist=inparams[src]['framelist']
297 #freqlist=inparams[src]['freqlist']
298 #freqlistraw=inparams[src]['freqlistraw']
299 #reffreqs=inparams[src]['reffreqs']
300 #spwids=inparams[src]['spwids']
302 framelist=inparams[vfid]['framelist']
303 freqlist=inparams[vfid]['freqlist']
304 freqlistraw=inparams[vfid]['freqlistraw']
305 reffreqs=inparams[vfid]['reffreqs']
306 spwids=inparams[vfid]['spwids']
308 clrecs=odict( )
309 # casalog.post("LOOP over multiple dir...")
310 labels = []
311 # loop for over for multiple directions (=multiple MJDs) for a given src
312 for i in range(len(dirs)): # this is currently only length of 1 since no multiple timestamps were used
313 # check errcode - error code would be there per flux per time (dir)
314 reterr=testerrs(errcodes[i],src)
315 if reterr == 2:
316 continue
318 #dirstring = [dirs[i]['refer'],qa.tos(dirs[i]['m0']),qa.tos(dirs[i]['m1'])]
319 if useephemdir:
320 dirstring = [dirs[i]['refer'],qa.tos(dirs[i]['m0']),qa.tos(dirs[i]['m1'])]
321 dirmsg = "Using the source position from the ephemeris table in the casa data directory: "
322 else:
323 #dirstring = [fieldref, str(fielddirs[fieldids[0]][0][0])+'rad', str(fielddirs[fieldids[0]][1][0])+'rad']
324 # extract field direction of first id of the selected field ids
325 #dirstring = [fieldref, "%.18frad" % (fielddirs[fieldids[0]][0][0]), "%.18frad" % (fielddirs[fieldids[0]][1][0])]
326 #dirstring = [fieldref, "%.18frad" % (fielddirs[vfid][0][0]), "%.18frad" % (fielddirs[vfid][1][0])]
327 # by now the direction should be in J2000
328 dirstring = ['J2000', "%.18frad" % (fielddirs[vfid][0][0]), "%.18frad" % (fielddirs[vfid][1][0])]
329 dirmsg = "Using the source position in the MS (or in the attached ephemeris table, if any): "
331 self._casalog.post(dirmsg+str(dirstring))
333 # casalog.post("componentlist construction")
334 # setup componentlists
335 # need to set per dir
336 # if scalebychan=F, len(freqs) corresponds to nspw selected
337 # Need to put in for-loop to create cl for each time stamp? or scan?
339 #clpath='/tmp/'
340 clpath='./'
341 # construct cl name (multiple spws in one componentlist table)
342 selreffreqs = [reffreqs[indx] for indx in spwids]
343 minfreq = min(selreffreqs)
344 maxfreq = max(selreffreqs)
345 freqlabel = '%.3f-%.3fGHz' % (minfreq/1.e9, maxfreq/1.e9)
346 #tmlabel = '%.1fd' % (tc/86400.)
347 mjd=inparams[vfid]['mjds'][0]
348 tmlabel = '%.1fd' % (mjd)
349 #clabel = src+'_spw'+str(spwids[j])+'_'+freqlabel+'_'+tmlabel
350 #clabel0 = src+'_'+freqlabel+'_'+tmlabel
351 #pid=str(os.getpid())
352 #clname = clpath+clabel0+"_"+pid+'.cl'
353 clabel0 = src+'_'+freqlabel+'_'+tmlabel+"_"+os.path.basename(self.vis)
354 clname = clpath+clabel0+'.cl'
355 #debug
356 self._casalog.post("Create componentlist: %s for vfid=%s" % (clname,vfid),'DEBUG1')
358 if (os.path.exists(clname)):
359 shutil.rmtree(clname)
361 # casalog.post("Logging/output dict")
362 iiflux=[]
363 iinfreq=[]
364 # for logging/output dict - pack each freq list for each spw
365 freqsforlog=[]
366 for j in range(len(freqlist)): # loop over nspw
367 freqlabel = '%.3fGHz' % (reffreqs[int(spwids[j])]/1.e9)
368 clabel=src+'_spw'+str(spwids[j])+'_'+freqlabel+'_'+tmlabel
369 # casalog.post("addcomponent...for spw=",spwids[j]," i=",i," flux=",fluxes[j][i])
370 if scalebychan:
371 index= 2.0
372 sptype = 'spectral index'
373 else:
374 index= 0.0
375 sptype = 'constant'
376 # adjust to change in returned flux shape. An extra [] now seems to be gone. 2012-09-27
377 iflux=fluxes[j][i][0]
378 # casalog.post( "addcomponent with flux=%s at frequency=%s" %
379 # (iflux,str(reffreqs[int(spwids[j])]/1.e9)+'GHz'))
381 # i - time stamps = 0 for now, j = a freq range
382 infreq=freqlist[j][0][0] if type(freqlist[j][0])==list else freqlist[j][0]
384 # for a single CL with multple spws
385 if j ==0: infreq0=infreq
387 # --- version for 'multi-spws in a single comp single CL'
388 # accumulate all frequencies from all spws to a list
389 #if not scalebychan:
390 #freqlist[j] for j spw should contain a list with two freqs (min and max of
391 # of the spw and the same fluxes should be set for the freqs so that cl.setspectrum
392 # with tabular would not interpolate inside each spw
393 ##iinfreq.extend(freqlist[j])
394 # nch = len(freqlistraw[j])
395 # iinfreq.extend(freqlistraw[j])
396 # assigning the same flux to the two freqs
397 ##iiflux.extend([iflux,iflux])
398 # for ich in range(nch):
399 # iiflux.extend([iflux])
400 #else:
401 # iiflux.extend(fluxes[j][i])
402 # if type(freqlist[j][0])==list and len(freqlist[j][0])>1:
403 # for fr in freqlist[j]:
404 # iinfreq.append((fr[1]+fr[0])/2)
405 # else:
406 # if type(freqlist[j])==list:
407 # iinfreq.append((freqlist[j][1]+freqlist[j][0])/2)
408 # else:
409 # iinfreq.append(freqlist[j])
410 # freqsforlog.append(iinfreq)
411 # -----------------------------------------------------------------------------------
413 # casalog.post("Logging/output dict addcomp")
414 self._casalog.post("addcomponent with flux=%s at frequency=%s" %\
415 # (fluxes[j][i][0],str(reffreqs[int(spwids[j])]/1.e9)+'GHz'), 'INFO1')
416 (iflux,str(reffreqs[int(spwids[j])]/1.e9)+'GHz'), 'INFO1')
418 #
419 mycl.addcomponent(iflux,fluxunit='Jy', polarization="Stokes", dir=dirstring,
420 shape='disk', majoraxis=str(sizes[i][0])+'arcsec', minoraxis=str(sizes[i][1])+'arcsec',
421 positionangle=str(sizes[i][2])+'deg', freq=[framelist[0],str(infreq)+'Hz'],
422 spectrumtype=sptype, index=index, label=clabel)
424 if scalebychan:
425 # use tabular form to set flux for each channel
427 # This may be redundant
428 if type(fluxes[j][i]) ==list and len(fluxes[j][i])> 1:
429 if type(freqlist[j][0])==list and len(freqlist[j][0])>1:
430 freqs=[]
431 for fr in freqlist[j]:
432 freqs.append((fr[1]+fr[0])/2)
433 clind = mycl.length() - 1
434 mycl.setspectrum(which=clind, type='tabular', tabularfreqs=freqs, tabularflux=fluxes[j][0],
435 tabularframe=framelist[j])
436 ### === per spw process end here
438 # - version that put all spws into a component -
439 # set a single component freq set at the the lower edge of first spw
440 #mycl.addcomponent(iiflux[0],fluxunit='Jy', polarization="Stokes", dir=dirstring,
441 # shape='disk', majoraxis=str(sizes[i][0])+'arcsec', minoraxis=str(sizes[i][1])+'arcsec',
442 # positionangle=str(sizes[i][2])+'deg', freq=[framelist[0],str(infreq0)+'Hz'],
443 # spectrumtype=sptype, index=index, label=clabel)
444 # if it's list of fluxes try to put in tabular form
445 #if type(fluxes[j][i]) ==list and len(fluxes[j][i])> 1:
446 # if len(iiflux)> 1:
447 # casalog.post("framelist[j]=",framelist[j])
448 #if type(freqlist[j][0])==list and len(freqlist[j][0])>1:
449 # freqs=[]
450 # for fr in freqlist[j]:
451 # freqs.append((fr[1]+fr[0])/2)
452 #else:
453 # freqs=freqlist[j]
454 #clind = mycl.length() - 1
455 # mycl.setspectrum(which=clind, type='tabular', tabularfreqs=freqs, tabularflux=fluxes[j][0],
456 #tabularframe=framelist[j])
458 # mycl.setspectrum(which=clind, type='tabular', tabularfreqs=iinfreq, tabularflux=iiflux,
459 # tabularframe=framelist[0])
460 #mycl.rename(clname)
461 #put in a record for log output
462 #clrecs[clabel] = mycl.torecord()
464 clrecs[clabel0] = mycl.torecord()
465 clrecs[clabel0]['allfluxinfo']=fluxes
466 clrecs[clabel0]['allfreqinfo']=freqsforlog
467 clrecs[clabel0]['spwids']=spwids
468 mycl.rename(clname) # - A CL per field
469 # casalog.post("clrecs=" + clrecs)
470 # casalog.post("clname=" + clname)
471 mycl.close(False) # False for not to send a warning message
472 mycl.done()
474 # if scratch=F check if the virtual model already exist
475 # for the field if it is clear it.
476 # if j==0 and (not usescratch):
477 # mytb.open(self.vis, nomodify=False)
478 # kwds = mytb.getkeywords()
479 # modelkwd='definedmodel_field_'+str(vfid)
480 # if modelkwd in kwds:
481 # clmodname=kwds[modelkwd]
482 # mytb.removekeyword(clmodname)
483 # mytb.removekeyword(modelkwd)
484 # mytb.close()
486 mycb.open(self.vis,addcorr=False,addmodel=False)
487 mycb.delmod(otf=True,field=str(vfid),spw=list(spwids), scr=False)
488 mycb.close()
490 # finally, put the componentlist as model
491 #tqlstr=''
492 #if intent!='':
493 # tqlstr='any(STATE_ID=='+str(stateids.tolist())+')'
494 # Currently still need to do per spw (see the comment below...)
495 # assuming CL contains nrow = nspw components
496 mycl.open(clname)
497 clinrec = mycl.torecord()
498 mycl.close(False)
499 mycl.done()
500 ncomp = clinrec['nelements']
501 if ncomp != len(spwids):
502 raise Exception("Inconsistency in generated componentlist...Please submit a bug report.")
503 for icomp in range(ncomp):
504 #self.im.selectvis(spw=spwids[icomp],field=field,observation=observation,time=timerange,intent=intent)
505 ok = self.im.selectvis(spw=spwids[icomp],field=vfid,observation=observation,time=timerange,intent=intent)
506 # for MMS, particular subms may not contain the particular spw, so skip for such a case
507 if ok:
508 newclinrec = {}
509 newclinrec['nelements']=1
510 newclinrec['component0']=clinrec['component'+str(icomp)]
511 mycl.fromrecord(newclinrec)
512 cdir = os.getcwd()
513 if os.access(cdir,os.W_OK):
514 tmpclpath = cdir+"/"
515 elif os.access("/tmp",os.W_OK):
516 tmpclpath = "/tmp/"
517 #tmpclpath=tmpclpath+"_tmp_setjyCLfile_"+pid
518 tmpclpath=tmpclpath+os.path.basename(self.vis)+"_tmp_setjyCLfile"
519 if os.path.exists(tmpclpath):
520 shutil.rmtree(tmpclpath)
521 mycl.rename(tmpclpath)
522 mycl.close(False)
523 self.im.ft(complist=tmpclpath, phasecentertime=tc)
526 # --------------------------------------------------------------------------------------------
527 # Currently this does not work properly for some case. Need ft can handle multi-component CL
528 # with the way to map which component apply to which spw
529 # Unable when such functionality is implemented in im.ft()
530 #self.im.selectvis(spw=spwids,field=field,observation=observation,time=timerange,intent=intent)
531 #self.im.ft(complist=clname) <= e.g. im.ft(comprec) where comrec = {'spw-mapping':[...],comprecs}..
532 # --------------------------------------------------------------------------------------------
534 #debug: set locally saved 2010-version component list instead
535 #cl2010='mod_setjy_spw0_Titan_230.543GHz55674.1d.cl'
536 # casalog.post("Using complist=" + cl2010)
537 #self.im.ft(complist=cl2010)
539 #if cleanupcomps:
540 # shutil.rmtree(clname)
541 #self.clnamelist.append(clname)
542 self.clnamelist.append(clname)
544 msg="Using channel dependent " if scalebychan else "Using spw dependent "
546 if clrecs=={}:
547 raise Exception("No componentlist is generated.")
548 #self._casalog.post(clrecs)
549 self._casalog.post(msg+" flux densities")
550 #self._reportoLog(clrecs,self._casalog)
551 self._reportoLog2(clrecs,self._casalog)
552 #self._updateHistory(clrecs,self.vis)
553 self._updateHistory2(clrecs,self.vis)
554 # dictionary of dictionary for each field
555 retdict[vfid]=clrecs
556 # ==== end of for loop over fields
557 #output=self._makeRetFluxDict(retdict)
558 # casalog.post("makeRetFluxDict2")
559 output=self._makeRetFluxDict2(retdict)
560 #return retval
561 # casalog.post("done setSolarSetJY")
562 return output
564 def getclnamelist(self):
565 return self.clnamelist
567 def _reportoLog(self,clrecs,casalog):
568 """
569 send model parameters to log
570 """
571 # casalog.post("clrecs=" + clrecs)
572 for ky in clrecs.keys():
573 # expect only one component for each
574 nelm = clrecs[ky]['nelements']
575 for icomp in range(nelm):
576 comp = clrecs[ky]['component'+str(icomp)]
577 complab = comp['label']
578 #srcn = ky.split('_')[0]
579 #ispw = ky.split('_')[1]
580 srcn = complab.split('_')[0]
581 ispw = complab.split('_')[1]
582 if icomp==0:
583 msg=" direction set in the componentlist: RA=%s rad, Dec%s rad" %\
584 (float('%.6g' % comp['shape']['direction']['m0']['value']),
585 float('%.6g' % comp['shape']['direction']['m1']['value']))
586 casalog.post(msg,'INFO2')
588 msg=" %s: %s Flux:[I=%s,Q=%s,U=%s,V=%s] +/- [I=%s,Q=%s,U=%s,V=%s] Jy" %\
589 (srcn, ispw, float('%.5g' % comp['flux']['value'][0]),
590 float('%.5g' % comp['flux']['value'][1]),
591 float('%.5g' % comp['flux']['value'][2]),
592 float('%.5g' % comp['flux']['value'][3]),
593 float('%.5g' % comp['flux']['error'][0]),
594 float('%.5g' % comp['flux']['error'][1]),
595 float('%.5g' % comp['flux']['error'][2]),
596 float('%.5g' % comp['flux']['error'][3]))
597 casalog.post(msg, 'INFO')
599 def _reportoLog2(self,clrecs,casalog):
600 """
601 send model parameters to log
602 (This version handles for a single componentlist with possible multiple componet)
603 Assume componentlist has a name src_spw_xxxx or src_spw
604 """
605 nelm = -1
606 ncl = len(clrecs.keys())
607 # casalog.post("clrecs=" + clrecs)
608 if ncl == 1:
609 clname = list(clrecs.keys())[0]
610 comprec = clrecs[clname]
611 if 'nelements' in comprec:
612 nelm = comprec['nelements']
614 if ncl != 1 or nelm == -1:
615 casalog.post("Cannot process the generated componentlist, possibly a bug.\
616 Please submit a bug report", "SEVERE")
618 for icomp in range(nelm):
619 comp = comprec['component'+str(icomp)]
620 complab = comp['label']
621 srcn = complab.split('_')[0]
622 if icomp==0:
623 msg=" direction set in the componentlist: RA=%s rad, Dec%s rad" %\
624 (float('%.6g' % comp['shape']['direction']['m0']['value']),
625 float('%.6g' % comp['shape']['direction']['m1']['value']))
626 casalog.post(msg,'INFO2')
627 freq0=''
628 if 'spectrum' in comp:
629 freqv = float('%.5g' % comp['spectrum']['frequency']['m0']['value'])
630 freq0 = str(freqv)+comp['spectrum']['frequency']['m0']['unit']
631 if nelm==1:
632 for i in range(len(comprec['spwids'])):
633 ispw=comprec['spwids'][i]
634 iflux = comprec['allfluxinfo'][i][0][0]
635 # currently only I flux is reported and no flux error is reported
636 msg=" %s: spw%s Flux:[I=%s,Q=%s,U=%s,V=%s] +/- [I=%s,Q=%s,U=%s,V=%s] Jy @ %s" %\
637 (srcn, ispw, float('%.5g' % iflux),0.0,0.0,0.0,
638 0.0,0.0,0.0,0.0, freq0)
639 casalog.post(msg, 'INFO')
640 else:
641 ispw = complab.split('_')[1]
642 #if 'spectrum'in comp:
643 # freqv = float('%.5g' % comp['spectrum']['frequency']['m0']['value'])
644 # freq0 = str(freqv)+comp['spectrum']['frequency']['m0']['unit']
645 msg=" %s: %s Flux:[I=%s,Q=%s,U=%s,V=%s] +/- [I=%s,Q=%s,U=%s,V=%s] Jy @ %s" %\
646 (srcn, ispw, float('%.5g' % comp['flux']['value'][0]),
647 float('%.5g' % comp['flux']['value'][1]),
648 float('%.5g' % comp['flux']['value'][2]),
649 float('%.5g' % comp['flux']['value'][3]),
650 float('%.5g' % comp['flux']['error'][0]),
651 float('%.5g' % comp['flux']['error'][1]),
652 float('%.5g' % comp['flux']['error'][2]),
653 float('%.5g' % comp['flux']['error'][3]), freq0)
654 casalog.post(msg, 'INFO')
656 def _updateHistory(self,clrecs,vis):
657 """
658 Update history table when setSolarObjectJy is run
659 """
660 #
661 mytb = table( )
663 mytb.open(os.path.join(vis,'HISTORY'),nomodify=False)
664 nrow = mytb.nrows()
665 lasttime=mytb.getcol('TIME')[nrow-1]
666 rown=nrow
667 #
668 for ky in clrecs.keys():
669 # expect only one component for each
670 comp = clrecs[ky]['component0']
671 srcn = ky.split('_')[0]
672 ispw = ky.split('_')[1]
673 mytb.addrows(1)
674 appl='setjy'
675 msg = (" %s: spw %s Flux:[I=%s,Q=%s,U=%s,V=%s] +/- [I=%s,Q=%s,U=%s,V=%s] Jy" %
676 (srcn, ispw, comp['flux']['value'][0], comp['flux']['value'][1],
677 comp['flux']['value'][2], comp['flux']['value'][3],
678 comp['flux']['error'][0], comp['flux']['error'][1],comp['flux']['error'][2],
679 comp['flux']['error'][3]))
681 origin='setjy'
682 priority='INFO'
683 time=lasttime
684 emptystrarr=numpy.array([''])
685 mytb.putcell('APP_PARAMS', rown, [''])
686 mytb.putcell('CLI_COMMAND',rown, [''])
687 mytb.putcell('APPLICATION',rown, appl)
688 mytb.putcell('MESSAGE',rown, msg)
689 mytb.putcell('OBSERVATION_ID',rown, -1)
690 mytb.putcell('ORIGIN',rown,origin)
691 mytb.putcell('PRIORITY',rown,priority)
692 mytb.putcell('TIME',rown, time)
693 rown += 1
694 mytb.close()
696 def _updateHistory2(self,clrecs,vis):
697 """
698 Update history table when setSolarObjectJy is run
699 (mulitple comopents in one version)
700 """
701 #
702 mytb = table( )
704 mytb.open(os.path.join(vis,'HISTORY'),nomodify=False)
705 nrow = mytb.nrows()
706 lasttime=mytb.getcol('TIME')[nrow-1]
707 rown=nrow
708 #
709 clname = list(clrecs.keys())[0]
711 #for ky in clrecs.keys():
712 srcn = clname.split('_')[0]
713 comprec = clrecs[clname]
714 nelm = comprec['nelements']
715 for icomp in range(nelm):
716 comp = comprec['component'+str(icomp)]
717 complab = comp['label']
718 origin='setjy'
719 priority='INFO'
720 time=lasttime
721 appl='setjy'
722 emptystrarr=numpy.array([''])
723 if nelm==1:
724 for i in range(len(comprec['spwids'])):
725 mytb.addrows(1)
726 ispw=comprec['spwids'][i]
727 iflux = comprec['allfluxinfo'][i][0][0]
728 msg = (" %s: spw %s Flux:[I=%s,Q=%s,U=%s,V=%s] +/- [I=%s,Q=%s,U=%s,V=%s] Jy" %
729 (srcn, ispw, iflux, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
730 mytb.putcell('APP_PARAMS', rown, [''])
731 mytb.putcell('CLI_COMMAND',rown, [''])
732 mytb.putcell('APPLICATION',rown, appl)
733 mytb.putcell('MESSAGE',rown, msg)
734 mytb.putcell('OBSERVATION_ID',rown, -1)
735 mytb.putcell('ORIGIN',rown,origin)
736 mytb.putcell('PRIORITY',rown,priority)
737 mytb.putcell('TIME',rown, time)
738 rown += 1
740 else:
741 ispw = complab.split('_')[1]
742 mytb.addrows(1)
743 msg = (" %s: spw %s Flux:[I=%s,Q=%s,U=%s,V=%s] +/- [I=%s,Q=%s,U=%s,V=%s] Jy" %
744 (srcn, ispw, comp['flux']['value'][0], comp['flux']['value'][1],
745 comp['flux']['value'][2], comp['flux']['value'][3],
746 comp['flux']['error'][0], comp['flux']['error'][1],comp['flux']['error'][2],
747 comp['flux']['error'][3]))
749 mytb.putcell('APP_PARAMS', rown, [''])
750 mytb.putcell('CLI_COMMAND',rown, [''])
751 mytb.putcell('APPLICATION',rown, appl)
752 mytb.putcell('MESSAGE',rown, msg)
753 mytb.putcell('OBSERVATION_ID',rown, -1)
754 mytb.putcell('ORIGIN',rown,origin)
755 mytb.putcell('PRIORITY',rown,priority)
756 mytb.putcell('TIME',rown, time)
757 rown += 1
758 mytb.close()
760 def _makeRetFluxDict(self, flxdict):
761 """
762 re-arrange the calculated flux density info for returned dict.
763 """
764 # flxdict should contains a ditionary per field
765 retflxdict={}
766 for fid in flxdict.keys():
767 comp=flxdict[fid]
768 tmpdict={}
769 for ky in comp.keys():
770 srcn = ky.split('_')[0]
771 ispw = ky.split('_')[1].strip('spw')
773 tmpdict[ispw]={}
774 tmpdict[ispw]['fluxd']=comp[ky]['component0']['flux']['value']
775 #tmpdict[ispw]['fluxderr']=comp[ky]['component0']['flux']['error']
776 tmpdict['fieldName']=srcn
777 retflxdict[str(fid)]=tmpdict
778 retflxdict['format']=\
779 '{field Id: {spw Id: {fluxd:[I,Q,U,V] in %s}, \'fieldName\':field name}}' % \
780 comp[ky]['component0']['flux']['unit']
781 return retflxdict
783 def _makeRetFluxDict2(self, flxdict):
784 """
785 re-arrange the calculated flux density info for returned dict.
786 (for each field id as key, single dict)
787 """
788 # flxdict should contains a ditionary per field
789 retflxdict={}
790 for fid in flxdict.keys():
791 clrecs=flxdict[fid]
792 clname=list(clrecs.keys())[0]
793 srcn = clname.split('_')[0]
794 comprec=clrecs[clname]
795 nelm = comprec['nelements']
796 tmpdict={}
797 for icomp in range(nelm):
798 comp = comprec['component'+str(icomp)]
799 complab = comp['label']
800 if nelm==1:
801 for i in range(len(comprec['spwids'])):
802 ispw = str(comprec['spwids'][i])
803 iflux = comprec['allfluxinfo'][i][0][0]
804 tmpdict[ispw]={}
805 tmpdict[ispw]['fluxd']=[iflux,0.0,0.0,0.0]
806 else:
807 ispw = complab.split('_')[1].strip('spw')
808 tmpdict[ispw]={}
809 #tmpdict[ispw]['fluxd']=comp[ky]['component0']['flux']['value']
810 tmpdict[ispw]['fluxd']=comp['flux']['value']
811 tmpdict['fieldName']=srcn
812 retflxdict[str(fid)]=tmpdict
813 retflxdict['format']=\
814 '{field Id: {spw Id: {fluxd:[I,Q,U,V] in %s}, \'fieldName\':field name}}' % \
815 comprec['component0']['flux']['unit']
816 return retflxdict
819def testerrs(errcode,srcname):
820 """
821 Check error codes from solar_system_fd
822 return = 0 all ok
823 return = 1 partly ok
824 return = 2 all bad - should not proceed to set component
825 """
826 errcount = 0
827 if type(errcode)!=list:
828 errcode=[errcode]
829 for ec in errcode:
830 if ec != 0:
831 errcount += 1
832 if ec == 1:
833 default_casalog.post("The model for %s is not supported" % srcname, 'WARN')
834 elif ec == 2:
835 default_casalog.post("Unsupported frequency range",'WARN')
836 elif ec == 3:
837 default_casalog.post("Tb model not found",'WARN')
838 elif ec == 4:
839 default_casalog.post("The ephemeris table is not found or the time is out of range",'WARN')
840 if errcount == len(errcode):
841 return 2
842 if errcount != 0:
843 return 1
844 else:
845 return 0