Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_simobserve.py: 1%
835 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 18:48 +0000
1import os
2import re
3import pylab as pl
5from casatools import ctsys, quanta, imager
6from casatasks import casalog
7from .simutil import *
8from .simutil import is_array_type
10qa = quanta()
12def simobserve(
13 project=None,
14 skymodel=None, inbright=None, indirection=None, incell=None,
15 incenter=None, inwidth=None, # innchan=None,
16 complist=None, compwidth=None, comp_nchan=1,
17 setpointings=None,
18 ptgfile=None, integration=None, direction=None, mapsize=None,
19 maptype=None, pointingspacing=None, caldirection=None, calflux=None,
20 # observe=None,
21 obsmode=None,
22 refdate=None, hourangle=None,
23 totaltime=None, antennalist=None,
24 sdantlist=None,
25 sdant=None,
26 outframe=None,
27 thermalnoise=None,
28 user_pwv=None, t_ground=None, t_sky=None, tau0=None, seed=None,
29 leakage=None,
30 graphics=None,
31 verbose=None,
32 overwrite=None
33 ):
35 # Collect a list of parameter values to save inputs
36 in_params = locals()
38 try:
40 #########################
41 # some hardcoded variables
42 pbcoeff = 1.13 ## PB defined as pbcoeff*lambda/d
43 nyquist = 0.5/pbcoeff ## Nyquist spacing = PB*nyquist
45 gridratio_int = 1./pl.sqrt(3) # times lambda/d
46 gridratio_tp = 1./3
48 relmargin = .5 # number of PB between edge of model and ptg centers
49 scanlength = 1 # number of integrations per scan
52 # RI TODO for inbright=unchanged, need to scale input image to jy/pix
53 # according to actual units in the input image
55 casalog.origin('simobserve')
56 if verbose: casalog.filter(level="DEBUG2")
58 # create the utility object
59 # this is the dir of the observation (could be "")
60 util = simutil(direction)
61 if verbose: util.verbose = True
62 msg = util.msg
64 # it was requested to make the user interface "observe" for what
65 # is sm.observe and sm.predict.
66 # interally the code is clearer if we stick with predict so
67 predict = obsmode.startswith('i') or obsmode.startswith('s')
68 if predict:
69 uvmode = obsmode.startswith('i')
70 if not uvmode: antennalist = sdantlist
71 elif sdantlist != "":
72 if antennalist == "":
73 uvmode = False
74 antennalist = sdantlist
75 else:
76 #uvmode = True
77 #msg("Both antennalist and sdantlist are defined. sdantlist will be ignored",priority="warn")
78 emsg = "Both antennalist and sdantlist are defined. Define one of them."
79 raise ValueError(emsg)
80 else:
81 uvmode = True
82 # uvmode = (sdant < 0) #when flexible default values come available
85 # put output in directory called "project"
86 fileroot = project
87 if not os.path.exists(fileroot):
88 os.mkdir(fileroot)
91 # filename parsing of cfg file here so that the project filenames
92 # can contain the cfg
93 repodir = ctsys.resolve("alma/simmos")
96 # convert "alma;0.4arcsec" to an actual configuration
97 # can only be done after reading skymodel, so here, we just string parse
98 if str.upper(antennalist[0:4]) == "ALMA":
99 foo=antennalist.replace(';','_')
100 elif len(antennalist) > 0:
101 foo=antennalist
102 else:
103 msg("The name of antenna list (antennalist/sdantlist) is not specified",priority="error")
105 if foo:
106 foo=foo.replace(".cfg","")
107 sfoo=foo.split('/')
108 if len(sfoo)>1:
109 foo=sfoo[-1]
110 project=project+"."+foo
113 if not overwrite:
114 if (predict and uvmode and os.path.exists(fileroot+"/"+
115 project+".ms")):
116 emsg = fileroot+"/"+project+".ms exists but overwrite=F"
117 raise RuntimeError(emsg)
118 if (predict and (not uvmode) and os.path.exists(fileroot+"/"+
119 project+".sd.ms")):
120 emsg = fileroot+"/"+project+".sd.ms exists but overwrite=F"
121 raise RuntimeError(emsg)
124 if is_array_type(skymodel):
125 skymodel = skymodel[0]
126 skymodel = skymodel.replace('$project',project)
128 if is_array_type(complist):
129 complist = complist[0]
131 if((not os.path.exists(skymodel)) and (not os.path.exists(complist))):
132 if len(skymodel)>0:
133 msg("Your skymodel '"+skymodel+"' could not be found.",
134 priority="warn")
135 if len(complist)>0:
136 msg("Your complist '"+complist+"' could not be found.",
137 priority="warn")
138 if len(skymodel)==0 and len(complist)==0:
139 msg("At least one of skymodel or complist must be set.",
140 priority="error")
142 else:
143 msg("No sky input found."+
144 " At least one of skymodel or complist must exist.",
145 priority="error")
147 ### WORKAROUND for wrong flux in COMP TP simulations (CAS-5095)
148 if (obsmode.startswith("s") and os.path.exists(complist)):
149 emsg = "Single dish simulation has a flux recovery issue when using a components list as an input.\nPlease generate compskymodel image first by obsmode='' and use the image as the skymodel input.\nSorry for the inconvenience."
150 raise RuntimeError(emsg)
151 ### End of WORKAROUND
153 grscreen = False
154 grfile = False
155 if graphics == "both":
156 grscreen = True
157 grfile = True
158 if graphics == "screen":
159 grscreen = True
160 if graphics == "file":
161 grfile = True
163 ##################################################################
164 # set up skymodel image
166 if os.path.exists(skymodel):
167 components_only = False
168 # create a new skymodel called skymodel,
169 # or if it's already there, called newmodel
170 default_model = project + ".skymodel"
171 if skymodel == default_model:
172 newmodel = fileroot + "/" + project + ".newmodel"
173 else:
174 newmodel = fileroot + "/" + default_model
175 if os.path.exists(newmodel):
176 if overwrite:
177 shutil.rmtree(newmodel)
178 else:
179 emsg = (newmodel+" exists -- "+
180 "please delete it, change skymodel, or set overwrite=T")
181 raise RuntimeError(emsg)
183 # modifymodel just collects info if skymodel==newmodel
184 innchan = -1
185 returnpars = util.modifymodel(skymodel,newmodel,
186 inbright,indirection,incell,
187 incenter,inwidth,innchan,
188 flatimage=False)
189 if not returnpars:
190 raise RuntimeError('Failure in modifymodel. It returned: {}'.
191 format(returnpars))
193 (model_refdir,model_cell,model_size,
194 model_nchan,model_specrefval,model_specrefpix,model_width,
195 model_stokes) = returnpars
197 modelflat = fileroot + "/" + project + ".skymodel.flat"
198 if os.path.exists(modelflat) and (not predict):
199 # if we're not predicting, then we want to use the previously
200 # created modelflat, because it may have components added
201 msg("flat sky model "+modelflat+
202 " exists, predict not requested",priority="warn")
203 msg(" working from existing model image - "+
204 "please delete it if you wish to overwrite.",
205 priority="warn")
206 else:
207 # create and add components into modelflat with util.flatimage()
208 util.flatimage(newmodel,complist=complist,verbose=verbose)
209 # we want skymodel.flat image to be called that no matter what
210 # the skymodel image is called, since it's used in analysis
211 if modelflat != newmodel+".flat":
212 if os.path.exists(modelflat):
213 shutil.rmtree(modelflat)
214 shutil.move(newmodel+".flat",modelflat)
216 casalog.origin('simobserve')
218 # set startfeq and bandwidth in util object after modifymodel
219 bandwidth = qa.mul(qa.quantity(model_nchan),
220 qa.quantity(model_width))
221 util.bandwidth = bandwidth
223 else:
224 components_only = True
225 msg("component-only simulation",priority="info")
226 # calculate model parameters from the component list:
228 compdirs = []
229 cl.done()
230 cl.open(complist)
232 for i in range(cl.length()):
233 compdirs.append(util.dir_m2s(cl.getrefdir(i)))
235 model_refdir, coffs = util.average_direction(compdirs)
236 model_specrefval = cl.getspectrum(0)['frequency']['m0']
238 if util.isquantity(compwidth,halt=False):
239 model_width = compwidth
240 msg("compwidth set: setting model bandwidth to input",
241 priority="info")
242 else:
243 model_width = "2GHz"
244 msg("compwidth unset: setting bandwidth to 2GHz",
245 priority="warn")
247 model_nchan = comp_nchan
248 # channelize component-only MS
249 # currently assuming equal width and center as frequency reference
250 model_specrefpix = 0.5*(comp_nchan-1)
251 msg("scaling model bandwidth by model_nchan",priority="info")
252 model_width = qa.div(model_width,model_nchan)
253 model_stokes = "I"
255 cmax = 0.0014 # ~5 arcsec
256 for i in range(coffs.shape[1]):
257 xc = pl.absolute(coffs[0,i]) # offsets in deg
258 yc = pl.absolute(coffs[1,i])
259 if xc > cmax:
260 cmax = xc
261 if yc > cmax:
262 cmax = yc
264 model_size = ["%fdeg" % (3*cmax), "%fdeg" % (3*cmax)]
267 # for cases either if there is a skymodel or are only components,
268 # if user has not input a map size (for setpointings), use model_size
269 if len(mapsize) == 0:
270 mapsize = model_size
271 if verbose: msg("setting map size to "+str(model_size),
272 origin='simobserve')
273 else:
274 if is_array_type(mapsize):
275 if len(mapsize[0]) == 0:
276 mapsize = model_size
277 if verbose: msg("setting map size to "+str(model_size),
278 origin="simobserve")
280 if components_only:
281 if is_array_type(mapsize):
282 map_asec = qa.convert(mapsize[0],"arcsec")['value']
283 else:
284 map_asec = qa.convert(mapsize,"arcsec")['value']
285 if is_array_type(model_size):
286 mod_asec = qa.convert(model_size[0],"arcsec")['value']
287 else:
288 mod_asec = qa.convert(model_size,"arcsec")['value']
289 if map_asec>mod_asec:
290 model_size=["%farcsec" % map_asec,"%farcsec" % map_asec]
295 ##################################################################
296 # read antenna file here to get Primary Beam, and estimate psfsize
297 aveant = -1
298 stnx = [] # for later, to know if we read an array in or not
299 pb = 0. # primary beam
301 # convert "alma;0.4arcsec" to an actual configuration
302 if str.upper(antennalist[0:4]) == "ALMA":
304 resparsed=False
305 # test for cycle 1
306 q = re.compile('.*CYCLE.?1.?;(.*)')
307 qq = q.match(antennalist.upper())
308 if qq:
309 z = qq.groups()
310 tail=z[0]
311 tail=tail.lower()
312 if util.isquantity(tail,halt=False):
313 resl = qa.convert(tail,"arcsec")['value']
314 if os.path.exists(repodir):
315 confnum = (1.044 - 6.733 * pl.log10(resl * qa.convert(model_specrefval,"GHz")['value'] / 345.))
316 confnum = max(1,min(6,confnum))
317 conf = str(int(round(confnum)))
318 antennalist = os.path.join(repodir,"alma.cycle1." + conf + ".cfg")
319 msg("converted resolution to antennalist "+antennalist)
320 resparsed=True
321 else:
322 msg("failed to find antenna configuration repository"+
323 " at "+repodir,priority="error")
324 if not resparsed:
325 q = re.compile('.*CYCLE.?2.?;(.*)')
326 qq = q.match(antennalist.upper())
327 if qq:
328 z = qq.groups()
329 tail=z[0]
330 tail=tail.lower()
331 if util.isquantity(tail,halt=False):
332 resl = qa.convert(tail,"arcsec")['value']
333 if os.path.exists(repodir):
334 confnum = 10. ** (0.91 - 0.74 * (resl * qa.convert(model_specrefval,"GHz")['value']/345.))
335 confnum = max(1,min(7,confnum))
336 conf = str(int(round(confnum)))
337 antennalist = os.path.join(repodir,"alma.cycle2." +
338 conf + ".cfg")
339 msg("converted resolution to antennalist "+
340 antennalist)
341 resparsed=True
342 else:
343 msg("failed to find antenna configuration repository at "+repodir,priority="error")
345 if not resparsed: # assume FS
346 tail = antennalist[5:]
347 if util.isquantity(tail,halt=False):
348 resl = qa.convert(tail,"arcsec")['value']
349 if os.path.exists(repodir):
350 confnum = (2.867-pl.log10(resl*1000*qa.convert(model_specrefval,"GHz")['value']/672.))/0.0721
351 confnum = max(1,min(28,confnum))
352 conf = str(int(round(confnum)))
353 if len(conf) < 2: conf = '0' + conf
354 antennalist = os.path.join(repodir,"alma.out"
355 + conf + ".cfg")
356 msg("converted resolution to antennalist "+antennalist)
357 else:
358 msg("failed to find antenna configuration repository at "+repodir,priority="error")
361 # Search order is fileroot/ -> specified path -> repository
362 if len(antennalist) > 0:
363 if os.path.exists(fileroot+"/"+antennalist):
364 antennalist = fileroot + "/" + antennalist
365 elif not os.path.exists(antennalist) and \
366 os.path.exists(os.path.join(repodir,antennalist)):
367 antennalist = os.path.join(repodir,antennalist)
368 # Now make sure the antennalist exists
369 if not os.path.exists(antennalist):
370 emsg = "Couldn't find antennalist: %s" % antennalist
371 raise RuntimeError(emsg)
372 elif predict or components_only:
373 # antennalist is required when predicting or components only
374 util.msg("Must specify antennalist", priority="error")
376 # Read antennalist
377 if os.path.exists(antennalist):
378 stnx, stny, stnz, stnd, padnames, antnames, telescopename, posobs = util.readantenna(antennalist)
379 nant=len(stnx)
380 if nant == 1:
381 if predict and uvmode:
382 # observe="int" but antennalist is SD
383 util.msg("antennalist contains only 1 antenna",
384 priority="error")
385 uvmode = False
386 if not uvmode: #Single-dish
387 # KS TODO: what if not predicting but SD with multi-Ants
388 # in antennalist (e.g., aca.tp)? In that case, PB on plots and
389 # pointingspacing="??PB" will not be correct for heterogeneous list.
390 if sdant > nant-1 or sdant < -nant:
391 msg("antenna index %d is out of range. setting sdant=0"%sdant,priority="warn")
392 sdant = 0
393 stnx = [stnx[sdant]]
394 stny = [stny[sdant]]
395 stnz = [stnz[sdant]]
396 stnd = pl.array(stnd[sdant])
397 padnames = [padnames[sdant]]
398 antnames = [antnames[sdant]]
399 nant = 1
401 # (set back to simdata - there must be an automatic way to do this)
402 casalog.origin('simobserve')
404 aveant = stnd.mean()
405 # TODO use max ant = min PB instead?
406 pb = pbcoeff*0.29979/qa.convert(qa.quantity(model_specrefval),'GHz')['value']/aveant*3600.*180/pl.pi # arcsec
408 # PSF size
409 if uvmode:
410 # approx beam, to compare with model_cell:
411 psfsize = util.approxBeam(antennalist,qa.convert(qa.quantity(model_specrefval),'GHz')['value'])
412 else: # Single-dish
413 psfsize = pb
414 # check for model size
415 if not components_only:
416 minsize = min(qa.convert(model_size[0],'arcsec')['value'],
417 qa.convert(model_size[1],'arcsec')['value'])
418 if minsize < (2.5 * pb):
419 msg("skymodel should be larger than 2.5*primary beam. Your skymodel: %.3f arcsec < %.3f arcsec: 2.5*primary beam" % (minsize, 2.5*pb),priority="error")
420 del minsize
421 else:
422 emsg = "Can't find antennalist"
423 raise RuntimeError(emsg)
426 # now we have an estimate of the psf from the antenna configuration,
427 # so we can guess a model_cell for the case of component-only
428 # simulation,
429 if components_only:
430 # first set based on psfsize:
431 # needs high subsampling because small shifts in placement of
432 # components lead to large changes in the difference image.
433 model_cell = [ qa.quantity(str(psfsize/20)+"arcsec"),
434 qa.quantity(str(psfsize/20)+"arcsec") ]
436 # XXX if the user has set direction should we center the compskymodel there?
437 # if len(direction)>0: model_refdir = direction
439 # and can create a compskymodel image (tmp) and
440 # skymodel.flat which is what is needed for analysis.
442 if components_only:
443 newmodel = fileroot + "/" + project + ".compskymodel"
444 needmodel=True
446 modimsize=int((qa.convert(model_size[0],"arcsec")['value'])/(qa.convert(model_cell[0],"arcsec")['value']))
447# modimsize=max([modimsize,32])
448 newepoch,newlat,newlon = util.direction_splitter(model_refdir)
450 if os.path.exists(newmodel):
451 if overwrite:
452 shutil.rmtree(newmodel)
453 else:
454 needmodel=False
455 ia.open(newmodel)
456 oldshape=ia.shape()
457 if len(oldshape) != 2:
458 needmodel=True
459 else:
460 if oldshape[0] != modimsize or oldshape[1]==modimsize:
461 needmodel=True
462 oldcs=ia.coordsys()
463 ia.done()
464 olddir = (oldcs.referencevalue())['numeric']
465 if ( olddir[0] != qa.convert(newlat,oldcs.units()[0])['value'] or
466 olddir[1] != qa.convert(newlon,oldcs.units()[1])['value'] or
467 newepoch != oldcs.referencecode() ):
468 needmodel=True
469 oldcs.done()
470 del oldcs, olddir
471 if needmodel:
472 emsg = newmodel+" exists and is inconsistent with required size="+str(modimsize)+" and direction. Please set overwrite=True"
473 raise RuntimeError(emsg)
475 if needmodel:
476 csmodel = ia.newimagefromshape(newmodel,[modimsize,modimsize,1,1])
477 modelcsys = csmodel.coordsys()
478 modelshape = csmodel.shape()
479 cell0_rad=qa.convert(model_cell[0],'rad')['value']
480 cell1_rad=qa.convert(model_cell[1],'rad')['value']
481 modelcsys.setdirection(refpix=[modimsize/2,modimsize/2],
482 refval=[qa.tos(newlat),qa.tos(newlon)],
483 refcode=newepoch)
484 modelcsys.setincrement([-cell0_rad,cell1_rad],'direction')
485 modelcsys.setreferencevalue(type="spectral",value=qa.convert(model_specrefval,'Hz')['value'])
486 modelcsys.setreferencepixel(type="spectral",value=model_specrefpix)
487 modelcsys.setrestfrequency(model_specrefval)
488 modelcsys.setincrement(type="spectral",value=compwidth)
489 csmodel.setcoordsys(modelcsys.torecord())
491 modelcsys.done()
492 cl.open(complist)
493 csmodel.setbrightnessunit("Jy/pixel")
494 csmodel.modify(cl.torecord(),subtract=False)
495 cl.done()
496 csmodel.done()
497 # as noted, compskymodel doesn't need to exist, only skymodel.flat
498 # flatimage adds in components if complist!=None
499 #util.flatimage(newmodel,complist=complist,verbose=verbose)
500 util.flatimage(newmodel,verbose=verbose)
501 modelflat = fileroot + "/" + project + ".compskymodel.flat"
502# modelflat = fileroot + "/" + project + ".skymodel.flat"
503# if modelflat != newmodel+".flat":
504# if os.path.exists(modelflat):
505# shutil.rmtree(modelflat)
506# shutil.move(newmodel+".flat",modelflat)
509 # and finally, with model_cell set either from an actual skymodel,
510 # or from the antenna configuration in components_only case,
511 # we can check for the user that the psf is likely to be sampled enough:
512 cell_asec=qa.convert(model_cell[0],'arcsec')['value']
513 if psfsize < cell_asec:
514 emsg = "Sky model cell of "+str(cell_asec)+" asec is very large compared to highest resolution "+str(psfsize)+" asec - this will lead to blank or erroneous output. (Did you set incell?)"
515 shutil.rmtree(modelflat)
516 raise RuntimeError(emsg)
517 if psfsize < 2*cell_asec:
518 msg("Sky model cell of "+str(cell_asec)+" asec is large compared to highest resolution "+str(psfsize)+" asec. (Did you set incell?)",priority="warn")
520 # set this for future minimum image size
521 minimsize = 8* int(psfsize/cell_asec)
525 ##################################################################
526 # set up pointings
527 dir = model_refdir
528 dir0 = dir
529 if is_array_type(direction):
530 if len(direction) > 0:
531 if util.isdirection(direction[0],halt=False):
532 dir = direction
533 dir0 = direction[0]
534 else:
535 if util.isdirection(direction,halt=False):
536 dir = direction
537 dir0 = dir
538 util.direction = dir0
540 # if the integration time is a real time quantity
541 # test for weird units
542 if not util.isquantity(integration):
543 emsg = "integration time "+integration+" does not appear to represent a time interval (use 's','min'; not 'sec','m')"
544 raise RuntimeError(emsg)
546 if qa.quantity(integration)['unit'] != '':
547 if not qa.compare(integration,"1s"):
548 emsg = "integration time "+integration+" does not appear to represent a time interval ('s','min'; not 'sec','m')"
549 raise RuntimeError(emsg)
550 intsec = qa.convert(qa.quantity(integration),"s")['value']
551 else:
552 if len(integration)>0:
553 intsec = float(integration)
554 msg("interpreting integration time parameter as "+
555 str(intsec)+"s",priority="warn")
556 else:
557 intsec = 0
558 integration="%fs" %intsec
561 if setpointings:
562 util.msg("calculating map pointings centered at "+
563 str(dir0),origin='simobserve')
565 if len(pointingspacing) < 1:
566 if uvmode:
567 # ALMA OT uses lambda/d/sqrt(3)
568 pointingspacing = "%fPB" % (gridratio_int/pbcoeff)
569 else:
570 pointingspacing = "%fPB" % (gridratio_tp/pbcoeff)
571 if str.upper(pointingspacing)=="NYQUIST":
572 pointingspacing="%fPB" % nyquist
573 q = re.compile('(\d+.?\d+)\s*PB')
574 qq = q.match(pointingspacing.upper())
575 if qq:
576 z = qq.groups()
577 if pb <= 0:
578 emsg = "Can't calculate pointingspacing in terms of primary beam because antennalist is not specified"
579 raise RuntimeError(emsg)
580 pointingspacing = "%farcsec" % (float(z[0])*pb)
581 # todo make more robust to nonconforming z[0] strings
583 if verbose:
584 msg("pointing spacing in mosaic = "+
585 pointingspacing,origin='simobserve')
586 pointings = util.calc_pointings2(pointingspacing,
587 mapsize,maptype=maptype,
588 direction=dir, beam=pb)
589 nfld=len(pointings)
590 etime = qa.convert(qa.mul(qa.quantity(integration),scanlength),"s")['value']
591 # etime is an array of scan lengths - here they're all the same.
592 etime = etime + pl.zeros(nfld)
593 # totaltime might not allow all fields to be observed, or it might
594 # repeat
595 ptgfile = fileroot + "/" + project + ".ptg.txt"
596 else:
597 if is_array_type(ptgfile):
598 ptgfile = ptgfile[0]
599 ptgfile = ptgfile.replace('$project',project)
600 # precedence to ptg file outside the project dir
601 if os.path.exists(ptgfile):
602 shutil.copyfile(ptgfile,fileroot+"/"+project + ".ptg.txt")
603 ptgfile = fileroot + "/" + project + ".ptg.txt"
604 else:
605 if os.path.exists(fileroot+"/"+ptgfile):
606 ptgfile = fileroot + "/" + ptgfile
607 else:
608 emsg = "Can't find pointing file "+ptgfile
609 raise RuntimeError(emsg)
611 nfld, pointings, etime = util.read_pointings(ptgfile)
612 if max(etime) <= 0:
613 # integration is a string in input params
614 etime = intsec
615 # make etime into an array
616 etime = etime + pl.zeros(nfld)
617 # etimes determine stop/start i.e. the scan
618 # if a longer etime is in the file, it'll do multiple integrations
619 # per scan
620 # expects that the cal is separate, and this is just one round of the mosaic
621 # furthermore, the cal will use _integration_ from the inputs, and that
622 # needs to be less than the min etime:
623 if min(etime) < intsec:
624 integration = str(min(etime))+"s"
625 msg("Setting integration to "+integration+
626 " to match the shortest time in the pointing file.",
627 priority="warn")
628 intsec = min(etime)
631 # find imcenter - phase center
632 imcenter , offsets = util.median_direction(pointings)
633 epoch, ra, dec = util.direction_splitter(imcenter)
635 # model is centered at model_refdir, and has model_size;
636 mepoch, mra, mdec = util.direction_splitter(model_refdir)
637 # ra/mra should be in degrees, but just in case
638 ra=qa.convert(ra,'deg')
639 mra=qa.convert(mra,'deg')
640 dec=qa.convert(dec,'deg')
641 mdec=qa.convert(mdec,'deg')
643 # observation near ra=0:
644 if abs(mra['value'])<10 or abs(mra['value']-360.)<10 or abs(ra['value'])<10 or abs(mra['value']-360.)<10:
645 nearzero=True
646 else:
647 nearzero=False
649 # fix wraps
650 if nearzero: # put break at 180
651 if ra['value'] >= 180.:
652 ra['value'] = ra['value'] - 360.
653 if mra['value'] >= 180.:
654 mra['value'] = mra['value'] - 360.
655 else:
656 if ra['value'] >= 360.:
657 ra['value'] = ra['value'] - 360.
658 if mra['value'] >= 360.:
659 mra['value'] = mra['value'] - 360.
661 # shift is the offset from the model to imcenter
662 shift = [ (qa.convert(ra,'deg')['value'] -
663 qa.convert(mra,'deg')['value'])*pl.cos(qa.convert(dec,'rad')['value'] ),
664 (qa.convert(dec,'deg')['value'] - qa.convert(mdec,'deg')['value']) ]
665 if verbose:
666 msg("pointings are shifted relative to the model by %g,%g arcsec" % (shift[0]*3600,shift[1]*3600),origin='simobserve')
668 xmax = qa.convert(model_size[0],'deg')['value']*0.5
669 ymax = qa.convert(model_size[1],'deg')['value']*0.5
670 # add PB halfwidth (relmargin=0.5)
671 # for mosaics of small model images
672 xmax=xmax+pb*relmargin/3600
673 ymax=ymax+pb*relmargin/3600
674 overlap = False
675 # wrapang in median_direction should make offsets always small, not >360
676 for i in range(offsets.shape[1]):
677 xc = pl.absolute(offsets[0,i]+shift[0]) # offsets and shift are in degrees
678 yc = pl.absolute(offsets[1,i]+shift[1])
679 if xc < xmax and yc < ymax:
680 overlap = True
681 break
683 if setpointings:
684 if os.path.exists(ptgfile):
685 if overwrite:
686 os.remove(ptgfile)
687 else:
688 emsg = "pointing file "+ptgfile+" already exists and user does not want to overwrite"
689 raise RuntimeError(emsg)
690 util.write_pointings(ptgfile,pointings,etime.tolist())
692 msg("center = "+imcenter,origin='simobserve')
693 if nfld > 1 and verbose:
694 for idir in range(min(len(pointings),20)):
695 msg(" "+pointings[idir],origin='simobserve')
696 if nfld >= 20:
697 msg(" (printing only first 20 - see pointing file for full list)")
699 if not overlap:
700 emsg = "No overlap between model and pointings"
701 raise RuntimeError(emsg)
705 ##################################################################
706 # calibrator is not explicitly contained in the pointing file
707 # but interleaved with etime=intergration
708 util.isquantity(calflux)
709 calfluxjy = qa.convert(calflux,'Jy')['value']
710 # XML returns a list even for a string:
711 if is_array_type(caldirection): caldirection = caldirection[0]
712 if len(caldirection) < 4: caldirection = ""
713 if calfluxjy > 0 and caldirection != "" and uvmode:
714 docalibrator = True
715 if os.path.exists(fileroot+"/"+project+'.cal.cclist'):
716 shutil.rmtree(fileroot+"/"+project+'.cal.cclist')
717 util.isdirection(caldirection)
718 cl.done()
719 cl.addcomponent(flux=calfluxjy,dir=caldirection,
720 label="phase calibrator")
721 # set reference freq to center freq of model
722 cl.rename(fileroot+"/"+project+'.cal.cclist')
723 cl.done()
724 else:
725 docalibrator = False
731 ##################################################################
732 # create one figure for model and pointings - need antenna diam
733 # to determine primary beam
734 if grfile:
735 file = fileroot + "/" + project + ".skymodel.png"
736 else:
737 file = ""
739 if grscreen or grfile:
740 util.newfig(show=grscreen)
742# remove after we fix the scaling algorithm for the images
743 if components_only:
744 pl.plot()
745 # TODO add symbols at locations of components
746 pl.plot(coffs[0,]*3600,coffs[1,]*3600,'o',c="#dddd66")
747 pl.axis("equal")
748 else:
749 discard = util.statim(modelflat,plot=True,incell=model_cell)
751 lims = pl.xlim(),pl.ylim()
752 if pb <= 0 and verbose:
753 msg("unknown primary beam size for plot",priority="warn")
754 if max(max(lims)) > pb and not components_only:
755 plotcolor = 'w'
756 else:
757 plotcolor = 'k'
759 #if offsets.shape[1] > 16 or pb <= 0 or pb > pl.absolute(max(max(lims))):
760 if offsets.shape[1] > 19 or pb <= 0:
761 lims = pl.xlim(),pl.ylim()
762 pl.plot((offsets[0]+shift[0])*3600.,
763 (offsets[1]+shift[1])*3600.,
764 plotcolor+'+',markeredgewidth=1)
765 #if pb > 0 and pl.absolute(lims[0][0]) > pb:
766 if pb > 0:
767 plotpb(pb,pl.gca(),lims=lims,color=plotcolor)
768 else:
769 from matplotlib.patches import Circle
770 for i in range(offsets.shape[1]):
771 pl.gca().add_artist(Circle(
772 ((offsets[0,i]+shift[0])*3600,
773 (offsets[1,i]+shift[1])*3600),
774 radius=pb/2.,edgecolor=plotcolor,fill=False,
775 label='beam',transform=pl.gca().transData,clip_on=True))
777 xlim = max(abs(pl.array(lims[0])))
778 ylim = max(abs(pl.array(lims[1])))
779 # show entire pb: (statim doesn't by default)
780 pl.xlim([max([xlim,pb/2]),min([-xlim,-pb/2])])
781 pl.ylim([min([-ylim,-pb/2]),max([ylim,pb/2])])
782 pl.xlabel("resized model sky",fontsize="x-small")
783 util.endfig(show=grscreen,filename=file)
792 ##################################################################
793 # set up observatory, feeds, etc
794 quickpsf_current = False
796 if uvmode:
797 msfile = fileroot + "/" + project + '.ms'
798 else:
799 msfile = fileroot + "/" + project + '.sd.ms'
801 if predict:
802 # TODO check for frequency overlap here - if zero stop
803 # position overlap already checked above in pointing section
805 message = "preparing empty measurement set"
806 if verbose:
807 msg(" ",priority="info")
808 msg(message,origin="simobserve",priority="warn")
809 else:
810 msg(message,origin="simobserve")
812 nbands = 1;
813 fband = util.bandname(qa.convert(model_specrefval, 'GHz')['value'])
815 ############################################
816 # predict observation
818 # if someone has the old style refdate with the included, discard
819 q = re.compile('(\d*/\d+/\d+)([/:\d]*)')
820 qq = q.match(refdate)
821 if not qq:
822 emsg = "Invalid reference date "+refdate
823 raise RuntimeError(emsg)
824 else:
825 z = qq.groups()
826 refdate=z[0]
827 if len(z)>1:
828 if len(z[1])>1:
829 msg("Discarding time part of refdate, '"+z[1]+
830 "', in favor of hourangle parameter = "+hourangle,
831 origin='simobserve')
833 if hourangle=="transit":
834 haoffset=0.0
835 else:
836 haoffset="no"
837 # is this a time quantity?
838 if qa.isquantity(str(hourangle)+"h"):
839 if qa.compare(str(hourangle)+"h","s"):
840 haoffset=qa.convert(qa.quantity(str(hourangle)+
841 "h"),'s')['value']
842 if qa.isquantity(hourangle):
843 qha=qa.convert(hourangle,"s")
844 if qa.compare(qha,"s"):
845 haoffset=qa.convert(qha,'s')['value']
846 if haoffset=="no":
847 msg("Cannot interpret your hourangle parameter "+hourangle+
848 " as a time quantity e.g. '5h', 30min'",
849 origin="simobserve",priority="error")
850 else:
851 msg("You desire an hour angle of "+
852 str(haoffset/3600.)+" hours",origin="simobserve")
854 refdate=refdate+"/00:00:00"
855 usehourangle=True
858 intsec = qa.convert(qa.quantity(integration),"s")['value']
860 # totaltime as an integer for # times through the mosaic:
861 if not util.isquantity(totaltime,halt=False):
862 emsg = "totaltime "+totaltime+" does not appear to represent a time interval (use 's','min','h'; not 'sec','m','hr')"
863 raise RuntimeError(emsg)
865 if qa.quantity(totaltime)['value'] < 0.:
866 # casapy crashes for negative totaltime
867 emsg = "Negative totaltime is not allowed."
868 raise ValueError(emsg)
869 if qa.quantity(totaltime)['unit'] == '':
870 # assume it means number of maps, or # repetitions.
871 totalsec = sum(etime)
872 if docalibrator:
873 totalsec = totalsec + intsec # cal gets one int-time
874 totalsec = float(totaltime) * totalsec
875 msg("Total observing time = "+str(totalsec)+"s.")
876 else:
877 if not qa.compare(totaltime,"1s"):
878 emsg = "totaltime "+totaltime+" does not appear to represent a time interval (use 's','min','h'; not 'sec','m','hr')"
879 raise ValueError(emsg)
880 totalsec = qa.convert(qa.quantity(totaltime),'s')['value']
882 if os.path.exists(msfile) and not overwrite: #redundant check?
883 emsg = ("measurement set "+msfile+
884 " already exists and user does not wish to overwrite")
885 raise RuntimeError(emsg)
886 sm.open(msfile)
888 diam = stnd;
889 # WARNING: sm.setspwindow is not consistent with clean::center
890 # but the "start" is the center of the first channel:
891 model_start = qa.sub(model_specrefval,
892 qa.mul(model_width,model_specrefpix))
894 mounttype = 'alt-az'
895 if telescopename in ['DRAO', 'WSRT']:
896 mounttype = 'EQUATORIAL'
897 # Should ASKAP be BIZARRE or something else? It may be effectively equatorial.
899 sm.setconfig(telescopename=telescopename, x=stnx, y=stny, z=stnz,
900 dishdiameter=diam.tolist(),
901 mount=[mounttype], antname=antnames, padname=padnames,
902 coordsystem='global', referencelocation=posobs)
903 if str.upper(telescopename).find('VLA') >= 0:
904 sm.setspwindow(spwname=fband, freq=qa.tos(model_start),
905 deltafreq=qa.tos(model_width),
906 freqresolution=qa.tos(model_width),
907 nchannels=model_nchan, refcode=outframe,
908 stokes='RR LL')
909 sm.setfeed(mode='perfect R L',pol=[''])
910 else:
911 sm.setspwindow(spwname=fband, freq=qa.tos(model_start),
912 deltafreq=qa.tos(model_width),
913 freqresolution=qa.tos(model_width),
914 nchannels=model_nchan, refcode=outframe,
915 stokes='XX YY')
916 sm.setfeed(mode='perfect X Y',pol=[''])
918 if verbose:
919 msg(" spectral window set at %s" % qa.tos(model_specrefval),
920 origin='simobserve')
921 sm.setlimits(shadowlimit=0.01, elevationlimit='10deg')
922 if uvmode:
923 sm.setauto(0.0)
924 else: #Single-dish
925 # auto-correlation should be unity for single dish obs.
926 sm.setauto(1.0)
928 mereftime = me.epoch('UTC', refdate)
929 # integration is a scalar quantity, etime is a vector of seconds
930 sm.settimes(integrationtime=integration, usehourangle=usehourangle,
931 referencetime=mereftime)
933 for k in range(0,nfld):
934 src = project + '_%d' % k
935 sm.setfield(sourcename=src, sourcedirection=pointings[k],
936 calcode="OBJ", distance='0m')
937 if k == 0:
938 sourcefieldlist = src
939 else:
940 sourcefieldlist = sourcefieldlist + ',' + src
941 if docalibrator:
942 sm.setfield(sourcename="phase calibrator",
943 sourcedirection=caldirection,calcode='C',
944 distance='0m')
946 # time required to observe all planned scanes in etime array:
947 totalscansec = sum(etime)
948 kfld = 0
950 if totalsec < totalscansec:
951 msg("Not all pointings in the mosaic will be observed - check mosaic setup and exposure time parameters!",priority="warn")
952 ###
953 #casalog.post("you need at least %16.12e sec but you have %16.12e sec (%f < %f = %s)" % (totalscansec, totalsec, totalsec, totalscansec, str(totalsec<totalscansec)))
954 ###
956 # sm.observemany
957 srces = []
958 starttimes = []
959 stoptimes = []
960 dirs = []
962 if usehourangle:
963 sttime = -totalsec/2.0
964 else:
965 sttime = 0. # leave start at the reftime
966 sttime=sttime+haoffset
967 scanstart=sttime
969 # can before sources
970 if docalibrator:
971 endtime = sttime + qa.convert(integration,'s')['value']
972 sm.observe(sourcename="phase calibrator", spwname=fband,
973 starttime=qa.quantity(sttime, "s"),
974 stoptime=qa.quantity(endtime, "s"),
975 state_obs_mode="CALIBRATE_PHASE.ON_SOURCE",
976 state_sig=True,project=project);
977 sttime = endtime
979 while (sttime-scanstart) < totalsec: # the last scan could exceed totaltime
980 endtime = sttime + etime[kfld]
981 src = project + '_%d' % kfld
982 srces.append(src)
983 starttimes.append(str(sttime)+"s")
984 stoptimes.append(str(endtime)+"s")
985 dirs.append(pointings[kfld])
986 kfld = kfld + 1
987 # advance start time - XX someday slew goes here
988 sttime = endtime
990 if kfld == nfld:
991 if docalibrator:
992 endtime = sttime + qa.convert(integration,'s')['value']
994 # need to observe cal singly to get new row in obs table
995 # so first observemany the on-source pointing(s)
996 sm.observemany(sourcenames=srces,
997 spwname=fband,
998 starttimes=starttimes,
999 stoptimes=stoptimes,
1000 project=project)
1001 # and clear the list
1002 srces = []
1003 starttimes = []
1004 stoptimes = []
1005 dirs = []
1006 sm.observe(sourcename="phase calibrator", spwname=fband,
1007 starttime=qa.quantity(sttime, "s"),
1008 stoptime=qa.quantity(endtime, "s"),
1009 state_obs_mode="CALIBRATE_PHASE.ON_SOURCE",
1010 state_sig=True,project=project);
1011 kfld = kfld + 1
1012 sttime = endtime
1013# if kfld > nfld: kfld = 0
1014 if kfld > nfld-1: kfld = 0
1015 # if directions is unset, NewMSSimulator::observemany
1017 # looks up the direction in the field table.
1018 if not docalibrator:
1019 sm.observemany(sourcenames=srces,
1020 spwname=fband,
1021 starttimes=starttimes,
1022 stoptimes=stoptimes,
1023 project=project)
1025 sm.setdata(fieldid=list(range(0,nfld)))
1026 if uvmode or components_only: #Interferometer only
1027 sm.setvp(dovp=True,usedefaultvp=False)
1028 # only use mosaic gridding for Het arrays for now -
1029 # the standard gridding with a VPSkyJones is less susceptible
1030 # to issues if the image is too small which can happen a lot
1031 # in Simulation.
1032 if len(pl.unique(diam))>1:
1033 if nfld>1:
1034 sm.setoptions(ftmachine="mosaic")
1035 else:
1036 msg("Heterogeneous array only supported for mosaics (nfld>1), and make sure that your image is larger than the primary beam or results may be unstable",priority="error")
1037 else:
1038 # checks have to be manual since there's no way to
1039 # get the "diam" out of PBMath AFAIK
1040 if telescopename=="ALMA":
1041 if (diam[0]<10)|(diam[0]>13):
1042 msg("Diameter = %f is inconsistent with telescope=ALMA in the configuration file. *12m ALMA PB will be used*. Use observatory='ACA' in the configuration file to use a 7m diameter."%diam[0],priority="warn")
1043 n=len(diam)
1044 diam=12.+pl.zeros(n)
1045 elif telescopename=="ACA":
1046 if (diam[0]<6)|(diam[0]>7.5):
1047 msg("Diameter = %f is inconsistent with telescope=ALMA in the configuration file. *7m ACA PB will be used*"%diam[0],priority="warn")
1048 n=len(diam)
1049 diam=7.+pl.zeros(n)
1050 else:
1051 msg("Note: diameters in configuration file will not be used - PB for "+telescopename+" will be used",priority="info")
1054 msg("done setting up observations (blank visibilities)",
1055 origin='simobserve')
1056 if verbose: sm.summary()
1058 # do actual calculation of visibilities:
1060 if not uvmode: #Single-dish
1061 sm.setoptions(gridfunction='pb',ftmachine="sd",location=posobs)
1063 if not components_only:
1064 if docalibrator:
1065 if len(complist) <=0:
1066 complist=fileroot+"/"+project+'.cal.cclist'
1067 else:
1068 # XXX will 2 cl work?
1069 complist=complist+","+fileroot+"/"+project+'.cal.cclist'
1071 if len(complist) > 1:
1072 message = "predicting from "+newmodel+" and "+complist
1073 if verbose:
1074 msg(" ",priority="info")
1075 msg(message,priority="warn",origin="simobserve")
1076 else:
1077 msg(message,origin="simobserve")
1078 else:
1079 message = "predicting from "+newmodel
1080 if verbose:
1081 msg(" ",priority="info")
1082 msg(message,priority="warn",origin="simobserve")
1083 else:
1084 msg(message,origin="simobserve")
1085 sm.predict(imagename=newmodel,complist=complist)
1086 else: # if we're doing only components
1087 # XXX will 2 cl work?
1088 if docalibrator:
1089 complist=complist+","+fileroot+"/"+project+'.cal.cclist'
1090 if verbose:
1091 msg("predicting from "+complist,priority="warn",
1092 origin="simobserve")
1093 else:
1094 msg("predicting from "+complist,origin="simobserve")
1095 sm.predict(complist=complist)
1097 sm.done()
1098 msg('generation of measurement set '+msfile+' complete',
1099 origin="simobserve")
1101 # rest freqs are hardcoded to the first freq in the spw in core
1102 tb.open(msfile+"/SPECTRAL_WINDOW/",nomodify=False)
1103 restfreq=tb.getcol("REF_FREQUENCY")
1104 for i in range(len(restfreq)):
1105 restfreq[i]=qa.convert(qa.quantity(model_specrefval),'Hz')['value']
1106 tb.putcol("REF_FREQUENCY",restfreq)
1107 tb.done()
1108 tb.open(msfile+"/SOURCE/",nomodify=False)
1109 restfreq=tb.getcol("REST_FREQUENCY")
1110 for i in range(len(restfreq)):
1111 restfreq[i]=qa.convert(qa.quantity(model_specrefval),'Hz')['value']
1112 tb.putcol("REST_FREQUENCY",restfreq)
1113 tb.done()
1119 ############################################
1120 # create figure
1121 if grfile:
1122 file = fileroot + "/" + project + ".observe.png"
1123 else:
1124 file = ""
1126 # update psfsize using uv coverage instead of maxbase above
1127 if os.path.exists(msfile):
1128 # psfsize was set from the antenna posns before, but uv is better
1129 tb.open(msfile)
1130 rawdata = tb.getcol("UVW")
1131 tb.done()
1132 # TODO make this use the 90% baseline as in aU.getBaselineStats
1133 maxbase = max([max(rawdata[0,]),max(rawdata[1,])]) # in m
1134 if maxbase > 0.:
1135 psfsize = 0.3/qa.convert(qa.quantity(model_specrefval),'GHz')['value']/maxbase*3600.*180/pl.pi # lambda/b converted to arcsec
1136 if not uvmode:
1137 msg("An single observation is requested but CROSS-correlation in "+msfile,priority="error")
1138 else: #Single-dish (zero-spacing)
1139 psfsize = pb
1140 if uvmode:
1141 msg("An interferometer observation is requested but only AUTO-correlation in "+msfile,priority="error")
1142 minimsize = 8* int(psfsize/cell_asec)
1143 else:
1144 msg("Couldn't find "+msfile,priority="error")
1146 if uvmode:
1147 multi = [2,2,1]
1148 else:
1149 multi = 0
1151 if (grscreen or grfile):
1152 util.newfig(multi=multi,show=grscreen)
1153 util.ephemeris(refdate,
1154 direction=util.direction,
1155 telescope=telescopename,
1156 ms=msfile,
1157 usehourangle=usehourangle,
1158 cofa=posobs)
1159 casalog.origin('simobserve')
1160 if uvmode:
1161 util.nextfig()
1162 util.plotants(stnx, stny, stnz, stnd, padnames)
1164 # uv coverage
1165 util.nextfig()
1166 klam_m = 300/qa.convert(model_specrefval,'GHz')['value']
1167 pl.box()
1168 pl.plot(rawdata[0,]/klam_m,rawdata[1,]/klam_m,'b,')
1169 pl.plot(-rawdata[0,]/klam_m,-rawdata[1,]/klam_m,'b,')
1170 ax = pl.gca()
1171 ax.yaxis.LABELPAD = -4
1172 pl.xlabel('u[klambda]',fontsize='x-small')
1173 pl.ylabel('v[klambda]',fontsize='x-small')
1174 pl.axis('equal')
1176 # show dirty beam from observed uv coverage
1177 util.nextfig()
1178 imt = imager()
1179 imt.open(msfile)
1180 # TODO spectral parms
1181 msg("using default model cell "+str(model_cell[0])+
1182 " for PSF calculation",origin='simobserve')
1183 imt.defineimage(cellx=str(model_cell[0]["value"])+
1184 str(model_cell[0]["unit"]),
1185 nx=int(max([minimsize,128])))
1186 # TODO trigger im.setoptions(ftmachine="mosaic")
1187 if os.path.exists(fileroot+"/"+project+".quick.psf"):
1188 shutil.rmtree(fileroot+"/"+project+".quick.psf")
1190 # if obs is unknown, casalog will send a warning to screen
1191 # temporarily(?) suppress that
1192 if not telescopename in me.obslist():
1193 casalog.filter("ERROR")
1194 imt.approximatepsf(psf=fileroot+"/"+project+".quick.psf")
1195 if not telescopename in me.obslist():
1196 casalog.filter() # set back to default level.
1198 quickpsf_current = True
1199 beam = imt.fitpsf(psf=fileroot+"/"+project+".quick.psf")
1200 imt.done()
1201 ia.open(fileroot+"/"+project+".quick.psf")
1202 beamcs = ia.coordsys()
1203 beam_array = ia.getchunk(axes=[beamcs.findcoordinate("spectral")['pixel'][0],beamcs.findcoordinate("stokes")['pixel'][0]],dropdeg=True)
1204 nn = beam_array.shape
1205 xextent = nn[0]*cell_asec*0.5
1206 xextent = [xextent,-xextent]
1207 yextent = nn[1]*cell_asec*0.5
1208 yextent = [-yextent,yextent]
1209 flipped_array = beam_array.transpose()
1210 ttrans_array = flipped_array.tolist()
1211 ttrans_array.reverse()
1212 pl.imshow(ttrans_array,
1213 interpolation='bilinear',
1214 cmap=pl.cm.jet,
1215 extent=xextent+yextent,
1216 origin="lower")
1217 pl.title(project+".quick.psf",fontsize="x-small")
1218 b = qa.convert(beam[1],'arcsec')['value']
1219 pl.xlim([-3*b,3*b])
1220 pl.ylim([-3*b,3*b])
1221 ax = pl.gca()
1222 pl.text(0.05,0.95,"bmaj=%7.1e\nbmin=%7.1e" % (beam[1]['value'],beam[2]['value']),transform = ax.transAxes,bbox=dict(facecolor='white', alpha=0.7),size="x-small",verticalalignment="top")
1223 ia.done()
1224 util.endfig(show=grscreen,filename=file)
1226 elif thermalnoise != "" or leakage > 0:
1227 # Not predicting this time but corrupting. Get obsmode from ms.
1228 if not os.path.exists(msfile):
1229 msg("Couldn't find "+msfile,priority="error")
1230 uvmode = (not util.ismstp(msfile,halt=False))
1233 ######################################################################
1234 # noisify
1236 noise_any = False
1237 msroot = fileroot + "/" + project # if leakage, can just copy from this project
1239 if thermalnoise != "":
1240 knowntelescopes = ["ALMASD", "ALMA", "ACA", "SMA", "EVLA", "VLA"]
1242 noise_any = True
1244 noisymsroot = msroot + ".noisy"
1245 if not uvmode: #Single-dish
1246 msroot += ".sd"
1247 noisymsroot += ".sd"
1249 # Cosmic background radiation temperature in K.
1250 t_cmb = 2.725
1253 # check for interferometric ms:
1254 if not os.path.exists(msroot+".ms"):
1255 msg("Couldn't find "+msroot+".ms",priority="error")
1257 # MS exists
1258 message = 'copying '+msroot+'.ms to ' + \
1259 noisymsroot+'.ms and adding thermal noise'
1260 if verbose:
1261 msg(" ",priority="info")
1262 msg(message,origin="noise",priority="warn")
1263 else:
1264 msg(message,origin="noise")
1266 if os.path.exists(noisymsroot+".ms"):
1267 shutil.rmtree(noisymsroot+".ms")
1268 shutil.copytree(msfile,noisymsroot+".ms")
1269 if sm.name() != '':
1270 emsg = "table persistence error on %s" % sm.name()
1271 raise RuntimeError(emsg)
1273 # if not predicted this time, get telescopename from ms
1274 if not predict:
1275 tb.open(noisymsroot+".ms/OBSERVATION")
1276 n = tb.getcol("TELESCOPE_NAME")
1277 telescopename = n[0]
1278 # todo add check that entire column is the same
1279 tb.done()
1280 msg("telescopename read from "+noisymsroot+".ms: "+
1281 telescopename)
1283 if telescopename not in knowntelescopes:
1284 msg("thermal noise only works properly for ALMA/ACA, (E)VLA, and SMA",origin="simobserve",priority="warn")
1285 eta_p, eta_s, eta_b, eta_t, eta_q, t_rx = util.noisetemp(telescope=telescopename,freq=model_specrefval)
1287 # antenna efficiency
1288 eta_a = eta_p * eta_s * eta_b * eta_t
1289 if verbose:
1290 msg('antenna efficiency = '+str(eta_a), origin="simobserve")
1291 msg('spillover efficiency = '+str(eta_s), origin="simobserve")
1292 msg('correlator efficiency = '+str(eta_q), origin="simobserve")
1293 # sensitivity constant
1294 scoeff = -1 #Force setting the default value, 1./sqrt(2.0)
1295 if not uvmode: #Single-dish
1296 scoeff = 1.0
1297 if verbose: msg('sensitivity constant = '+str(scoeff),
1298 origin="simobserve")
1300 sm.openfromms(noisymsroot+".ms") # an existing MS
1301 sm.setdata(fieldid=[]) # force to get all fields
1302 sm.setseed(seed)
1303 if thermalnoise == "tsys-manual":
1304 if verbose:
1305 message = "sm.setnoise(spillefficiency="+str(eta_s)+\
1306 ",correfficiency="+str(eta_q)+\
1307 ",antefficiency="+str(eta_a)+\
1308 ",trx="+str(t_rx)+",tau="+str(tau0)+\
1309 ",tatmos="+str(t_sky)+",tground="+str(t_ground)+\
1310 ",tcmb="+str(t_cmb)
1311 if not uvmode: message += ",senscoeff="+str(scoeff)
1312 message += ",mode='tsys-manual')"
1313 msg(message);
1314 msg("** this may be slow if your MS is finely sampled in time ** ",priority="warn")
1315 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q,
1316 antefficiency=eta_a,trx=t_rx,
1317 tau=tau0,tatmos=t_sky,tground=t_ground,tcmb=t_cmb,
1318 mode="tsys-manual",senscoeff=scoeff)
1319 else:
1320 if verbose:
1321 message = "sm.setnoise(spillefficiency="+str(eta_s)+\
1322 ",correfficiency="+str(eta_q)+\
1323 ",antefficiency="+str(eta_a)+\
1324 ",trx="+str(t_rx)+",tground="+str(t_ground)+\
1325 ",tcmb="+str(t_cmb)
1326 if not uvmode: message += ",senscoeff="+str(scoeff)
1327 message += ",mode='tsys-atm'"+\
1328 ",pground='560mbar',altitude='5000m'"+\
1329 ",waterheight='2km',relhum=20,pwv="+str(user_pwv)+"mm)"
1330 msg(message);
1331 msg("** this may be slow if your MS is finely sampled in time ** ",priority="warn")
1332 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q,
1333 antefficiency=eta_a,trx=t_rx,
1334 tground=t_ground,tcmb=t_cmb,pwv=str(user_pwv)+"mm",
1335 mode="tsys-atm",table=noisymsroot,senscoeff=scoeff)
1336 # don't set table, that way it won't save to disk
1337 # mode="calculate",table=noisymsroot)
1339 sm.corrupt();
1340 sm.done();
1343 msroot = noisymsroot
1344 if verbose: msg("done corrupting with thermal noise",origin="noise")
1347 if leakage > 0:
1348 noise_any = True
1349 # TODO: need to handle SD name when leakage is available
1350 if msroot == fileroot+"/"+project:
1351 noisymsroot = fileroot + "/" + project + ".noisy"
1352 else:
1353 noisymsroot = fileroot + "/" + project + ".noisier"
1354 if not uvmode: #Single-dish
1355 msg("Can't corrupt SD data with polarization leakage",
1356 priority="warn")
1357 if os.path.exists(msfile):
1358 msg('copying '+msfile+' to ' +
1359 noisymsroot+'.ms and adding polarization leakage',
1360 origin="noise",priority="warn")
1361 if os.path.exists(noisymsroot+".ms"):
1362 shutil.rmtree(noisymsroot+".ms")
1363 shutil.copytree(msfile,noisymsroot+".ms")
1364 if sm.name() != '':
1365 emsg = "table persistence error on %s" % sm.name()
1366 raise RuntimeError(emsg)
1368 sm.openfromms(noisymsroot+".ms") # an existing MS
1369 sm.setdata(fieldid=[]) # force to get all fields
1370 sm.setleakage(amplitude=leakage,table=noisymsroot+".cal")
1371 sm.corrupt();
1372 sm.done();
1376 # cleanup - delete newmodel, newmodel.flat etc
1377# if os.path.exists(modelflat):
1378# shutil.rmtree(modelflat)
1379 if os.path.exists(modelflat+".regrid"):
1380 shutil.rmtree(modelflat+".regrid")
1381 if os.path.exists(fileroot+"/"+project+".noisy.T.cal"):
1382 shutil.rmtree(fileroot+"/"+project+".noisy.T.cal")
1383 if os.path.exists(fileroot+"/"+project+".noisy.sd.T.cal"):
1384 shutil.rmtree(fileroot+"/"+project+".noisy.sd.T.cal")
1386 finally:
1387 finalize_tools()
1389##### Helper functions to plot primary beam
1390def plotpb(pb,axes,lims=None,color='k'):
1391 # This beam is automatically scaled when you zoom in/out but
1392 # not anchored in plot area. We'll wait for Matplotlib 0.99
1393 # for that function.
1394 #major=major
1395 #minor=minor
1396 #rangle=rangle
1397 #bwidth=max(major*pl.cos(rangle),minor*pl.sin(rangle))*1.1
1398 #bheight=max(major*pl.sin(rangle),minor*pl.cos(rangle))*1.1
1399 from matplotlib.patches import Rectangle, Circle #,Ellipse
1400 try:
1401 from matplotlib.offsetbox import AnchoredOffsetbox, AuxTransformBox
1402 box = AuxTransformBox(axes.transData)
1403 box.set_alpha(0.7)
1404 circ = Circle((pb,pb),radius=pb/2.,color=color,fill=False,\
1405 label='primary beam',linewidth=2.0)
1406 box.add_artist(circ)
1407 pblegend = AnchoredOffsetbox(loc=3,pad=0.2,borderpad=0.,\
1408 child=box,prop=None,frameon=False)#,frameon=True)
1409 pblegend.set_alpha(0.7)
1410 axes.add_artist(pblegend)
1411 except:
1412 casalog.post("Using old matplotlib substituting with circle")
1413 # work around for old matplotlib
1414 boxsize = pb*1.1
1415 if not lims: lims = axes.get_xlim(),axes.get_ylim()
1416 incx = 1
1417 incy = 1
1418 if axes.xaxis_inverted(): incx = -1
1419 if axes.yaxis_inverted(): incy = -1
1420 #ecx = lims[0][0] + bwidth/2.*incx
1421 #ecy = lims[1][0] + bheight/2.*incy
1422 ccx = lims[0][0] + boxsize/2.*incx
1423 ccy = lims[1][0] + boxsize/2.*incy
1425 #box = Rectangle((lims[0][0],lims[1][0]),incx*bwidth,incy*bheight,
1426 box = Rectangle((lims[0][0],lims[1][0]),incx*boxsize,incy*boxsize,
1427 alpha=0.7,facecolor='w',
1428 transform=axes.transData) #Axes
1429 #beam = Ellipse((ecx,ecy),major,minor,angle=rangle,
1430 beam = Circle((ccx,ccy), radius=pb/2.,
1431 edgecolor='k',fill=False,
1432 label='beam',transform=axes.transData)
1433 #props = {'pad': 3, 'edgecolor': 'k', 'linewidth':2, 'facecolor': 'w', 'alpha': 0.5}
1434 #pl.matplotlib.patches.bbox_artist(beam,axes.figure.canvas.get_renderer(),props=props)
1435 axes.add_artist(box)
1436 axes.add_artist(beam)
1438def finalize_tools():
1439 if ia.isopen(): ia.close()
1440 sm.close()
1441 #cl.close() # crashes casa