Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/setjy_helper.py: 4%

1# setjy helper functions 

2import os 

3import sys 

4import shutil 

5import numpy 

6 

7from casatasks import casalog as default_casalog 

8from casatools import quanta, ms, table, componentlist, measures, calibrater, msmetadata 

9from collections import OrderedDict as odict 

10 

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 

18 

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  

30 

31 qa = quanta() 

32 myms = ms( ) 

33 mytb = table( ) 

34 mycl = componentlist( ) 

35 myme = measures( ) 

36 mycb = calibrater( ) 

37 mymsmd = msmetadata( ) 

38 

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 

44 

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 

52 

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

62 

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

70 

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 

101 

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 

106 

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

115 

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 

131 

132 

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] 

140 

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 " 

151 

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) 

162 

163 # get a correct direction measure 

164 mymsmd.open(self.vis) 

165 dirdic = mymsmd.phasecenter(int(fid)) 

166 mymsmd.close() 

167 

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') 

174 

175 fielddirs[fid]=(numpy.array([[dirdic['m0']['value']],[dirdic['m1']['value']]])) 

176 

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) 

183 

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] 

194 

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))]) 

212 

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 

219 

220 inparams[fid]['freqlist']=freqlist 

221 inparams[fid]['freqlistraw']=freqlistraw 

222 inparams[fid]['framelist']=framelist 

223 inparams[fid]['reffreqs']=reffreqs 

224 myms.close() 

225 

226 

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 

235 

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']) 

249 

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") 

260 

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') 

281 

282 # packed fluxes for all spws  

283 fluxes.append(subfluxes) 

284 #casalog.post( "fluxes=" + fluxes) 

285 # fluxes has fluxesi[nspw][nf][nt] 

286 

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'] 

301 

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'] 

307 

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 

317 

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): " 

330 

331 self._casalog.post(dirmsg+str(dirstring)) 

332 

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? 

338 

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') 

357 

358 if (os.path.exists(clname)): 

359 shutil.rmtree(clname) 

360 

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')) 

380 

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] 

383 

384 # for a single CL with multple spws 

385 if j ==0: infreq0=infreq 

386 

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

412 

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') 

417 

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) 

423 

424 if scalebychan: 

425 # use tabular form to set flux for each channel 

426 

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 

437 

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]) 

457 

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

463 

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

473 

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

485 

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

489 

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) 

524 

525 

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

533 

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) 

538 

539 #if cleanupcomps:  

540 # shutil.rmtree(clname) 

541 #self.clnamelist.append(clname) 

542 self.clnamelist.append(clname) 

543 

544 msg="Using channel dependent " if scalebychan else "Using spw dependent " 

545 

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 

563 

564 def getclnamelist(self): 

565 return self.clnamelist 

566 

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') 

587 

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') 

598 

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'] 

613 

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") 

617 

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') 

655 

656 def _updateHistory(self,clrecs,vis): 

657 """ 

658 Update history table when setSolarObjectJy is run 

659 """ 

660 #  

661 mytb = table( ) 

662 

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])) 

680 

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

695 

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

703 

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] 

710 

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 

739 

740 else: 

741 ispw = complab.split('_')[1] 

742 mytb.addrows(1)