|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 from taskinit import * 00002 from simutil import * 00003 import os 00004 import re 00005 import pylab as pl 00006 import pdb 00007 00008 def simobserve( 00009 project=None, 00010 skymodel=None, inbright=None, indirection=None, incell=None, 00011 incenter=None, inwidth=None, # innchan=None, 00012 complist=None, compwidth=None, 00013 setpointings=None, 00014 ptgfile=None, integration=None, direction=None, mapsize=None, 00015 maptype=None, pointingspacing=None, caldirection=None, calflux=None, 00016 # observe=None, 00017 obsmode=None, 00018 refdate=None, hourangle=None, 00019 totaltime=None, antennalist=None, 00020 sdantlist=None, 00021 sdant=None, 00022 thermalnoise=None, 00023 user_pwv=None, t_ground=None, t_sky=None, tau0=None, seed=None, 00024 leakage=None, 00025 graphics=None, 00026 verbose=None, 00027 overwrite=None, 00028 async=False): 00029 00030 # Collect a list of parameter values to save inputs 00031 in_params = locals() 00032 00033 import re 00034 00035 try: 00036 00037 # RI TODO for inbright=unchanged, need to scale input image to jy/pix 00038 # according to actual units in the input image 00039 00040 casalog.origin('simobserve') 00041 if verbose: casalog.filter(level="DEBUG2") 00042 00043 a = inspect.stack() 00044 stacklevel = 0 00045 for k in range(len(a)): 00046 if (string.find(a[k][1], 'ipython console') > 0): 00047 stacklevel = k 00048 myf = sys._getframe(stacklevel).f_globals 00049 00050 # create the utility object: 00051 util = simutil(direction) # this is the dir of the observation - could be "" 00052 if verbose: util.verbose = True 00053 msg = util.msg 00054 00055 # it was requested to make the user interface "observe" for what 00056 # is sm.observe and sm.predict. 00057 # interally the code is clearer if we stick with predict so 00058 predict = obsmode.startswith('i') or obsmode.startswith('s') 00059 if predict: 00060 uvmode = obsmode.startswith('i') 00061 if not uvmode: antennalist = sdantlist 00062 elif sdantlist != "": 00063 if antennalist == "": 00064 uvmode = False 00065 antennalist = sdantlist 00066 else: 00067 #uvmode = True 00068 #msg("Both antennalist and sdantlist are defined. sdantlist will be ignored",priority="warn") 00069 msg("Both antennalist and sdantlist are defined. Define one of them.",priority="error") 00070 return False 00071 else: 00072 uvmode = True 00073 # uvmode = (sdant < 0) #when flexible default values come available 00074 00075 00076 # put output in directory called "project" 00077 fileroot = project 00078 if not os.path.exists(fileroot): 00079 os.mkdir(fileroot) 00080 00081 00082 # filename parsing of cfg file here so that the project filenames 00083 # can contain the cfg 00084 repodir = os.getenv("CASAPATH").split(' ')[0] + "/data/alma/simmos/" 00085 00086 # convert "alma;0.4arcsec" to an actual configuration 00087 # can only be done after reading skymodel, so here, we just string parse 00088 if str.upper(antennalist[0:4]) == "ALMA": 00089 foo=antennalist.replace(';','_') 00090 elif len(antennalist) > 0: 00091 foo=antennalist 00092 else: 00093 msg("The name of antenna list (antennalist/sdantlist) is not specified",priority="error") 00094 00095 if foo: 00096 foo=foo.replace(".cfg","") 00097 sfoo=foo.split('/') 00098 if len(sfoo)>1: 00099 foo=sfoo[-1] 00100 project=project+"."+foo 00101 00102 00103 if not overwrite: 00104 if (predict and uvmode and os.path.exists(fileroot+"/"+project+".ms")): 00105 msg(fileroot+"/"+project+".ms exists but overwrite=F",priority="error") 00106 return False 00107 if (predict and (not uvmode) and os.path.exists(fileroot+"/"+project+".sd.ms")): 00108 msg(fileroot+"/"+project+".sd.ms exists but overwrite=F",priority="error") 00109 return False 00110 00111 00112 saveinputs = myf['saveinputs'] 00113 saveinputs('simobserve',fileroot+"/"+project+".simobserve.last", 00114 myparams=in_params) 00115 00116 00117 00118 # some hardcoded variables that may be reintroduced in future development 00119 relmargin = .5 # number of PB between edge of model and pointing centers 00120 scanlength = 1 # number of integrations per scan 00121 00122 if type(skymodel) == type([]): 00123 skymodel = skymodel[0] 00124 skymodel = skymodel.replace('$project',project) 00125 00126 if type(complist) == type([]): 00127 complist = complist[0] 00128 00129 if((not os.path.exists(skymodel)) and (not os.path.exists(complist))): 00130 msg("No sky input found. At least one of skymodel or complist must be set.",priority="error") 00131 return False 00132 00133 00134 grscreen = False 00135 grfile = False 00136 if graphics == "both": 00137 grscreen = True 00138 grfile = True 00139 if graphics == "screen": 00140 grscreen = True 00141 if graphics == "file": 00142 grfile = True 00143 00144 ################################################################## 00145 # set up skymodel image 00146 00147 if os.path.exists(skymodel): 00148 components_only = False 00149 # create a new skymodel called skymodel, or if its already there, called newmodel 00150 default_model = project + ".skymodel" 00151 if skymodel == default_model: 00152 newmodel = fileroot + "/" + project + ".newmodel" 00153 else: 00154 newmodel = fileroot + "/" + default_model 00155 if os.path.exists(newmodel): 00156 if overwrite: 00157 shutil.rmtree(newmodel) 00158 else: 00159 msg(newmodel+" exists -- please delete it, change skymodel, or set overwrite=T",priority="error") 00160 return False 00161 00162 # modifymodel just collects info if skymodel==newmodel 00163 innchan = -1 00164 returnpars = util.modifymodel(skymodel,newmodel, 00165 inbright,indirection,incell, 00166 incenter,inwidth,innchan, 00167 flatimage=False) 00168 if not returnpars: 00169 return False 00170 00171 (model_refdir,model_cell,model_size, 00172 model_nchan,model_center,model_width, 00173 model_stokes) = returnpars 00174 00175 modelflat = fileroot + "/" + project + ".skymodel.flat" 00176 if os.path.exists(modelflat) and (not predict): 00177 # if we're not predicting, then we want to use the previously 00178 # created modelflat, because it may have components added 00179 msg("flat sky model "+modelflat+" exists, predict not requested",priority="warn") 00180 msg(" working from existing model image - please delete it if you wish to overwrite.",priority="warn") 00181 else: 00182 # create and add components into modelflat with util.flatimage() 00183 util.flatimage(newmodel,complist=complist,verbose=verbose) 00184 # we want the skymodel.flat image to be called that no matter what 00185 # the skymodel image is called, since that's what used in analysis 00186 if modelflat != newmodel+".flat": 00187 if os.path.exists(modelflat): 00188 shutil.rmtree(modelflat) 00189 shutil.move(newmodel+".flat",modelflat) 00190 00191 casalog.origin('simobserve') 00192 00193 # set startfeq and bandwidth in util object after modifymodel 00194 bandwidth = qa.mul(qa.quantity(model_nchan),qa.quantity(model_width)) 00195 util.bandwidth = bandwidth 00196 00197 else: 00198 components_only = True 00199 # calculate model parameters from the component list: 00200 00201 compdirs = [] 00202 cl.done() 00203 cl.open(complist) 00204 00205 for i in range(cl.length()): 00206 compdirs.append(util.dir_m2s(cl.getrefdir(i))) 00207 00208 model_refdir, coffs = util.average_direction(compdirs) 00209 model_center = cl.getspectrum(0)['frequency']['m0'] 00210 # components don't yet support spectrum 00211 if util.isquantity(compwidth,halt=False): 00212 model_width = compwidth 00213 else: 00214 model_width = "2GHz" 00215 msg("component-only simulation, compwidth unset: setting bandwidth to 2GHz",priority="warn") 00216 00217 model_nchan = 1 00218 model_stokes = "I" 00219 00220 cmax = 0.0014 # ~5 arcsec 00221 for i in range(coffs.shape[1]): 00222 xc = pl.absolute(coffs[0,i]) # offsets in deg 00223 yc = pl.absolute(coffs[1,i]) 00224 if xc > cmax: 00225 cmax = xc 00226 if yc > cmax: 00227 cmax = yc 00228 00229 model_size = ["%fdeg" % (3*cmax), "%fdeg" % (3*cmax)] 00230 00231 00232 # for cases either if there is a skymodel or if there are only components, 00233 # if the user has not input a map size (for setpointings), then use model_size 00234 if len(mapsize) == 0: 00235 mapsize = model_size 00236 if verbose: msg("setting map size to "+str(model_size)) 00237 else: 00238 if type(mapsize) == type([]): 00239 if len(mapsize[0]) == 0: 00240 mapsize = model_size 00241 if verbose: msg("setting map size to "+str(model_size)) 00242 00243 if components_only: 00244 if type(mapsize) == type([]): 00245 map_asec = qa.convert(mapsize[0],"arcsec")['value'] 00246 else: 00247 map_asec = qa.convert(mapsize,"arcsec")['value'] 00248 if type(model_size) == type([]): 00249 mod_asec = qa.convert(model_size[0],"arcsec")['value'] 00250 else: 00251 mod_asec = qa.convert(model_size,"arcsec")['value'] 00252 if map_asec>mod_asec: 00253 model_size=["%farcsec" % map_asec,"%farcsec" % map_asec] 00254 00255 00256 00257 00258 ################################################################## 00259 # read antenna file here to get Primary Beam, and estimate maxbase and psfsize 00260 aveant = -1 00261 stnx = [] # for later, to know if we read an array in or not 00262 pb = 0. # primary beam 00263 00264 # convert "alma;0.4arcsec" to an actual configuration 00265 if str.upper(antennalist[0:4]) == "ALMA": 00266 00267 # test for cycle 1 00268 q = re.compile('.*CYCLE.?1.?;(.*)') 00269 qq = q.match(antennalist.upper()) 00270 if qq: 00271 z = qq.groups() 00272 tail=z[0] 00273 tail=tail.lower() 00274 if util.isquantity(tail,halt=False): 00275 resl = qa.convert(tail,"arcsec")['value'] 00276 if os.path.exists(repodir): 00277 confnum = (1.044-6.733*pl.log10(resl*qa.convert(model_center,"GHz")['value']/345.)) 00278 confnum = max(1,min(6,confnum)) 00279 conf = str(int(round(confnum))) 00280 antennalist = repodir + "alma_cycle1_" + conf + ".cfg" 00281 msg("converted resolution to antennalist "+antennalist) 00282 else: 00283 msg("failed to find antenna configuration repository at "+repodir,priority="warn") 00284 00285 else: # assume FS 00286 tail = antennalist[5:] 00287 if util.isquantity(tail,halt=False): 00288 resl = qa.convert(tail,"arcsec")['value'] 00289 if os.path.exists(repodir): 00290 confnum = (2.867-pl.log10(resl*1000*qa.convert(model_center,"GHz")['value']/672.))/0.0721 00291 confnum = max(1,min(28,confnum)) 00292 conf = str(int(round(confnum))) 00293 if len(conf) < 2: conf = '0' + conf 00294 antennalist = repodir + "alma.out" + conf + ".cfg" 00295 msg("converted resolution to antennalist "+antennalist) 00296 else: 00297 msg("failed to find antenna configuration repository at "+repodir,priority="warn") 00298 00299 00300 # Search order is fileroot/ -> specified path -> repository 00301 if len(antennalist) > 0: 00302 if os.path.exists(fileroot+"/"+antennalist): 00303 antennalist = fileroot + "/" + antennalist 00304 elif not os.path.exists(antennalist) and \ 00305 os.path.exists(repodir+antennalist): 00306 antennalist = repodir + antennalist 00307 # Now make sure the antennalist exists 00308 if not os.path.exists(antennalist): 00309 util.msg("Couldn't find antennalist: %s" % antennalist, priority="error") 00310 return False 00311 elif predict or components_only: 00312 # antennalist is required when predicting or components only 00313 util.msg("Must specify antennalist", priority="error") 00314 00315 # Read antennalist 00316 if os.path.exists(antennalist): 00317 stnx, stny, stnz, stnd, padnames, nant, telescopename = util.readantenna(antennalist) 00318 if len(stnx) == 1: 00319 if predict and uvmode: 00320 # observe="int" but antennalist is SD 00321 util.msg("antennalist contains only 1 antenna", priority="error") 00322 uvmode = False 00323 antnames = [] 00324 if not uvmode: #Single-dish 00325 # KS TODO: what if not predicting but SD with multi-Ants 00326 # in antennalist (e.g., aca.tp)? In that case, PB on plots and 00327 # pointingspacing="??PB" will not be correct for heterogeneous list. 00328 if sdant > nant-1 or sdant < -nant: 00329 msg("antenna index %d is out of range. setting sdant=0"%sdant,priority="warn") 00330 sdant = 0 00331 stnx = [stnx[sdant]] 00332 stny = [stny[sdant]] 00333 stnz = [stnz[sdant]] 00334 stnd = pl.array(stnd[sdant]) 00335 padnames = [padnames[sdant]] 00336 nant = 1 00337 00338 for k in xrange(0,nant): antnames.append('A%02d'%k) 00339 aveant = stnd.mean() 00340 # TODO use max ant = min PB instead? 00341 # (set back to simdata - there must be an automatic way to do this) 00342 casalog.origin('simobserve') 00343 pb = 1.2*0.29979/qa.convert(qa.quantity(model_center),'GHz')['value']/aveant*3600.*180/pl.pi # arcsec 00344 00345 # PSF size 00346 if uvmode: 00347 # approx max baseline, to compare with model_cell: 00348 cx=pl.mean(stnx) 00349 cy=pl.mean(stny) 00350 cz=pl.mean(stnz) 00351 lat,lon = util.itrf2loc(stnx,stny,stnz,cx,cy,cz) 00352 maxbase=max(lat)-min(lat) # in meters 00353 maxbase2=max(lon)-min(lon) 00354 if maxbase2>maxbase: 00355 maxbase=maxbase2 00356 # estimate the psf size from the minimum spatial scale 00357 psfsize = 0.3/qa.convert(qa.quantity(model_center),'GHz')['value']/maxbase*3600.*180/pl.pi # lambda/b converted to arcsec 00358 del cx, cy, cz, lat, lon, maxbase, maxbase2 00359 else: # Single-dish 00360 psfsize = pb 00361 # check for model size 00362 if not components_only: 00363 minsize = min(qa.convert(model_size[0],'arcsec')['value'],\ 00364 qa.convert(model_size[1],'arcsec')['value']) 00365 if minsize < 2*pb: 00366 msg("skymodel should be larger than 2*primary beam. Your skymodel: %.3f arcsec < %.3f arcsec: 2*primary beam" % (minsize, 2*pb),priority="error") 00367 del minsize 00368 else: 00369 msg("Can't find antennalist",priority="error") 00370 return False 00371 00372 00373 # now we have an estimate of the psf from the antenna configuration, 00374 # so we can guess a model_cell for the case of component-only 00375 # simulation, 00376 if components_only: 00377 # first set based on psfsize: 00378 # needs high subsampling because small shifts in placement of 00379 # components lead to large changes in the difference image. 00380 model_cell = [ str(psfsize/20)+"arcsec", str(psfsize/20)+"arcsec" ] 00381 00382 # XXX if the user has set direction should we center the compskymodel there? 00383 # if len(direction)>0: model_refdir = direction 00384 00385 # and can create a compskymodel image (tmp) and 00386 # skymodel.flat which is what is needed for analysis. 00387 00388 if components_only: 00389 newmodel = fileroot + "/" + project + ".compskymodel" 00390 needmodel=True 00391 00392 modimsize=int((qa.convert(model_size[0],"arcsec")['value'])/(qa.convert(model_cell[0],"arcsec")['value'])) 00393 # modimsize=max([modimsize,32]) 00394 newepoch,newlat,newlon = util.direction_splitter(model_refdir) 00395 00396 if os.path.exists(newmodel): 00397 if overwrite: 00398 shutil.rmtree(newmodel) 00399 else: 00400 needmodel=False 00401 ia.open(newmodel) 00402 oldshape=ia.shape() 00403 if len(oldshape) != 2: 00404 needmodel=True 00405 else: 00406 if oldshape[0] != modimsize or oldshape[1]==modimsize: 00407 needmodel=True 00408 oldcs=ia.coordsys() 00409 ia.done() 00410 olddir = (oldcs.referencevalue())['numeric'] 00411 if ( olddir[0] != qa.convert(newlat,oldcs.units()[0])['value'] or 00412 olddir[1] != qa.convert(newlon,oldcs.units()[1])['value'] or 00413 newepoch != oldcs.referencecode() ): 00414 needmodel=True 00415 oldcs.done() 00416 del oldcs, olddir 00417 if needmodel: 00418 msg(newmodel+" exists and is inconsistent with required size="+str(modimsize)+" and direction. Please set overwrite=True",priority="error") 00419 return False 00420 00421 if needmodel: 00422 csmodel = ia.newimagefromshape(newmodel,[modimsize,modimsize,1,1]) 00423 modelcsys = csmodel.coordsys() 00424 modelshape = csmodel.shape() 00425 cell0_rad=qa.convert(model_cell[0],'rad')['value'] 00426 cell1_rad=qa.convert(model_cell[1],'rad')['value'] 00427 modelcsys.setdirection(refpix=[modimsize/2,modimsize/2], 00428 refval=[qa.tos(newlat),qa.tos(newlon)], 00429 refcode=newepoch) 00430 modelcsys.setincrement([-cell0_rad,cell1_rad],'direction') 00431 modelcsys.setreferencevalue(type="spectral",value=qa.tos(model_center)) 00432 modelcsys.setrestfrequency(qa.tos(model_center)) 00433 modelcsys.setincrement(type="spectral",value=compwidth) 00434 csmodel.setcoordsys(modelcsys.torecord()) 00435 modelcsys.done() 00436 cl.open(complist) 00437 csmodel.setbrightnessunit("Jy/pixel") 00438 csmodel.modify(cl.torecord(),subtract=False) 00439 cl.done() 00440 csmodel.done() 00441 # as noted, compskymodel doesn't need to exist, only skymodel.flat 00442 # flatimage adds in components if complist!=None 00443 #util.flatimage(newmodel,complist=complist,verbose=verbose) 00444 util.flatimage(newmodel,verbose=verbose) 00445 modelflat = fileroot + "/" + project + ".compskymodel.flat" 00446 # modelflat = fileroot + "/" + project + ".skymodel.flat" 00447 # if modelflat != newmodel+".flat": 00448 # if os.path.exists(modelflat): 00449 # shutil.rmtree(modelflat) 00450 # shutil.move(newmodel+".flat",modelflat) 00451 00452 00453 # and finally, with model_cell set either from an actual skymodel, 00454 # or from the antenna configuration in components_only case, 00455 # we can check for the user that the psf is likely to be sampled enough: 00456 cell_asec=qa.convert(model_cell[0],'arcsec')['value'] 00457 if psfsize < cell_asec: 00458 msg("Sky model cell of "+str(cell_asec)+" asec is very large compared to highest resolution "+str(psfsize)+" asec - this will lead to blank or erroneous output. (Did you set incell?)",priority="error") 00459 shutil.rmtree(modelflat) 00460 return False 00461 if psfsize < 2*cell_asec: 00462 msg("Sky model cell of "+str(cell_asec)+" asec is large compared to highest resolution "+str(psfsize)+" asec. (Did you set incell?)",priority="warn") 00463 00464 # set this for future minimum image size 00465 minimsize = 8* int(psfsize/cell_asec) 00466 00467 00468 00469 ################################################################## 00470 # set up pointings 00471 dir = model_refdir 00472 dir0 = dir 00473 if type(direction) == type([]): 00474 if len(direction) > 0: 00475 if util.isdirection(direction[0],halt=False): 00476 dir = direction 00477 dir0 = direction[0] 00478 else: 00479 if util.isdirection(direction,halt=False): 00480 dir = direction 00481 dir0 = dir 00482 util.direction = dir0 00483 00484 # if the integration time is a real time quantity 00485 # test for weird units 00486 if not util.isquantity(integration): 00487 msg("integration time "+integration+" does not appear to represent a time interval (use 's','min'; not 'sec','m')",priority="error") 00488 return False 00489 00490 if qa.quantity(integration)['unit'] != '': 00491 if not qa.compare(integration,"1s"): 00492 msg("integration time "+integration+" does not appear to represent a time interval ('s','min'; not 'sec','m')",priority="error") 00493 return False 00494 intsec = qa.convert(qa.quantity(integration),"s")['value'] 00495 else: 00496 if len(integration)>0: 00497 intsec = float(integration) 00498 msg("interpreting integration time parameter as "+str(intsec)+"s",priority="warn") 00499 else: 00500 intsec = 0 00501 integration="%fs" %intsec 00502 00503 00504 if setpointings: 00505 if verbose: util.msg("calculating map pointings centered at "+str(dir0)) 00506 00507 if len(pointingspacing) < 1: 00508 pointingspacing = "Nyquist" 00509 if str.upper(pointingspacing)=="NYQUIST": 00510 pointingspacing="0.48113PB" 00511 q = re.compile('(\d+.?\d+)\s*PB') 00512 qq = q.match(pointingspacing.upper()) 00513 if qq: 00514 z = qq.groups() 00515 if pb <= 0: 00516 util.msg("Can't calculate pointingspacing in terms of primary beam because antennalist is not specified",priority="error") 00517 return False 00518 pointingspacing = "%farcsec" % (float(z[0])*pb) 00519 # todo make more robust to nonconforming z[0] strings 00520 if verbose: 00521 msg("pointing spacing in mosaic = "+pointingspacing) 00522 pointings = util.calc_pointings2(pointingspacing,mapsize,maptype=maptype, direction=dir, beam=pb) 00523 nfld=len(pointings) 00524 etime = qa.convert(qa.mul(qa.quantity(integration),scanlength),"s")['value'] 00525 # etime is an array of scan lengths - here they're all the same. 00526 etime = etime + pl.zeros(nfld) 00527 # totaltime might not allow all fields to be observed, or it might 00528 # repeat 00529 ptgfile = fileroot + "/" + project + ".ptg.txt" 00530 else: 00531 if type(ptgfile) == type([]): 00532 ptgfile = ptgfile[0] 00533 ptgfile = ptgfile.replace('$project',project) 00534 # precedence to ptg file outside the project dir 00535 if os.path.exists(ptgfile): 00536 shutil.copyfile(ptgfile,fileroot+"/"+project + ".ptg.txt") 00537 ptgfile = fileroot + "/" + project + ".ptg.txt" 00538 else: 00539 if os.path.exists(fileroot+"/"+ptgfile): 00540 ptgfile = fileroot + "/" + ptgfile 00541 else: 00542 util.msg("Can't find pointing file "+ptgfile,priority="error") 00543 return False 00544 00545 nfld, pointings, etime = util.read_pointings(ptgfile) 00546 if max(etime) <= 0: 00547 # integration is a string in input params 00548 etime = intsec 00549 # make etime into an array 00550 etime = etime + pl.zeros(nfld) 00551 # etimes determine stop/start i.e. the scan 00552 # if a longer etime is in the file, it'll do multiple integrations 00553 # per scan 00554 # expects that the cal is separate, and this is just one round of the mosaic 00555 # furthermore, the cal will use _integration_ from the inputs, and that 00556 # needs to be less than the min etime: 00557 if min(etime) < intsec: 00558 integration = str(min(etime))+"s" 00559 msg("Setting integration to "+integration+" to match the shortest time in the pointing file.",priority="warn") 00560 intsec = min(etime) 00561 00562 00563 # find imcenter - phase center 00564 imcenter , offsets = util.median_direction(pointings) 00565 epoch, ra, dec = util.direction_splitter(imcenter) 00566 00567 # model is centered at model_refdir, and has model_size; 00568 mepoch, mra, mdec = util.direction_splitter(model_refdir) 00569 # ra/mra should be in degrees, but just in case 00570 ra=qa.convert(ra,'deg') 00571 mra=qa.convert(mra,'deg') 00572 dec=qa.convert(dec,'deg') 00573 mdec=qa.convert(mdec,'deg') 00574 00575 # observation near ra=0: 00576 if abs(mra['value'])<10 or abs(mra['value']-360.)<10 or abs(ra['value'])<10 or abs(mra['value']-360.)<10: 00577 nearzero=True 00578 else: 00579 nearzero=False 00580 00581 # fix wraps 00582 if nearzero: # put break at 180 00583 if ra['value'] >= 180.: 00584 ra['value'] = ra['value'] - 360. 00585 if mra['value'] >= 180.: 00586 mra['value'] = mra['value'] - 360. 00587 else: 00588 if ra['value'] >= 360.: 00589 ra['value'] = ra['value'] - 360. 00590 if mra['value'] >= 360.: 00591 mra['value'] = mra['value'] - 360. 00592 00593 # shift is the offset from the model to imcenter 00594 shift = [ (qa.convert(ra,'deg')['value'] - 00595 qa.convert(mra,'deg')['value'])*pl.cos(qa.convert(dec,'rad')['value'] ), 00596 (qa.convert(dec,'deg')['value'] - qa.convert(mdec,'deg')['value']) ] 00597 if verbose: 00598 msg("pointings are shifted relative to the model by %g,%g arcsec" % (shift[0]*3600,shift[1]*3600)) 00599 00600 xmax = qa.convert(model_size[0],'deg')['value']*0.5 00601 ymax = qa.convert(model_size[1],'deg')['value']*0.5 00602 overlap = False 00603 # wrapang in median_direction should make offsets always small, not >360 00604 for i in xrange(offsets.shape[1]): 00605 xc = pl.absolute(offsets[0,i]+shift[0]) # offsets and shift are in degrees 00606 yc = pl.absolute(offsets[1,i]+shift[1]) 00607 if xc < xmax and yc < ymax: 00608 overlap = True 00609 break 00610 00611 if setpointings: 00612 if os.path.exists(ptgfile): 00613 if overwrite: 00614 os.remove(ptgfile) 00615 else: 00616 util.msg("pointing file "+ptgfile+" already exists and user does not want to overwrite",priority="error") 00617 return False 00618 util.write_pointings(ptgfile,pointings,etime.tolist()) 00619 00620 msg("center = "+imcenter) 00621 if nfld > 1 and verbose: 00622 for idir in range(min(len(pointings),20)): 00623 msg(" "+pointings[idir]) 00624 if nfld >= 20: 00625 msg(" (printing only first 20 - see pointing file for full list)") 00626 00627 if not overlap: 00628 msg("No overlap between model and pointings",priority="error") 00629 return False 00630 00631 00632 00633 ################################################################## 00634 # calibrator is not explicitly contained in the pointing file 00635 # but interleaved with etime=intergration 00636 util.isquantity(calflux) 00637 calfluxjy = qa.convert(calflux,'Jy')['value'] 00638 # XML returns a list even for a string: 00639 if type(caldirection) == type([]): caldirection = caldirection[0] 00640 if len(caldirection) < 4: caldirection = "" 00641 if calfluxjy > 0 and caldirection != "" and uvmode: 00642 docalibrator = True 00643 if os.path.exists(fileroot+"/"+project+'.cal.cclist'): 00644 shutil.rmtree(fileroot+"/"+project+'.cal.cclist') 00645 util.isdirection(caldirection) 00646 cl.done() 00647 cl.addcomponent(flux=calfluxjy,dir=caldirection,label="phase calibrator") 00648 # set reference freq to center freq of model 00649 cl.rename(fileroot+"/"+project+'.cal.cclist') 00650 cl.done() 00651 else: 00652 docalibrator = False 00653 00654 00655 00656 00657 00658 ################################################################## 00659 # create one figure for model and pointings - need antenna diam 00660 # to determine primary beam 00661 if grfile: 00662 file = fileroot + "/" + project + ".skymodel.png" 00663 else: 00664 file = "" 00665 00666 if grscreen or grfile: 00667 util.newfig(show=grscreen) 00668 00669 # remove after we fix the scaling algorithm for the images 00670 if components_only: 00671 pl.plot() 00672 # TODO add symbols at locations of components 00673 pl.plot(coffs[0,]*3600,coffs[1,]*3600,'o',c="#dddd66") 00674 pl.axis("equal") 00675 else: 00676 discard = util.statim(modelflat,plot=True,incell=model_cell) 00677 00678 lims = pl.xlim(),pl.ylim() 00679 if pb <= 0 and verbose: 00680 msg("unknown primary beam size for plot",priority="warn") 00681 if max(max(lims)) > pb and not components_only: 00682 plotcolor = 'w' 00683 else: 00684 plotcolor = 'k' 00685 00686 #if offsets.shape[1] > 16 or pb <= 0 or pb > pl.absolute(max(max(lims))): 00687 if offsets.shape[1] > 19 or pb <= 0: 00688 lims = pl.xlim(),pl.ylim() 00689 pl.plot((offsets[0]+shift[0])*3600.,(offsets[1]+shift[1])*3600., 00690 plotcolor+'+',markeredgewidth=1) 00691 #if pb > 0 and pl.absolute(lims[0][0]) > pb: 00692 if pb > 0: 00693 plotpb(pb,pl.gca(),lims=lims,color=plotcolor) 00694 else: 00695 from matplotlib.patches import Circle 00696 for i in xrange(offsets.shape[1]): 00697 pl.gca().add_artist(Circle( 00698 ((offsets[0,i]+shift[0])*3600, 00699 (offsets[1,i]+shift[1])*3600), 00700 radius=pb/2.,edgecolor=plotcolor,fill=False, 00701 label='beam',transform=pl.gca().transData,clip_on=True)) 00702 00703 xlim = max(abs(pl.array(lims[0]))) 00704 ylim = max(abs(pl.array(lims[1]))) 00705 # show entire pb: (statim doesn't by default) 00706 pl.xlim([max([xlim,pb/2]),min([-xlim,-pb/2])]) 00707 pl.ylim([min([-ylim,-pb/2]),max([ylim,pb/2])]) 00708 pl.xlabel("resized model sky",fontsize="x-small") 00709 util.endfig(show=grscreen,filename=file) 00710 00711 00712 00713 00714 00715 00716 00717 00718 ################################################################## 00719 # set up observatory, feeds, etc 00720 quickpsf_current = False 00721 00722 if uvmode: 00723 msfile = fileroot + "/" + project + '.ms' 00724 else: 00725 msfile = fileroot + "/" + project + '.sd.ms' 00726 00727 if predict: 00728 # TODO check for frequency overlap here - if zero stop 00729 # position overlap already checked above in pointing section 00730 00731 message = "preparing empty measurement set" 00732 if verbose: 00733 msg(message,origin="simobserve",priority="warn") 00734 else: 00735 msg(message,origin="simobserve") 00736 00737 nbands = 1; 00738 fband = util.bandname(qa.convert(model_center, 'GHz')['value']) 00739 00740 ############################################ 00741 # predict observation 00742 00743 # if someone has the old style refdate with the included, discard 00744 q = re.compile('(\d*/\d+/\d+)([/:\d]*)') 00745 qq = q.match(refdate) 00746 if not qq: 00747 msg("Invalid reference date "+refdate,priority="error") 00748 return False 00749 else: 00750 z = qq.groups() 00751 refdate=z[0] 00752 if len(z)>1: 00753 msg("Discarding time part of refdate, "+z[1]+", in favor of hourangle parameter = "+hourangle) 00754 00755 if hourangle=="transit": 00756 haoffset=0.0 00757 else: 00758 haoffset=qa.convert(qa.quantity(hourangle),'s')['value'] 00759 00760 refdate=refdate+"/00:00:00" 00761 usehourangle=True 00762 00763 00764 intsec = qa.convert(qa.quantity(integration),"s")['value'] 00765 00766 # totaltime as an integer for # times through the mosaic: 00767 if not util.isquantity(totaltime): 00768 msg("total time "+totaltime+" does not appear to represent a time interval (use 's','min','h'; not 'sec','m','hr')",priority="error") 00769 return False 00770 00771 if qa.quantity(totaltime)['value'] < 0.: 00772 # casapy crashes for negative totaltime 00773 msg("Negative totaltime is not allowed.",priority="error") 00774 return False 00775 if qa.quantity(totaltime)['unit'] == '': 00776 # assume it means number of maps, or # repetitions. 00777 totalsec = sum(etime) 00778 if docalibrator: 00779 totalsec = totalsec + intsec # cal gets one int-time 00780 totalsec = float(totaltime) * totalsec 00781 msg("Total observing time = "+str(totalsec)+"s.",priority="warn") 00782 else: 00783 if not qa.compare(totaltime,"1s"): 00784 msg("total time "+totaltime+" does not appear to represent a time interval (use 's','min','h'; not 'sec','m','hr')",priority="error") 00785 return False 00786 totalsec = qa.convert(qa.quantity(totaltime),'s')['value'] 00787 00788 if os.path.exists(msfile) and not overwrite: #redundant check? 00789 util.msg("measurement set "+msfile+" already exists and user does not wish to overwrite",priority="error") 00790 return False 00791 sm.open(msfile) 00792 posobs = me.observatory(telescopename) 00793 diam = stnd; 00794 # WARNING: sm.setspwindow is not consistent with clean::center 00795 #model_start=qa.sub(model_center,qa.mul(model_width,0.5*model_nchan)) 00796 # but the "start" is the center of the first channel: 00797 model_start = qa.sub(model_center,qa.mul(model_width,0.5*(model_nchan-1))) 00798 00799 mounttype = 'alt-az' 00800 if telescopename in ['DRAO', 'WSRT']: 00801 mounttype = 'EQUATORIAL' 00802 # Should ASKAP be BIZARRE or something else? It may be effectively equatorial. 00803 00804 sm.setconfig(telescopename=telescopename, x=stnx, y=stny, z=stnz, 00805 dishdiameter=diam.tolist(), 00806 mount=[mounttype], antname=antnames, padname=padnames, 00807 coordsystem='global', referencelocation=posobs) 00808 if str.upper(telescopename).find('VLA') > 0: 00809 sm.setspwindow(spwname=fband, freq=qa.tos(model_start), 00810 deltafreq=qa.tos(model_width), 00811 freqresolution=qa.tos(model_width), 00812 nchannels=model_nchan, refcode="LSRK", 00813 stokes='RR LL') 00814 sm.setfeed(mode='perfect R L',pol=['']) 00815 else: 00816 sm.setspwindow(spwname=fband, freq=qa.tos(model_start), 00817 deltafreq=qa.tos(model_width), 00818 freqresolution=qa.tos(model_width), 00819 nchannels=model_nchan, refcode="LSRK", 00820 stokes='XX YY') 00821 sm.setfeed(mode='perfect X Y',pol=['']) 00822 00823 if verbose: msg(" spectral window set at %s" % qa.tos(model_center)) 00824 sm.setlimits(shadowlimit=0.01, elevationlimit='10deg') 00825 if uvmode: 00826 sm.setauto(0.0) 00827 else: #Single-dish 00828 # auto-correlation should be unity for single dish obs. 00829 sm.setauto(1.0) 00830 00831 for k in xrange(0,nfld): 00832 src = project + '_%d' % k 00833 sm.setfield(sourcename=src, sourcedirection=pointings[k], 00834 calcode="OBJ", distance='0m') 00835 if k == 0: 00836 sourcefieldlist = src 00837 else: 00838 sourcefieldlist = sourcefieldlist + ',' + src 00839 if docalibrator: 00840 sm.setfield(sourcename="phase calibrator", 00841 sourcedirection=caldirection,calcode='C', 00842 distance='0m') 00843 00844 mereftime = me.epoch('TAI', refdate) 00845 # integration is a scalar quantity, etime is a vector of seconds 00846 sm.settimes(integrationtime=integration, usehourangle=usehourangle, 00847 referencetime=mereftime) 00848 00849 # time required to observe all planned scanes in etime array: 00850 totalscansec = sum(etime) 00851 kfld = 0 00852 00853 if totalsec < totalscansec: 00854 msg("Not all pointings in the mosaic will be observed - check mosaic setup and exposure time parameters!",priority="warn") 00855 ### 00856 #print "you need at least %16.12e sec but you have %16.12e sec (%f < %f = %s)" % (totalscansec, totalsec, totalsec, totalscansec, str(totalsec<totalscansec)) 00857 ### 00858 00859 # sm.observemany 00860 srces = [] 00861 starttimes = [] 00862 stoptimes = [] 00863 dirs = [] 00864 00865 if usehourangle: 00866 sttime = -totalsec/2.0 00867 else: 00868 sttime = 0. # leave start at the reftime 00869 sttime=sttime+haoffset 00870 scanstart=sttime 00871 00872 # can before sources 00873 if docalibrator: 00874 endtime = sttime + qa.convert(integration,'s')['value'] 00875 sm.observe(sourcename="phase calibrator", spwname=fband, 00876 starttime=qa.quantity(sttime, "s"), 00877 stoptime=qa.quantity(endtime, "s"), 00878 state_obs_mode="CALIBRATE_PHASE.ON_SOURCE",state_sig=True, 00879 project=project); 00880 sttime = endtime 00881 00882 while (sttime-scanstart) < totalsec: # the last scan could exceed totaltime 00883 endtime = sttime + etime[kfld] 00884 src = project + '_%d' % kfld 00885 srces.append(src) 00886 starttimes.append(str(sttime)+"s") 00887 stoptimes.append(str(endtime)+"s") 00888 dirs.append(pointings[kfld]) 00889 kfld = kfld + 1 00890 # advance start time - XX someday slew goes here 00891 sttime = endtime 00892 00893 if kfld == nfld: 00894 if docalibrator: 00895 endtime = sttime + qa.convert(integration,'s')['value'] 00896 00897 # need to observe cal singly to get new row in obs table, so 00898 # first observemany the on-source pointing(s) 00899 sm.observemany(sourcenames=srces,spwname=fband,starttimes=starttimes,stoptimes=stoptimes,project=project) 00900 # and clear the list 00901 srces = [] 00902 starttimes = [] 00903 stoptimes = [] 00904 dirs = [] 00905 sm.observe(sourcename="phase calibrator", spwname=fband, 00906 starttime=qa.quantity(sttime, "s"), 00907 stoptime=qa.quantity(endtime, "s"), 00908 state_obs_mode="CALIBRATE_PHASE.ON_SOURCE",state_sig=True, 00909 project=project); 00910 kfld = kfld + 1 00911 sttime = endtime 00912 # if kfld > nfld: kfld = 0 00913 if kfld > nfld-1: kfld = 0 00914 # if directions is unset, NewMSSimulator::observemany 00915 00916 # looks up the direction in the field table. 00917 if not docalibrator: 00918 sm.observemany(sourcenames=srces,spwname=fband,starttimes=starttimes,stoptimes=stoptimes,project=project) 00919 00920 sm.setdata(fieldid=range(0,nfld)) 00921 sm.setvp() 00922 00923 msg("done setting up observations (blank visibilities)") 00924 if verbose: sm.summary() 00925 00926 # do actual calculation of visibilities: 00927 00928 if not uvmode: #Single-dish 00929 sm.setoptions(gridfunction='pb', ftmachine="sd", location=posobs) 00930 00931 if not components_only: 00932 if docalibrator: 00933 if len(complist) <=0: 00934 complist=fileroot+"/"+project+'.cal.cclist' 00935 else: 00936 # XXX will 2 cl work? 00937 complist=complist+","+fileroot+"/"+project+'.cal.cclist' 00938 00939 if len(complist) > 1: 00940 message = "predicting from "+newmodel+" and "+complist 00941 if verbose: 00942 msg(message,priority="warn",origin="simobserve") 00943 else: 00944 msg(message,origin="simobserve") 00945 else: 00946 message = "predicting from "+newmodel 00947 if verbose: 00948 msg(message,priority="warn",origin="simobserve") 00949 else: 00950 msg(message,origin="simobserve") 00951 sm.predict(imagename=newmodel,complist=complist) 00952 else: # if we're doing only components 00953 # XXX will 2 cl work? 00954 if docalibrator: 00955 complist=complist+","+fileroot+"/"+project+'.cal.cclist' 00956 if verbose: 00957 msg("predicting from "+complist,priority="warn",origin="simobserve") 00958 else: 00959 msg("predicting from "+complist,origin="simobserve") 00960 sm.predict(complist=complist) 00961 00962 sm.done() 00963 msg('generation of measurement set '+msfile+' complete',origin="simobserve") 00964 00965 00966 ############################################ 00967 # create figure 00968 if grfile: 00969 file = fileroot + "/" + project + ".observe.png" 00970 else: 00971 file = "" 00972 00973 # update psfsize using uv coverage instead of maxbase above 00974 if os.path.exists(msfile): 00975 # psfsize was set from the antenna posns before, but uv is better 00976 tb.open(msfile) 00977 rawdata = tb.getcol("UVW") 00978 tb.done() 00979 maxbase = max([max(rawdata[0,]),max(rawdata[1,])]) # in m 00980 if maxbase > 0.: 00981 psfsize = 0.3/qa.convert(qa.quantity(model_center),'GHz')['value']/maxbase*3600.*180/pl.pi # lambda/b converted to arcsec 00982 if not uvmode: 00983 msg("An single observation is requested but CROSS-correlation in "+msfile,priority="error") 00984 else: #Single-dish (zero-spacing) 00985 psfsize = pb 00986 if uvmode: 00987 msg("An interferometer observation is requested but only AUTO-correlation in "+msfile,priority="error") 00988 minimsize = 8* int(psfsize/cell_asec) 00989 else: 00990 msg("Couldn't find "+msfile,priority="error") 00991 00992 if uvmode: 00993 multi = [2,2,1] 00994 else: 00995 multi = 0 00996 00997 00998 if (grscreen or grfile): 00999 util.newfig(multi=multi,show=grscreen) 01000 util.ephemeris(refdate,direction=util.direction,telescope=telescopename,ms=msfile,usehourangle=usehourangle) 01001 casalog.origin('simobserve') 01002 if uvmode: 01003 util.nextfig() 01004 util.plotants(stnx, stny, stnz, stnd, padnames) 01005 01006 # uv coverage 01007 util.nextfig() 01008 klam_m = 300/qa.convert(model_center,'GHz')['value'] 01009 pl.box() 01010 pl.plot(rawdata[0,]/klam_m,rawdata[1,]/klam_m,'b,') 01011 pl.plot(-rawdata[0,]/klam_m,-rawdata[1,]/klam_m,'b,') 01012 ax = pl.gca() 01013 ax.yaxis.LABELPAD = -4 01014 pl.xlabel('u[klambda]',fontsize='x-small') 01015 pl.ylabel('v[klambda]',fontsize='x-small') 01016 pl.axis('equal') 01017 01018 # show dirty beam from observed uv coverage 01019 util.nextfig() 01020 im.open(msfile) 01021 # TODO spectral parms 01022 msg("using default model cell "+qa.tos(model_cell[0])+" for PSF calculation",priority="warn") 01023 im.defineimage(cellx=qa.tos(model_cell[0]),nx=int(max([minimsize,128]))) 01024 if os.path.exists(fileroot+"/"+project+".quick.psf"): 01025 shutil.rmtree(fileroot+"/"+project+".quick.psf") 01026 im.approximatepsf(psf=fileroot+"/"+project+".quick.psf") 01027 quickpsf_current = True 01028 beam = im.fitpsf(psf=fileroot+"/"+project+".quick.psf") 01029 im.done() 01030 ia.open(fileroot+"/"+project+".quick.psf") 01031 beamcs = ia.coordsys() 01032 beam_array = ia.getchunk(axes=[beamcs.findcoordinate("spectral")[1][0],beamcs.findcoordinate("stokes")[1][0]],dropdeg=True) 01033 nn = beam_array.shape 01034 xextent = nn[0]*cell_asec*0.5 01035 xextent = [xextent,-xextent] 01036 yextent = nn[1]*cell_asec*0.5 01037 yextent = [-yextent,yextent] 01038 flipped_array = beam_array.transpose() 01039 ttrans_array = flipped_array.tolist() 01040 ttrans_array.reverse() 01041 pl.imshow(ttrans_array,interpolation='bilinear',cmap=pl.cm.jet,extent=xextent+yextent,origin="bottom") 01042 pl.title(project+".quick.psf",fontsize="x-small") 01043 b = qa.convert(beam[1],'arcsec')['value'] 01044 pl.xlim([-3*b,3*b]) 01045 pl.ylim([-3*b,3*b]) 01046 ax = pl.gca() 01047 pl.text(0.05,0.95,"bmaj=%7.1e\nbmin=%7.1e" % (beam[1]['value'],beam[2]['value']),transform = ax.transAxes,bbox=dict(facecolor='white', alpha=0.7),size="x-small",verticalalignment="top") 01048 ia.done() 01049 util.endfig(show=grscreen,filename=file) 01050 01051 elif thermalnoise != "" or leakage > 0: 01052 # Not predicting this time but corrupting. Get obsmode from ms. 01053 if not os.path.exists(msfile): 01054 msg("Couldn't find "+msfile,priority="error") 01055 uvmode = (not util.ismstp(msfile,halt=False)) 01056 01057 01058 ###################################################################### 01059 # noisify 01060 01061 noise_any = False 01062 msroot = fileroot + "/" + project # if leakage, can just copy from this project 01063 01064 if thermalnoise != "": 01065 knowntelescopes = ["ALMA", "ACA", "SMA", "EVLA", "VLA"] 01066 01067 noise_any = True 01068 01069 noisymsroot = msroot + ".noisy" 01070 if not uvmode: #Single-dish 01071 msroot += ".sd" 01072 noisymsroot += ".sd" 01073 01074 # Cosmic background radiation temperature in K. 01075 t_cmb = 2.725 01076 01077 01078 # check for interferometric ms: 01079 if not os.path.exists(msroot+".ms"): 01080 msg("Couldn't find "+msroot+".ms",priority="error") 01081 01082 # MS exists 01083 message = 'copying '+msroot+'.ms to ' + \ 01084 noisymsroot+'.ms and adding thermal noise' 01085 if verbose: 01086 msg(message,origin="noise",priority="warn") 01087 else: 01088 msg(message,origin="noise") 01089 01090 if os.path.exists(noisymsroot+".ms"): 01091 shutil.rmtree(noisymsroot+".ms") 01092 shutil.copytree(msfile,noisymsroot+".ms") 01093 if sm.name() != '': 01094 msg("table persistence error on %s" % sm.name(),priority="error") 01095 return False 01096 01097 # if not predicted this time, get telescopename from ms 01098 if not predict: 01099 tb.open(noisymsroot+".ms/OBSERVATION") 01100 n = tb.getcol("TELESCOPE_NAME") 01101 telescopename = n[0] 01102 # todo add check that entire column is the same 01103 tb.done() 01104 msg("telescopename read from "+noisymsroot+".ms: "+telescopename) 01105 01106 if telescopename not in knowntelescopes: 01107 msg("thermal noise only works properly for ALMA/ACA or EVLA",origin="noise",priority="warn") 01108 eta_p, eta_s, eta_b, eta_t, eta_q, t_rx = util.noisetemp(telescope=telescopename,freq=model_center) 01109 01110 # antenna efficiency 01111 eta_a = eta_p * eta_s * eta_b * eta_t 01112 if verbose: 01113 msg('antenna efficiency = '+str(eta_a), origin="noise") 01114 msg('spillover efficiency = '+str(eta_s), origin="noise") 01115 msg('correlator efficiency = '+str(eta_q), origin="noise") 01116 # sensitivity constant 01117 scoeff = -1 #Force setting the default value, 1./sqrt(2.0) 01118 if not uvmode: #Single-dish 01119 scoeff = 1.0 01120 if verbose: msg('sensitivity constant = '+str(scoeff), origin="noise") 01121 01122 sm.openfromms(noisymsroot+".ms") # an existing MS 01123 sm.setdata(fieldid=[]) # force to get all fields 01124 sm.setseed(seed) 01125 if thermalnoise == "tsys-manual": 01126 if verbose: 01127 message = "sm.setnoise(spillefficiency="+str(eta_s)+\ 01128 ",correfficiency="+str(eta_q)+\ 01129 ",antefficiency="+str(eta_a)+\ 01130 ",trx="+str(t_rx)+",tau="+str(tau0)+\ 01131 ",tatmos="+str(t_sky)+",tground="+str(t_ground)+\ 01132 ",tcmb="+str(t_cmb) 01133 if not uvmode: message += ",senscoeff="+str(scoeff) 01134 message += ",mode='tsys-manual')" 01135 msg(message); 01136 msg("** this may be slow if your MS is finely sampled in time ** ",priority="warn") 01137 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q, 01138 antefficiency=eta_a,trx=t_rx, 01139 tau=tau0,tatmos=t_sky,tground=t_ground,tcmb=t_cmb, 01140 mode="tsys-manual",senscoeff=scoeff) 01141 else: 01142 if verbose: 01143 message = "sm.setnoise(spillefficiency="+str(eta_s)+\ 01144 ",correfficiency="+str(eta_q)+\ 01145 ",antefficiency="+str(eta_a)+\ 01146 ",trx="+str(t_rx)+",tground="+str(t_ground)+\ 01147 ",tcmb="+str(t_cmb) 01148 if not uvmode: message += ",senscoeff="+str(scoeff) 01149 message += ",mode='tsys-atm'"+\ 01150 ",pground='560mbar',altitude='5000m'"+\ 01151 ",waterheight='2km',relhum=20,pwv="+str(user_pwv)+"mm)" 01152 msg(message); 01153 msg("** this may be slow if your MS is finely sampled in time ** ",priority="warn") 01154 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q, 01155 antefficiency=eta_a,trx=t_rx, 01156 tground=t_ground,tcmb=t_cmb,pwv=str(user_pwv)+"mm", 01157 mode="tsys-atm",table=noisymsroot,senscoeff=scoeff) 01158 # don't set table, that way it won't save to disk 01159 # mode="calculate",table=noisymsroot) 01160 01161 sm.corrupt(); 01162 sm.done(); 01163 01164 01165 msroot = noisymsroot 01166 if verbose: msg("done corrupting with thermal noise",origin="noise") 01167 01168 01169 if leakage > 0: 01170 noise_any = True 01171 # TODO: need to handle SD name when leakage is available 01172 if msroot == fileroot+"/"+project: 01173 noisymsroot = fileroot + "/" + project + ".noisy" 01174 else: 01175 noisymsroot = fileroot + "/" + project + ".noisier" 01176 if not uvmode: #Single-dish 01177 msg("Can't corrupt SD data with polarization leakage",priority="warn") 01178 if os.path.exists(msfile): 01179 msg('copying '+msfile+' to ' + 01180 noisymsroot+'.ms and adding polarization leakage', 01181 origin="noise",priority="warn") 01182 if os.path.exists(noisymsroot+".ms"): 01183 shutil.rmtree(noisymsroot+".ms") 01184 shutil.copytree(msfile,noisymsroot+".ms") 01185 if sm.name() != '': 01186 msg("table persistence error on %s" % sm.name(),priority="error") 01187 return False 01188 01189 sm.openfromms(noisymsroot+".ms") # an existing MS 01190 sm.setdata(fieldid=[]) # force to get all fields 01191 sm.setleakage(amplitude=leakage,table=noisymsroot+".cal") 01192 sm.corrupt(); 01193 sm.done(); 01194 01195 01196 01197 # cleanup - delete newmodel, newmodel.flat etc 01198 # if os.path.exists(modelflat): 01199 # shutil.rmtree(modelflat) 01200 if os.path.exists(modelflat+".regrid"): 01201 shutil.rmtree(modelflat+".regrid") 01202 if os.path.exists(fileroot+"/"+project+".noisy.T.cal"): 01203 shutil.rmtree(fileroot+"/"+project+".noisy.T.cal") 01204 if os.path.exists(fileroot+"/"+project+".noisy.sd.T.cal"): 01205 shutil.rmtree(fileroot+"/"+project+".noisy.sd.T.cal") 01206 01207 except TypeError, e: 01208 finalize_tools() 01209 #msg("simobserve -- TypeError: %s" % e, priority="error") 01210 casalog.post("simobserve -- TypeError: %s" % e, priority="ERROR") 01211 raise TypeError, e 01212 return False 01213 except ValueError, e: 01214 finalize_tools() 01215 #print "task_simobserve -- OptionError: ", e 01216 #msg("simobserve -- OptionError: %s" % e, priority="error") 01217 casalog.post("simobserve -- OptionError: %s" % e, priority="ERROR") 01218 raise ValueError, e 01219 return False 01220 except Exception, instance: 01221 finalize_tools() 01222 #print '***Error***',instance 01223 #msg("simobserve -- Exception: %s" % instance, priority="error") 01224 casalog.post("simobserve -- Exception: %s" % instance, priority="ERROR") 01225 raise Exception, instance 01226 return False 01227 return True 01228 01229 ##### Helper functions to plot primary beam 01230 def plotpb(pb,axes,lims=None,color='k'): 01231 # This beam is automatically scaled when you zoom in/out but 01232 # not anchored in plot area. We'll wait for Matplotlib 0.99 01233 # for that function. 01234 #major=major 01235 #minor=minor 01236 #rangle=rangle 01237 #bwidth=max(major*pl.cos(rangle),minor*pl.sin(rangle))*1.1 01238 #bheight=max(major*pl.sin(rangle),minor*pl.cos(rangle))*1.1 01239 from matplotlib.patches import Rectangle, Circle #,Ellipse 01240 try: 01241 from matplotlib.offsetbox import AnchoredOffsetbox, AuxTransformBox 01242 box = AuxTransformBox(axes.transData) 01243 box.set_alpha(0.7) 01244 circ = Circle((pb,pb),radius=pb/2.,color=color,fill=False,\ 01245 label='primary beam',linewidth=2.0) 01246 box.add_artist(circ) 01247 pblegend = AnchoredOffsetbox(loc=3,pad=0.2,borderpad=0.,\ 01248 child=box,prop=None,frameon=False)#,frameon=True) 01249 pblegend.set_alpha(0.7) 01250 axes.add_artist(pblegend) 01251 except: 01252 print "Using old matplotlib substituting with circle" 01253 # work around for old matplotlib 01254 boxsize = pb*1.1 01255 if not lims: lims = axes.get_xlim(),axes.get_ylim() 01256 incx = 1 01257 incy = 1 01258 if axes.xaxis_inverted(): incx = -1 01259 if axes.yaxis_inverted(): incy = -1 01260 #ecx = lims[0][0] + bwidth/2.*incx 01261 #ecy = lims[1][0] + bheight/2.*incy 01262 ccx = lims[0][0] + boxsize/2.*incx 01263 ccy = lims[1][0] + boxsize/2.*incy 01264 01265 #box = Rectangle((lims[0][0],lims[1][0]),incx*bwidth,incy*bheight, 01266 box = Rectangle((lims[0][0],lims[1][0]),incx*boxsize,incy*boxsize, 01267 alpha=0.7,facecolor='w', 01268 transform=axes.transData) #Axes 01269 #beam = Ellipse((ecx,ecy),major,minor,angle=rangle, 01270 beam = Circle((ccx,ccy), radius=pb/2., 01271 edgecolor='k',fill=False, 01272 label='beam',transform=axes.transData) 01273 #props = {'pad': 3, 'edgecolor': 'k', 'linewidth':2, 'facecolor': 'w', 'alpha': 0.5} 01274 #pl.matplotlib.patches.bbox_artist(beam,axes.figure.canvas.get_renderer(),props=props) 01275 axes.add_artist(box) 01276 axes.add_artist(beam) 01277 01278 def finalize_tools(): 01279 if ia.isopen(): ia.close() 01280 sm.close() 01281 #cl.close() # crashes casa
1.8.0
