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

977 statements  

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

1import os 

2import shutil 

3import re 

4#import pdb 

5 

6from casatools import ctsys 

7from casatasks import concat, imregrid, immath, sdimaging, impbcor, simobserve, simanalyze, feather, casalog 

8from .simutil import * 

9from . import sdbeamutil 

10 

11def simalma( 

12 project=None, 

13 dryrun=None, 

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

15 incenter=None, inwidth=None, 

16 complist=None, compwidth=None, 

17 ######## 

18 setpointings=None, 

19 ptgfile=None, 

20 integration=None, direction=None, mapsize=None, 

21 antennalist=None, 

22 hourangle=None, 

23 totaltime=None, 

24 ### 

25 tpnant = None, 

26 tptime = None, 

27 ### 

28 pwv=None, 

29 image=None, 

30 imsize=None, imdirection=None,cell=None, 

31 niter=None, threshold=None, 

32 graphics=None, 

33 verbose=None, 

34 overwrite=None 

35 ): 

36 

37 # Collect a list of parameter values to save inputs 

38 in_params = locals() 

39 

40 #------------------------- 

41 # Create the utility object 

42 myutil = simutil(direction) 

43 if verbose: myutil.verbose = True 

44 msg = myutil.msg 

45 

46 try: 

47 

48 ########################### 

49 # preliminaries 

50 ########################### 

51 

52 # Predefined parameters  

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

54 nyquist = 0.5/1.13 ## Nyquist spacing = PB*nyquist 

55 maptype_int = 'ALMA' 

56 maptype_tp = 'square' 

57 

58 # time ratios for 12extended, 12compact,7m, and TP: 

59 default_timeratio = [1,0.5,2,4] 

60 

61 # pbgridratio_tp = 0.25 # -> would be gridratio_tp = 1/3.4 

62 # the grid spacing is defined in terms of lambda/d times these factors 

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

64 gridratio_int = 1./pl.sqrt(3) # 0.5 is nyquist  

65 gridratio_tp = 1./3 

66 

67 # number of 12m primary beams to pad the total power image during 

68 # the gridding stage (i.e. even larger pad than the padding  

69 # added for the observation).  

70 tppad_npb = 2. 

71 

72 # weight of 7m data relative to 12m data 

73 weightratio_7_12 = 0.34 

74 

75 

76 # the scale factor to correct expected Gauss PB size to empirical simPB 

77 simpb_factor = 0.96 

78 caldirection = "" 

79 calflux = "0Jy" 

80 tpantid = 0 

81 t_ground = 270. 

82 if pwv > 0: 

83 thermalnoise = "tsys-atm" 

84 else: 

85 thermalnoise = "" 

86 leakage = 0. 

87 weighting = "briggs" 

88 

89 antlist_tp_default="aca.tp.cfg" 

90 

91 #---------------------------------------- 

92 # Save outputs in a directory called "project" 

93 fileroot = project 

94 # simalma is not supposed to run multiple times. 

95 if os.path.exists(fileroot): 

96 infomsg = "Project directory, '%s', already exists." % fileroot 

97 if overwrite: 

98 casalog.post(infomsg) 

99 casalog.post("Removing old project directory '%s'" % fileroot) 

100 shutil.rmtree(fileroot) 

101 else: 

102 raise RuntimeError(infomsg) 

103 

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

105 os.mkdir(fileroot) 

106 

107 concatname=fileroot+"/"+project+".concat" 

108 

109 #------------------------- 

110 # set up reporting to file, terminal, logger 

111 casalog.origin('simalma') 

112 if verbose: 

113 casalog.filter(level="DEBUG2") 

114 v_priority = "INFO" 

115 else: 

116 v_priority = None 

117 

118 simobserr = "simalma caught an exception in task simobserve" 

119 simanaerr = "simalma caught an exception in task simanalyze" 

120 

121 # open report file 

122 myutil.reportfile=fileroot+"/"+project+".simalma.report.txt" 

123 myutil.openreport() 

124 

125 

126 

127 #---------------------------------------- 

128 # Get globals to call saveinputs() 

129 casalog.post("saveinputs not available in casatasks, skipping saving simalma inputs", priority='WARN') 

130 

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

132 # can contain the cfg 

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

134 

135 #-------------------------- 

136 # format mapsize 

137 if not is_array_type(mapsize): 

138 mapsize = [mapsize, mapsize] 

139 elif len(mapsize) < 2: 

140 mapsize = [mapsize[0], mapsize[0]] 

141 

142 #--------------------------- 

143 # Operation flags 

144 addnoise = (thermalnoise != '') 

145 # Rectangle setup mode 

146 multiptg = (not setpointings) \ 

147 or (is_array_type(direction) and len(direction) > 1) 

148 rectmode = (not multiptg) 

149 

150 # Use full model image as a mapsize = ["", ""] 

151 fullsize = (len(mapsize[0]) == 0) or (len(mapsize[1]) == 0) 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 ########################### 

162 # analyze input sky model 

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

164 msg("="*60,priority="INFO") 

165 msg("Checking the sky model",priority="INFO",origin="simalma") 

166 msg(" ",priority="INFO") 

167 

168 #---------------------------- 

169 # Either skymodel or complist should exist 

170 if is_array_type(skymodel): 

171 if len(skymodel)>1: 

172 msg("You have given more than one skymodel - simalma will only use the first one, "+skymodel[0],priority="INFO") 

173 skymodel = skymodel[0] 

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

175 

176 if is_array_type(complist): 

177 if len(complist)>1: 

178 msg("You have given more than one componentlist - simalma will only use the first one, "+componentlist[0],priority="INFO") 

179 complist = complist[0] 

180 

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

182 if len(skymodel)>0: 

183 msg("Your skymodel '"+skymodel+"' could not be found.",priority="warn") 

184 if len(complist)>0: 

185 msg("Your complist '"+complist+"' could not be found.",priority="warn") 

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

187 msg("At least one of skymodel or complist must be set.",priority="error") 

188 

189 else: 

190 msg("No sky input found. At least one of skymodel or complist must exist.",priority="error") 

191 

192 

193 #------------------------- 

194 # Get model_size and model_center 

195 # TODO: check if outmodel==inmodel works (just collect info)  

196 if os.path.exists(skymodel): 

197 components_only=False 

198 outmodel = fileroot+"/"+project+"temp.skymodel" 

199 model_vals = myutil.modifymodel(skymodel, outmodel, inbright, 

200 indirection, incell, incenter, 

201 inwidth, -1, False) 

202 shutil.rmtree(outmodel) 

203 model_refdir = model_vals[0] 

204 model_cell = model_vals[1] 

205 model_size = model_vals[2] 

206 model_nchan = model_vals[3] 

207 model_center = model_vals[4] 

208 model_width = model_vals[5] 

209 del model_vals 

210 msg("You will be simulating from sky model image "+skymodel,priority="info") 

211 msg(" pixel(cell) size = "+ qa.tos(model_cell[0]),priority="info") 

212 msg(" image size = "+ qa.tos(model_size[0]),priority="info") 

213 msg(" reference direction = "+model_refdir,priority="info") 

214 if model_nchan>1: 

215 msg(" # channels = "+ qa.tos(model_nchan),priority="info") 

216 msg(" channel width = "+ qa.tos(model_width),priority="info") 

217 else: 

218 msg(" single channel / continuum image, with bandwidth = "+qa.tos(model_width),priority="info") 

219 

220 if os.path.exists(complist): 

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

222 msg("You will also be simulating the components in "+complist,priority="info") 

223 msg(" These will get added to a copy of the skymodel image.",priority="info") 

224 else: 

225 # XXX TODO make sure components AND image work here 

226 components_only=True 

227 compdirs = [] 

228 cl.open(complist) 

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

230 compdirs.append(myutil.dir_m2s(cl.getrefdir(i))) 

231 

232 model_refdir, coffs = myutil.average_direction(compdirs) 

233 model_center = cl.getspectrum(0)['frequency']['m0'] 

234 cmax = 0.0014 # ~5 arcsec 

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

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

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

238 cmax = max(cmax, xc) 

239 cmax = max(cmax, yc) 

240 model_size = ["%fdeg" % (2*cmax), "%fdeg" % (2*cmax)] 

241 cl.done() 

242 model_cell = None 

243 model_nchan = 1 

244 del compdirs, coffs, xc, yc, cmax 

245 msg("You will be simulating only from component list "+complist,priority="info") 

246 msg("Based on the spatial distribution of components, a sky model image will be generated with these parameters:",priority="info") 

247 msg(" image size = "+ qa.tos(model_size[0]),priority="info") 

248 msg(" reference direction = "+model_refdir,priority="info") 

249 msg("and each call to simobserve will chose a skymodel cell size of 1/20 the expected PSF FWHM (to accurately locate components).",priority="info") 

250 msg("Simulation from components-only produces a single channel / continuum observation",priority="info") 

251 msg(" with bandwidth = "+compwidth,priority="info") 

252 

253 

254 #----------------------------------- 

255 # determine mapsize - copied code from simobserve 

256 # if the user has not input a map size, then use model_size 

257 if len(mapsize) == 0: 

258 mapsize = model_size 

259 msg("setting map size to "+str(model_size)) 

260 else: 

261 if type(mapsize) == type([]): 

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

263 mapsize = model_size 

264 msg("setting map size to "+str(model_size)) 

265 

266 if components_only: 

267 # if map is greater tham model (defined from components)  

268 # then use mapsize as modelsize 

269 if type(mapsize) == type([]): 

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

271 else: 

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

273 if type(model_size) == type([]): 

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

275 else: 

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

277 if map_asec>mod_asec: 

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

279 

280 

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

282 msg("You will be mapping an area of size "+qa.tos(mapsize[0])+','+qa.tos(mapsize[1])) 

283 

284 

285 

286 

287 

288 

289 

290 

291 ########################### 

292 # Calculate 12-m PB and grid spacing for int and tp 

293 ########################### 

294 

295 Dant = 12. 

296 wave_length = 0.2997924/ qa.convert(qa.quantity(model_center),'GHz')['value'] 

297 # (wave length)/D_ant in arcsec 

298 lambda_D = wave_length / Dant * 3600. * 180 / pl.pi 

299 PB12 = qa.quantity(lambda_D*pbcoeff, "arcsec") 

300 # Correction factor for PB in simulation 

301 # (PS of simulated image is somehow smaller than PB12) 

302 PB12sim = qa.mul(PB12, simpb_factor) 

303 

304 msg(" primary beam size: %s" % (qa.tos(PB12)),priority="info") 

305 

306 # Pointing spacing of observations - define in terms of lambda/d 

307 # to avoid ambiguities in primary beam definition.  

308 # it would be best for simobserve to accept the shorthand "LD" 

309 # but until it does that, we'll divide by pbcoeff and use "PB": 

310 # this is a fragile solution since it depends on pbcoeff being 

311 # the same in simalma and simobserve. 

312 pointingspacing = str(gridratio_int/pbcoeff)+"PB" 

313 ptgspacing_tp = qa.tos(qa.mul(PB12sim, (gridratio_tp/pbcoeff))) 

314 

315 

316 

317 

318 

319 

320 

321 

322 

323 

324 ########################### 

325 # analyze antennalist(s)  

326 ########################### 

327 

328 if type(antennalist)==type(" "): 

329 antennalist=[antennalist] 

330 

331 # number of arrays/configurations to run: 

332 nconfigs=len(antennalist) 

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

334 msg("="*60,priority="info") 

335 msg("You are requesting %i configurations:" % nconfigs,origin="simalma",priority="info") 

336 

337 configtypes=[] # ALMA, ACA, or ALMASD (latter is special case) 

338 cycles=[] 

339 resols=[] 

340 tp_inconfiglist=False 

341 

342 # check filename consistency, get resolution, etc  

343 for configfile in antennalist: 

344 

345 isalma=0 

346 

347 # antennalist should contain ALMA or ACA 

348 if configfile.find(";")>=0: 

349 telescopename="BYRES" 

350 configparms=configfile.split(";") 

351 res=configparms[-1] 

352 res_arcsec=-1 

353 if myutil.isquantity(res): 

354 if qa.compare(res,"arcsec"): 

355 res_arcsec=qa.convert(res,"arcsec")['value'] 

356 

357 if res_arcsec>0: 

358 resols.append(res_arcsec) 

359 else: 

360 msg("simalma cannot interpret the antennalist entry '"+configfile+"' as desired array and resolution e.g. ALMA;0.5arcsec",priority="ERROR") 

361 

362 else: 

363 # we can only verify explicit files right now, not  

364 # configs specified as "ALMA;0.5arcsec" 

365 # 

366 configfile_short = (configfile.split("/"))[-1] 

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

368 if len(configfile) > 0: 

369 if os.path.exists(fileroot+"/"+configfile): 

370 configfile = fileroot + "/" + configfile 

371 elif not os.path.exists(configfile) and \ 

372 os.path.exists(os.path.join(repodir,configfile)): 

373 configfile = os.path.join(repodir,configfile) 

374 # Now make sure the configfile exists 

375 if not os.path.exists(configfile): 

376 msg("Couldn't find configfile: %s" % configfile, priority="error") 

377 stnx, stny, stnz, stnd, padnames, antnames, telescopename, posobs = myutil.readantenna(configfile) 

378 if telescopename=="ALMASD": 

379 resols.append(qa.convert(PB12,'arcsec')['value']) 

380 else: 

381 psfsize = myutil.approxBeam(configfile,qa.convert(qa.quantity(model_center),'GHz')['value']) 

382 resols.append(psfsize) # FWHM in arcsec 

383 

384 q = re.compile('CYCLE\d?\d') 

385 isCycle = q.search(configfile_short.upper()) 

386 if isCycle: 

387 whatCycle = isCycle.group()[-1] # will break if cycle>9 

388 cycles.append(whatCycle) 

389 else: 

390 cycles.append(-1) # -1 is unknown; defaults to full ALMA 

391 

392 if configfile_short.upper().find("ALMA") >= 0: 

393 if telescopename=="ALMA" or telescopename=="BYRES": 

394 configtypes.append("ALMA") 

395 isalma=isalma+1 

396 else: 

397 msg("Inconsistent antennalist entry '"+configfile_short+"' has ALMA in the name but not set as the observatory",priority="error") 

398 if configfile_short.upper().find("ACA") >= 0: 

399 if telescopename=="ACA" or telescopename=="BYRES": 

400 configtypes.append("ACA") 

401 isalma=isalma+1 

402 elif telescopename=="ALMASD" or telescopename=="BYRES": 

403 tp_inconfigfile=True 

404 configtypes.append("ALMASD") 

405 isalma=isalma+1 

406 else: 

407 msg("Inconsistent antennalist entry '"+configfile_short+"' has ACA in the name but the observatory is not ACA or ALMASD",priority="error") 

408 #if configfile.upper().find("CYCLE0") 

409 

410 

411 if isalma==0: 

412 s="simalma can't accept antennalist entry '"+configfile_short+"' because it is neither ALMA nor ACA (in the name and the observatory in the file)" 

413 msg(s,origin="simalma",priority="error") 

414# raise ValueError, s # not ness - msg2 raises the exception 

415 if isalma==2: 

416 s="simalma doesn't understand your antennalist entry '"+configfile_short+"'" 

417 msg(s,origin="simalma",priority="error") 

418# raise ValueError,s 

419 

420 #----------------------------- 

421 # total power parameter: 

422 tptime_min=0. 

423 if myutil.isquantity(tptime,halt=False): 

424 if qa.compare(tptime,'s'): 

425 tptime_min=qa.convert(tptime,'min')['value'] 

426 else: 

427 msg("Can't interpret tptime='" 

428 +tptime+"' as a time quantity e.g. '3h'", 

429 priority="error") 

430 else: 

431 msg("Can't interpret tptime='" 

432 +tptime+"' as a time quantity e.g. '3h'", 

433 priority="error") 

434 

435 

436 #----------------------------- 

437 # is there a mix of cycle0,1,2,etc,full? 

438 whatCycle=cycles[0] 

439 cyclesInconsistent=False 

440 for i in range(nconfigs): 

441 if cycles[i]!=whatCycle: 

442 cyclesInconsistent=True 

443 

444 

445 #----------------------------- 

446 # exposure time parameter totaltime 

447 totaltime_min=[] 

448 if len(totaltime)==1: 

449 # scalar input - use defaults 

450 totaltime_min0=0. 

451 if not myutil.isquantity(totaltime[0],halt=False): 

452 raise ValueError("Can't interpret totaltime parameter '"+totaltime[0]+"' as a time quantity - example quantities: '1h', '20min', '3600sec'") 

453 if qa.compare(totaltime[0],'s'): 

454 totaltime_min0=qa.convert(totaltime[0],'min')['value'] 

455 else: 

456 raise ValueError("Can't convert totaltime parameter '"+totaltime[0]+"' to minutes - example quantities: '1h', '20min','3600sec'") 

457 totaltime_min=pl.zeros(nconfigs) 

458 # sort by res'l - TP could still be on here 

459 resols=pl.array(resols) 

460 intorder=resols.argsort() 

461 # assume the scalar user input refers to the highest res 12m 

462 totaltime_min[intorder[0]]=totaltime_min0 

463 for j in intorder[1:]: 

464 if configtypes[j]=='ALMA': 

465 # lower res 12m 

466 totaltime_min[j]=totaltime_min0*default_timeratio[1] 

467 elif configtypes[j]=='ACA': 

468 # 7m 

469 totaltime_min[j]=totaltime_min0*default_timeratio[2] 

470 elif configtypes[j]=='ALMASD': 

471 # tp 

472 totaltime_min[j]=totaltime_min0*default_timeratio[3] 

473 else: 

474 raise ValueError("configuration types = "+str(configtypes)) 

475 else: 

476 for time in totaltime: 

477 if not myutil.isquantity(time,halt=False): 

478 raise ValueError("Can't interpret totaltime vector element '"+time+"' as a time quantity - example quantities: '1h', '20min', '3600sec'") 

479 if qa.compare(time,'s'): 

480 time_min=qa.convert(time,'min')['value'] 

481 else: 

482 raise ValueError("Can't convert totaltime vector element '"+time+"' to minutes - example quantities: '1h', '20min','3600sec'") 

483 totaltime_min.append(time_min) 

484 

485 if len(totaltime_min)!=len(antennalist): 

486 raise ValueError("totaltime must either be the same length vector as antennalist or a scalar") 

487 

488 

489 

490 #----------------------------- 

491 # print out what's requested.  

492 antlist_tp=None 

493 for i in range(nconfigs): 

494 configfile=antennalist[i] 

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

496 msg(" "+configfile+":",priority="info") 

497 

498 if configtypes[i]=="ALMA": 

499 msg(" 12m ALMA array",priority="info") 

500 elif configtypes[i]=="ACA": 

501 msg(" 7m ACA array",priority="info") 

502 elif configtypes[i]=="ALMASD": 

503 msg(" 12m total power observation",priority="info") 

504 if tpnant>0: 

505 msg(" This antennalist entry will be ignored in favor of the tpnant requesting %d total power antennas" % tpnant,priority="info") 

506 antlist_tp=antlist_tp_default 

507 if tptime_min>0: 

508 msg(" The total map integration time will be %s minutes as per the tptime parameter." % tptime_min,priority="info") 

509 else: 

510 tptime_min=totaltime_min[i] 

511 msg(" The total map integration time will be %s minutes as per the totaltime parameter." % tptime_min,priority="info") 

512 else: 

513 # Note: assume the user won't put in >1 TP antlist. 

514 msg(" Note: it is preferred to specify total power with the tpnant,tptime parameters so that one can also specify the number of antennas.",priority="info") 

515 msg(" Specifying the total power array as one of the antenna configurations will only allow you to use a single antenna for the simulation",priority="info") 

516 msg(" simalma will proceed with a single total power antenna observing for %d minutes." % totaltime_min[i],priority="info") 

517 tpnant=1 

518 tptime_min=totaltime_min[i] 

519 antlist_tp=configfile 

520 

521 # print cycle and warn if mixed 

522 if cycles[i]>'-1': 

523 msg(" This is a cycle "+cycles[i]+" configuration",priority="info") 

524 else: 

525 msg(" This is a full ALMA configuration",priority="info") 

526 if cyclesInconsistent: 

527 msg(" WARNING: Your choices of configurations mix different cycles and/or Full ALMA. Assuming you know what you want.",priority="info") 

528 

529 

530 if configtypes[i]=='ALMASD': 

531 # imsize check for PB: 

532 if not components_only: 

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

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

535 PB12sec=qa.convert(PB12,"arcsec")['value'] 

536 if minsize < 2.5*PB12sec: 

537 msg(" WARNING: For TP imaging, skymodel should be larger than 2.5*primary beam. Your skymodel: %.3f arcsec < %.3f arcsec: 2.5*primary beam" % (minsize, 2.5*PB12sec),priority="warn") 

538 msg(" Total power imaging when the source is this small is not necessary, and may fail with an error.",priority="warn") 

539 del minsize 

540 else: 

541 # casalog.post resolution for INT, and warn if model_cell too large 

542 msg(" approximate synthesized beam FWHM = %f arcsec" % resols[i],priority="info") 

543 msg(" (at zenith; the actual beam will depend on declination, hourangle, and uv coverage)",priority="info") 

544 

545 if is_array_type(model_cell): 

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

547 else: 

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

549 if cell_asec > 0.2*resols[i]: 

550 if cell_asec >= resols[i]: 

551 # XXX if not dryrun raise exception here 

552 msg(" ERROR: your sky model cell %f arcsec is too large to simulate this beam. Simulation will fail" % cell_asec,priority="info") 

553 else: 

554 msg(" WARNING: your sky model cell %f arcsec does not sample this beam well. Simulation may be unreliable or clean may fail." % cell_asec,priority="info") 

555 

556 msg(" Observation time = %f min" % totaltime_min[i],priority="info") 

557 if totaltime_min[i]>360: 

558 msg(" WARNING: this is longer than 6hr - simalma won't (yet) split the observation into multiple nights, so your results may be unrealistic",priority="info") 

559 

560 

561 if not tp_inconfiglist and tpnant>0: 

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

563 if tptime_min>0: 

564 msg("You are also requesting Total Power observations:",priority="info") 

565 msg(" %d antennas for %f minutes" % (tpnant,tptime_min),priority="info") 

566 # imsize check for PB: 

567 if not components_only: 

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

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

570 PB12sec=qa.convert(PB12,"arcsec")['value'] 

571 if minsize < 2.5*PB12sec: 

572 msg(" WARNING: For TP imaging, skymodel should be larger than 2.5*primary beam. Your skymodel: %.3f arcsec < %.3f arcsec: 2.5*primary beam" % (minsize, 2.5*PB12sec),priority="warn") 

573 msg(" Total power imaging when the source is this small is not necessary, and may fail with an error.",priority="warn") 

574 

575 else: 

576 msg("You have requested %d total power antennas (tpnant), but no finite integration (tptime) -- check your inputs; no Total Power will be simulated." % tpnant,priority="info") 

577 

578 ### WORKAROUND for wrong flux in COMP TP simulations 

579 if (tpnant > 0) and os.path.exists(complist): 

580 idx_min = pl.where(resols == min(resols))[0] 

581 idx = idx_min[0] if len(idx_min) > 0 else 0 

582 dummy_proj = "gen_skymodel" 

583 errmsg = "You requested Single dish simulation with components list.\n" 

584 errmsg += "Single dish simulation has flux recovery issue "+\ 

585 "when using a components list as an input.\n" 

586 errmsg += "Please generate compskymodel image first by task "+\ 

587 "simobserve and use the image as the skymodel input. " 

588 errmsg += "Sorry for the inconvenience.\n\n" 

589 errmsg += "How to workaround the issue:\n" 

590 errmsg += "1. Generate skymodel image by simobserve\n" 

591 errmsg += ("\tsimobserve(project='%s', complist='%s', compwidth='%s', "\ 

592 % (dummy_proj, complist, compwidth)) 

593 if os.path.exists(skymodel): 

594 skysuffix = '.skymodel' 

595 errmsg += ( "skymodel='%s', inbright='%s', indirection='%s', " \ 

596 % (skymodel, inbright, indirection)) 

597 errmsg += ( "incell='%s', incenter='%s', inwidth='%s', " \ 

598 % (skymodel, inbright, indirection, incell, incenter, inwidth) ) 

599 else: 

600 skysuffix = '.compskymodel' 

601 errmsg += ("setpointings=True, obsmode='', antennalist='%s', thermalnoise='')\n" \ 

602 % antennalist[idx]) 

603 errmsg += "2. Use the generated skymodel image in project directory as an input of simalma.\n" 

604 errmsg += ("\tsimalma(project='%s', skymodel='%s/%s', complist='', ....)" % \ 

605 (project, dummy_proj, \ 

606 get_data_prefix(antennalist[idx], dummy_proj)+skysuffix)) 

607 msg(errmsg,priority="error") 

608 ### End of WORKAROUND 

609 

610 # remove tp from configlist 

611 antennalist=pl.array(antennalist) 

612 configtypes=pl.array(configtypes) 

613 totaltime_min=pl.array(totaltime_min) 

614 resols=pl.array(resols) 

615 z=pl.where(configtypes!='ALMASD')[0] 

616 antennalist=antennalist[z] 

617 configtypes=configtypes[z] 

618 totaltime_min=totaltime_min[z] 

619 resols=resols[z] 

620 nconfigs=len(antennalist) 

621 if nconfigs < 1: 

622 msg("No interferometer configuration is requested. At least one interferometer configuration should be selected.", \ 

623 origin="simalma", priority="error") 

624 

625 

626# TODO check model_size against mapsize - separately after this? 

627 

628 

629 

630 

631 

632 

633 

634 

635 

636 

637 

638 

639 

640 

641# ########################### 

642# # Resolve prefixes of simulation data (as defined in  

643# # simobserve and simanalyze) 

644#  

645# pref_bl = get_data_prefix(antennalist, project) 

646# pref_aca = get_data_prefix(antlist_aca, project) 

647 # Resolve output names (as defined in simobserve and simanalyze) 

648# if addnoise: 

649# msname_bl = pref_bl+".noisy.ms" 

650# msname_aca = pref_aca+".noisy.ms" 

651# msname_tp = pref_tp+".noisy.sd.ms" 

652# #imagename_aca = pref_aca+".noisy.image" 

653# else: 

654# msname_bl = pref_bl+".ms" 

655# msname_aca = pref_aca+".ms" 

656# msname_tp = pref_tp+".sd.ms" 

657# #imagename_aca = pref_aca+".image" 

658#  

659# imagename_tp = project+".sd.image" 

660# imagename_int = project+".concat.image" 

661# msname_concat = project+".concat.ms" 

662# 

663 combimage = project+".feather.image" 

664 simana_file = project+".simanalyze.last" 

665 

666 

667 

668 

669 ############################################################ 

670 # run simobserve  

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

672 msg("="*60,priority="info") 

673 msg("simalma calls simobserve to simulate each component:",priority="info",origin="simalma") 

674 

675 # sort by res'l - save the ptgfile from the max res'l config  

676 # for use in TP pointing definition. 

677 resols=pl.array(resols) 

678 intorder=resols.argsort() 

679 pref_bl=get_data_prefix(antennalist[intorder[0]],project) 

680 ptgfile_bl=fileroot+"/"+pref_bl+".ptg.txt" 

681 

682 step = 0 

683 

684 for i in range(nconfigs): 

685 if configtypes[i]=="ALMA": 

686 s="12m ALMA array" 

687 elif configtypes[i]=='ACA': 

688 s="7m ACA array" 

689 else: 

690 s="12m total power map" 

691 

692 if configtypes[i]!='ALMASD': 

693 step += 1 

694 

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

696 msg("-"*60, origin="simalma", priority="warn") 

697 msg(("Step %d: simulating " % step)+s, origin="simalma", priority="warn") 

698 msg("-"*60, origin="simalma", priority="warn") 

699 

700 # filename prefixes 

701 pref = get_data_prefix(antennalist[i], project) 

702 

703 obsmode_int = 'int' 

704 ptgfile_int = fileroot+"/"+pref+".ptg.txt" 

705 

706 task_param = {} 

707 task_param['project'] = project 

708 task_param['skymodel'] = skymodel 

709 task_param['inbright'] = inbright 

710 task_param['indirection'] = indirection 

711 task_param['incell'] = incell 

712 task_param['incenter'] = incenter 

713 task_param['inwidth'] = inwidth 

714 task_param['complist'] = complist 

715 task_param['compwidth'] = compwidth 

716 task_param['setpointings'] = setpointings 

717 task_param['ptgfile'] = ptgfile 

718 task_param['integration'] = integration 

719 task_param['direction'] = direction 

720 task_param['mapsize'] = mapsize 

721 task_param['maptype'] = maptype_int 

722 # this is approriate for 7m or 12m since its in terms of PB 

723 task_param['pointingspacing'] = pointingspacing 

724 task_param['caldirection'] = caldirection 

725 task_param['calflux'] = calflux 

726 task_param['obsmode'] = obsmode_int 

727 task_param['hourangle'] = hourangle 

728 task_param['totaltime'] = "%fmin" % totaltime_min[i] 

729 task_param['antennalist'] = antennalist[i] 

730 task_param['sdantlist'] = "" 

731 task_param['sdant'] = 0 

732 task_param['thermalnoise'] = thermalnoise 

733 task_param['user_pwv'] = pwv 

734 task_param['t_ground'] = t_ground 

735 task_param['leakage'] = leakage 

736 task_param['graphics'] = graphics 

737 task_param['verbose'] = verbose 

738 task_param['overwrite'] = overwrite 

739 

740 msg(get_taskstr('simobserve', task_param), priority="info") 

741 if not dryrun: 

742 try: 

743 simobserve(**task_param) 

744 del task_param 

745 except: 

746 raise RuntimeError(simobserr) 

747 finally: 

748 casalog.origin('simalma') 

749 

750 qimgsize_tp = None 

751 

752 if tpnant>0 and tptime_min<=0: 

753 raise ValueError("You requested total power (tpnant=%d) but did not specify a valid nonzero tptime" % tpnant) 

754 

755 

756 mslist_tp = [] 

757 if tptime_min>0: 

758 ######################################################## 

759 # ACA-TP simulation - always do this last 

760 obsmode_sd = "sd" 

761 step += 1 

762 

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

764 msg("-"*60, origin="simalma", priority="warn") 

765 msg(("Step %d: simulating Total Power" % step), origin="simalma", priority="warn") 

766 msg("-"*60, origin="simalma", priority="warn") 

767 

768 if antlist_tp is None: 

769 antlist_tp = antlist_tp_default 

770 

771 pref_tp = get_data_prefix(antlist_tp, project) 

772 if addnoise: 

773 msname_tp = pref_tp+".noisy.sd.ms" 

774 else: 

775 msname_tp = pref_tp+".sd.ms" 

776 

777 

778 ########################### 

779 # Resolve pointing directions of ACA-TP. 

780 # 

781 # Pointing directions of TP simulation is defined as follows: 

782 # 

783 # [I] if ALMA-12m maps a rectangle region (rectmode=T), 

784 # TP maps slightly larger region than ALMA-12m by adding 1 PB to 

785 # mapsize (pointing extent of ALMA-12m). 

786 # 

787 # [II] if a list of pointing deirections are specified for the  

788 # ALMA-12m observation (multiptg=T), TP pointings are defined as 

789 # rectangle areas of 2PB x 2PB centered at each pointing direction 

790 # of ALMA-12m. However, in some cases, it is more efficient to 

791 # just map a rectangle region that covers all ALMA-12m pointings. 

792 # In such case, ACA-TP maps a rectangle region whose extent is 2PB 

793 # larger than the extent of all ALMA-12m pointings. 

794 if rectmode: 

795 # Add 1PB to mapsize 

796 if fullsize: 

797 mapx = qa.add(PB12,model_size[0]) # in the unit same as PB 

798 mapy = qa.add(PB12,model_size[1]) # in the unit same as PB 

799 mapsize_tp = [qa.tos(mapx), qa.tos(mapy)] 

800 msg("The full skymodel is being mapped by ALMA 12-m and ACA 7-m arrays. The total power antenna observes 1PB larger extent.", origin="simalma", priority='warn') 

801 else: 

802 # mapsize is defined. Add 1 PB to mapsize. 

803 mapx = qa.add(qa.quantity(mapsize[0]), PB12) 

804 mapy = qa.add(qa.quantity(mapsize[1]), PB12) 

805 mapsize_tp = [qa.tos(mapx), qa.tos(mapy)] 

806 msg("Only part of the skymodel is being mapped by ALMA 12-m and ACA 7-m arrays. The total power antenna observes 1PB larger extent.", priority='warn') 

807 else: 

808 # multi-pointing mode 

809 npts, pointings, time = myutil.read_pointings(ptgfile_bl) 

810 center, offsets = myutil.average_direction(pointings) 

811 del time 

812 qx = qa.quantity(max(offsets[0])-min(offsets[0]),"deg") 

813 qy = qa.quantity(max(offsets[1])-min(offsets[1]),"deg") 

814 # map extent to cover all pointings + 2PB  

815 mapx = qa.add(qa.mul(PB12,2.),qx) # in the unit same as PB 

816 mapy = qa.add(qa.mul(PB12,2.),qy) # in the unit same as PB 

817 mapsize_tp = [qa.tos(mapx), qa.tos(mapy)] 

818 # number of pointings to map vicinity of each pointings 

819 qptgspc_tp = qa.quantity(ptgspacing_tp) 

820 dirs_multi_tp = myutil.calc_pointings2(qptgspc_tp, 

821 qa.tos(qa.mul(PB12,2.)), 

822 "square", pointings[0]) 

823 npts_multi = npts * len(dirs_multi_tp) 

824 

825 msg("Number of pointings to map vicinity of each direction = %d" % npts_multi) 

826 del qptgspc_tp, dirs_multi_tp 

827 

828 # save imsize for TP image generation 

829# qimgsize_tp = [mapx, mapy] 

830 # for imaging later, we need to pad even more! 

831 mapx_im = qa.add(mapx, qa.mul(PB12,(2*tppad_npb))) 

832 mapy_im = qa.add(mapy, qa.mul(PB12,(2*tppad_npb))) 

833 qimgsize_tp = [mapx_im, mapy_im] 

834 

835 

836 qptgspc_tp = qa.quantity(ptgspacing_tp) 

837 pbunit = PB12['unit'] 

838 # number of pointings to map pointing region 

839 # TODO: use calc pointings for consistent calculation 

840 npts_rect = int(qa.convert(mapx, pbunit)['value'] \ 

841 / qa.convert(qptgspc_tp, pbunit)['value']) \ 

842 * int(qa.convert(mapy, pbunit)['value'] \ 

843 / qa.convert(qptgspc_tp, pbunit)['value']) 

844 msg("Number of pointings to map a rect region = %d" % npts_rect, priority="DEBUG2") 

845 

846 if rectmode: 

847 dir_tp = direction 

848 npts_tp = npts_rect 

849 msg("Rectangle mode: The total power antenna observes 1PB larger region compared to ALMA 12-m and ACA 7-m arrays", priority='info') 

850 else: 

851 if npts_multi < npts_rect: 

852 # Map 2PB x 2PB extent centered at each pointing direction 

853 # need to get a list of pointings 

854 dir_tp = [] 

855 locsize = qa.mul(2, PB12) 

856 for dir in pointings: 

857 dir_tp += myutil.calc_pointings2(qa.tos(qptgspc_tp), 

858 qa.tos(locsize), 

859 "square", dir) 

860 

861 mapsize_tp = ["", ""] 

862 #npts_tp = npts_multi 

863 npts_tp = len(dir_tp) 

864 msg("Multi-pointing mode: The total power antenna observes +-1PB of each point", origin="simalma", priority='warn') 

865 else: 

866 # Map a region that covers all directions 

867 dir_tp = center 

868 npts_tp = npts_rect 

869 msg("Multi-pointing mode: The total power antenna maps a region that covers all pointings", origin="simalma", priority='warn') 

870 msg("- Center of poinings: %s" % center, origin="simalma", priority='warn') 

871 msg("- Map size: [%s, %s]" % (mapsize_tp[0], mapsize_tp[1]), origin="simalma", priority='warn') 

872 

873 

874 # Scale integration time of TP (assure >= 1 visit per direction) 

875 tottime_tp = "%dmin" % tptime_min 

876 integration_tp = integration 

877 ndump = int(qa.convert(tottime_tp, 's')['value'] 

878 / qa.convert(integration, 's')['value']) 

879 msg("Max number of dump in %s (integration %s): %d" % \ 

880 (tottime_tp, integration, ndump), origin="simalma") 

881 

882 if ndump < npts_tp: 

883 t_scale = float(ndump)/float(npts_tp) 

884 integration_tp = qa.tos(qa.mul(integration, t_scale)) 

885 msg("Integration time is scaled to cover all pointings in observation time.", origin="simalma", priority='warn') 

886 msg("-> Scaled total power integration time: %s" % integration_tp, origin="simalma", priority='warn') 

887 ## Sometimes necessary to avoid the effect of round-off error 

888 #iunit = qa.quantity(integration_tp)['unit'] 

889 #intsec = qa.convert(integration_tp,"s") 

890 #totsec = intsec['value']*npts_tp#+0.000000001) 

891 ##tottime_tp = qa.tos(qa.convert(qa.quantity(totsec, "s"), iunit)) 

892 #tottime_tp = qa.tos(qa.quantity(totsec, "s")) 

893 

894 task_param = {} 

895 task_param['project'] = project 

896 task_param['skymodel'] = skymodel 

897 task_param['inbright'] = inbright 

898 task_param['indirection'] = indirection 

899 task_param['incell'] = incell 

900 task_param['incenter'] = incenter 

901 task_param['inwidth'] = inwidth 

902 task_param['complist'] = complist 

903 task_param['compwidth'] = compwidth 

904 task_param['setpointings'] = True 

905 task_param['ptgfile'] = '$project.ptg.txt' 

906 task_param['integration'] = integration_tp 

907 task_param['direction'] = dir_tp 

908 task_param['mapsize'] = mapsize_tp 

909 task_param['maptype'] = maptype_tp 

910 task_param['pointingspacing'] = ptgspacing_tp 

911 task_param['caldirection'] = caldirection 

912 task_param['calflux'] = calflux 

913 task_param['obsmode'] = obsmode_sd 

914 task_param['hourangle'] = hourangle 

915 task_param['totaltime'] = "%dmin" % tptime_min 

916 task_param['antennalist'] = "" 

917 task_param['sdantlist'] = antlist_tp 

918 task_param['sdant'] = tpantid 

919 task_param['thermalnoise'] = thermalnoise 

920 task_param['user_pwv'] = pwv 

921 task_param['t_ground'] = t_ground 

922 task_param['leakage'] = leakage 

923 task_param['graphics'] = graphics 

924 task_param['verbose'] = verbose 

925 task_param['overwrite'] = overwrite 

926 

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

928 msg("Simulating %d TP antennas" % tpnant, priority="info") 

929 for iant in range(tpnant): 

930 task_param['sdant'] = iant 

931 task_param['seed'] = int(pl.random()*100000) 

932 msg(" ",priority=v_priority) 

933 msg("Running TP simulation with sdant = %d" % task_param['sdant'], priority=v_priority) 

934 msg(get_taskstr('simobserve', task_param), priority="info") 

935 if not dryrun: 

936 try: 

937 simobserve(**task_param) 

938 except: 

939 raise RuntimeError(simobserr) 

940 finally: 

941 casalog.origin('simalma') 

942 if tpnant == 1: 

943 mslist_tp = [msname_tp] 

944 else: 

945 # copy MSes 

946 orig_tp = fileroot + "/" + msname_tp 

947 suffix = (".Ant%d" % iant) 

948 msg(" ",priority=v_priority) 

949 msg("Renaming '%s' to '%s'" % (orig_tp, orig_tp+suffix), priority=v_priority) 

950 if not dryrun: 

951 shutil.move(orig_tp, orig_tp+suffix) 

952 

953 mslist_tp.append(msname_tp+suffix) 

954 # noiseless MS 

955 if addnoise: 

956 orig_tp = orig_tp.replace(".noisy", "") 

957 msg("Renaming '%s' to '%s'" % (orig_tp, orig_tp+suffix), priority=v_priority) 

958 if not dryrun: 

959 shutil.move(orig_tp, orig_tp+suffix) 

960 

961 del task_param 

962 

963 

964 

965 

966 

967 

968 

969 

970 

971 

972 

973 ################################################################ 

974 # Imaging 

975 if dryrun: image=True # why not? 

976 

977 # for 4.2 print more info 

978 v_priority="info" 

979 

980 if image: 

981 modelimage = "" 

982 imagename_tp = project+".sd.image" 

983 if tptime_min > 0: 

984 ######################################################## 

985 # Image ACA-TP 

986 step += 1 

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

988 msg("-"*60, origin="simalma", priority="info") 

989 msg("Step %d: generating a total power image. " % step, origin="simalma", priority="info") 

990 msg(" WARNING: Optimal gridding parameters are being analyzed by ALMA and may change.",priority="warn",origin="simalma") 

991 msg("-"*60, origin="simalma", priority="info") 

992 

993 vis_tp = [] 

994 for msname_tp in mslist_tp: 

995 if dryrun: 

996 vis_tp.append(fileroot+"/"+msname_tp) 

997 elif os.path.exists(fileroot+"/"+msname_tp): 

998 vis_tp.append(fileroot+"/"+msname_tp) 

999 else: 

1000 msg("Total power MS '%s' is not found" \ 

1001 % msname_tp, origin="simalma", priority="error") 

1002 

1003 tp_kernel = 'SF' 

1004 #tp_kernel = 'GJINC' 

1005 

1006 # Define imsize to cover TP map region 

1007 msg(" ",priority=v_priority) 

1008 msg("Defining image size to cover map region of total power simulation", priority=v_priority) 

1009 msg("-> The total power map size: [%s, %s]" % \ 

1010 (qa.tos(qimgsize_tp[0]), qa.tos(qimgsize_tp[1])), \ 

1011 priority=v_priority) 

1012 if tp_kernel.upper() == 'SF': 

1013 beamsamp=9 

1014 pb_asec = sdbeamutil.primaryBeamArcsec(qa.tos(qa.convert(qa.quantity(model_center),'GHz')),12.0,0.75,10.0) 

1015 qcell=qa.quantity(pb_asec/beamsamp, 'arcsec') 

1016 cell_tp = [qa.tos(qcell), qa.tos(qcell)] 

1017 msg("Using fixed cell size for SF grid: %s" % cell_tp[0], \ 

1018 priority=v_priority) 

1019 elif cell != '': 

1020 # user-defined cell size 

1021 msg("The user defined cell size: %s" % cell, \ 

1022 priority=v_priority) 

1023 cell_tp = [qa.tos(cell), qa.tos(cell)] 

1024 else: 

1025 if model_cell == None: 

1026 # components only simulation 

1027 compmodel = fileroot+"/"+pref_bl+".compskymodel" 

1028 msg("getting the cell size of input compskymodel", \ 

1029 priority=v_priority) 

1030 if not os.path.exists(compmodel): 

1031 msg("Could not find the skymodel, '%s'" % \ 

1032 compmodel, priority='error') 

1033 # modifymodel just collects info if outmodel==inmodel 

1034 model_vals = myutil.modifymodel(compmodel,compmodel, 

1035 "","","","","",-1, 

1036 flatimage=False) 

1037 model_cell = model_vals[1] 

1038 model_size = model_vals[2] 

1039 

1040 cell_tp = [qa.tos(model_cell[0]), qa.tos(model_cell[1])] 

1041 

1042 ##################################################### 

1043 

1044 imsize_tp = calc_imsize(mapsize=qimgsize_tp, cell=cell_tp) 

1045 

1046 msg(" ",priority=v_priority) 

1047 msg("-> The number of pixels needed to cover the map region: [%d, %d]" % \ 

1048 (imsize_tp[0], imsize_tp[1]), \ 

1049 priority=v_priority) 

1050 

1051 msg("Compare with interferometer image area and adopt the larger one:", \ 

1052 priority=v_priority) 

1053 # Compare with imsize of BL (note: imsize is an intArray) 

1054 imsize_bl = [] 

1055 if is_array_type(imsize) and imsize[0] > 0: 

1056 if len(imsize) > 1: 

1057 imsize_bl = imsize[0:2] 

1058 else: 

1059 imsize_bl = [imsize[0], imsize[0]] 

1060 

1061 if tp_kernel.upper() == 'SF': 

1062 if len(imsize_bl) > 0: 

1063 if cell != '': 

1064 # user defined cell size for INT 

1065 imarea = [qa.tos(qa.mul(cell, imsize[0])), 

1066 qa.tos(qa.mul(cell, imsize[1]))] 

1067 else: 

1068 # using model_cell for INT 

1069 imarea = [qa.tos(qa.mul(model_cell[0], imsize[0])), 

1070 qa.tos(qa.mul(model_cell[1], imsize[1]))] 

1071 tmpimsize = calc_imsize(mapsize=imarea, cell=cell_tp) 

1072 else: 

1073 msg("estimating imsize from input sky model.", \ 

1074 priority=v_priority) 

1075 tmpimsize = calc_imsize(mapsize=model_size, cell=cell_tp) 

1076 msg("-> TP imsize to cover interferometrer image area: [%d, %d]" % \ 

1077 (tmpimsize[0], tmpimsize[1]), \ 

1078 priority=v_priority) 

1079 elif len(imsize_bl) > 0: 

1080 # User has defined imsize 

1081 tmpimsize = imsize_bl 

1082 msg("-> Interferometer imsize (user defined): [%d, %d]" % \ 

1083 (tmpimsize[0], tmpimsize[1]), \ 

1084 priority=v_priority) 

1085 else: 

1086 # the same as input model (calculate from model_size) 

1087 msg("estimating imsize of interferometer from input sky model.", \ 

1088 priority=v_priority) 

1089 tmpimsize = calc_imsize(mapsize=model_size, cell=cell_tp) 

1090 msg("-> Estimated interferometer imsize (sky model): [%d, %d]" % \ 

1091 (tmpimsize[0], tmpimsize[1]), \ 

1092 priority=v_priority) 

1093 

1094 

1095 imsize_tp = [max(imsize_tp[0], tmpimsize[0]), \ 

1096 max(imsize_tp[1], tmpimsize[1])] 

1097 

1098 msg("The image pixel size of TP: [%d, %d]" % \ 

1099 (imsize_tp[0], imsize_tp[1]), \ 

1100 priority=v_priority) 

1101 

1102 # Generate TP image 

1103 msg(" ",priority=v_priority) 

1104 temp_out = fileroot+"/"+imagename_tp + '0' 

1105 task_param = dict(infiles=vis_tp, outfile=temp_out, imsize=imsize, 

1106 cell=cell_tp, phasecenter=model_refdir, 

1107 mode='channel', nchan= model_nchan) 

1108 if tp_kernel.upper() == 'SF': 

1109 msg("Generating TP image using 'SF' kernel.",\ 

1110 priority=v_priority) 

1111 # Parameters for sdimaging 

1112 task_param['gridfunction'] = 'sf' 

1113 task_param['convsupport'] = 6 

1114 else: 

1115 msg("Generating TP image using 'GJinc' kernel.",\ 

1116 priority=v_priority) 

1117 gfac = 2.52 # b in Mangum et al. (2007) 

1118 jfac = 1.55 # c in Mangum et al. (2007) 

1119 convfac = 1.8 # The conversion factor to get HWHM of kernel roughly equal to qhwhm 

1120 kernelfac = 0.7 # ratio of (kernel HWHM)/(TP pointingspacing) 

1121 #qfwhm = PB12 # FWHM of GJinc kernel. 

1122 #gwidth = qa.tos(qa.mul(gfac/3.,qfwhm)) 

1123 #jwidth = qa.tos(qa.mul(jfac/3.,qfwhm)) 

1124 qhwhm = qa.mul(qptgspc_tp, kernelfac) # hwhm of GJinc kernel 

1125 gwidth = qa.tos(qa.mul(qhwhm, convfac)) 

1126 jwidth = qa.tos(qa.mul(jfac/gfac/pl.log(2.),gwidth)) 

1127 #casalog.post("Kernel parameter: [qhwhm, gwidth, jwidth] = [%s, %s, %s]" % (qa.tos(qhwhm), gwidth, jwidth)) 

1128 # Parameters for sdimaging 

1129 task_param['gridfunction'] = 'gjinc' 

1130 task_param['gwidth'] = gwidth 

1131 task_param['jwidth'] = jwidth 

1132 

1133 casalog.post("saveinputs not available in casatasks, skipping saving sdimaging inputs", priority='WARN') 

1134 

1135 msg("Having set up the gridding parameters, the sdimaging task is called to actually creage the image:",priority=v_priority) 

1136 msg(get_taskstr('sdimaging', task_param), priority="info") 

1137 

1138 if not dryrun: 

1139 sdimaging(**task_param) 

1140 #del task_param 

1141 # TODO: scale TP image 

1142 

1143 # Scale TP image 

1144 if dryrun: #emulate beam calc in sdimaging 

1145 bu = sdbeamutil.TheoreticalBeam() 

1146 bu.set_antenna("12m", "0.75m") 

1147 bu.set_sampling([ptgspacing_tp, ptgspacing_tp], "0.0deg") 

1148 bu.set_image_param(task_param['cell'], model_center,task_param['gridfunction'], 

1149 task_param['convsupport'] if 'convsupport' in task_param else -1, 

1150 -1, 

1151 task_param['gwidth'] if 'gwidth' in task_param else -1, 

1152 task_param['jwidth'] if 'jwidth' in task_param else -1, 

1153 is_alma=True) 

1154 #bu.summary() 

1155 imbeam = bu.get_beamsize_image() 

1156 else: 

1157 ia.open(temp_out) 

1158 imbeam = ia.restoringbeam() 

1159 ia.close() 

1160 bmsize = imbeam['major'] 

1161 beam_area_ratio = qa.getvalue(qa.convert(bmsize, 'arcsec'))**2 \ 

1162 / qa.getvalue(qa.convert(PB12sim, 'arcsec'))**2 

1163 msg(" ",priority=v_priority) 

1164 msg("Scaling TP image intensity by beam area before and after gridding: %f" % beam_area_ratio) 

1165 task_param = dict(imagename=temp_out, mode='evalexpr', 

1166 expr=("IM0*%f" % (beam_area_ratio)), 

1167 outfile = fileroot+"/"+imagename_tp) 

1168 casalog.post("saveinputs not available in casatasks, skipping saving inmath inputs", priority='WARN') 

1169 

1170 msg(get_taskstr('immath', task_param), priority="info") 

1171 if not dryrun: 

1172 immath(**task_param) 

1173 del task_param 

1174 

1175 # Analyze TP image 

1176 tpskymodel = fileroot+"/"+pref_tp+".skymodel" 

1177 if components_only: 

1178 tpskymodel = fileroot+"/"+pref_tp+".compskymodel" 

1179 

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

1181 msg("Analyzing TP image", priority="info") 

1182 

1183 task_param = {} 

1184 task_param['project'] = project 

1185 task_param['image'] = False 

1186 task_param['imagename'] = fileroot+"/"+imagename_tp 

1187 task_param['skymodel'] = tpskymodel 

1188 task_param['analyze'] = True 

1189 task_param['showuv'] = False 

1190 task_param['showpsf'] = False 

1191 task_param['showconvolved'] = True 

1192 task_param['graphics'] = graphics 

1193 task_param['verbose'] = verbose 

1194 task_param['overwrite'] = overwrite 

1195 task_param['dryrun'] = dryrun 

1196 task_param['logfile'] = myutil.reportfile 

1197 

1198 msg(get_taskstr('simanalyze', task_param), priority="info") 

1199 

1200 try: 

1201 myutil.closereport() 

1202 simanalyze(**task_param) 

1203 del task_param 

1204 myutil.openreport() 

1205 except: 

1206 raise RuntimeError(simanaerr) 

1207 finally: 

1208 casalog.origin('simalma') 

1209 

1210 # Back up simanalyze.last file 

1211 if os.path.exists(fileroot+"/"+simana_file): 

1212 simana_new = imagename_tp.replace(".image",".simanalyze.last") 

1213 msg("Back up input parameter set to '%s'" % simana_new, \ 

1214 priority=v_priority) 

1215 shutil.move(fileroot+"/"+simana_file, fileroot+"/"+simana_new) 

1216 

1217 if not os.path.exists(fileroot+"/"+imagename_tp) and not dryrun: 

1218 msg("TP image '%s' is not found" \ 

1219 % imagename_tp, priority="error") 

1220 #modelimage = imagename_aca 

1221 

1222 

1223 

1224 

1225 

1226 

1227 ############################################################ 

1228 # Image each INT array separately 

1229 nconfig=len(antennalist) 

1230 

1231 for i in range(nconfig): 

1232 step += 1 

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

1234 msg("-"*60, origin="simalma", priority="info") 

1235 msg(("Step %d: imaging and analyzing " % step)+antennalist[i], origin="simalma", priority="info") 

1236 msg(" This step is optional, but useful to assess the result from just one configuration.",priority="warn",origin="simalma") 

1237 msg(" WARNING: The example clean shown here uses no mask, may diverge, and almost certainly is not optimal.",priority="warn",origin="simalma") 

1238 msg(" Users are HIGHLY recommended to use interactive clean masking (in simanalyze or directly in clean)",priority="warn",origin="simalma") 

1239 msg(" Auto-masking is under development for use in the ALMA pipeline and will be included here in a future release",priority="warn",origin="simalma") 

1240 msg("-"*60, origin="simalma", priority="info") 

1241 

1242 pref=get_data_prefix(antennalist[i], project) 

1243 if addnoise: 

1244 msname = pref+".noisy.ms" 

1245 imagename_int=pref+".noisy.image" 

1246 else: 

1247 msname= pref+".ms" 

1248 imagename_int=pref+".image" 

1249 

1250 if dryrun: 

1251 vis_int = msname 

1252 else: 

1253 if os.path.exists(fileroot+"/"+msname): 

1254 vis_int = msname 

1255 else: 

1256 msg("Could not find MS to image, '%s'" \ 

1257 % msname, origin="simalma", priority="error") 

1258 

1259# i think simanalye is fixed 20130826 

1260# # TMP fix: get correct skymodel file so that simanalyze picks it 

1261# if acaratio > 0: 

1262# if os.path.exists(tpskymodel): 

1263# shutil.move(tpskymodel,tpskymodel+".save") 

1264# else: 

1265# msg("TP skymodel '%s' is not found" \ 

1266# % tpskymodel, origin="simalma", priority="error") 

1267# 

1268# if os.path.exists(acaskymodel+".save"): 

1269# shutil.move(acaskymodel+".save",acaskymodel) 

1270# else: 

1271# msg("ACA skymodel '%s' is not found" \ 

1272# % acaskymodel+".save", origin="simalma", priority="error") 

1273 

1274 # dryrun requires feeding cell and imsize from here. 

1275 task_param = {} 

1276 task_param['project'] = project 

1277 task_param['image'] = image 

1278 task_param['vis'] = vis_int 

1279 task_param['modelimage'] = "" 

1280 if cell: 

1281 task_param['cell'] = cell 

1282 else: 

1283 task_param['cell'] = qa.tos(model_cell[0]) 

1284 if imsize[0]>0: 

1285 task_param['imsize'] = imsize 

1286 else: 

1287 task_param['imsize'] = int(qa.div(model_size[0],model_cell[0])['value']) 

1288 if len(imdirection)>0: 

1289 task_param['imdirection'] = imdirection 

1290 else: 

1291 task_param['imdirection'] = model_refdir 

1292 task_param['niter'] = niter 

1293 task_param['threshold'] = threshold 

1294 task_param['weighting'] = weighting 

1295 task_param['mask'] = [] 

1296 task_param['outertaper'] = [] 

1297 task_param['stokes'] = 'I' 

1298 task_param['analyze'] = True 

1299 task_param['graphics'] = graphics 

1300 task_param['verbose'] = verbose 

1301 task_param['overwrite'] = overwrite 

1302 task_param['dryrun'] = dryrun 

1303 task_param['logfile'] = myutil.reportfile 

1304 

1305 msg(get_taskstr('simanalyze', task_param), priority="info") 

1306 

1307 try: 

1308 myutil.closereport() 

1309 simanalyze(**task_param) 

1310 del task_param 

1311 myutil.openreport() 

1312 except: 

1313 raise RuntimeError(simanaerr) 

1314 finally: 

1315 casalog.origin('simalma') 

1316 

1317 

1318 

1319 if nconfig>1: 

1320 ############################################################ 

1321 # concat 

1322 step += 1 

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

1324 msg("-"*60, origin="simalma", priority="info") 

1325 msg("Step %d: concatenating interferometric visibilities." % step, origin="simalma", priority="info") 

1326 msg("-"*60, origin="simalma", priority="info") 

1327 

1328 weights=pl.zeros(nconfig)+1 

1329 z=pl.where(configtypes=='ACA')[0] 

1330 if len(z)>0: 

1331 weights[z]=weightratio_7_12 

1332 

1333 mslist=[] 

1334 for i in range(nconfig): 

1335 pref=get_data_prefix(antennalist[i],project) 

1336 if addnoise: 

1337 msname = fileroot + "/" + pref+".noisy.ms" 

1338 else: 

1339 msname = fileroot + "/" + pref+".ms" 

1340 mslist.append(msname) 

1341 

1342 msg("concat(vis="+str(mslist)+",concatvis="+concatname+".ms,visweightscale="+str(weights.tolist()),priority="info") 

1343 if not dryrun: 

1344 try: 

1345 concat(vis=mslist,concatvis=concatname+".ms",visweightscale=weights) 

1346 except: 

1347 raise RuntimeError(simanaerr) 

1348 finally: 

1349 casalog.origin('simalma') 

1350 

1351 

1352 

1353 

1354 ############################################################ 

1355 # Image ALMA-BL + ACA-7m  

1356 step += 1 

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

1358 msg("-"*60, origin="simalma", priority="info") 

1359 msg(("Step %d: imaging and analyzing " % step)+concatname+".ms", origin="simalma", priority="info") 

1360 msg(" WARNING: The example clean shown here uses no mask, may diverge, and almost certainly is not optimal.",priority="warn",origin="simalma") 

1361 msg(" Users are HIGHLY recommended to use interactive clean masking (in simanalyze or directly in clean)",priority="warn",origin="simalma") 

1362 msg(" Auto-masking is under development for use in the ALMA pipeline and will be included here in a future release",priority="warn",origin="simalma") 

1363 msg("-"*60, origin="simalma", priority="info") 

1364 

1365 

1366 if dryrun: 

1367 vis_int = concatname+".ms" 

1368 else: 

1369 if os.path.exists(fileroot+"/"+concatname+".ms"): 

1370 vis_int = fileroot+"/"+concatname+".ms" 

1371 elif os.path.exists(concatname+".ms"): 

1372 vis_int = concatname+".ms" 

1373 else: msg("Could not find MS to image, "+concatname+".ms", origin="simalma", priority="error") 

1374 

1375 task_param = {} 

1376 task_param['project'] = project 

1377 task_param['image'] = image 

1378 task_param['vis'] = vis_int 

1379 task_param['modelimage'] = "" 

1380 if cell: 

1381 task_param['cell'] = cell 

1382 else: 

1383 task_param['cell'] = model_cell 

1384 if imsize[0]>0: 

1385 task_param['imsize'] = imsize 

1386 else: 

1387 task_param['imsize'] = int(qa.div(model_size[0],model_cell[0])['value']) 

1388 if len(imdirection)>0: 

1389 task_param['imdirection'] = imdirection 

1390 else: 

1391 task_param['imdirection'] = model_refdir 

1392 task_param['niter'] = niter 

1393 task_param['threshold'] = threshold 

1394 task_param['weighting'] = weighting 

1395 task_param['mask'] = [] 

1396 task_param['outertaper'] = [] 

1397 task_param['stokes'] = 'I' 

1398 task_param['analyze'] = True 

1399 task_param['graphics'] = graphics 

1400 task_param['verbose'] = verbose 

1401 task_param['overwrite'] = overwrite 

1402 task_param['dryrun'] = dryrun 

1403 task_param['logfile'] = myutil.reportfile 

1404 

1405 msg(get_taskstr('simanalyze', task_param), priority="info") 

1406 imagename_int=os.path.basename(concatname.rstrip("/"))+".image" 

1407 

1408 try: 

1409 myutil.closereport() 

1410 simanalyze(**task_param) 

1411 del task_param 

1412 myutil.openreport() 

1413 except: 

1414 raise RuntimeError(simanaerr) 

1415 finally: 

1416 casalog.origin('simalma') 

1417 

1418 

1419 

1420 

1421 

1422 

1423 

1424 

1425 if tptime_min > 0: 

1426 ######################################################## 

1427 # Combine TP + INT image 

1428 step += 1 

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

1430 msg("-"*60, origin="simalma", priority="info") 

1431 msg("Step %d: combining a total power and synthesis image. " % step, origin="simalma", priority="info") 

1432 msg(" WARNING: feathering the two images is only one way to combine them. ",priority="warn",origin="simalma") 

1433 msg(" Using the total power image as a model in cleaning the interferometric visibilities may work better in some circumstances.",priority="warn",origin="simalma") 

1434 msg("-"*60, origin="simalma", priority="info") 

1435 

1436 if os.path.exists(fileroot+"/"+imagename_int) or dryrun: 

1437 highimage0 = fileroot+"/"+imagename_int 

1438 else: 

1439 msg("The synthesized image '%s' is not found" \ 

1440 % imagename_int, origin="simalma", priority="error") 

1441 if (not os.path.exists(fileroot+"/"+imagename_tp)) and (not dryrun): 

1442 msg("ACA is requested but total power image '%s' is not found" \ 

1443 % imagename_tp, origin="simalma", priority="error") 

1444 #lowimage = fileroot+"/"+imagename_tp 

1445 

1446 # Need to manipulate TP image here 

1447 outimage0 = fileroot+"/" + combimage+"0" 

1448 outimage = fileroot+"/" + combimage 

1449 # CAS-13369 -- imclean call in simanalyze is now imtclean 

1450 #pbcov = highimage0.rstrip("image") + "flux.pbcoverage" 

1451 pbcov = highimage0.rstrip("image") + "pb" 

1452 regridimg = fileroot + "/" + imagename_tp + ".regrid" 

1453 scaledimg = fileroot + "/" + imagename_tp + ".pbscaled" 

1454 lowimage = scaledimg 

1455 

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

1457 msg("Regrid total power image to interferometric image grid:",priority=v_priority) 

1458 # regrid TP image 

1459 msg("inttemplate = imregrid(imagename = '"+highimage0+"', template='get')",priority=v_priority) 

1460 if not dryrun: 

1461 inttemplate = imregrid(imagename = highimage0, template='get') 

1462 msg("imregrid(imagename = '"+fileroot+"/"+imagename_tp+ 

1463 "',interpolation='cubic',template = inttemplate, output = '"+regridimg+"')",priority="info") 

1464 if not dryrun: 

1465 imregrid(imagename = fileroot+"/"+imagename_tp, 

1466 interpolation="cubic", 

1467 template = inttemplate, output = regridimg) 

1468 # multiply SD image with INT PB coverage 

1469 if not os.path.exists(pbcov): 

1470 msg("The flux image '%s' is not found" \ 

1471 % pbcov, origin="simalma", priority="error") 

1472 

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

1474 msg("Multiply total power image by interferometric sensitivity map:",priority=v_priority) 

1475# msg("immath(imagename=['"+regridimg+"','"+pbcov+"'],expr='IM1*IM0',outfile='"+scaledimg+"')",priority="info") 

1476# if not dryrun: 

1477# immath(imagename=[regridimg, pbcov], 

1478# expr='IM1*IM0',outfile=scaledimg) 

1479# 

1480# msg(" ",priority="info")  

1481# msg("Set total power image beam and brightness unit:",priority=v_priority)  

1482# msg("ia.open('"+fileroot+"/"+imagename_tp+"')",priority="info") 

1483# msg("beam_tp = ia.restoringbeam()",priority="info") 

1484# msg("bunit_tp = ia.brightnessunit()",priority="info") 

1485# msg("ia.close()",priority="info") 

1486# msg("ia.open('"+scaledimg+"')",priority="info") 

1487# msg("ia.setrestoringbeam(beam=beam_tp)",priority="info") 

1488# msg("ia.setbrightnessunit(bunit_tp)",priority="info") 

1489# msg("ia.close()",priority="info") 

1490# 

1491# if not dryrun: 

1492# pdb.set_trace() 

1493# # restore TP beam and brightness unit 

1494# ia.open(fileroot+"/"+imagename_tp) 

1495# beam_tp = ia.restoringbeam() 

1496# bunit_tp = ia.brightnessunit() 

1497# ia.close() 

1498# ia.open(scaledimg) 

1499# ia.setrestoringbeam(beam=beam_tp) 

1500# ia.setbrightnessunit(bunit_tp) 

1501# ia.close() 

1502 

1503 msg("impbcor('"+regridimg+"', '"+pbcov+"', outfile='"+scaledimg+"',mode='multiply')",priority="info") 

1504 if not dryrun: 

1505 impbcor(regridimg, pbcov, outfile=scaledimg,mode='multiply') 

1506 

1507 # de-pbcor the INT image 

1508 # not needed now that imtclean is used and .image from tclean is flat sky by default, see CAS-13370 

1509 #highimage = fileroot+"/"+imagename_int+".pbscaled" 

1510 #immath(imagename=[highimage0, pbcov], 

1511 # expr='IM1/IM0',outfile=highimage)  

1512# msg(" ",priority="info") 

1513# msg("Multiply interferometric image by sensitivity map (un-pbcor):",priority="info") 

1514# msg("immath(imagename=['"+highimage0+"','"+pbcov+"'],expr='IM1*IM0',outfile='"+highimage+"')",priority="info") 

1515# msg("Restore interferometric beam and brightness unit:",priority="info") 

1516# msg("ia.open('"+highimage0+"')",priority="info") 

1517# msg("beam_int = ia.restoringbeam()",priority="info") 

1518# msg("bunit_int = ia.brightnessunit()",priority="info") 

1519# msg("ia.close()",priority="info") 

1520# msg("ia.open('"+highimage+"')",priority="info") 

1521# msg("ia.setrestoringbeam(beam=beam_int)",priority="info") 

1522# msg("ia.setbrightnessunit(bunit_int)",priority="info") 

1523# msg("ia.close()",priority="info") 

1524# 

1525# if not dryrun: 

1526# immath(imagename=[highimage0, pbcov], 

1527# expr='IM1*IM0',outfile=highimage) 

1528# # restore INT beam and brightness unit 

1529# ia.open(highimage0) 

1530# beam_int = ia.restoringbeam() 

1531# bunit_int = ia.brightnessunit() 

1532# ia.close() 

1533# ia.open(highimage) 

1534# ia.setrestoringbeam(beam=beam_int) 

1535# ia.setbrightnessunit(bunit_int) 

1536# ia.close() 

1537# 

1538# msg("impbcor('"+highimage0+"', '"+pbcov+"', outfile='"+highimage+"',mode='multiply')",priority="info") 

1539# if not dryrun: 

1540# impbcor(highimage0, pbcov, outfile=highimage,mode='multiply') 

1541 

1542 

1543 

1544 

1545 

1546 # Feathering 

1547 task_param = {} 

1548 task_param['imagename'] = outimage0 

1549 task_param['highres'] = highimage0 

1550 task_param['lowres'] = lowimage 

1551 

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

1553 msg(get_taskstr('feather', task_param), priority="info") 

1554 try: 

1555 casalog.post("saveinputs not available in casatasks, skipping saving feather inputs", priority='WARN') 

1556 

1557 if not dryrun: 

1558 feather(**task_param) 

1559 del task_param 

1560 

1561 # This seems not necessary anymore. 

1562 ## transfer mask - feather should really do this 

1563 #ia.open(outimage0) 

1564 #ia.maskhandler('copy',[highimage+":mask0",'mask0']) 

1565 #ia.maskhandler('set','mask0') 

1566 #ia.done() 

1567 except Exception as exc: 

1568 raise Exception("simalma caught an exception in task feather: {}". 

1569 format(exc)) 

1570 finally: 

1571 if not dryrun: shutil.rmtree(regridimg) 

1572 #shutil.rmtree(scaledimg) 

1573 casalog.origin('simalma') 

1574 

1575 

1576 # re-pbcor the result 

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

1578 msg("Re-apply the primary beam correction to the feathered result:",priority=v_priority) 

1579# msg("immath(imagename=['"+outimage0+"','"+pbcov+"'],expr='IM0/IM1',outfile='"+outimage+"')",priority="info") 

1580# if not dryrun: 

1581# immath(imagename=[outimage0, pbcov], 

1582# expr='IM0/IM1',outfile=outimage) 

1583 

1584 msg("impbcor('"+outimage0+"', '"+pbcov+"', outfile='"+outimage+"')",priority="info") 

1585 if not dryrun: 

1586 impbcor(outimage0, pbcov, outfile=outimage) 

1587 

1588 

1589 

1590 

1591 

1592 ######################################################## 

1593 # Generate Summary Plot 

1594 grscreen = False 

1595 grfile = False 

1596 if not dryrun: 

1597 if graphics == "both": 

1598 grscreen = True 

1599 grfile = True 

1600 if graphics == "screen": 

1601 grscreen = True 

1602 if graphics == "file": 

1603 grfile = True 

1604 

1605 if grscreen or grfile: 

1606 if grfile: 

1607 file = fileroot + "/" + project + ".combine.png" 

1608 else: 

1609 file = "" 

1610 

1611 # check for image pathes 

1612 if os.path.exists(skymodel): 

1613 flatsky = pref_bl + ".skymodel.flat" 

1614 else: 

1615 flatsky = pref_bl + ".compskymodel.flat" 

1616 if not os.path.exists(fileroot+"/"+flatsky): 

1617 raise RuntimeError("Coud not find a skymodel image '%s'" % flatsky) 

1618 

1619 if not os.path.exists(fileroot+"/"+combimage): 

1620 raise RuntimeError("Coud not find the combined image '%s'" % combimage) 

1621 

1622 if not os.path.exists(fileroot+"/"+imagename_int): 

1623 raise RuntimeError("Coud not find the synthesized image '%s'" % imagename_int) 

1624 

1625 if not os.path.exists(fileroot+"/"+imagename_tp): 

1626 raise RuntimeError("Coud not find the total power image '%s'" % (imagename_tp)) 

1627 # Now the actual plotting 

1628 disprange = None 

1629 myutil.newfig(multi=[2,2,1],show=grscreen) 

1630 # skymodel 

1631 #discard = myutil.statim(fileroot+"/"+flatsky,disprange=disprange) 

1632 

1633 #disprange = [] 

1634 # generate flat synthesized (7m+12m) image 

1635 flatint = fileroot + "/" + imagename_int + ".flat" 

1636 myutil.flatimage(fileroot+"/"+imagename_int,verbose=verbose) 

1637 if not os.path.exists(flatint): 

1638 raise RuntimeError("Failed to generate '%s'" % (flatint)) 

1639 

1640 # generate convolved sky model image 

1641 myutil.convimage(fileroot+"/"+flatsky, flatint) 

1642 discard = myutil.statim(fileroot+"/"+flatsky+".regrid.conv",disprange=disprange) 

1643 shutil.rmtree(fileroot+"/"+flatsky+".regrid") 

1644 shutil.rmtree(fileroot+"/"+flatsky+".regrid.conv") 

1645 

1646 # total power image 

1647 flattp = fileroot + "/" + imagename_tp + ".flat" 

1648 myutil.flatimage(fileroot+"/"+imagename_tp,verbose=verbose) 

1649 #flattp = scaledimg + ".flat" 

1650 #myutil.flatimage(scaledimg,verbose=verbose) 

1651 if not os.path.exists(flattp): 

1652 raise RuntimeError("Failed to generate '%s'" % (flattp)) 

1653 myutil.nextfig() 

1654 discard = myutil.statim(flattp,disprange=disprange) 

1655 shutil.rmtree(flattp) 

1656 

1657 #disprange = [] 

1658 # display flat synthesized (7m+12m) image 

1659 myutil.nextfig() 

1660 discard = myutil.statim(flatint,disprange=disprange) 

1661 shutil.rmtree(flatint) 

1662 

1663 # combined image 

1664 flatcomb = fileroot + "/" + combimage + ".flat" 

1665 myutil.flatimage(fileroot+"/"+combimage,verbose=verbose) 

1666 if not os.path.exists(flatcomb): 

1667 raise RuntimeError("Failed to generate '%s'" % (flatcomb)) 

1668 myutil.nextfig() 

1669 discard = myutil.statim(flatcomb,disprange=disprange) 

1670 myutil.endfig(show=grscreen,filename=file) 

1671 shutil.rmtree(flatcomb) 

1672 

1673 if myutil.isreport(): 

1674 myutil.closereport() 

1675 finalize_tools() 

1676 

1677 finally: 

1678 finalize_tools() 

1679 if myutil.isreport(): 

1680 myutil.closereport() 

1681 

1682 

1683def finalize_tools(): 

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

1685 sm.close() 

1686 #cl.close() # crashes casa 

1687 

1688def get_data_prefix(cfgname, project=""): 

1689 if str.upper(cfgname[0:4]) == "ALMA": 

1690 foo=cfgname.replace(';','_') 

1691 else: 

1692 foo = cfgname 

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

1694 sfoo=foo.split('/') 

1695 if len(sfoo)>1: foo=sfoo[-1] 

1696 return project+"."+foo 

1697 

1698def calc_imsize(mapsize=None, cell=None): 

1699 if mapsize == None: 

1700 raise ValueError("mapsize is not defined") 

1701 if cell == None: 

1702 raise ValueError("cell is not defined") 

1703 # get a list of cell size 

1704 if is_array_type(cell): 

1705 if len(cell) < 2: 

1706 cell = [cell[0], cell[0]] 

1707 else: 

1708 cell = [cell, cell] 

1709 

1710 for qval in cell: 

1711 if not qa.compare(qval, "deg"): 

1712 raise TypeError("cell should be an angular size") 

1713 

1714 qcellx = qa.quantity(cell[0]) 

1715 qcelly = qa.quantity(cell[1]) 

1716 

1717 # get a list of map size 

1718 if is_array_type(mapsize): 

1719 if len(mapsize) < 2: 

1720 mapsize = [mapsize[0], mapsize[0]] 

1721 else: 

1722 mapsize = [mapsize, mapsize] 

1723 

1724 for qval in mapsize: 

1725 if not qa.compare(qval, "deg"): 

1726 raise TypeError("mapsize should be an angular size") 

1727 

1728 vsizex = qa.convert(mapsize[0], qcellx['unit'])['value'] 

1729 vsizey = qa.convert(mapsize[1], qcelly['unit'])['value'] 

1730 

1731 # Calculate the number of pixels to cover the map size 

1732 npixx = int(pl.ceil(vsizex/qcellx['value'])) 

1733 npixy = int(pl.ceil(vsizey/qcelly['value'])) 

1734 

1735 return [npixx, npixy] 

1736 

1737