|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 import os 00002 import shutil 00003 import re 00004 import glob 00005 import numpy 00006 00007 from taskinit import * 00008 from simutil import * 00009 from simobserve import simobserve 00010 from simanalyze import simanalyze 00011 00012 def simalma( 00013 project=None, 00014 skymodel=None, inbright=None, indirection=None, incell=None, 00015 incenter=None, inwidth=None, 00016 complist=None, compwidth=None, 00017 ######## 00018 setpointings=None, 00019 ptgfile=None, 00020 integration=None, direction=None, mapsize=None, 00021 #maptype=None, 00022 #pointingspacing=None, 00023 #caldirection=None, calflux=None, 00024 #observe=None, 00025 #refdate=None, 00026 antennalist=None, 00027 hourangle=None, 00028 totaltime=None, 00029 #sdantlist=None, sdant=None, 00030 ### 00031 acaratio = None, 00032 acaconfig = None, 00033 ### 00034 #thermalnoise=None, 00035 #user_pwv=None, t_ground=None, t_sky=None, tau0=None, 00036 pwv=None, 00037 #seed=None, 00038 #leakage=None, 00039 image=None, 00040 #vis=None, modelimage=None, 00041 imsize=None, imdirection=None,cell=None, 00042 niter=None, threshold=None, 00043 #weighting=None, mask=None, outertaper=None, stokes=None, 00044 #analyze=None, 00045 #showarray=None, showuv=None, showpsf=None, showmodel=None, 00046 #showconvolved=None, showclean=None, showresidual=None, showdifference=None, 00047 #showfidelity=None, 00048 graphics=None, 00049 verbose=None, 00050 overwrite=None, 00051 async=False): 00052 00053 # Collect a list of parameter values to save inputs 00054 in_params = locals() 00055 00056 try: 00057 casalog.origin('simalma') 00058 if verbose: 00059 casalog.filter(level="DEBUG2") 00060 v_priority = "WARN" 00061 else: 00062 v_priority = "INFO" 00063 00064 simobserr = "simalma caught an exception in task simobserve" 00065 simanaerr = "simalma caught an exception in task simanalyze" 00066 00067 # Get globals to call saveinputs() 00068 a = inspect.stack() 00069 stacklevel = 0 00070 for k in range(len(a)): 00071 if (string.find(a[k][1], 'ipython console') > 0): 00072 stacklevel = k 00073 myf = sys._getframe(stacklevel).f_globals 00074 00075 # Parameter verifications 00076 # antennalist should be one of ALMA config. 00077 if antennalist.upper().find("ALMA") < 0: 00078 raise ValueError, "antennalist should be ALMA configuration in simalma." 00079 00080 # Save outputs in a directory called "project" 00081 fileroot = project 00082 # simalma is not supposed to run multiple times. 00083 if os.path.exists(fileroot): 00084 infomsg = "Project directory, '%s', already exists." % fileroot 00085 if overwrite: 00086 casalog.post(infomsg) 00087 casalog.post("Removing old project directory '%s'" % fileroot) 00088 shutil.rmtree(fileroot) 00089 else: 00090 raise Exception, infomsg 00091 00092 if not os.path.exists(fileroot): 00093 os.mkdir(fileroot) 00094 00095 # Save input parameters of simalma 00096 saveinputs = myf['saveinputs'] 00097 saveinputs('simalma',fileroot+"/"+project+".simalma.last", 00098 myparams=in_params) 00099 00100 # Create the utility object 00101 myutil = simutil(direction) 00102 if verbose: myutil.verbose = True 00103 msg = myutil.msg 00104 00105 # Preset parameters in simalma 00106 nyquist = 0.48113 ## Nyquist 00107 maptype_bl = 'ALMA' 00108 maptype_aca = 'ALMA' 00109 maptype_tp = 'square' 00110 pbgridratio_tp = 0.36 00111 refdate = '2013/05/21' 00112 caldirection = "" 00113 calflux = "0Jy" 00114 tpantid = 0 00115 t_ground = 270. 00116 if pwv > 0: 00117 thermalnoise = "tsys-atm" 00118 else: 00119 thermalnoise = "" 00120 leakage = 0. 00121 weighting = "briggs" 00122 00123 # format mapsize 00124 if not is_array_type(mapsize): 00125 mapsize = [mapsize, mapsize] 00126 elif len(mapsize) < 2: 00127 mapsize = [mapsize[0], mapsize[0]] 00128 00129 # Operation flags 00130 addnoise = (thermalnoise != '') 00131 # Rectangle setup mode 00132 multiptg = (not setpointings) \ 00133 or (is_array_type(direction) and len(direction) > 1) 00134 rectmode = (not multiptg) 00135 00136 # Use full model image as a mapsize = ["", ""] 00137 fullsize = (len(mapsize[0]) == 0) or (len(mapsize[1]) == 0) 00138 00139 # Test for cycle 1 00140 #q = re.compile('.*CYCLE.?1.?;(.*)') 00141 q = re.compile('.*CYCLE.?1.*') 00142 isC1 = q.match(antennalist.upper()) 00143 if isC1: 00144 msg("Cycle-1 ALMA simulation", origin="simalma", priority="warn") 00145 else: 00146 msg("Full Science ALMA simulation", origin="simalma", priority="warn") 00147 00148 # antennalist of ACA and TP 00149 antlist_tp = "aca.tp.cfg" 00150 00151 if not acaconfig: # default 00152 if isC1: 00153 acaconfig = "cycle1" 00154 else: 00155 acaconfig = "i" 00156 00157 warnmsg = "You are simulating Cycle1 observation but requested Full Science configuration for ACA. Assuming you know what you want and using config '%s'." 00158 if acaconfig.upper().startswith("I"): 00159 antlist_aca = "aca.i.cfg" 00160 if isC1: 00161 msg(warnmsg % antlist_aca, origin="simalma", priority="warn") 00162 elif acaconfig.upper().startswith("N"): 00163 antlist_aca = "aca.ns.cfg" 00164 if isC1: 00165 msg(warnmsg % antlist_aca, origin="simalma", priority="warn") 00166 else: 00167 antlist_aca = "aca_cycle1.cfg" 00168 00169 # Resolve prefixes of simulation data (as defined in 00170 # simobserve and simanalyze) 00171 pref_bl = get_data_prefix(antennalist, project) 00172 pref_aca = get_data_prefix(antlist_aca, project) 00173 pref_tp = get_data_prefix(antlist_tp, project) 00174 # Resolve output names (as defined in simobserve and simanalyze) 00175 ptgfile_bl = fileroot+"/"+pref_bl+".ptg.txt" 00176 if addnoise: 00177 msname_bl = pref_bl+".noisy.ms" 00178 msname_aca = pref_aca+".noisy.ms" 00179 msname_tp = pref_tp+".noisy.sd.ms" 00180 imagename_aca = pref_aca+".noisy.image" 00181 else: 00182 msname_bl = pref_bl+".ms" 00183 msname_aca = pref_aca+".ms" 00184 msname_tp = pref_tp+".sd.ms" 00185 imagename_aca = pref_aca+".image" 00186 00187 simana_file = project+".simanalyze.last" 00188 00189 # Either skymodel or complist should exists 00190 if is_array_type(skymodel): 00191 skymodel = skymodel[0] 00192 skymodel = skymodel.replace('$project',pref_bl) 00193 00194 if is_array_type(complist): 00195 complist = complist[0] 00196 00197 if((not os.path.exists(skymodel)) and (not os.path.exists(complist))): 00198 raise ValueError, "No sky input found. At least one of skymodel or complist must be set." 00199 00200 ########################### 00201 # Get model_size and model_center 00202 if os.path.exists(skymodel): 00203 outmodel = fileroot+"/"+project+"temp.skymodel" 00204 model_vals = myutil.modifymodel(skymodel, outmodel, inbright, 00205 indirection, incell, incenter, 00206 inwidth, -1, False) 00207 shutil.rmtree(outmodel) 00208 model_cell = model_vals[1] 00209 model_size = model_vals[2] 00210 model_center = model_vals[4] 00211 del model_vals 00212 else: 00213 compdirs = [] 00214 cl.open(complist) 00215 for i in range(cl.length()): 00216 compdirs.append(myutil.dir_m2s(cl.getrefdir(i))) 00217 00218 model_refdir, coffs = myutil.average_direction(compdirs) 00219 model_center = cl.getspectrum(0)['frequency']['m0'] 00220 cmax = 0.0014 # ~5 arcsec 00221 for i in range(coffs.shape[1]): 00222 xc = numpy.absolute(coffs[0,i]) # offsets in deg 00223 yc = numpy.absolute(coffs[1,i]) 00224 cmax = max(cmax, xc) 00225 cmax = max(cmax, yc) 00226 model_size = ["%fdeg" % (2*cmax), "%fdeg" % (2*cmax)] 00227 cl.done() 00228 model_cell = None 00229 del compdirs, model_refdir, coffs, xc, yc, cmax 00230 00231 # Calculate 12-m PB 00232 Dant = 12. 00233 pbval = 0.2997924 / qa.convert(qa.quantity(model_center),'GHz')['value'] \ 00234 / Dant * 3600. * 180 / numpy.pi # arcsec 00235 PB12 = qa.quantity(pbval*1.2, "arcsec") 00236 msg("PB size: %s" % (qa.tos(PB12)), origin="simalma", priority='DEBUG2') 00237 pointingspacing = str(nyquist)+"PB" 00238 00239 ############################################################ 00240 # ALMA-BL simulation 00241 step = 1 00242 msg("="*60, origin="simalma", priority="warn") 00243 msg(" Step %d: simulating ALMA 12-m array observation." % step, origin="simalma", priority="warn") 00244 msg("="*60, origin="simalma", priority="warn") 00245 00246 obsmode_int = 'int' 00247 # BL mapsize should be 1 PB smaller than skymodel when using ACA 00248 #if acaratio > 0 and rectmode and fullsize: 00249 # mapx = qa.sub(model_size[0], PB12ot) 00250 # mapy = qa.sub(model_size[1], PB12ot) 00251 # mapsize_bl = [qa.tos(mapx), qa.tos(mapy)] 00252 #else: 00253 # mapsize_bl = mapsize 00254 mapsize_bl = mapsize 00255 00256 taskstr = "simobserve(project='"+project+"', skymodel='"+skymodel+"', inbright='"+inbright+"', indirection='"+indirection+"', incell='"+incell+"', incenter='"+incenter+"', inwidth='"+inwidth+"', complist='"+complist+"', compwidth='"+compwidth+"', setpointings="+str(setpointings)+", ptgfile='"+ptgfile+"', integration='"+integration+"', direction='"+str(direction)+"', mapsize="+str(mapsize_bl)+", maptype='"+maptype_bl+"', pointingspacing='"+pointingspacing+"', caldirection='"+caldirection+"', calflux='"+calflux+"', obsmode='"+obsmode_int+"', refdate='"+refdate+"', hourangle='"+hourangle+"', totaltime='"+totaltime+"', antennalist='"+antennalist+"', sdantlist='', sdant="+str(0)+", thermalnoise='"+thermalnoise+"', user_pwv="+str(pwv)+", t_ground="+str(t_ground)+", leakage="+str(leakage)+", graphics='"+graphics+"', verbose="+str(verbose)+", overwrite="+str(overwrite)+")" 00257 msg("Executing: "+taskstr, origin="simalma", priority=v_priority) 00258 00259 try: 00260 simobserve(project=project, 00261 skymodel=skymodel, inbright=inbright, 00262 indirection=indirection, incell=incell, 00263 incenter=incenter, inwidth=inwidth, 00264 complist=complist, compwidth=compwidth, 00265 setpointings=setpointings, ptgfile=ptgfile, 00266 integration=integration, 00267 direction=direction, mapsize=mapsize_bl, 00268 maptype=maptype_bl, pointingspacing=pointingspacing, 00269 caldirection=caldirection, calflux=calflux, 00270 obsmode=obsmode_int, refdate=refdate, 00271 hourangle=hourangle, totaltime=totaltime, 00272 antennalist=antennalist, 00273 sdantlist="", sdant=0, 00274 thermalnoise=thermalnoise, user_pwv=pwv, 00275 t_ground=t_ground, #t_sky=t_sky, tau0=tau0, seed=seed, 00276 leakage=leakage, 00277 graphics=graphics, verbose=verbose, overwrite=overwrite)#, 00278 #async=False) 00279 except: 00280 raise Exception, simobserr 00281 finally: 00282 casalog.origin('simalma') 00283 00284 qimgsize_tp = None 00285 00286 if acaratio > 0: 00287 ######################################################## 00288 # ACA-7m simulation 00289 step += 1 00290 msg("="*60, origin="simalma", priority="warn") 00291 msg(" Step %d: simulating ACA 7-m array observation." % step, origin="simalma", priority="warn") 00292 msg("="*60, origin="simalma", priority="warn") 00293 00294 # Calculate total time for ACA and TP 00295 if qa.compare(totaltime,'s'): 00296 tottime_aca = qa.tos(qa.mul(totaltime, acaratio)) 00297 else: # number of visit (calc ACA tottime) 00298 npts, pointings, time = myutil.read_pointings(ptgfile_bl) 00299 if len(time) == npts: 00300 tottime_aca = qa.tos(qa.quantity(sum(time)*acaratio,'s')) 00301 else: 00302 tottime_aca = qa.tos(qa.mul(integration, ntps*acaratio)) 00303 del npts, pointings, time 00304 msg("Total observing time of ACA = %s" % tottime_aca, origin="simalma", priority='warn') 00305 00306 # Same pointings as BL 00307 #ptgfile_aca = ptgfile_bl 00308 00309 taskstr = "simobserve(project='"+project+"', skymodel='"+skymodel+"', inbright='"+inbright+"', indirection='"+indirection+"', incell='"+incell+"', incenter='"+incenter+"', inwidth='"+inwidth+"', complist='"+complist+"', compwidth='"+compwidth+"', setpointings="+str(setpointings)+", ptgfile='"+ptgfile+"', integration='"+integration+"', direction='"+str(direction)+"', mapsize="+str(mapsize_bl)+", maptype='"+maptype_aca+"', pointingspacing='"+pointingspacing+"', caldirection='"+caldirection+"', calflux='"+calflux+"', obsmode='"+obsmode_int+"', refdate='"+refdate+"', hourangle='"+hourangle+"', totaltime='"+tottime_aca+"', antennalist='"+antlist_aca+"', sdantlist='', sdant="+str(0)+", thermalnoise='"+thermalnoise+"', user_pwv="+str(pwv)+", t_ground="+str(t_ground)+", leakage="+str(leakage)+", graphics='"+graphics+"', verbose="+str(verbose)+", overwrite="+str(overwrite)+")" 00310 msg("Executing: "+taskstr, origin="simalma", priority=v_priority) 00311 00312 try: 00313 simobserve(project=project, 00314 skymodel=skymodel, inbright=inbright, 00315 indirection=indirection, incell=incell, 00316 incenter=incenter, inwidth=inwidth, 00317 complist=complist, compwidth=compwidth, 00318 #setpointings=False, ptgfile=ptgfile_aca, 00319 setpointings=setpointings, ptgfile=ptgfile, 00320 integration=integration, 00321 direction=direction, mapsize=mapsize_bl, 00322 maptype=maptype_aca, pointingspacing=pointingspacing, 00323 caldirection=caldirection, calflux=calflux, 00324 obsmode=obsmode_int, refdate=refdate, 00325 hourangle=hourangle, totaltime=tottime_aca, 00326 antennalist=antlist_aca, 00327 sdantlist="", sdant=0, 00328 thermalnoise=thermalnoise, user_pwv=pwv, 00329 t_ground=t_ground, #t_sky=t_sky, tau0=tau0, seed=seed, 00330 leakage=leakage, 00331 graphics=graphics, verbose=verbose, overwrite=overwrite)#, 00332 #async=False) 00333 except: 00334 raise Exception, simobserr 00335 finally: 00336 casalog.origin('simalma') 00337 00338 ######################################################## 00339 # ACA-TP simulation 00340 step += 1 00341 msg("="*60, origin="simalma", priority="warn") 00342 msg(" Step %d: simulating ACA total power observation." % step, origin="simalma", priority="warn") 00343 msg("="*60, origin="simalma", priority="warn") 00344 obsmode_sd = "sd" 00345 00346 # Resolve mapsize of TP. Add 1 PB to pointing extent of BL 00347 if rectmode: 00348 # Add 1PB to mapsize 00349 if fullsize: 00350 mapx = qa.add(PB12ot,model_size[0]) # in the unit same as PB 00351 mapy = qa.add(PB12ot,model_size[1]) # in the unit same as PB 00352 mapsize_tp = [qa.tos(mapx), qa.tos(mapy)] 00353 msg("Full skymodel mapped by ALMA 12-m and ACA 7-m arrays. The total power antenna observes 1PB larger extent.", origin="simalma", priority='warn') 00354 else: 00355 # mapsize is defined. Add 1 PB to mapsize. 00356 mapx = qa.add(qa.quantity(mapsize[0]), PB12ot) 00357 mapy = qa.add(qa.quantity(mapsize[1]), PB12ot) 00358 mapsize_tp = [qa.tos(mapx), qa.tos(mapy)] 00359 msg("A part of skymodel mapped by ALMA 12-m and ACA 7-m arrays. The total power antenna observes 1PB larger extent.", origin="simalma", priority='warn') 00360 else: 00361 # multi-pointing mode 00362 npts, pointings, time = myutil.read_pointings(ptgfile_bl) 00363 center, offsets = myutil.average_direction(pointings) 00364 del time 00365 qx = qa.quantity(max(offsets[0])-min(offsets[0]),"deg") 00366 qy = qa.quantity(max(offsets[1])-min(offsets[1]),"deg") 00367 mapx = qa.add(qa.mul(PB12ot,2.),qx) # in the unit same as PB 00368 mapy = qa.add(qa.mul(PB12ot,2.),qy) # in the unit same as PB 00369 mapsize_tp = [qa.tos(mapx), qa.tos(mapy)] 00370 # number of pointings to map vicinity of each pointings 00371 npts_multi = npts * int(2./pbgridratio_tp)**2 00372 msg("Number of pointings to map vicinity of each direction = %d" % npts_multi, origin="simalma", priority="DEBUG2") 00373 00374 # back-up imsize for TP image generation 00375 qimgsize_tp = [mapx, mapy] 00376 00377 grid_tp = qa.mul(PB12ot, pbgridratio_tp) 00378 pbunit = PB12ot['unit'] 00379 # number of pointings to map pointing region 00380 npts_rect = int(qa.convert(mapx, pbunit)['value'] \ 00381 / qa.convert(grid_tp, pbunit)['value']) \ 00382 * int(qa.convert(mapy, pbunit)['value'] \ 00383 / qa.convert(grid_tp, pbunit)['value']) 00384 msg("Number of pointings to map a rect region = %d" % npts_rect, origin="simalma", priority="DEBUG2") 00385 00386 if rectmode: 00387 dir_tp = direction 00388 npts_tp = npts_rect 00389 msg("Rectangle mode: The total power antenna observes 1PB larger region compared to ALMA 12-m and ACA 7-m arrays", origin="simalma", priority='warn') 00390 else: 00391 if npts_multi < npts_rect: 00392 # Map +-1PB extent of each direction 00393 # need to get a list of pointings 00394 dir_tp = [] 00395 locsize = qa.mul(2, PB12ot) 00396 for dir in pointings: 00397 dir_tp += myutil.calc_pointings2(qa.tos(grid_tp), 00398 qa.tos(locsize), 00399 "square", dir) 00400 00401 mapsize_tp = ["", ""] 00402 npts_tp = npts_multi 00403 msg("Multi-pointing mode: The total power antenna observes +-1PB of each point", origin="simalma", priority='warn') 00404 else: 00405 # Map a region that covers all directions 00406 dir_tp = center 00407 npts_tp = npts_rect 00408 msg("Multi-pointing mode: The total power antenna maps a region that covers all pointings", origin="simalma", priority='warn') 00409 msg("- Center of poinings: %s" % center, origin="simalma", priority='warn') 00410 msg("- Map size: [%s, %s]" % (mapsize_tp[0], mapsize_tp[1]), origin="simalma", priority='warn') 00411 00412 ptgspacing_tp = str(pbgridratio_tp*PB12ot['value']/PB12sim['value'])+"PB" 00413 00414 # Scale integration time of TP (assure >= 1 visit per direction) 00415 tottime_tp = tottime_aca 00416 integration_tp = integration 00417 ndump = int(qa.convert(tottime_tp, 's')['value'] 00418 / qa.convert(integration, 's')['value']) 00419 msg("Max number of dump in %s (integration %s): %d" % \ 00420 (tottime_tp, integration, ndump), origin="simalma", \ 00421 priority="DEBUG2") 00422 00423 if ndump < npts_tp: 00424 t_scale = float(ndump)/float(npts_tp) 00425 integration_tp = qa.tos(qa.mul(integration, t_scale)) 00426 msg("Integration time is scaled to cover all pointings in observation time.", origin="simalma", priority='warn') 00427 msg("- Scaled total power integration time: %s" % integration_tp, origin="simalma", priority='warn') 00428 ## Sometimes necessary to avoid the effect of round-off error 00429 #iunit = qa.quantity(integration_tp)['unit'] 00430 #intsec = qa.convert(integration_tp,"s") 00431 #totsec = intsec['value']*npts_tp#+0.000000001) 00432 ##tottime_tp = qa.tos(qa.convert(qa.quantity(totsec, "s"), iunit)) 00433 #tottime_tp = qa.tos(qa.quantity(totsec, "s")) 00434 00435 taskstr = "simobserve(project='"+project+"', skymodel='"+skymodel+"', inbright='"+inbright+"', indirection='"+indirection+"', incell='"+incell+"', incenter='"+incenter+"', inwidth='"+inwidth+"', complist='"+complist+"', compwidth='"+compwidth+"', setpointings="+str(True)+", ptgfile='$project.ptg.txt', integration='"+integration_tp+"', direction='"+str(dir_tp)+"', mapsize="+str(mapsize_tp)+", maptype='"+maptype_tp+"', pointingspacing='"+ptgspacing_tp+"', caldirection='"+caldirection+"', calflux='"+calflux+"', obsmode='"+obsmode_sd+"', refdate='"+refdate+"', hourangle='"+hourangle+"', totaltime='"+tottime_tp+"', antennalist='', sdantlist='"+antlist_tp+"', sdant="+str(tpantid)+", thermalnoise='"+thermalnoise+"', user_pwv="+str(pwv)+", t_ground="+str(t_ground)+", leakage="+str(leakage)+", graphics='"+graphics+"', verbose="+str(verbose)+", overwrite="+str(overwrite)+")" 00436 msg("Executing: "+taskstr, origin="simalma", priority=v_priority) 00437 00438 try: 00439 simobserve(project=project, 00440 skymodel=skymodel, inbright=inbright, 00441 indirection=indirection, incell=incell, 00442 incenter=incenter, inwidth=inwidth, 00443 complist=complist, compwidth=compwidth, 00444 setpointings=True, ptgfile='$project.ptg.txt', 00445 integration=integration_tp, 00446 direction=dir_tp, mapsize=mapsize_tp, 00447 maptype=maptype_tp, pointingspacing=ptgspacing_tp, 00448 caldirection=caldirection, calflux=calflux, 00449 obsmode=obsmode_sd, refdate=refdate, 00450 hourangle=hourangle, totaltime=tottime_tp, 00451 antennalist="", sdantlist=antlist_tp, sdant=tpantid, 00452 thermalnoise=thermalnoise, user_pwv=pwv, 00453 t_ground=t_ground, #t_sky=t_sky, tau0=tau0, seed=seed, 00454 leakage=leakage, 00455 graphics=graphics, verbose=verbose, overwrite=overwrite)#, 00456 #async=False) 00457 except: 00458 raise Exception, simobserr 00459 finally: 00460 casalog.origin('simalma') 00461 00462 ################################################################ 00463 # Imaging 00464 if image: 00465 modelimage = "" 00466 if acaratio > 0: 00467 ######################################################## 00468 # Image ACA-7m + ACA-TP 00469 step += 1 00470 msg("="*60, origin="simalma", priority="warn") 00471 msg(" Step %d: generating ACA 7-m array + total power image. " % step, origin="simalma", priority="warn") 00472 msg("="*60, origin="simalma", priority="warn") 00473 if os.path.exists(fileroot+"/"+msname_aca): 00474 vis_aca = msname_aca+"," 00475 else: 00476 msg("ACA is requested but ACA 7-m MS '%s' is not found" \ 00477 % msname_aca, origin="simalma", priority="error") 00478 if os.path.exists(fileroot+"/"+msname_tp): 00479 vis_aca += msname_tp 00480 else: 00481 msg("ACA is requested but total power MS '%s' is not found" \ 00482 % msname_tp, origin="simalma", priority="error") 00483 00484 cell_aca = cell 00485 00486 # Define imsize to cover TP map region 00487 #imsize_aca = 0 00488 msg("Defining image size of ACA to cover map region of total power simulation", origin="simalma", priority=v_priority) 00489 msg("- The total power map size: [%s, %s]" % \ 00490 (qa.tos(qimgsize_tp[0]), qa.tos(qimgsize_tp[1])), \ 00491 origin="simalma", priority=v_priority) 00492 if cell != '': 00493 # user-defined cell size 00494 msg("- The user defined cell size: %s" % cell, \ 00495 origin="simalma", priority=v_priority) 00496 imgcell = [cell, cell] 00497 else: 00498 if model_cell == None: 00499 # components only simulation 00500 compmodel = fileroot+"/"+pref_bl+".compskymodel" 00501 msg("getting the cell size of input compskymodel", \ 00502 origin="simalma", priority=v_priority) 00503 if not os.path.exists(compmodel): 00504 msg("Could not find the skymodel, '%s'" % \ 00505 compmodel, priority='error') 00506 # modifymodel just collects info if outmodel==inmodel 00507 model_vals = myutil.modifymodel(compmodel,compmodel, 00508 "","","","","",-1, 00509 flatimage=False) 00510 model_cell = model_vals[1] 00511 model_size = model_vals[2] 00512 00513 # skymodel (+ components list) simulation 00514 msg("- The cell size of input skymodel: [%s, %s]" % \ 00515 (qa.tos(model_cell[0]), qa.tos(model_cell[1])), \ 00516 origin="simalma", priority=v_priority) 00517 imgcell = model_cell 00518 00519 imsize_aca = calc_imsize(mapsize=qimgsize_tp, cell=imgcell) 00520 00521 msg("---> The number of pixels needed to cover the map region: [%d, %d]" % \ 00522 (imsize_aca[0], imsize_aca[1]), \ 00523 origin="simalma", priority=v_priority) 00524 00525 msg("Compare with BL imsize and adopt the larger one", \ 00526 origin="simalma", priority=v_priority) 00527 # Compare with imsize of BL (note: imsize is an intArray) 00528 if is_array_type(imsize) and imsize[0] > 0: 00529 # User has defined imsize 00530 if len(imsize) > 1: 00531 imsize_bl = imsize[0:2] 00532 else: 00533 imsize_bl = [imsize[0], imsize[0]] 00534 msg("---> BL imsize (user defined): [%d, %d]" % \ 00535 (imsize_bl[0], imsize_bl[1]), \ 00536 origin="simalma", priority=v_priority) 00537 else: 00538 # the same as input model (calculate from model_size) 00539 msg("estimating imsize of BL from input sky model.", \ 00540 origin="simalma", priority=v_priority) 00541 imsize_bl = calc_imsize(mapsize=model_size, cell=imgcell) 00542 msg("---> Estimated BL imsize (sky model): [%d, %d]" % \ 00543 (imsize_bl[0], imsize_bl[1]), \ 00544 origin="simalma", priority=v_priority) 00545 00546 imsize_aca = [max(imsize_aca[0], imsize_bl[0]), \ 00547 max(imsize_aca[1], imsize_bl[1])] 00548 00549 msg("The image size of ACA: [%d, %d]" % \ 00550 (imsize_aca[0], imsize_aca[1]), \ 00551 origin="simalma", priority=v_priority) 00552 00553 taskstr = "simanalyze(project='"+project+"', image="+str(image)+", vis='"+vis_aca+"', modelimage='', cell='"+str(cell_aca)+"', imsize="+str(imsize_aca)+", imdirection='"+imdirection+"', niter="+str(niter)+", threshold='"+threshold+"', weighting='"+weighting+"', mask="+str([])+", outertaper="+str([])+", stokes='I', analyze="+str(True)+", graphics='"+graphics+"', verbose="+str(verbose)+", overwrite="+str(overwrite)+")" 00554 msg("Executing: "+taskstr, origin="simalma", priority=v_priority) 00555 00556 try: 00557 simanalyze(project=project, image=image, 00558 vis=vis_aca, modelimage="", 00559 cell=cell_aca, imsize=imsize_aca, 00560 imdirection=imdirection, niter=niter, 00561 threshold=threshold, weighting=weighting, 00562 mask=[], outertaper=[], stokes='I', 00563 analyze=True, 00564 #showuv=None, showpsf=None, showmodel=None, 00565 #showconvolved=None, showclean=None, 00566 #showresidual=None, showdifference=None, 00567 #showfidelity=None, 00568 graphics=graphics, verbose=verbose, 00569 overwrite=overwrite)#, async=False) 00570 except: 00571 raise Exception, simanaerr 00572 finally: 00573 casalog.origin('simalma') 00574 00575 # Back up simanalyze.last file 00576 if os.path.exists(fileroot+"/"+simana_file): 00577 simana_new = imagename_aca.replace(".image",".simanalyze.last") 00578 msg("Back up input parameter set to '%s'" % simana_new, \ 00579 origin="simalma", priority=v_priority) 00580 shutil.move(fileroot+"/"+simana_file, fileroot+"/"+simana_new) 00581 00582 if not os.path.exists(fileroot+"/"+imagename_aca): 00583 msg("ACA model image '%s' is not found" \ 00584 % imagename_aca, origin="simalma", priority="error") 00585 modelimage = imagename_aca 00586 ############################################################ 00587 # Image ALMA-BL 00588 step += 1 00589 msg("="*60, origin="simalma", priority="warn") 00590 if acaratio > 0: 00591 msg(" Step %d: generating ALMA 12-m + ACA image." % step, origin="simalma", priority="warn") 00592 else: 00593 msg(" Step %d: generating ALMA 12-m array image." % step, origin="simalma", priority="warn") 00594 msg("="*60, origin="simalma", priority="warn") 00595 00596 if os.path.exists(fileroot+"/"+msname_bl): 00597 vis_bl = msname_bl 00598 else: 00599 msg("Could not find MS to image, '%s'" \ 00600 % msname_bl, origin="simalma", priority="error") 00601 00602 taskstr = "simanalyze(project='"+ project+"', image="+str(image)+", vis='"+ vis_bl+"', modelimage='"+ modelimage+"', cell='"+str(cell)+"', imsize="+str(imsize)+", imdirection='"+ imdirection+"', niter="+str(niter)+", threshold='"+ threshold+"', weighting='"+ weighting+"', mask="+str([])+", outertaper="+str([])+", stokes='I', analyze="+str(True)+", graphics='"+ graphics+"', verbose="+str(verbose)+", overwrite="+ str(overwrite)+")" 00603 msg("Executing: "+taskstr, origin="simalma", priority=v_priority) 00604 00605 try: 00606 simanalyze(project=project, image=image, 00607 vis=vis_bl, modelimage=modelimage, 00608 cell=cell, imsize=imsize, imdirection=imdirection, 00609 niter=niter, threshold=threshold, weighting=weighting, 00610 mask=[], outertaper=[], stokes='I', 00611 analyze=True, 00612 #showuv=None, showpsf=None, showmodel=None, 00613 #showconvolved=None, showclean=None, 00614 #showresidual=None, showdifference=None, 00615 #showfidelity=None, 00616 graphics=graphics, verbose=verbose, 00617 overwrite=overwrite)#, async=False) 00618 except: 00619 raise Exception, simanaerr 00620 finally: 00621 casalog.origin('simalma') 00622 00623 00624 except TypeError, e: 00625 finalize_tools() 00626 casalog.post("simalma -- TypeError: %s" % str(e), priority="ERROR") 00627 raise TypeError, e 00628 return False 00629 except ValueError, e: 00630 finalize_tools() 00631 casalog.post("simalma -- OptionError: %s" % str(e), priority="ERROR") 00632 raise ValueError, e 00633 return False 00634 except Exception, instance: 00635 finalize_tools() 00636 casalog.post("simalma -- Exception: %s" % str(instance), 00637 priority="ERROR") 00638 raise Exception, instance 00639 return False 00640 return True 00641 00642 def finalize_tools(): 00643 if ia.isopen(): ia.close() 00644 sm.close() 00645 #cl.close() # crashes casa 00646 00647 def get_data_prefix(cfgname, project=""): 00648 if str.upper(cfgname[0:4]) == "ALMA": 00649 foo=cfgname.replace(';','_') 00650 else: 00651 foo = cfgname 00652 foo=foo.replace(".cfg","") 00653 sfoo=foo.split('/') 00654 if len(sfoo)>1: foo=sfoo[-1] 00655 return project+"."+foo 00656 00657 def is_array_type(value): 00658 array_type = [list, tuple, numpy.ndarray] 00659 if type(value) in array_type: 00660 return True 00661 else: 00662 return False 00663 00664 def calc_imsize(mapsize=None, cell=None): 00665 if mapsize == None: 00666 raise ValueError, "mapsize is not defined" 00667 if cell == None: 00668 raise ValueError, "cell is not defined" 00669 # get a list of cell size 00670 if is_array_type(cell): 00671 if len(cell) < 2: 00672 cell = [cell[0], cell[0]] 00673 else: 00674 cell = [cell, cell] 00675 00676 for qval in cell: 00677 if not qa.compare(qval, "deg"): 00678 raise TypeError, "cell should be an angular size" 00679 00680 qcellx = qa.quantity(cell[0]) 00681 qcelly = qa.quantity(cell[1]) 00682 00683 # get a list of map size 00684 if is_array_type(mapsize): 00685 if len(mapsize) < 2: 00686 mapsize = [mapsize[0], mapsize[0]] 00687 else: 00688 mapsize = [mapsize, mapsize] 00689 00690 for qval in mapsize: 00691 if not qa.compare(qval, "deg"): 00692 raise TypeError, "mapsize should be an angular size" 00693 00694 vsizex = qa.convert(mapsize[0], qcellx['unit'])['value'] 00695 vsizey = qa.convert(mapsize[1], qcelly['unit'])['value'] 00696 00697 # Calculate the number of pixels to cover the map size 00698 npixx = int(numpy.ceil(vsizex/qcellx['value'])) 00699 npixy = int(numpy.ceil(vsizey/qcelly['value'])) 00700 00701 return [npixx, npixy] 00702
1.8.0
