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

835 statements  

« prev     ^ index     » next       coverage.py v7.6.4, created at 2024-11-01 07:19 +0000

1import os 

2import re 

3import pylab as pl 

4 

5from casatools import ctsys, quanta, imager 

6from casatasks import casalog 

7from .simutil import * 

8from .simutil import is_array_type 

9 

10qa = quanta() 

11 

12def simobserve( 

13 project=None, 

14 skymodel=None, inbright=None, indirection=None, incell=None, 

15 incenter=None, inwidth=None, # innchan=None, 

16 complist=None, compwidth=None, comp_nchan=1, 

17 setpointings=None, 

18 ptgfile=None, integration=None, direction=None, mapsize=None, 

19 maptype=None, pointingspacing=None, caldirection=None, calflux=None, 

20 # observe=None, 

21 obsmode=None, 

22 refdate=None, hourangle=None, 

23 totaltime=None, antennalist=None, 

24 sdantlist=None, 

25 sdant=None, 

26 outframe=None, 

27 thermalnoise=None, 

28 user_pwv=None, t_ground=None, t_sky=None, tau0=None, seed=None, 

29 leakage=None, 

30 graphics=None, 

31 verbose=None, 

32 overwrite=None 

33 ): 

34 

35 # Collect a list of parameter values to save inputs 

36 in_params = locals() 

37 

38 try: 

39 

40 ######################### 

41 # some hardcoded variables 

42 pbcoeff = 1.13 ## PB defined as pbcoeff*lambda/d 

43 nyquist = 0.5/pbcoeff ## Nyquist spacing = PB*nyquist 

44 

45 gridratio_int = 1./pl.sqrt(3) # times lambda/d 

46 gridratio_tp = 1./3 

47 

48 relmargin = .5 # number of PB between edge of model and ptg centers 

49 scanlength = 1 # number of integrations per scan 

50 

51 

52 # RI TODO for inbright=unchanged, need to scale input image to jy/pix 

53 # according to actual units in the input image 

54 

55 casalog.origin('simobserve') 

56 if verbose: casalog.filter(level="DEBUG2") 

57 

58 # create the utility object 

59 # this is the dir of the observation (could be "") 

60 util = simutil(direction) 

61 if verbose: util.verbose = True 

62 msg = util.msg 

63 

64 # it was requested to make the user interface "observe" for what  

65 # is sm.observe and sm.predict. 

66 # interally the code is clearer if we stick with predict so 

67 predict = obsmode.startswith('i') or obsmode.startswith('s') 

68 if predict: 

69 uvmode = obsmode.startswith('i') 

70 if not uvmode: antennalist = sdantlist 

71 elif sdantlist != "": 

72 if antennalist == "": 

73 uvmode = False 

74 antennalist = sdantlist 

75 else: 

76 #uvmode = True 

77 #msg("Both antennalist and sdantlist are defined. sdantlist will be ignored",priority="warn") 

78 emsg = "Both antennalist and sdantlist are defined. Define one of them." 

79 raise ValueError(emsg) 

80 else: 

81 uvmode = True 

82 # uvmode = (sdant < 0) #when flexible default values come available 

83 

84 

85 # put output in directory called "project" 

86 fileroot = project 

87 if not os.path.exists(fileroot): 

88 os.mkdir(fileroot) 

89 

90 

91 # filename parsing of cfg file here so that the project filenames  

92 # can contain the cfg 

93 repodir = ctsys.resolve("alma/simmos") 

94 

95 

96 # convert "alma;0.4arcsec" to an actual configuration 

97 # can only be done after reading skymodel, so here, we just string parse 

98 if str.upper(antennalist[0:4]) == "ALMA": 

99 foo=antennalist.replace(';','_') 

100 elif len(antennalist) > 0: 

101 foo=antennalist 

102 else: 

103 msg("The name of antenna list (antennalist/sdantlist) is not specified",priority="error") 

104 

105 if foo: 

106 foo=foo.replace(".cfg","") 

107 sfoo=foo.split('/') 

108 if len(sfoo)>1: 

109 foo=sfoo[-1] 

110 project=project+"."+foo 

111 

112 

113 if not overwrite: 

114 if (predict and uvmode and os.path.exists(fileroot+"/"+ 

115 project+".ms")): 

116 emsg = fileroot+"/"+project+".ms exists but overwrite=F" 

117 raise RuntimeError(emsg) 

118 if (predict and (not uvmode) and os.path.exists(fileroot+"/"+ 

119 project+".sd.ms")): 

120 emsg = fileroot+"/"+project+".sd.ms exists but overwrite=F" 

121 raise RuntimeError(emsg) 

122 

123 

124 if is_array_type(skymodel): 

125 skymodel = skymodel[0] 

126 skymodel = skymodel.replace('$project',project) 

127 

128 if is_array_type(complist): 

129 complist = complist[0] 

130 

131 if((not os.path.exists(skymodel)) and (not os.path.exists(complist))): 

132 if len(skymodel)>0: 

133 msg("Your skymodel '"+skymodel+"' could not be found.", 

134 priority="warn") 

135 if len(complist)>0: 

136 msg("Your complist '"+complist+"' could not be found.", 

137 priority="warn") 

138 if len(skymodel)==0 and len(complist)==0: 

139 msg("At least one of skymodel or complist must be set.", 

140 priority="error") 

141 

142 else: 

143 msg("No sky input found."+ 

144 " At least one of skymodel or complist must exist.", 

145 priority="error") 

146 

147 ### WORKAROUND for wrong flux in COMP TP simulations (CAS-5095) 

148 if (obsmode.startswith("s") and os.path.exists(complist)): 

149 emsg = "Single dish simulation has a flux recovery issue when using a components list as an input.\nPlease generate compskymodel image first by obsmode='' and use the image as the skymodel input.\nSorry for the inconvenience." 

150 raise RuntimeError(emsg) 

151 ### End of WORKAROUND 

152 

153 grscreen = False 

154 grfile = False 

155 if graphics == "both": 

156 grscreen = True 

157 grfile = True 

158 if graphics == "screen": 

159 grscreen = True 

160 if graphics == "file": 

161 grfile = True 

162 

163 ################################################################## 

164 # set up skymodel image 

165 

166 if os.path.exists(skymodel): 

167 components_only = False 

168 # create a new skymodel called skymodel,  

169 # or if it's already there, called newmodel 

170 default_model = project + ".skymodel" 

171 if skymodel == default_model: 

172 newmodel = fileroot + "/" + project + ".newmodel" 

173 else: 

174 newmodel = fileroot + "/" + default_model 

175 if os.path.exists(newmodel): 

176 if overwrite: 

177 shutil.rmtree(newmodel) 

178 else: 

179 emsg = (newmodel+" exists -- "+ 

180 "please delete it, change skymodel, or set overwrite=T") 

181 raise RuntimeError(emsg) 

182 

183 # modifymodel just collects info if skymodel==newmodel 

184 innchan = -1 

185 returnpars = util.modifymodel(skymodel,newmodel, 

186 inbright,indirection,incell, 

187 incenter,inwidth,innchan, 

188 flatimage=False) 

189 if not returnpars: 

190 raise RuntimeError('Failure in modifymodel. It returned: {}'. 

191 format(returnpars)) 

192 

193 (model_refdir,model_cell,model_size, 

194 model_nchan,model_specrefval,model_specrefpix,model_width, 

195 model_stokes) = returnpars 

196 

197 modelflat = fileroot + "/" + project + ".skymodel.flat" 

198 if os.path.exists(modelflat) and (not predict): 

199 # if we're not predicting, then we want to use the previously 

200 # created modelflat, because it may have components added  

201 msg("flat sky model "+modelflat+ 

202 " exists, predict not requested",priority="warn") 

203 msg(" working from existing model image - "+ 

204 "please delete it if you wish to overwrite.", 

205 priority="warn") 

206 else: 

207 # create and add components into modelflat with util.flatimage() 

208 util.flatimage(newmodel,complist=complist,verbose=verbose) 

209 # we want skymodel.flat image to be called that no matter what  

210 # the skymodel image is called, since it's used in analysis 

211 if modelflat != newmodel+".flat": 

212 if os.path.exists(modelflat): 

213 shutil.rmtree(modelflat) 

214 shutil.move(newmodel+".flat",modelflat) 

215 

216 casalog.origin('simobserve') 

217 

218 # set startfeq and bandwidth in util object after modifymodel 

219 bandwidth = qa.mul(qa.quantity(model_nchan), 

220 qa.quantity(model_width)) 

221 util.bandwidth = bandwidth 

222 

223 else: 

224 components_only = True 

225 msg("component-only simulation",priority="info") 

226 # calculate model parameters from the component list: 

227 

228 compdirs = [] 

229 cl.done() 

230 cl.open(complist) 

231 

232 for i in range(cl.length()): 

233 compdirs.append(util.dir_m2s(cl.getrefdir(i))) 

234 

235 model_refdir, coffs = util.average_direction(compdirs) 

236 model_specrefval = cl.getspectrum(0)['frequency']['m0'] 

237 

238 if util.isquantity(compwidth,halt=False): 

239 model_width = compwidth 

240 msg("compwidth set: setting model bandwidth to input", 

241 priority="info") 

242 else: 

243 model_width = "2GHz" 

244 msg("compwidth unset: setting bandwidth to 2GHz", 

245 priority="warn") 

246 

247 model_nchan = comp_nchan 

248 # channelize component-only MS  

249 # currently assuming equal width and center as frequency reference 

250 model_specrefpix = 0.5*(comp_nchan-1) 

251 msg("scaling model bandwidth by model_nchan",priority="info") 

252 model_width = qa.div(model_width,model_nchan) 

253 model_stokes = "I" 

254 

255 cmax = 0.0014 # ~5 arcsec 

256 for i in range(coffs.shape[1]): 

257 xc = pl.absolute(coffs[0,i]) # offsets in deg 

258 yc = pl.absolute(coffs[1,i]) 

259 if xc > cmax: 

260 cmax = xc 

261 if yc > cmax: 

262 cmax = yc 

263 

264 model_size = ["%fdeg" % (3*cmax), "%fdeg" % (3*cmax)] 

265 

266 

267 # for cases either if there is a skymodel or are only components, 

268 # if user has not input a map size (for setpointings), use model_size 

269 if len(mapsize) == 0: 

270 mapsize = model_size 

271 if verbose: msg("setting map size to "+str(model_size), 

272 origin='simobserve') 

273 else: 

274 if is_array_type(mapsize): 

275 if len(mapsize[0]) == 0: 

276 mapsize = model_size 

277 if verbose: msg("setting map size to "+str(model_size), 

278 origin="simobserve") 

279 

280 if components_only: 

281 if is_array_type(mapsize): 

282 map_asec = qa.convert(mapsize[0],"arcsec")['value'] 

283 else: 

284 map_asec = qa.convert(mapsize,"arcsec")['value'] 

285 if is_array_type(model_size): 

286 mod_asec = qa.convert(model_size[0],"arcsec")['value'] 

287 else: 

288 mod_asec = qa.convert(model_size,"arcsec")['value'] 

289 if map_asec>mod_asec: 

290 model_size=["%farcsec" % map_asec,"%farcsec" % map_asec] 

291 

292 

293 

294 

295 ################################################################## 

296 # read antenna file here to get Primary Beam, and estimate psfsize 

297 aveant = -1 

298 stnx = [] # for later, to know if we read an array in or not 

299 pb = 0. # primary beam 

300 

301 # convert "alma;0.4arcsec" to an actual configuration 

302 if str.upper(antennalist[0:4]) == "ALMA": 

303 

304 resparsed=False 

305 # test for cycle 1 

306 q = re.compile('.*CYCLE.?1.?;(.*)') 

307 qq = q.match(antennalist.upper()) 

308 if qq: 

309 z = qq.groups() 

310 tail=z[0] 

311 tail=tail.lower() 

312 if util.isquantity(tail,halt=False): 

313 resl = qa.convert(tail,"arcsec")['value'] 

314 if os.path.exists(repodir): 

315 confnum = (1.044 - 6.733 * pl.log10(resl * qa.convert(model_specrefval,"GHz")['value'] / 345.)) 

316 confnum = max(1,min(6,confnum)) 

317 conf = str(int(round(confnum))) 

318 antennalist = os.path.join(repodir,"alma.cycle1." + conf + ".cfg") 

319 msg("converted resolution to antennalist "+antennalist) 

320 resparsed=True 

321 else: 

322 msg("failed to find antenna configuration repository"+ 

323 " at "+repodir,priority="error") 

324 if not resparsed: 

325 q = re.compile('.*CYCLE.?2.?;(.*)') 

326 qq = q.match(antennalist.upper()) 

327 if qq: 

328 z = qq.groups() 

329 tail=z[0] 

330 tail=tail.lower() 

331 if util.isquantity(tail,halt=False): 

332 resl = qa.convert(tail,"arcsec")['value'] 

333 if os.path.exists(repodir): 

334 confnum = 10. ** (0.91 - 0.74 * (resl * qa.convert(model_specrefval,"GHz")['value']/345.)) 

335 confnum = max(1,min(7,confnum)) 

336 conf = str(int(round(confnum))) 

337 antennalist = os.path.join(repodir,"alma.cycle2." + 

338 conf + ".cfg") 

339 msg("converted resolution to antennalist "+ 

340 antennalist) 

341 resparsed=True 

342 else: 

343 msg("failed to find antenna configuration repository at "+repodir,priority="error") 

344 

345 if not resparsed: # assume FS 

346 tail = antennalist[5:] 

347 if util.isquantity(tail,halt=False): 

348 resl = qa.convert(tail,"arcsec")['value'] 

349 if os.path.exists(repodir): 

350 confnum = (2.867-pl.log10(resl*1000*qa.convert(model_specrefval,"GHz")['value']/672.))/0.0721 

351 confnum = max(1,min(28,confnum)) 

352 conf = str(int(round(confnum))) 

353 if len(conf) < 2: conf = '0' + conf 

354 antennalist = os.path.join(repodir,"alma.out" 

355 + conf + ".cfg") 

356 msg("converted resolution to antennalist "+antennalist) 

357 else: 

358 msg("failed to find antenna configuration repository at "+repodir,priority="error") 

359 

360 

361 # Search order is fileroot/ -> specified path -> repository 

362 if len(antennalist) > 0: 

363 if os.path.exists(fileroot+"/"+antennalist): 

364 antennalist = fileroot + "/" + antennalist 

365 elif not os.path.exists(antennalist) and \ 

366 os.path.exists(os.path.join(repodir,antennalist)): 

367 antennalist = os.path.join(repodir,antennalist) 

368 # Now make sure the antennalist exists 

369 if not os.path.exists(antennalist): 

370 emsg = "Couldn't find antennalist: %s" % antennalist 

371 raise RuntimeError(emsg) 

372 elif predict or components_only: 

373 # antennalist is required when predicting or components only 

374 util.msg("Must specify antennalist", priority="error") 

375 

376 # Read antennalist 

377 if os.path.exists(antennalist): 

378 stnx, stny, stnz, stnd, padnames, antnames, telescopename, posobs = util.readantenna(antennalist) 

379 nant=len(stnx) 

380 if nant == 1: 

381 if predict and uvmode: 

382 # observe="int" but antennalist is SD 

383 util.msg("antennalist contains only 1 antenna", 

384 priority="error") 

385 uvmode = False 

386 if not uvmode: #Single-dish 

387 # KS TODO: what if not predicting but SD with multi-Ants 

388 # in antennalist (e.g., aca.tp)? In that case, PB on plots and 

389 # pointingspacing="??PB" will not be correct for heterogeneous list. 

390 if sdant > nant-1 or sdant < -nant: 

391 msg("antenna index %d is out of range. setting sdant=0"%sdant,priority="warn") 

392 sdant = 0 

393 stnx = [stnx[sdant]] 

394 stny = [stny[sdant]] 

395 stnz = [stnz[sdant]] 

396 stnd = pl.array(stnd[sdant]) 

397 padnames = [padnames[sdant]] 

398 antnames = [antnames[sdant]] 

399 nant = 1 

400 

401 # (set back to simdata - there must be an automatic way to do this) 

402 casalog.origin('simobserve') 

403 

404 aveant = stnd.mean() 

405 # TODO use max ant = min PB instead? 

406 pb = pbcoeff*0.29979/qa.convert(qa.quantity(model_specrefval),'GHz')['value']/aveant*3600.*180/pl.pi # arcsec 

407 

408 # PSF size 

409 if uvmode: 

410 # approx beam, to compare with model_cell: 

411 psfsize = util.approxBeam(antennalist,qa.convert(qa.quantity(model_specrefval),'GHz')['value']) 

412 else: # Single-dish 

413 psfsize = pb 

414 # check for model size 

415 if not components_only: 

416 minsize = min(qa.convert(model_size[0],'arcsec')['value'], 

417 qa.convert(model_size[1],'arcsec')['value']) 

418 if minsize < (2.5 * pb): 

419 msg("skymodel should be larger than 2.5*primary beam. Your skymodel: %.3f arcsec < %.3f arcsec: 2.5*primary beam" % (minsize, 2.5*pb),priority="error") 

420 del minsize 

421 else: 

422 emsg = "Can't find antennalist" 

423 raise RuntimeError(emsg) 

424 

425 

426 # now we have an estimate of the psf from the antenna configuration,  

427 # so we can guess a model_cell for the case of component-only  

428 # simulation,  

429 if components_only: 

430 # first set based on psfsize: 

431 # needs high subsampling because small shifts in placement of  

432 # components lead to large changes in the difference image. 

433 model_cell = [ qa.quantity(str(psfsize/20)+"arcsec"), 

434 qa.quantity(str(psfsize/20)+"arcsec") ] 

435 

436 # XXX if the user has set direction should we center the compskymodel there? 

437 # if len(direction)>0: model_refdir = direction 

438 

439 # and can create a compskymodel image (tmp) and  

440 # skymodel.flat which is what is needed for analysis. 

441 

442 if components_only: 

443 newmodel = fileroot + "/" + project + ".compskymodel" 

444 needmodel=True 

445 

446 modimsize=int((qa.convert(model_size[0],"arcsec")['value'])/(qa.convert(model_cell[0],"arcsec")['value'])) 

447# modimsize=max([modimsize,32]) 

448 newepoch,newlat,newlon = util.direction_splitter(model_refdir) 

449 

450 if os.path.exists(newmodel): 

451 if overwrite: 

452 shutil.rmtree(newmodel) 

453 else: 

454 needmodel=False 

455 ia.open(newmodel) 

456 oldshape=ia.shape() 

457 if len(oldshape) != 2: 

458 needmodel=True 

459 else: 

460 if oldshape[0] != modimsize or oldshape[1]==modimsize: 

461 needmodel=True 

462 oldcs=ia.coordsys() 

463 ia.done() 

464 olddir = (oldcs.referencevalue())['numeric'] 

465 if ( olddir[0] != qa.convert(newlat,oldcs.units()[0])['value'] or 

466 olddir[1] != qa.convert(newlon,oldcs.units()[1])['value'] or 

467 newepoch != oldcs.referencecode() ): 

468 needmodel=True 

469 oldcs.done() 

470 del oldcs, olddir 

471 if needmodel: 

472 emsg = newmodel+" exists and is inconsistent with required size="+str(modimsize)+" and direction. Please set overwrite=True" 

473 raise RuntimeError(emsg) 

474 

475 if needmodel: 

476 csmodel = ia.newimagefromshape(newmodel,[modimsize,modimsize,1,1]) 

477 modelcsys = csmodel.coordsys() 

478 modelshape = csmodel.shape() 

479 cell0_rad=qa.convert(model_cell[0],'rad')['value'] 

480 cell1_rad=qa.convert(model_cell[1],'rad')['value'] 

481 modelcsys.setdirection(refpix=[modimsize/2,modimsize/2], 

482 refval=[qa.tos(newlat),qa.tos(newlon)], 

483 refcode=newepoch) 

484 modelcsys.setincrement([-cell0_rad,cell1_rad],'direction') 

485 modelcsys.setreferencevalue(type="spectral",value=qa.convert(model_specrefval,'Hz')['value']) 

486 modelcsys.setreferencepixel(type="spectral",value=model_specrefpix) 

487 modelcsys.setrestfrequency(model_specrefval) 

488 modelcsys.setincrement(type="spectral",value=compwidth) 

489 csmodel.setcoordsys(modelcsys.torecord()) 

490 

491 modelcsys.done() 

492 cl.open(complist) 

493 csmodel.setbrightnessunit("Jy/pixel") 

494 csmodel.modify(cl.torecord(),subtract=False) 

495 cl.done() 

496 csmodel.done() 

497 # as noted, compskymodel doesn't need to exist, only skymodel.flat 

498 # flatimage adds in components if complist!=None 

499 #util.flatimage(newmodel,complist=complist,verbose=verbose) 

500 util.flatimage(newmodel,verbose=verbose) 

501 modelflat = fileroot + "/" + project + ".compskymodel.flat" 

502# modelflat = fileroot + "/" + project + ".skymodel.flat" 

503# if modelflat != newmodel+".flat": 

504# if os.path.exists(modelflat): 

505# shutil.rmtree(modelflat) 

506# shutil.move(newmodel+".flat",modelflat) 

507 

508 

509 # and finally, with model_cell set either from an actual skymodel, 

510 # or from the antenna configuration in components_only case, 

511 # we can check for the user that the psf is likely to be sampled enough: 

512 cell_asec=qa.convert(model_cell[0],'arcsec')['value'] 

513 if psfsize < cell_asec: 

514 emsg = "Sky model cell of "+str(cell_asec)+" asec is very large compared to highest resolution "+str(psfsize)+" asec - this will lead to blank or erroneous output. (Did you set incell?)" 

515 shutil.rmtree(modelflat) 

516 raise RuntimeError(emsg) 

517 if psfsize < 2*cell_asec: 

518 msg("Sky model cell of "+str(cell_asec)+" asec is large compared to highest resolution "+str(psfsize)+" asec. (Did you set incell?)",priority="warn") 

519 

520 # set this for future minimum image size 

521 minimsize = 8* int(psfsize/cell_asec) 

522 

523 

524 

525 ################################################################## 

526 # set up pointings 

527 dir = model_refdir 

528 dir0 = dir 

529 if is_array_type(direction): 

530 if len(direction) > 0: 

531 if util.isdirection(direction[0],halt=False): 

532 dir = direction 

533 dir0 = direction[0] 

534 else: 

535 if util.isdirection(direction,halt=False): 

536 dir = direction 

537 dir0 = dir 

538 util.direction = dir0 

539 

540 # if the integration time is a real time quantity 

541 # test for weird units 

542 if not util.isquantity(integration): 

543 emsg = "integration time "+integration+" does not appear to represent a time interval (use 's','min'; not 'sec','m')" 

544 raise RuntimeError(emsg) 

545 

546 if qa.quantity(integration)['unit'] != '': 

547 if not qa.compare(integration,"1s"): 

548 emsg = "integration time "+integration+" does not appear to represent a time interval ('s','min'; not 'sec','m')" 

549 raise RuntimeError(emsg) 

550 intsec = qa.convert(qa.quantity(integration),"s")['value'] 

551 else: 

552 if len(integration)>0: 

553 intsec = float(integration) 

554 msg("interpreting integration time parameter as "+ 

555 str(intsec)+"s",priority="warn") 

556 else: 

557 intsec = 0 

558 integration="%fs" %intsec 

559 

560 

561 if setpointings: 

562 util.msg("calculating map pointings centered at "+ 

563 str(dir0),origin='simobserve') 

564 

565 if len(pointingspacing) < 1: 

566 if uvmode: 

567 # ALMA OT uses lambda/d/sqrt(3) 

568 pointingspacing = "%fPB" % (gridratio_int/pbcoeff) 

569 else: 

570 pointingspacing = "%fPB" % (gridratio_tp/pbcoeff) 

571 if str.upper(pointingspacing)=="NYQUIST": 

572 pointingspacing="%fPB" % nyquist 

573 q = re.compile('(\d+.?\d+)\s*PB') 

574 qq = q.match(pointingspacing.upper()) 

575 if qq: 

576 z = qq.groups() 

577 if pb <= 0: 

578 emsg = "Can't calculate pointingspacing in terms of primary beam because antennalist is not specified" 

579 raise RuntimeError(emsg) 

580 pointingspacing = "%farcsec" % (float(z[0])*pb) 

581 # todo make more robust to nonconforming z[0] strings 

582 

583 if verbose: 

584 msg("pointing spacing in mosaic = "+ 

585 pointingspacing,origin='simobserve') 

586 pointings = util.calc_pointings2(pointingspacing, 

587 mapsize,maptype=maptype, 

588 direction=dir, beam=pb) 

589 nfld=len(pointings) 

590 etime = qa.convert(qa.mul(qa.quantity(integration),scanlength),"s")['value'] 

591 # etime is an array of scan lengths - here they're all the same. 

592 etime = etime + pl.zeros(nfld) 

593 # totaltime might not allow all fields to be observed, or it might 

594 # repeat 

595 ptgfile = fileroot + "/" + project + ".ptg.txt" 

596 else: 

597 if is_array_type(ptgfile): 

598 ptgfile = ptgfile[0] 

599 ptgfile = ptgfile.replace('$project',project) 

600 # precedence to ptg file outside the project dir 

601 if os.path.exists(ptgfile): 

602 shutil.copyfile(ptgfile,fileroot+"/"+project + ".ptg.txt") 

603 ptgfile = fileroot + "/" + project + ".ptg.txt" 

604 else: 

605 if os.path.exists(fileroot+"/"+ptgfile): 

606 ptgfile = fileroot + "/" + ptgfile 

607 else: 

608 emsg = "Can't find pointing file "+ptgfile 

609 raise RuntimeError(emsg) 

610 

611 nfld, pointings, etime = util.read_pointings(ptgfile) 

612 if max(etime) <= 0: 

613 # integration is a string in input params 

614 etime = intsec 

615 # make etime into an array 

616 etime = etime + pl.zeros(nfld) 

617 # etimes determine stop/start i.e. the scan 

618 # if a longer etime is in the file, it'll do multiple integrations 

619 # per scan 

620 # expects that the cal is separate, and this is just one round of the mosaic 

621 # furthermore, the cal will use _integration_ from the inputs, and that 

622 # needs to be less than the min etime: 

623 if min(etime) < intsec: 

624 integration = str(min(etime))+"s" 

625 msg("Setting integration to "+integration+ 

626 " to match the shortest time in the pointing file.", 

627 priority="warn") 

628 intsec = min(etime) 

629 

630 

631 # find imcenter - phase center 

632 imcenter , offsets = util.median_direction(pointings) 

633 epoch, ra, dec = util.direction_splitter(imcenter) 

634 

635 # model is centered at model_refdir, and has model_size; 

636 mepoch, mra, mdec = util.direction_splitter(model_refdir) 

637 # ra/mra should be in degrees, but just in case 

638 ra=qa.convert(ra,'deg') 

639 mra=qa.convert(mra,'deg') 

640 dec=qa.convert(dec,'deg') 

641 mdec=qa.convert(mdec,'deg') 

642 

643 # observation near ra=0: 

644 if abs(mra['value'])<10 or abs(mra['value']-360.)<10 or abs(ra['value'])<10 or abs(mra['value']-360.)<10: 

645 nearzero=True 

646 else: 

647 nearzero=False 

648 

649 # fix wraps 

650 if nearzero: # put break at 180 

651 if ra['value'] >= 180.: 

652 ra['value'] = ra['value'] - 360. 

653 if mra['value'] >= 180.: 

654 mra['value'] = mra['value'] - 360. 

655 else: 

656 if ra['value'] >= 360.: 

657 ra['value'] = ra['value'] - 360. 

658 if mra['value'] >= 360.: 

659 mra['value'] = mra['value'] - 360. 

660 

661 # shift is the offset from the model to imcenter 

662 shift = [ (qa.convert(ra,'deg')['value'] - 

663 qa.convert(mra,'deg')['value'])*pl.cos(qa.convert(dec,'rad')['value'] ), 

664 (qa.convert(dec,'deg')['value'] - qa.convert(mdec,'deg')['value']) ] 

665 if verbose: 

666 msg("pointings are shifted relative to the model by %g,%g arcsec" % (shift[0]*3600,shift[1]*3600),origin='simobserve') 

667 

668 xmax = qa.convert(model_size[0],'deg')['value']*0.5 

669 ymax = qa.convert(model_size[1],'deg')['value']*0.5 

670 # add PB halfwidth (relmargin=0.5)  

671 # for mosaics of small model images 

672 xmax=xmax+pb*relmargin/3600 

673 ymax=ymax+pb*relmargin/3600 

674 overlap = False 

675 # wrapang in median_direction should make offsets always small, not >360 

676 for i in range(offsets.shape[1]): 

677 xc = pl.absolute(offsets[0,i]+shift[0]) # offsets and shift are in degrees 

678 yc = pl.absolute(offsets[1,i]+shift[1]) 

679 if xc < xmax and yc < ymax: 

680 overlap = True 

681 break 

682 

683 if setpointings: 

684 if os.path.exists(ptgfile): 

685 if overwrite: 

686 os.remove(ptgfile) 

687 else: 

688 emsg = "pointing file "+ptgfile+" already exists and user does not want to overwrite" 

689 raise RuntimeError(emsg) 

690 util.write_pointings(ptgfile,pointings,etime.tolist()) 

691 

692 msg("center = "+imcenter,origin='simobserve') 

693 if nfld > 1 and verbose: 

694 for idir in range(min(len(pointings),20)): 

695 msg(" "+pointings[idir],origin='simobserve') 

696 if nfld >= 20: 

697 msg(" (printing only first 20 - see pointing file for full list)") 

698 

699 if not overlap: 

700 emsg = "No overlap between model and pointings" 

701 raise RuntimeError(emsg) 

702 

703 

704 

705 ################################################################## 

706 # calibrator is not explicitly contained in the pointing file 

707 # but interleaved with etime=intergration 

708 util.isquantity(calflux) 

709 calfluxjy = qa.convert(calflux,'Jy')['value'] 

710 # XML returns a list even for a string: 

711 if is_array_type(caldirection): caldirection = caldirection[0] 

712 if len(caldirection) < 4: caldirection = "" 

713 if calfluxjy > 0 and caldirection != "" and uvmode: 

714 docalibrator = True 

715 if os.path.exists(fileroot+"/"+project+'.cal.cclist'): 

716 shutil.rmtree(fileroot+"/"+project+'.cal.cclist') 

717 util.isdirection(caldirection) 

718 cl.done() 

719 cl.addcomponent(flux=calfluxjy,dir=caldirection, 

720 label="phase calibrator") 

721 # set reference freq to center freq of model 

722 cl.rename(fileroot+"/"+project+'.cal.cclist') 

723 cl.done() 

724 else: 

725 docalibrator = False 

726 

727 

728 

729 

730 

731 ################################################################## 

732 # create one figure for model and pointings - need antenna diam  

733 # to determine primary beam 

734 if grfile: 

735 file = fileroot + "/" + project + ".skymodel.png" 

736 else: 

737 file = "" 

738 

739 if grscreen or grfile: 

740 util.newfig(show=grscreen) 

741 

742# remove after we fix the scaling algorithm for the images 

743 if components_only: 

744 pl.plot() 

745 # TODO add symbols at locations of components 

746 pl.plot(coffs[0,]*3600,coffs[1,]*3600,'o',c="#dddd66") 

747 pl.axis("equal") 

748 else: 

749 discard = util.statim(modelflat,plot=True,incell=model_cell) 

750 

751 lims = pl.xlim(),pl.ylim() 

752 if pb <= 0 and verbose: 

753 msg("unknown primary beam size for plot",priority="warn") 

754 if max(max(lims)) > pb and not components_only: 

755 plotcolor = 'w' 

756 else: 

757 plotcolor = 'k' 

758 

759 #if offsets.shape[1] > 16 or pb <= 0 or pb > pl.absolute(max(max(lims))): 

760 if offsets.shape[1] > 19 or pb <= 0: 

761 lims = pl.xlim(),pl.ylim() 

762 pl.plot((offsets[0]+shift[0])*3600., 

763 (offsets[1]+shift[1])*3600., 

764 plotcolor+'+',markeredgewidth=1) 

765 #if pb > 0 and pl.absolute(lims[0][0]) > pb: 

766 if pb > 0: 

767 plotpb(pb,pl.gca(),lims=lims,color=plotcolor) 

768 else: 

769 from matplotlib.patches import Circle 

770 for i in range(offsets.shape[1]): 

771 pl.gca().add_artist(Circle( 

772 ((offsets[0,i]+shift[0])*3600, 

773 (offsets[1,i]+shift[1])*3600), 

774 radius=pb/2.,edgecolor=plotcolor,fill=False, 

775 label='beam',transform=pl.gca().transData,clip_on=True)) 

776 

777 xlim = max(abs(pl.array(lims[0]))) 

778 ylim = max(abs(pl.array(lims[1]))) 

779 # show entire pb: (statim doesn't by default) 

780 pl.xlim([max([xlim,pb/2]),min([-xlim,-pb/2])]) 

781 pl.ylim([min([-ylim,-pb/2]),max([ylim,pb/2])]) 

782 pl.xlabel("resized model sky",fontsize="x-small") 

783 util.endfig(show=grscreen,filename=file) 

784 

785 

786 

787 

788 

789 

790 

791 

792 ################################################################## 

793 # set up observatory, feeds, etc 

794 quickpsf_current = False 

795 

796 if uvmode: 

797 msfile = fileroot + "/" + project + '.ms' 

798 else: 

799 msfile = fileroot + "/" + project + '.sd.ms' 

800 

801 if predict: 

802 # TODO check for frequency overlap here - if zero stop 

803 # position overlap already checked above in pointing section 

804 

805 message = "preparing empty measurement set" 

806 if verbose: 

807 msg(" ",priority="info") 

808 msg(message,origin="simobserve",priority="warn") 

809 else: 

810 msg(message,origin="simobserve") 

811 

812 nbands = 1; 

813 fband = util.bandname(qa.convert(model_specrefval, 'GHz')['value']) 

814 

815 ############################################ 

816 # predict observation 

817 

818 # if someone has the old style refdate with the included, discard 

819 q = re.compile('(\d*/\d+/\d+)([/:\d]*)') 

820 qq = q.match(refdate) 

821 if not qq: 

822 emsg = "Invalid reference date "+refdate 

823 raise RuntimeError(emsg) 

824 else: 

825 z = qq.groups() 

826 refdate=z[0] 

827 if len(z)>1: 

828 if len(z[1])>1: 

829 msg("Discarding time part of refdate, '"+z[1]+ 

830 "', in favor of hourangle parameter = "+hourangle, 

831 origin='simobserve') 

832 

833 if hourangle=="transit": 

834 haoffset=0.0 

835 else: 

836 haoffset="no" 

837 # is this a time quantity? 

838 if qa.isquantity(str(hourangle)+"h"): 

839 if qa.compare(str(hourangle)+"h","s"): 

840 haoffset=qa.convert(qa.quantity(str(hourangle)+ 

841 "h"),'s')['value'] 

842 if qa.isquantity(hourangle): 

843 qha=qa.convert(hourangle,"s") 

844 if qa.compare(qha,"s"): 

845 haoffset=qa.convert(qha,'s')['value'] 

846 if haoffset=="no": 

847 msg("Cannot interpret your hourangle parameter "+hourangle+ 

848 " as a time quantity e.g. '5h', 30min'", 

849 origin="simobserve",priority="error") 

850 else: 

851 msg("You desire an hour angle of "+ 

852 str(haoffset/3600.)+" hours",origin="simobserve") 

853 

854 refdate=refdate+"/00:00:00" 

855 usehourangle=True 

856 

857 

858 intsec = qa.convert(qa.quantity(integration),"s")['value'] 

859 

860 # totaltime as an integer for # times through the mosaic: 

861 if not util.isquantity(totaltime,halt=False): 

862 emsg = "totaltime "+totaltime+" does not appear to represent a time interval (use 's','min','h'; not 'sec','m','hr')" 

863 raise RuntimeError(emsg) 

864 

865 if qa.quantity(totaltime)['value'] < 0.: 

866 # casapy crashes for negative totaltime 

867 emsg = "Negative totaltime is not allowed." 

868 raise ValueError(emsg) 

869 if qa.quantity(totaltime)['unit'] == '': 

870 # assume it means number of maps, or # repetitions. 

871 totalsec = sum(etime) 

872 if docalibrator: 

873 totalsec = totalsec + intsec # cal gets one int-time 

874 totalsec = float(totaltime) * totalsec 

875 msg("Total observing time = "+str(totalsec)+"s.") 

876 else: 

877 if not qa.compare(totaltime,"1s"): 

878 emsg = "totaltime "+totaltime+" does not appear to represent a time interval (use 's','min','h'; not 'sec','m','hr')" 

879 raise ValueError(emsg) 

880 totalsec = qa.convert(qa.quantity(totaltime),'s')['value'] 

881 

882 if os.path.exists(msfile) and not overwrite: #redundant check? 

883 emsg = ("measurement set "+msfile+ 

884 " already exists and user does not wish to overwrite") 

885 raise RuntimeError(emsg) 

886 sm.open(msfile) 

887 

888 diam = stnd; 

889 # WARNING: sm.setspwindow is not consistent with clean::center 

890 # but the "start" is the center of the first channel: 

891 model_start = qa.sub(model_specrefval, 

892 qa.mul(model_width,model_specrefpix)) 

893 

894 mounttype = 'alt-az' 

895 if telescopename in ['DRAO', 'WSRT']: 

896 mounttype = 'EQUATORIAL' 

897 # Should ASKAP be BIZARRE or something else? It may be effectively equatorial. 

898 

899 sm.setconfig(telescopename=telescopename, x=stnx, y=stny, z=stnz, 

900 dishdiameter=diam.tolist(), 

901 mount=[mounttype], antname=antnames, padname=padnames, 

902 coordsystem='global', referencelocation=posobs) 

903 if str.upper(telescopename).find('VLA') >= 0: 

904 sm.setspwindow(spwname=fband, freq=qa.tos(model_start), 

905 deltafreq=qa.tos(model_width), 

906 freqresolution=qa.tos(model_width), 

907 nchannels=model_nchan, refcode=outframe, 

908 stokes='RR LL') 

909 sm.setfeed(mode='perfect R L',pol=['']) 

910 else: 

911 sm.setspwindow(spwname=fband, freq=qa.tos(model_start), 

912 deltafreq=qa.tos(model_width), 

913 freqresolution=qa.tos(model_width), 

914 nchannels=model_nchan, refcode=outframe, 

915 stokes='XX YY') 

916 sm.setfeed(mode='perfect X Y',pol=['']) 

917 

918 if verbose: 

919 msg(" spectral window set at %s" % qa.tos(model_specrefval), 

920 origin='simobserve') 

921 sm.setlimits(shadowlimit=0.01, elevationlimit='10deg') 

922 if uvmode: 

923 sm.setauto(0.0) 

924 else: #Single-dish 

925 # auto-correlation should be unity for single dish obs. 

926 sm.setauto(1.0) 

927 

928 mereftime = me.epoch('UTC', refdate) 

929 # integration is a scalar quantity, etime is a vector of seconds 

930 sm.settimes(integrationtime=integration, usehourangle=usehourangle, 

931 referencetime=mereftime) 

932 

933 for k in range(0,nfld): 

934 src = project + '_%d' % k 

935 sm.setfield(sourcename=src, sourcedirection=pointings[k], 

936 calcode="OBJ", distance='0m') 

937 if k == 0: 

938 sourcefieldlist = src 

939 else: 

940 sourcefieldlist = sourcefieldlist + ',' + src 

941 if docalibrator: 

942 sm.setfield(sourcename="phase calibrator", 

943 sourcedirection=caldirection,calcode='C', 

944 distance='0m') 

945 

946 # time required to observe all planned scanes in etime array: 

947 totalscansec = sum(etime) 

948 kfld = 0 

949 

950 if totalsec < totalscansec: 

951 msg("Not all pointings in the mosaic will be observed - check mosaic setup and exposure time parameters!",priority="warn") 

952 ### 

953 #casalog.post("you need at least %16.12e sec but you have %16.12e sec (%f < %f = %s)" % (totalscansec, totalsec, totalsec, totalscansec, str(totalsec<totalscansec))) 

954 ### 

955 

956 # sm.observemany 

957 srces = [] 

958 starttimes = [] 

959 stoptimes = [] 

960 dirs = [] 

961 

962 if usehourangle: 

963 sttime = -totalsec/2.0 

964 else: 

965 sttime = 0. # leave start at the reftime 

966 sttime=sttime+haoffset 

967 scanstart=sttime 

968 

969 # can before sources 

970 if docalibrator: 

971 endtime = sttime + qa.convert(integration,'s')['value'] 

972 sm.observe(sourcename="phase calibrator", spwname=fband, 

973 starttime=qa.quantity(sttime, "s"), 

974 stoptime=qa.quantity(endtime, "s"), 

975 state_obs_mode="CALIBRATE_PHASE.ON_SOURCE", 

976 state_sig=True,project=project); 

977 sttime = endtime 

978 

979 while (sttime-scanstart) < totalsec: # the last scan could exceed totaltime 

980 endtime = sttime + etime[kfld] 

981 src = project + '_%d' % kfld 

982 srces.append(src) 

983 starttimes.append(str(sttime)+"s") 

984 stoptimes.append(str(endtime)+"s") 

985 dirs.append(pointings[kfld]) 

986 kfld = kfld + 1 

987 # advance start time - XX someday slew goes here 

988 sttime = endtime 

989 

990 if kfld == nfld: 

991 if docalibrator: 

992 endtime = sttime + qa.convert(integration,'s')['value'] 

993 

994 # need to observe cal singly to get new row in obs table 

995 # so first observemany the on-source pointing(s) 

996 sm.observemany(sourcenames=srces, 

997 spwname=fband, 

998 starttimes=starttimes, 

999 stoptimes=stoptimes, 

1000 project=project) 

1001 # and clear the list 

1002 srces = [] 

1003 starttimes = [] 

1004 stoptimes = [] 

1005 dirs = [] 

1006 sm.observe(sourcename="phase calibrator", spwname=fband, 

1007 starttime=qa.quantity(sttime, "s"), 

1008 stoptime=qa.quantity(endtime, "s"), 

1009 state_obs_mode="CALIBRATE_PHASE.ON_SOURCE", 

1010 state_sig=True,project=project); 

1011 kfld = kfld + 1 

1012 sttime = endtime 

1013# if kfld > nfld: kfld = 0 

1014 if kfld > nfld-1: kfld = 0 

1015 # if directions is unset, NewMSSimulator::observemany 

1016 

1017 # looks up the direction in the field table. 

1018 if not docalibrator: 

1019 sm.observemany(sourcenames=srces, 

1020 spwname=fband, 

1021 starttimes=starttimes, 

1022 stoptimes=stoptimes, 

1023 project=project) 

1024 

1025 sm.setdata(fieldid=list(range(0,nfld))) 

1026 if uvmode or components_only: #Interferometer only 

1027 sm.setvp(dovp=True,usedefaultvp=False) 

1028 # only use mosaic gridding for Het arrays for now -  

1029 # the standard gridding with a VPSkyJones is less susceptible 

1030 # to issues if the image is too small which can happen a lot 

1031 # in Simulation. 

1032 if len(pl.unique(diam))>1: 

1033 if nfld>1: 

1034 sm.setoptions(ftmachine="mosaic") 

1035 else: 

1036 msg("Heterogeneous array only supported for mosaics (nfld>1), and make sure that your image is larger than the primary beam or results may be unstable",priority="error") 

1037 else: 

1038 # checks have to be manual since there's no way to  

1039 # get the "diam" out of PBMath AFAIK 

1040 if telescopename=="ALMA": 

1041 if (diam[0]<10)|(diam[0]>13): 

1042 msg("Diameter = %f is inconsistent with telescope=ALMA in the configuration file. *12m ALMA PB will be used*. Use observatory='ACA' in the configuration file to use a 7m diameter."%diam[0],priority="warn") 

1043 n=len(diam) 

1044 diam=12.+pl.zeros(n) 

1045 elif telescopename=="ACA": 

1046 if (diam[0]<6)|(diam[0]>7.5): 

1047 msg("Diameter = %f is inconsistent with telescope=ALMA in the configuration file. *7m ACA PB will be used*"%diam[0],priority="warn") 

1048 n=len(diam) 

1049 diam=7.+pl.zeros(n) 

1050 else: 

1051 msg("Note: diameters in configuration file will not be used - PB for "+telescopename+" will be used",priority="info") 

1052 

1053 

1054 msg("done setting up observations (blank visibilities)", 

1055 origin='simobserve') 

1056 if verbose: sm.summary() 

1057 

1058 # do actual calculation of visibilities: 

1059 

1060 if not uvmode: #Single-dish 

1061 sm.setoptions(gridfunction='pb',ftmachine="sd",location=posobs) 

1062 

1063 if not components_only: 

1064 if docalibrator: 

1065 if len(complist) <=0: 

1066 complist=fileroot+"/"+project+'.cal.cclist' 

1067 else: 

1068 # XXX will 2 cl work? 

1069 complist=complist+","+fileroot+"/"+project+'.cal.cclist' 

1070 

1071 if len(complist) > 1: 

1072 message = "predicting from "+newmodel+" and "+complist 

1073 if verbose: 

1074 msg(" ",priority="info") 

1075 msg(message,priority="warn",origin="simobserve") 

1076 else: 

1077 msg(message,origin="simobserve") 

1078 else: 

1079 message = "predicting from "+newmodel 

1080 if verbose: 

1081 msg(" ",priority="info") 

1082 msg(message,priority="warn",origin="simobserve") 

1083 else: 

1084 msg(message,origin="simobserve") 

1085 sm.predict(imagename=newmodel,complist=complist) 

1086 else: # if we're doing only components 

1087 # XXX will 2 cl work? 

1088 if docalibrator: 

1089 complist=complist+","+fileroot+"/"+project+'.cal.cclist' 

1090 if verbose: 

1091 msg("predicting from "+complist,priority="warn", 

1092 origin="simobserve") 

1093 else: 

1094 msg("predicting from "+complist,origin="simobserve") 

1095 sm.predict(complist=complist) 

1096 

1097 sm.done() 

1098 msg('generation of measurement set '+msfile+' complete', 

1099 origin="simobserve") 

1100 

1101 # rest freqs are hardcoded to the first freq in the spw in core 

1102 tb.open(msfile+"/SPECTRAL_WINDOW/",nomodify=False) 

1103 restfreq=tb.getcol("REF_FREQUENCY") 

1104 for i in range(len(restfreq)): 

1105 restfreq[i]=qa.convert(qa.quantity(model_specrefval),'Hz')['value'] 

1106 tb.putcol("REF_FREQUENCY",restfreq) 

1107 tb.done() 

1108 tb.open(msfile+"/SOURCE/",nomodify=False) 

1109 restfreq=tb.getcol("REST_FREQUENCY") 

1110 for i in range(len(restfreq)): 

1111 restfreq[i]=qa.convert(qa.quantity(model_specrefval),'Hz')['value'] 

1112 tb.putcol("REST_FREQUENCY",restfreq) 

1113 tb.done() 

1114 

1115 

1116 

1117 

1118 

1119 ############################################ 

1120 # create figure  

1121 if grfile: 

1122 file = fileroot + "/" + project + ".observe.png" 

1123 else: 

1124 file = "" 

1125 

1126 # update psfsize using uv coverage instead of maxbase above 

1127 if os.path.exists(msfile): 

1128 # psfsize was set from the antenna posns before, but uv is better 

1129 tb.open(msfile) 

1130 rawdata = tb.getcol("UVW") 

1131 tb.done() 

1132 # TODO make this use the 90% baseline as in aU.getBaselineStats 

1133 maxbase = max([max(rawdata[0,]),max(rawdata[1,])]) # in m 

1134 if maxbase > 0.: 

1135 psfsize = 0.3/qa.convert(qa.quantity(model_specrefval),'GHz')['value']/maxbase*3600.*180/pl.pi # lambda/b converted to arcsec 

1136 if not uvmode: 

1137 msg("An single observation is requested but CROSS-correlation in "+msfile,priority="error") 

1138 else: #Single-dish (zero-spacing) 

1139 psfsize = pb 

1140 if uvmode: 

1141 msg("An interferometer observation is requested but only AUTO-correlation in "+msfile,priority="error") 

1142 minimsize = 8* int(psfsize/cell_asec) 

1143 else: 

1144 msg("Couldn't find "+msfile,priority="error") 

1145 

1146 if uvmode: 

1147 multi = [2,2,1] 

1148 else: 

1149 multi = 0 

1150 

1151 if (grscreen or grfile): 

1152 util.newfig(multi=multi,show=grscreen) 

1153 util.ephemeris(refdate, 

1154 direction=util.direction, 

1155 telescope=telescopename, 

1156 ms=msfile, 

1157 usehourangle=usehourangle, 

1158 cofa=posobs) 

1159 casalog.origin('simobserve') 

1160 if uvmode: 

1161 util.nextfig() 

1162 util.plotants(stnx, stny, stnz, stnd, padnames) 

1163 

1164 # uv coverage 

1165 util.nextfig() 

1166 klam_m = 300/qa.convert(model_specrefval,'GHz')['value'] 

1167 pl.box() 

1168 pl.plot(rawdata[0,]/klam_m,rawdata[1,]/klam_m,'b,') 

1169 pl.plot(-rawdata[0,]/klam_m,-rawdata[1,]/klam_m,'b,') 

1170 ax = pl.gca() 

1171 ax.yaxis.LABELPAD = -4 

1172 pl.xlabel('u[klambda]',fontsize='x-small') 

1173 pl.ylabel('v[klambda]',fontsize='x-small') 

1174 pl.axis('equal') 

1175 

1176 # show dirty beam from observed uv coverage 

1177 util.nextfig() 

1178 imt = imager() 

1179 imt.open(msfile) 

1180 # TODO spectral parms 

1181 msg("using default model cell "+str(model_cell[0])+ 

1182 " for PSF calculation",origin='simobserve') 

1183 imt.defineimage(cellx=str(model_cell[0]["value"])+ 

1184 str(model_cell[0]["unit"]), 

1185 nx=int(max([minimsize,128]))) 

1186 # TODO trigger im.setoptions(ftmachine="mosaic") 

1187 if os.path.exists(fileroot+"/"+project+".quick.psf"): 

1188 shutil.rmtree(fileroot+"/"+project+".quick.psf") 

1189 

1190 # if obs is unknown, casalog will send a warning to screen 

1191 # temporarily(?) suppress that 

1192 if not telescopename in me.obslist(): 

1193 casalog.filter("ERROR") 

1194 imt.approximatepsf(psf=fileroot+"/"+project+".quick.psf") 

1195 if not telescopename in me.obslist(): 

1196 casalog.filter() # set back to default level. 

1197 

1198 quickpsf_current = True 

1199 beam = imt.fitpsf(psf=fileroot+"/"+project+".quick.psf") 

1200 imt.done() 

1201 ia.open(fileroot+"/"+project+".quick.psf") 

1202 beamcs = ia.coordsys() 

1203 beam_array = ia.getchunk(axes=[beamcs.findcoordinate("spectral")['pixel'][0],beamcs.findcoordinate("stokes")['pixel'][0]],dropdeg=True) 

1204 nn = beam_array.shape 

1205 xextent = nn[0]*cell_asec*0.5 

1206 xextent = [xextent,-xextent] 

1207 yextent = nn[1]*cell_asec*0.5 

1208 yextent = [-yextent,yextent] 

1209 flipped_array = beam_array.transpose() 

1210 ttrans_array = flipped_array.tolist() 

1211 ttrans_array.reverse() 

1212 pl.imshow(ttrans_array, 

1213 interpolation='bilinear', 

1214 cmap=pl.cm.jet, 

1215 extent=xextent+yextent, 

1216 origin="lower") 

1217 pl.title(project+".quick.psf",fontsize="x-small") 

1218 b = qa.convert(beam[1],'arcsec')['value'] 

1219 pl.xlim([-3*b,3*b]) 

1220 pl.ylim([-3*b,3*b]) 

1221 ax = pl.gca() 

1222 pl.text(0.05,0.95,"bmaj=%7.1e\nbmin=%7.1e" % (beam[1]['value'],beam[2]['value']),transform = ax.transAxes,bbox=dict(facecolor='white', alpha=0.7),size="x-small",verticalalignment="top") 

1223 ia.done() 

1224 util.endfig(show=grscreen,filename=file) 

1225 

1226 elif thermalnoise != "" or leakage > 0: 

1227 # Not predicting this time but corrupting. Get obsmode from ms. 

1228 if not os.path.exists(msfile): 

1229 msg("Couldn't find "+msfile,priority="error") 

1230 uvmode = (not util.ismstp(msfile,halt=False)) 

1231 

1232 

1233 ###################################################################### 

1234 # noisify 

1235 

1236 noise_any = False 

1237 msroot = fileroot + "/" + project # if leakage, can just copy from this project 

1238 

1239 if thermalnoise != "": 

1240 knowntelescopes = ["ALMASD", "ALMA", "ACA", "SMA", "EVLA", "VLA"] 

1241 

1242 noise_any = True 

1243 

1244 noisymsroot = msroot + ".noisy" 

1245 if not uvmode: #Single-dish 

1246 msroot += ".sd" 

1247 noisymsroot += ".sd" 

1248 

1249 # Cosmic background radiation temperature in K. 

1250 t_cmb = 2.725 

1251 

1252 

1253 # check for interferometric ms: 

1254 if not os.path.exists(msroot+".ms"): 

1255 msg("Couldn't find "+msroot+".ms",priority="error") 

1256 

1257 # MS exists 

1258 message = 'copying '+msroot+'.ms to ' + \ 

1259 noisymsroot+'.ms and adding thermal noise' 

1260 if verbose: 

1261 msg(" ",priority="info") 

1262 msg(message,origin="noise",priority="warn") 

1263 else: 

1264 msg(message,origin="noise") 

1265 

1266 if os.path.exists(noisymsroot+".ms"): 

1267 shutil.rmtree(noisymsroot+".ms") 

1268 shutil.copytree(msfile,noisymsroot+".ms") 

1269 if sm.name() != '': 

1270 emsg = "table persistence error on %s" % sm.name() 

1271 raise RuntimeError(emsg) 

1272 

1273 # if not predicted this time, get telescopename from ms 

1274 if not predict: 

1275 tb.open(noisymsroot+".ms/OBSERVATION") 

1276 n = tb.getcol("TELESCOPE_NAME") 

1277 telescopename = n[0] 

1278 # todo add check that entire column is the same 

1279 tb.done() 

1280 msg("telescopename read from "+noisymsroot+".ms: "+ 

1281 telescopename) 

1282 

1283 if telescopename not in knowntelescopes: 

1284 msg("thermal noise only works properly for ALMA/ACA, (E)VLA, and SMA",origin="simobserve",priority="warn") 

1285 eta_p, eta_s, eta_b, eta_t, eta_q, t_rx = util.noisetemp(telescope=telescopename,freq=model_specrefval) 

1286 

1287 # antenna efficiency 

1288 eta_a = eta_p * eta_s * eta_b * eta_t 

1289 if verbose: 

1290 msg('antenna efficiency = '+str(eta_a), origin="simobserve") 

1291 msg('spillover efficiency = '+str(eta_s), origin="simobserve") 

1292 msg('correlator efficiency = '+str(eta_q), origin="simobserve") 

1293 # sensitivity constant 

1294 scoeff = -1 #Force setting the default value, 1./sqrt(2.0) 

1295 if not uvmode: #Single-dish 

1296 scoeff = 1.0 

1297 if verbose: msg('sensitivity constant = '+str(scoeff), 

1298 origin="simobserve") 

1299 

1300 sm.openfromms(noisymsroot+".ms") # an existing MS 

1301 sm.setdata(fieldid=[]) # force to get all fields 

1302 sm.setseed(seed) 

1303 if thermalnoise == "tsys-manual": 

1304 if verbose: 

1305 message = "sm.setnoise(spillefficiency="+str(eta_s)+\ 

1306 ",correfficiency="+str(eta_q)+\ 

1307 ",antefficiency="+str(eta_a)+\ 

1308 ",trx="+str(t_rx)+",tau="+str(tau0)+\ 

1309 ",tatmos="+str(t_sky)+",tground="+str(t_ground)+\ 

1310 ",tcmb="+str(t_cmb) 

1311 if not uvmode: message += ",senscoeff="+str(scoeff) 

1312 message += ",mode='tsys-manual')" 

1313 msg(message); 

1314 msg("** this may be slow if your MS is finely sampled in time ** ",priority="warn") 

1315 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q, 

1316 antefficiency=eta_a,trx=t_rx, 

1317 tau=tau0,tatmos=t_sky,tground=t_ground,tcmb=t_cmb, 

1318 mode="tsys-manual",senscoeff=scoeff) 

1319 else: 

1320 if verbose: 

1321 message = "sm.setnoise(spillefficiency="+str(eta_s)+\ 

1322 ",correfficiency="+str(eta_q)+\ 

1323 ",antefficiency="+str(eta_a)+\ 

1324 ",trx="+str(t_rx)+",tground="+str(t_ground)+\ 

1325 ",tcmb="+str(t_cmb) 

1326 if not uvmode: message += ",senscoeff="+str(scoeff) 

1327 message += ",mode='tsys-atm'"+\ 

1328 ",pground='560mbar',altitude='5000m'"+\ 

1329 ",waterheight='2km',relhum=20,pwv="+str(user_pwv)+"mm)" 

1330 msg(message); 

1331 msg("** this may be slow if your MS is finely sampled in time ** ",priority="warn") 

1332 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q, 

1333 antefficiency=eta_a,trx=t_rx, 

1334 tground=t_ground,tcmb=t_cmb,pwv=str(user_pwv)+"mm", 

1335 mode="tsys-atm",table=noisymsroot,senscoeff=scoeff) 

1336 # don't set table, that way it won't save to disk 

1337 # mode="calculate",table=noisymsroot) 

1338 

1339 sm.corrupt(); 

1340 sm.done(); 

1341 

1342 

1343 msroot = noisymsroot 

1344 if verbose: msg("done corrupting with thermal noise",origin="noise") 

1345 

1346 

1347 if leakage > 0: 

1348 noise_any = True 

1349 # TODO: need to handle SD name when leakage is available 

1350 if msroot == fileroot+"/"+project: 

1351 noisymsroot = fileroot + "/" + project + ".noisy" 

1352 else: 

1353 noisymsroot = fileroot + "/" + project + ".noisier" 

1354 if not uvmode: #Single-dish 

1355 msg("Can't corrupt SD data with polarization leakage", 

1356 priority="warn") 

1357 if os.path.exists(msfile): 

1358 msg('copying '+msfile+' to ' + 

1359 noisymsroot+'.ms and adding polarization leakage', 

1360 origin="noise",priority="warn") 

1361 if os.path.exists(noisymsroot+".ms"): 

1362 shutil.rmtree(noisymsroot+".ms") 

1363 shutil.copytree(msfile,noisymsroot+".ms") 

1364 if sm.name() != '': 

1365 emsg = "table persistence error on %s" % sm.name() 

1366 raise RuntimeError(emsg) 

1367 

1368 sm.openfromms(noisymsroot+".ms") # an existing MS 

1369 sm.setdata(fieldid=[]) # force to get all fields 

1370 sm.setleakage(amplitude=leakage,table=noisymsroot+".cal") 

1371 sm.corrupt(); 

1372 sm.done(); 

1373 

1374 

1375 

1376 # cleanup - delete newmodel, newmodel.flat etc 

1377# if os.path.exists(modelflat): 

1378# shutil.rmtree(modelflat) 

1379 if os.path.exists(modelflat+".regrid"): 

1380 shutil.rmtree(modelflat+".regrid") 

1381 if os.path.exists(fileroot+"/"+project+".noisy.T.cal"): 

1382 shutil.rmtree(fileroot+"/"+project+".noisy.T.cal") 

1383 if os.path.exists(fileroot+"/"+project+".noisy.sd.T.cal"): 

1384 shutil.rmtree(fileroot+"/"+project+".noisy.sd.T.cal") 

1385 

1386 finally: 

1387 finalize_tools() 

1388 

1389##### Helper functions to plot primary beam 

1390def plotpb(pb,axes,lims=None,color='k'): 

1391 # This beam is automatically scaled when you zoom in/out but 

1392 # not anchored in plot area. We'll wait for Matplotlib 0.99 

1393 # for that function. 

1394 #major=major 

1395 #minor=minor 

1396 #rangle=rangle 

1397 #bwidth=max(major*pl.cos(rangle),minor*pl.sin(rangle))*1.1 

1398 #bheight=max(major*pl.sin(rangle),minor*pl.cos(rangle))*1.1 

1399 from matplotlib.patches import Rectangle, Circle #,Ellipse 

1400 try: 

1401 from matplotlib.offsetbox import AnchoredOffsetbox, AuxTransformBox 

1402 box = AuxTransformBox(axes.transData) 

1403 box.set_alpha(0.7) 

1404 circ = Circle((pb,pb),radius=pb/2.,color=color,fill=False,\ 

1405 label='primary beam',linewidth=2.0) 

1406 box.add_artist(circ) 

1407 pblegend = AnchoredOffsetbox(loc=3,pad=0.2,borderpad=0.,\ 

1408 child=box,prop=None,frameon=False)#,frameon=True) 

1409 pblegend.set_alpha(0.7) 

1410 axes.add_artist(pblegend) 

1411 except: 

1412 casalog.post("Using old matplotlib substituting with circle") 

1413 # work around for old matplotlib 

1414 boxsize = pb*1.1 

1415 if not lims: lims = axes.get_xlim(),axes.get_ylim() 

1416 incx = 1 

1417 incy = 1 

1418 if axes.xaxis_inverted(): incx = -1 

1419 if axes.yaxis_inverted(): incy = -1 

1420 #ecx = lims[0][0] + bwidth/2.*incx 

1421 #ecy = lims[1][0] + bheight/2.*incy 

1422 ccx = lims[0][0] + boxsize/2.*incx 

1423 ccy = lims[1][0] + boxsize/2.*incy 

1424 

1425 #box = Rectangle((lims[0][0],lims[1][0]),incx*bwidth,incy*bheight, 

1426 box = Rectangle((lims[0][0],lims[1][0]),incx*boxsize,incy*boxsize, 

1427 alpha=0.7,facecolor='w', 

1428 transform=axes.transData) #Axes 

1429 #beam = Ellipse((ecx,ecy),major,minor,angle=rangle, 

1430 beam = Circle((ccx,ccy), radius=pb/2., 

1431 edgecolor='k',fill=False, 

1432 label='beam',transform=axes.transData) 

1433 #props = {'pad': 3, 'edgecolor': 'k', 'linewidth':2, 'facecolor': 'w', 'alpha': 0.5} 

1434 #pl.matplotlib.patches.bbox_artist(beam,axes.figure.canvas.get_renderer(),props=props) 

1435 axes.add_artist(box) 

1436 axes.add_artist(beam) 

1437 

1438def finalize_tools(): 

1439 if ia.isopen(): ia.close() 

1440 sm.close() 

1441 #cl.close() # crashes casa