|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 # geodesy and pointing and other helper functions that are useful 00002 # to be available outside of the simdata task 00003 import casac 00004 import os 00005 import shutil 00006 import commands 00007 #import pdb 00008 # all I really need is casalog, but how to get it:? 00009 from taskinit import * 00010 import pylab as pl 00011 00012 # qa doesn't hold state. 00013 #qatool = casac.homefinder.find_home_by_name('quantaHome') 00014 #qa = qatool.create() 00015 00016 im, cb, ms, tb, fg, me, ia, po, sm, cl, cs, rg, sl, dc, vp, msmd, fi, fn = gentools() 00017 00018 class compositenumber: 00019 def __init__(self, maxval=100): 00020 self.generate(maxval) 00021 def generate(self,maxval): 00022 n2 = int(log(float(maxval))/log(2.0) + 1) +1 00023 n3 = int(log(float(maxval))/log(3.0) + 1) +1 00024 n5 = int(log(float(maxval))/log(5.0) + 1) +1 00025 itsnumbers=pl.zeros(n2*n3*n5) 00026 n = 0 00027 for i2 in range(n2): 00028 for i3 in range(n3): 00029 for i5 in range(n5): 00030 composite=( 2.**i2 * 3.**i3 * 5.**i5 ) 00031 itsnumbers[n] = composite 00032 #print i2,i3,i5,composite 00033 n=n+1 00034 itsnumbers.sort() 00035 maxi=0 00036 while maxi<(n2*n3*n5) and itsnumbers[maxi]<=maxval: maxi=maxi+1 00037 self.itsnumbers=pl.int64(itsnumbers[0:maxi]) 00038 def list(self): 00039 print self.itsnumbers 00040 def nextlarger(self,x): 00041 if x>max(self.itsnumbers): self.generate(2*x) 00042 xi=0 00043 n=self.itsnumbers.__len__() 00044 while xi<n and self.itsnumbers[xi]<x: xi=xi+1 00045 return self.itsnumbers[xi] 00046 00047 00048 00049 00050 00051 class simutil: 00052 def __init__(self, direction="", 00053 centerfreq=qa.quantity("245GHz"), 00054 bandwidth=qa.quantity("1GHz"), 00055 totaltime=qa.quantity("0h"), 00056 verbose=False): 00057 self.direction=direction 00058 self.verbose=verbose 00059 self.centerfreq=centerfreq 00060 self.bandwidth=bandwidth 00061 self.totaltime=totaltime 00062 self.pmulti=0 # rows, cols, currsubplot 00063 00064 00065 def newfig(self,multi=0,show=True): # new graphics window/file 00066 if show: 00067 pl.ion() # creates a fig if ness 00068 else: 00069 pl.ioff() 00070 pl.clf() 00071 00072 if multi!=0: 00073 if type(multi)!=type([]): 00074 self.msg("internal error setting multi-panel figure with multi="+str(multi),priority="warn") 00075 if len(multi)!=3: 00076 self.msg("internal error setting multi-panel figure with multi="+str(multi),priority="warn") 00077 self.pmulti=multi 00078 pl.subplot(multi[0],multi[1],multi[2]) 00079 pl.subplots_adjust(left=0.05,right=0.98,bottom=0.09,top=0.95,hspace=0.2,wspace=0.2) 00080 00081 00082 ########################################################### 00083 00084 def nextfig(self): # advance subwindow 00085 ax=pl.gca() 00086 l=ax.get_xticklabels() 00087 pl.setp(l,fontsize="x-small") 00088 l=ax.get_yticklabels() 00089 pl.setp(l,fontsize="x-small") 00090 if self.pmulti!=0: 00091 self.pmulti[2] += 1 00092 multi=self.pmulti 00093 if multi[2] <= multi[0]*multi[1]: 00094 pl.subplot(multi[0],multi[1],multi[2]) 00095 # consider pl.draw() here - may be slow 00096 00097 ########################################################### 00098 00099 def endfig(self,show=True,filename=""): # set margins to smaller, save to file if required 00100 ax=pl.gca() 00101 l=ax.get_xticklabels() 00102 pl.setp(l,fontsize="x-small") 00103 l=ax.get_yticklabels() 00104 pl.setp(l,fontsize="x-small") 00105 if show: 00106 pl.draw() 00107 if len(filename)>0: 00108 pl.savefig(filename) 00109 pl.ioff() 00110 00111 00112 00113 ########################################################### 00114 00115 def msg(self, s, origin=None, priority=None): 00116 # ansi color codes: 00117 # Foreground colors 00118 # 30 Black 00119 # 31 Red 00120 # 32 Green 00121 # 33 Yellow 00122 # 34 Blue 00123 # 35 Magenta 00124 # 36 Cyan 00125 # 37 White 00126 toterm=False 00127 if priority==None: 00128 clr="\x1b[32m" 00129 priority="INFO" 00130 else: 00131 priority=priority.upper() 00132 if priority=="WARN": 00133 clr="\x1b[35m" 00134 toterm=True 00135 priority="INFO" # otherwise casalog will spew to term also. 00136 else: 00137 if priority=="ERROR": 00138 clr="\x1b[31m" 00139 toterm=False # casalog spews severe to term already 00140 else: 00141 if not (priority=="DEBUG" or priority[:-1]=="DEBUG"): 00142 priority="INFO" 00143 bw="\x1b[0m" 00144 if origin==None: 00145 origin="simutil" 00146 if toterm: 00147 print clr+"["+origin+"] "+bw+s 00148 if priority=="ERROR": 00149 #return False 00150 raise Exception, s 00151 else: 00152 casalog.post(s,priority=priority,origin=origin) 00153 00154 00155 00156 ########################################################### 00157 00158 def isquantity(self,s,halt=True): 00159 if type(s)!=type([]): 00160 t=[s] 00161 else: 00162 t=s 00163 for t0 in t: 00164 if not (len(t0)>0 and qa.isquantity(t0)): 00165 if halt: 00166 self.msg("can't interpret '"+str(t0)+"' as a CASA quantity",priority="error") 00167 return False 00168 return True 00169 00170 ########################################################### 00171 00172 def isdirection(self,s,halt=True): 00173 if type(s)==type([]): 00174 t=s[0] 00175 else: 00176 t=s 00177 try: 00178 x=self.direction_splitter(s) 00179 y=me.direction(x[0],x[1],x[2]) 00180 except: 00181 if halt: 00182 self.msg("can't interpret '"+str(s)+"' as a direction",priority="error") 00183 return False 00184 if not me.measure(y,y['refer']): 00185 if halt: 00186 self.msg("can't interpret '"+str(s)+"' as a direction",priority="error") 00187 return False 00188 return True 00189 00190 ########################################################### 00191 00192 def ismstp(self,s,halt=False): 00193 try: 00194 istp = False 00195 # check if the ms is tp data or not. 00196 tb.open(s+'/ANTENNA') 00197 antname = tb.getcol('NAME') 00198 tb.close() 00199 if antname[0].find('TP') > -1: 00200 istp = True 00201 elif len(antname) == 1: 00202 istp = True 00203 else: 00204 # need complete testing of UVW 00205 tb.open(s) 00206 uvw = tb.getcol("UVW") 00207 tb.close() 00208 if uvw.all() == 0.: 00209 istp = True 00210 except: 00211 if halt: 00212 self.msg("can't understand the file '"+str(s)+"'",priority="error") 00213 return False 00214 if not istp: 00215 if halt: 00216 self.msg("input file '"+str(s)+"' is not a totalpower ms",priority="error") 00217 return False 00218 return True 00219 00220 00221 00222 ########################################################### 00223 # plot an image (optionally), and calculate its statistics 00224 00225 def statim(self,image,plot=True,incell=None,disprange=None,bar=True,showstats=True): 00226 pix=self.cellsize(image) # cell positive by convention 00227 pixarea=abs(qa.convert(pix[0],'arcsec')['value']* 00228 qa.convert(pix[1],'arcsec')['value']) 00229 ia.open(image) 00230 imunit=ia.brightnessunit() 00231 if imunit == 'Jy/beam': 00232 bm=ia.restoringbeam() 00233 if len(bm)>0: 00234 toJyarcsec=1./pixarea 00235 else: 00236 toJyarcsec=1. 00237 toJypix=toJyarcsec*pixarea 00238 elif imunit == 'Jy/pixel': 00239 toJyarcsec=1./pixarea 00240 toJypix=1. 00241 else: 00242 self.msg("%s: unknown units" % image,origin="statim") 00243 toJyarcsec=1. 00244 toJypix=1. 00245 stats=ia.statistics(robust=True,verbose=False,list=False) 00246 im_min=stats['min']*toJypix 00247 plarr=pl.zeros(1) 00248 badim=False 00249 if type(im_min)==type([]) or type(im_min)==type(plarr): 00250 if len(im_min)<1: 00251 badim=True 00252 im_min=0. 00253 im_max=stats['max']*toJypix 00254 if type(im_max)==type([]) or type(im_max)==type(plarr): 00255 if len(im_max)<1: 00256 badim=True 00257 im_max=1. 00258 imsize=ia.shape()[0:2] 00259 reg1=rg.box([0,0],[imsize[0]*.25,imsize[1]*.25]) 00260 stats=ia.statistics(region=reg1,verbose=False,list=False) 00261 im_rms=stats['rms']*toJypix 00262 if type(im_rms)==type([]) or type(im_rms)==type(plarr): 00263 if len(im_rms)==0: 00264 badim=True 00265 im_rms=0. 00266 data_array=ia.getchunk([-1,-1,1,1],[-1,-1,1,1],[1],[],True,True,False) 00267 data_array=pl.array(data_array) 00268 tdata_array=pl.transpose(data_array) 00269 ttrans_array=tdata_array.tolist() 00270 ttrans_array.reverse() 00271 if (plot): 00272 pixsize=[qa.convert(pix[0],'arcsec')['value'],qa.convert(pix[1],'arcsec')['value']] 00273 xextent=imsize[0]*abs(pixsize[0])*0.5 00274 yextent=imsize[1]*abs(pixsize[1])*0.5 00275 if self.verbose: 00276 self.msg("plotting %fx%f\" im with %fx%f\" pix" % 00277 (xextent,yextent,pixsize[0],pixsize[1]),origin="statim") 00278 xextent=[xextent,-xextent] 00279 yextent=[-yextent,yextent] 00280 # remove top .5% of pixels: 00281 nbin=200 00282 if badim: 00283 highvalue=im_max 00284 lowvalue=im_min 00285 else: 00286 imhist=ia.histograms(cumu=True,nbins=nbin,list=False)#['histout'] 00287 ii=0 00288 lowcounts=imhist['counts'][ii] 00289 while imhist['counts'][ii]<0.005*lowcounts and ii<nbin: 00290 ii=ii+1 00291 lowvalue=imhist['values'][ii] 00292 ii=nbin-1 00293 highcounts=imhist['counts'][ii] 00294 while imhist['counts'][ii]>0.995*highcounts and ii>0 and 0.995*highcounts>lowcounts: 00295 ii=ii-1 00296 highvalue=imhist['values'][ii] 00297 if disprange != None: 00298 if type(disprange)==type([]): 00299 if len(disprange)>0: 00300 highvalue=disprange[-1] 00301 if len(disprange)>1: 00302 lowvalue=disprange[0] 00303 if len(disprange)>2: 00304 throw("internal error disprange="+str(disprange)+" has too many elements") 00305 else: # if passed an empty list [], return low.high 00306 disprange.append(lowvalue) 00307 disprange.append(highvalue) 00308 else: 00309 highvalue=disprange # assume if scalar passed its the max 00310 00311 if plot: 00312 img=pl.imshow(ttrans_array,interpolation='bilinear',cmap=pl.cm.jet,extent=xextent+yextent,vmax=highvalue,vmin=lowvalue) 00313 ax=pl.gca() 00314 #l=ax.get_xticklabels() 00315 #pl.setp(l,fontsize="x-small") 00316 #l=ax.get_yticklabels() 00317 #pl.setp(l,fontsize="x-small") 00318 foo=image.split("/") 00319 if len(foo)==1: 00320 imagestrip=image 00321 else: 00322 imagestrip=foo[1] 00323 pl.title(imagestrip,fontsize="x-small") 00324 if showstats: 00325 pl.text(0.05,0.95,"min=%7.1e\nmax=%7.1e\nRMS=%7.1e\n%s" % (im_min,im_max,im_rms,imunit),transform = ax.transAxes,bbox=dict(facecolor='white', alpha=0.7),size="x-small",verticalalignment="top") 00326 if bar: 00327 cb=pl.colorbar(pad=0) 00328 cl = pl.getp(cb.ax,'yticklabels') 00329 pl.setp(cl,fontsize='x-small') 00330 ia.done() 00331 return im_min,im_max,im_rms,imunit 00332 00333 00334 00335 00336 00337 00338 00339 ########################################################### 00340 00341 def calc_pointings2(self, spacing, size, maptype="hex", direction=None, relmargin=0.5, beam=0.): 00342 """ 00343 If direction is a list, simply returns direction and the number of 00344 pointings in it. 00345 00346 Otherwise, returns a hexagonally packed list of pointings separated by 00347 spacing and fitting inside an area specified by direction and mapsize, 00348 as well as the number of pointings. The hexagonal packing starts with a 00349 horizontal row centered on direction, and the other rows alternate 00350 being horizontally offset by a half spacing. 00351 """ 00352 # make size 2-dimensional and ensure it is quantity 00353 if type(size) != type([]): 00354 size=[size,size] 00355 if len(size) <2: 00356 size=[size[0],size[0]] 00357 self.isquantity(size) 00358 00359 # parse and check direction 00360 if direction==None: 00361 # if no direction is specified, use the object's direction 00362 direction=self.direction 00363 else: 00364 # if one is specified, use it to set the object's direction 00365 self.direction=direction 00366 self.isdirection(direction) 00367 00368 # direction is always a list of strings (defined by .xml) 00369 if type(direction)==type([]): 00370 if len(direction) > 1: 00371 if self.verbose: self.msg("you are inputing the precise pointings in 'direction' - if you want to calculate a mosaic, give a single direction string and set maptype",priority="warn") 00372 return direction 00373 else: direction=direction[0] 00374 00375 00376 # haveing eliminated other options, we need to calculate: 00377 epoch, centx, centy = self.direction_splitter() 00378 00379 shorttype=str.upper(maptype[0:3]) 00380 # if not shorttype=="HEX": 00381 # self.msg("can't calculate map of maptype "+maptype,priority="error") 00382 if shorttype == "HEX": 00383 # this is hexagonal grid - Kana will add other types here 00384 self.isquantity(spacing) 00385 spacing = qa.quantity(spacing) 00386 yspacing = qa.mul(0.866025404, spacing) 00387 00388 xsize=qa.quantity(size[0]) 00389 ysize=qa.quantity(size[1]) 00390 00391 nrows = 1+ int(pl.floor(qa.convert(qa.div(ysize, yspacing), '')['value'] 00392 - 2.309401077 * relmargin)) 00393 00394 availcols = 1 + qa.convert(qa.div(xsize, spacing), 00395 '')['value'] - 2.0 * relmargin 00396 ncols = int(pl.floor(availcols)) 00397 00398 # By making the even rows shifted spacing/2 ahead, and possibly shorter, 00399 # the top and bottom rows (nrows odd), are guaranteed to be short. 00400 if availcols - ncols >= 0.5 and nrows>1: # O O O 00401 evencols = ncols # O O O 00402 ncolstomin = 0.5 * (ncols - 0.5) 00403 else: 00404 evencols = ncols - 1 # O O 00405 ncolstomin = 0.5 * (ncols - 1) # O O O 00406 pointings = [] 00407 00408 # Start from the top because in the Southern hemisphere it sets first. 00409 y = qa.add(centy, qa.mul(0.5 * (nrows - 1), yspacing)) 00410 for row in xrange(0, nrows): # xrange stops early. 00411 xspacing = qa.mul(1.0 / pl.cos(qa.convert(y, 'rad')['value']),spacing) 00412 ystr = qa.formxxx(y, format='dms',prec=5) 00413 00414 if row % 2: # Odd 00415 xmin = qa.sub(centx, qa.mul(ncolstomin, xspacing)) 00416 stopcolp1 = ncols 00417 else: # Even (including 0) 00418 xmin = qa.sub(centx, qa.mul(ncolstomin - 0.5, 00419 xspacing)) 00420 stopcolp1 = evencols 00421 for col in xrange(0, stopcolp1): # xrange stops early. 00422 x = qa.formxxx(qa.add(xmin, qa.mul(col, xspacing)), 00423 format='hms',prec=5) 00424 pointings.append("%s%s %s" % (epoch, x, ystr)) 00425 y = qa.sub(y, yspacing) 00426 elif shorttype =="SQU": 00427 # lattice gridding 00428 self.isquantity(spacing) 00429 spacing = qa.quantity(spacing) 00430 yspacing = spacing 00431 00432 xsize=qa.quantity(size[0]) 00433 ysize=qa.quantity(size[1]) 00434 00435 nrows = 1+ int(pl.floor(qa.convert(qa.div(ysize, yspacing), '')['value'] 00436 - 2.0 * relmargin)) 00437 00438 availcols = 1 + qa.convert(qa.div(xsize, spacing), '')['value'] \ 00439 - 2.0 * relmargin 00440 ncols = int(pl.floor(availcols)) 00441 00442 00443 ncolstomin = 0.5 * (ncols - 1) # O O O O 00444 pointings = [] # O O O O 00445 00446 # Start from the top because in the Southern hemisphere it sets first. 00447 y = qa.add(centy, qa.mul(0.5 * (nrows - 1), yspacing)) 00448 for row in xrange(0, nrows): # xrange stops early. 00449 xspacing = qa.mul(1.0 / pl.cos(qa.convert(y, 'rad')['value']),spacing) 00450 ystr = qa.formxxx(y, format='dms',prec=5) 00451 00452 xmin = qa.sub(centx, qa.mul(ncolstomin, xspacing)) 00453 stopcolp1 = ncols 00454 00455 for col in xrange(0, stopcolp1): # xrange stops early. 00456 x = qa.formxxx(qa.add(xmin, qa.mul(col, xspacing)), 00457 format='hms',prec=5) 00458 pointings.append("%s%s %s" % (epoch, x, ystr)) 00459 y = qa.sub(y, yspacing) 00460 if shorttype == "ALM": 00461 # OT algorithm 00462 self.isquantity(spacing) 00463 spacing = qa.quantity(spacing) 00464 xsize = qa.quantity(size[0]) 00465 ysize = qa.quantity(size[1]) 00466 00467 spacing_asec=qa.convert(spacing,'arcsec')['value'] 00468 xsize_asec=qa.convert(xsize,'arcsec')['value'] 00469 ysize_asec=qa.convert(ysize,'arcsec')['value'] 00470 angle = 0. # deg 00471 00472 if str.upper(maptype[0:8]) == 'ALMA2012': 00473 x,y = self.getTrianglePoints(xsize_asec, ysize_asec, angle, spacing_asec) 00474 else: 00475 if beam<=0: 00476 beam=spacing_asec*1.2*sqrt(3) # ASSUMES Nyquist and arcsec 00477 x,y = self.getTriangularTiling(xsize_asec, ysize_asec, angle, spacing_asec, beam) 00478 00479 pointings = [] 00480 nx=len(x) 00481 for i in range(nx): 00482 # Start from the top because in the Southern hemisphere it sets first. 00483 y1=qa.add(centy, str(y[nx-i-1])+"arcsec") 00484 ycos=pl.cos(qa.convert(y1,"rad")['value']) 00485 ystr = qa.formxxx(y1, format='dms',prec=5) 00486 xstr = qa.formxxx(qa.add(centx, str(x[nx-i-1]/ycos)+"arcsec"), format='hms',prec=5) 00487 pointings.append("%s%s %s" % (epoch, xstr, ystr)) 00488 00489 00490 # if could not fit any pointings, then return single pointing 00491 if(len(pointings)==0): 00492 pointings.append(direction) 00493 00494 self.msg("using %i generated pointing(s)" % len(pointings)) 00495 self.pointings=pointings 00496 return pointings 00497 00498 00499 00500 00501 00502 00503 00504 00505 00506 00507 ########################################################### 00508 00509 def read_pointings(self, filename): 00510 """ 00511 read pointing list from file containing epoch, ra, dec, 00512 and scan time (optional,in sec). 00513 Parameter: 00514 filename: (str) the name of input file 00515 00516 The input file (ASCII) should contain at least 3 fields separated 00517 by a space which specify positions with epoch, ra and dec (in dms 00518 or hms). 00519 The optional field, time, shoud be a list of decimal numbers 00520 which specifies integration time at each position in second. 00521 The lines which start with '#' is ignored and can be used 00522 as comment lines. 00523 00524 Example of an input file: 00525 #Epoch RA DEC TIME(optional) 00526 J2000 23h59m28.10 -019d52m12.35 10.0 00527 J2000 23h59m32.35 -019d52m12.35 10.0 00528 J2000 23h59m36.61 -019d52m12.35 60.0 00529 00530 """ 00531 f=open(filename) 00532 line= ' ' 00533 time=[] 00534 pointings=[] 00535 00536 # add option of different epoch in a header line like read_antenna? 00537 00538 while (len(line)>0): 00539 try: 00540 line=f.readline() 00541 if not line.startswith('#'): 00542 ### ignoring line that has less than 3 elements 00543 if(len(line.split()) >2): 00544 splitline=line.split() 00545 epoch=splitline[0] 00546 ra0=splitline[1] 00547 de0=splitline[2] 00548 if len(splitline)>3: 00549 time.append(float(splitline[3])) 00550 else: 00551 time.append(0.) 00552 xstr = qa.formxxx(qa.quantity(ra0), format='hms',prec=5) 00553 ystr = qa.formxxx(qa.quantity(de0), format='dms',prec=5) 00554 pointings.append("%s %s %s" % (epoch,xstr,ystr)) 00555 except: 00556 break 00557 f.close() 00558 00559 # need an error check here if zero valid pointings were read 00560 if len(pointings) < 1: 00561 s="No valid lines found in pointing file" 00562 self.msg(s,priority="error") 00563 self.msg("read in %i pointing(s) from file" % len(pointings)) 00564 self.pointings=pointings 00565 #self.direction=pointings 00566 00567 return len(pointings), pointings, time 00568 00569 00570 00571 00572 00573 00574 00575 ########################################################### 00576 00577 def write_pointings(self, filename,pointings,time=1.): 00578 """ 00579 write pointing list to file containing epoch, ra, dec, 00580 and scan time (optional,in sec). 00581 00582 Example of an output file: 00583 #Epoch RA DEC TIME(optional) 00584 J2000 23h59m28.10 -019d52m12.35 10.0 00585 J2000 23h59m32.35 -019d52m12.35 10.0 00586 J2000 23h59m36.61 -019d52m12.35 60.0 00587 00588 """ 00589 f=open(filename,"write") 00590 f.write('#Epoch RA DEC TIME[sec]\n') 00591 if type(pointings)!=type([]): 00592 pointings=[pointings] 00593 npos=len(pointings) 00594 if type(time)!=type([]): 00595 time=[time] 00596 if len(time)==1: 00597 time=list(time[0] for x in range(npos)) 00598 00599 for i in range(npos): 00600 #epoch,ra,dec=direction_splitter(pointings[i]) 00601 #xstr = qa.formxxx(qa.quantity(ra), format='hms') 00602 #ystr = qa.formxxx(qa.quantity(dec), format='dms') 00603 #line = "%s %s %s" % (epoch,xstr,ystr) 00604 #self.isdirection(line) # extra check 00605 #f.write(line+" "+str(time[i])+"\n") 00606 f.write(pointings[i]+" "+str(time[i])+"\n") 00607 00608 f.close() 00609 return 00610 00611 00612 00613 ########################################################### 00614 00615 def average_direction(self, directions=None): 00616 # RI TODO make deal with list of measures as well as list of strings 00617 """ 00618 Returns the average of directions as a string, and relative offsets 00619 """ 00620 if directions==None: 00621 directions=self.direction 00622 epoch0, x, y = self.direction_splitter(directions[0]) 00623 i = 1 00624 avgx = 0.0 00625 avgy = 0.0 00626 for drn in directions: 00627 epoch, x, y = self.direction_splitter(drn) 00628 # in principle direction_splitter returns directions in degrees, 00629 # but can we be sure? 00630 x=qa.convert(x,'deg') 00631 y=qa.convert(y,'deg') 00632 x = x['value'] 00633 y = y['value'] 00634 if epoch != epoch0: # Paranoia 00635 print "[simutil] WARN: precession not handled by average_direction()" 00636 x = self.wrapang(x, avgx, 360.0) 00637 avgx += (x - avgx) / i 00638 avgy += (y - avgy) / i 00639 i += 1 00640 offsets=pl.zeros([2,i-1]) 00641 i=0 00642 cosdec=pl.cos(avgy*pl.pi/180.) 00643 for drn in directions: 00644 epoch, x, y = self.direction_splitter(drn) 00645 x=qa.convert(x,'deg') 00646 y=qa.convert(y,'deg') 00647 x = x['value'] 00648 y = y['value'] 00649 x = self.wrapang(x, avgx, 360.0) 00650 offsets[:,i]=[(x-avgx)*cosdec,y-avgy] # apply cosdec to make offsets on sky 00651 i+=1 00652 avgx = qa.toangle('%17.12fdeg' % avgx) 00653 avgy = qa.toangle('%17.12fdeg' % avgy) 00654 avgx = qa.formxxx(avgx, format='hms',prec=5) 00655 avgy = qa.formxxx(avgy, format='dms',prec=5) 00656 return "%s%s %s" % (epoch0, avgx, avgy), offsets 00657 00658 00659 00660 ########################################################### 00661 00662 def median_direction(self, directions=None): 00663 # RI TODO make deal with list of measures as well as list of strings 00664 """ 00665 Returns the median of directions as a string, and relative offsets 00666 """ 00667 if directions==None: 00668 directions=self.direction 00669 epoch0, x, y = self.direction_splitter(directions[0]) 00670 i = 1 00671 avgx = 0.0 00672 avgy = 0.0 00673 xx=[] 00674 yy=[] 00675 for drn in directions: 00676 epoch, x, y = self.direction_splitter(drn) 00677 # in principle direction_splitter returns directions in degrees, 00678 # but can we be sure? 00679 x=qa.convert(x,'deg') 00680 y=qa.convert(y,'deg') 00681 x = x['value'] 00682 y = y['value'] 00683 if epoch != epoch0: # Paranoia 00684 print "[simutil] WARN: precession not handled by average_direction()" 00685 x = self.wrapang(x, avgx, 360.0) 00686 xx.append(x) 00687 yy.append(y) 00688 i += 1 00689 avgx = pl.median(xx) 00690 avgy = pl.median(yy) 00691 offsets=pl.zeros([2,i-1]) 00692 i=0 00693 cosdec=pl.cos(avgy*pl.pi/180.) 00694 for drn in directions: 00695 epoch, x, y = self.direction_splitter(drn) 00696 x=qa.convert(x,'deg') 00697 y=qa.convert(y,'deg') 00698 x = x['value'] 00699 y = y['value'] 00700 x = self.wrapang(x, avgx, 360.0) 00701 offsets[:,i]=[(x-avgx)*cosdec,y-avgy] # apply cosdec to make offsets on sky 00702 i+=1 00703 avgx = qa.toangle('%17.12fdeg' % avgx) 00704 avgy = qa.toangle('%17.12fdeg' % avgy) 00705 avgx = qa.formxxx(avgx, format='hms',prec=5) 00706 avgy = qa.formxxx(avgy, format='dms',prec=5) 00707 return "%s%s %s" % (epoch0, avgx, avgy), offsets 00708 00709 00710 00711 ########################################################### 00712 00713 def direction_splitter(self, direction=None): 00714 """ 00715 Given a direction, return its epoch, x, and y parts. Epoch will be '' 00716 if absent, or '%s ' % epoch if present. x and y will be angle qa's in 00717 degrees. 00718 """ 00719 import re 00720 if direction == None: 00721 direction=self.direction 00722 if type(direction) == type([]): 00723 direction=self.average_direction(direction)[0] 00724 dirl = direction.split() 00725 if len(dirl) == 3: 00726 epoch = dirl[0] + ' ' 00727 else: 00728 epoch = '' 00729 # x, y = map(qa.toangle, dirl[-2:]) 00730 x=qa.toangle(dirl[1]) 00731 # qa is stupid when it comes to dms vs hms, and assumes hms with colons and dms for periods. 00732 decstr=dirl[2] 00733 # parse with regex to get three numbers and reinstall them as dms 00734 q=re.compile('([+-]?\d+).(\d+).(\d+\.?\d*)') 00735 qq=q.match(decstr) 00736 z=qq.groups() 00737 decstr=z[0]+'d' 00738 if len(z)>1: 00739 decstr=decstr+z[1]+'m' 00740 if len(z)>2: 00741 decstr=decstr+z[2]+'s' 00742 y=qa.toangle(decstr) 00743 00744 return epoch, qa.convert(x, 'deg'), qa.convert(y, 'deg') 00745 00746 00747 ########################################################### 00748 00749 def dir_s2m(self, direction=None): 00750 """ 00751 Given a direction as a string 'refcode lon lat', return it as qa measure. 00752 """ 00753 if direction == None: 00754 direction=self.direction 00755 if type(direction) == type([]): 00756 direction=self.average_direction(direction)[0] 00757 dirl = direction.split() 00758 if len(dirl) == 3: 00759 refcode = dirl[0] + ' ' 00760 else: 00761 refcode = 'J2000' 00762 if self.verbose: self.msg("assuming J2000 for "+direction,origin="simutil.s2m") 00763 x, y = map(qa.quantity, dirl[-2:]) 00764 if x['unit'] == '': x['unit']='deg' 00765 if y['unit'] == '': y['unit']='deg' 00766 return me.direction(refcode,qa.toangle(x),qa.toangle(y)) 00767 00768 00769 ########################################################### 00770 00771 def dir_m2s(self, dir): 00772 """ 00773 Given a direction as a measure, return it as astring 'refcode lon lat'. 00774 """ 00775 if dir['type'] != 'direction': 00776 self.msg("converting direction measure",priority="error",origin="simutil.m2s") 00777 return False 00778 ystr = qa.formxxx(dir['m1'], format='dms',prec=5) 00779 xstr = qa.formxxx(dir['m0'], format='hms',prec=5) 00780 return "%s %s %s" % (dir['refer'], xstr, ystr) 00781 00782 ########################################################### 00783 00784 def wrapang(self, ang, target, period = 360.0): 00785 """ 00786 Returns ang wrapped so that it is within +-period/2 of target. 00787 """ 00788 dang = ang - target 00789 period = pl.absolute(period) 00790 halfperiod = 0.5 * period 00791 if pl.absolute(dang) > halfperiod: 00792 nwraps = pl.floor(0.5 + float(dang) / period) 00793 ang -= nwraps * period 00794 return ang 00795 00796 00797 00798 00799 00800 00801 00802 00803 00804 00805 ########################################################### 00806 #========================== tsys ========================== 00807 00808 def noisetemp(self, telescope=None, freq=None, 00809 diam=None, epsilon=None): 00810 00811 """ 00812 Noise temperature and efficiencies for several telescopes: 00813 ALMA, ACA, EVLA, VLA, and SMA 00814 Input: telescope name, frequency as a quantity string "300GHz", 00815 dish diameter (optional - knows diameters for arrays above) 00816 epsilon = rms surface accuracy in microns (also optional - 00817 this method contains the spec values for each telescope) 00818 Output: eta_p phase efficieny (from Ruze formula), 00819 eta_s spill (main beam) efficiency, 00820 eta_b geometrical blockage efficiency, 00821 eta_t taper efficiency, 00822 eta_q correlator efficiency including quantization, 00823 t_rx reciever temperature. 00824 antenna efficiency = eta_p * eta_s * eta_b * eta_t 00825 Notes: VLA correlator efficiency includes waveguide loss 00826 EVLA correlator efficiency is probably optimistic at 0.88 00827 """ 00828 00829 if telescope==None: telescope=self.telescopename 00830 telescope=str.upper(telescope) 00831 00832 obs =['ALMA','ACA','EVLA','VLA','SMA'] 00833 d =[ 12. ,7., 25. , 25. , 6. ] 00834 ds =[ 0.75, 0.75, 0.364, 0.364,0.35] # subreflector size for ACA? 00835 eps =[ 25., 20., 300, 300 ,15. ] # antenna surface accuracy 00836 00837 cq =[ 0.845, 0.845, 0.88, 0.79, 0.86] # correlator eff 00838 # ALMA includes quantization eff of 0.96 00839 # VLA includes additional waveguide loss from correlator loss of 0.809 00840 # EVLA is probably optimistic 00841 00842 # things hardcoded in ALMA etimecalculator 00843 # t_ground=270. 00844 # t_cmb=2.73 00845 # eta_q*eta_corr = 0.88*.961 00846 # eta_ap = 0.72*eta_ruze 00847 00848 if obs.count(telescope)>0: 00849 iobs=obs.index(telescope) 00850 else: 00851 if self.verbose: self.msg("I don't know much about "+telescope+" so I'm going to use ALMA specs") 00852 iobs=0 # ALMA is the default ;) 00853 00854 if diam==None: diam=d[iobs] 00855 diam_subreflector=ds[iobs] 00856 if self.verbose: self.msg("subreflector diameter="+str(diam_subreflector),origin="noisetemp") 00857 00858 # blockage efficiency. 00859 eta_b = 1.-(diam_subreflector/diam)**2 00860 00861 # spillover efficiency. 00862 eta_s = 0.95 # these are ALMA values 00863 # taper efficiency. 00864 #eta_t = 0.86 # these are ALMA values 00865 eta_t = 0.819 # 20100914 OT value 00866 eta_t = 0.72 00867 00868 # Ruze phase efficiency. 00869 if epsilon==None: epsilon = eps[iobs] # microns RMS 00870 if freq==None: 00871 freq_ghz=qa.convert(qa.quantity(self.centerfreq),'GHz') 00872 bw_ghz=qa.convert(qa.quantity(self.bandwidth),'GHz') 00873 else: 00874 freq_ghz=qa.convert(qa.quantity(freq),'GHz') 00875 eta_p = pl.exp(-(4.0*3.1415926535*epsilon*freq_ghz.get("value")/2.99792458e5)**2) 00876 if self.verbose: self.msg("ruze phase efficiency for surface accuracy of "+str(epsilon)+"um = " + str(eta_p) + " at "+str(freq),origin="noisetemp") 00877 00878 # antenna efficiency 00879 # eta_a = eta_p*eta_s*eta_b*eta_t 00880 00881 # correlator quantization efficiency. 00882 eta_q = cq[iobs] 00883 00884 # Receiver radiation temperature in K. 00885 if telescope=='ALMA' or telescope=='ACA': 00886 # ALMA-40.00.00.00-001-A-SPE.pdf 00887 # http://www.eso.org/sci/facilities/alma/system/frontend/ 00888 00889 # limits instead of centers, go to higher band in gaps 00890 f0=[31.3,45,84,116,163,211,275,373,500,720] 00891 # cycle 1 OT values 7/12 00892 t0=[ 17, 30, 45, 51, 65, 55, 75, 196, 100, 230] 00893 00894 flim=[31.3,950] 00895 if self.verbose: self.msg("using ALMA/ACA Rx specs",origin="noisetemp") 00896 else: 00897 if telescope=='EVLA': 00898 # 201009114 from rick perley: 00899 # f0=[1.5,3,6,10,15,23,33,45] 00900 t0=[10.,15,12,15,10,12,15,28] 00901 # limits 00902 f0=[1,2,4,8,12,18,26.5,40,50] 00903 00904 flim=[0.8,50] 00905 if self.verbose: self.msg("using EVLA Rx specs",origin="noisetemp") 00906 else: 00907 if telescope=='VLA': 00908 # http://www.vla.nrao.edu/genpub/overview/ 00909 # f0=[0.0735,0.32, 1.5, 4.75, 8.4, 14.9, 23, 45 ] 00910 # t0=[5000, 165, 56, 44, 34, 110, 110, 110] 00911 # exclude P band for now 00912 # f0=[0.32, 1.5, 4.75, 8.4, 14.9, 23, 45 ] 00913 # limits 00914 # http://www.vla.nrao.edu/genpub/overview/ 00915 f0=[0.30,0.34,1.73,5,8.8,15.4,24,50] 00916 t0=[165, 56, 44, 34, 110, 110, 110] 00917 flim=[0.305,50] 00918 if self.verbose: self.msg("using old VLA Rx specs",origin="noisetemp") 00919 else: 00920 if telescope=='SMA': 00921 # f0=[212.,310.,383.,660.] 00922 # limits 00923 f0=[180,250,320,420,720] 00924 t0=[67, 116, 134, 500] 00925 flim=[180.,720] 00926 else: 00927 self.msg("I don't know about the "+telescope+" receivers, using 200K",priority="warn",origin="noisetemp") 00928 f0=[10,900] 00929 t0=[200,200] 00930 flim=[0,5000] 00931 00932 00933 obsfreq=freq_ghz.get("value") 00934 # z=pl.where(abs(obsfreq-pl.array(f0)) == min(abs(obsfreq-pl.array(f0)))) 00935 # t_rx=t0[z[0]] 00936 z=0 00937 while(f0[z]<obsfreq and z<len(t0)): 00938 z+=1 00939 t_rx=t0[z-1] 00940 00941 if obsfreq<flim[0]: 00942 self.msg("observing freqency is lower than expected for "+telescope,priority="warn",origin="noise") 00943 self.msg("proceeding with extrapolated receiver temp="+str(t_rx),priority="warn",origin="noise") 00944 if obsfreq>flim[1]: 00945 self.msg("observing freqency is higher than expected for "+telescope,priority="warn",origin="noise") 00946 self.msg("proceeding with extrapolated receiver temp="+str(t_rx),priority="warn",origin="noise") 00947 if obsfreq<=flim[1] and obsfreq>=flim[0]: 00948 self.msg("interpolated receiver temp="+str(t_rx),origin="noise") 00949 00950 return eta_p, eta_s, eta_b, eta_t, eta_q, t_rx 00951 00952 00953 00954 00955 00956 00957 def sensitivity(self, freq, bandwidth, etime, elevation, pwv=None, 00958 telescope=None, diam=None, nant=None, 00959 antennalist=None, 00960 doimnoise=None, 00961 integration=None,debug=None, 00962 method="tsys-atm",tau0=None,t_sky=None): 00963 00964 import glob 00965 tmpname="tmp"+str(os.getpid()) 00966 i=0 00967 while i<10 and len(glob.glob(tmpname+"*"))>0: 00968 tmpname="tmp"+str(os.getpid())+str(i) 00969 i=i+1 00970 if i>=9: 00971 xx=glob.glob(tmpname+"*") 00972 for k in range(len(xx)): 00973 if os.path.isdir(xx[k]): 00974 cu.removetable(xx[k]) 00975 else: 00976 os.remove(xx[k]) 00977 00978 msfile=tmpname+".ms" 00979 sm.open(msfile) 00980 00981 rxtype=0 # 2SB 00982 if antennalist==None: 00983 if telescope==None: 00984 self.msg("Telescope name has not been set.",priority="error") 00985 return False 00986 if diam==None: 00987 self.msg("Antenna diameter has not been set.",priority="error") 00988 return False 00989 if nant==None: 00990 self.msg("Number of antennas has not been set.",priority="error") 00991 return False 00992 00993 posobs=me.observatory(telescope) 00994 obs=me.measure(posobs,'WGS84') 00995 obslat=qa.convert(obs['m1'],'deg')['value'] 00996 obslon=qa.convert(obs['m0'],'deg')['value'] 00997 obsalt=qa.convert(obs['m2'],'m')['value'] 00998 stnx,stny,stnz = self.locxyz2itrf(obslat,obslon,0,0,obsalt) 00999 antnames="A00" 01000 01001 else: 01002 if str.upper(antennalist[0:4])=="ALMA": 01003 tail=antennalist[5:] 01004 if self.isquantity(tail,halt=False): 01005 resl=qa.convert(tail,"arcsec")['value'] 01006 repodir=os.getenv("CASAPATH").split(' ')[0]+"/data/alma/simmos/" 01007 if os.path.exists(repodir): 01008 confnum=(2.867-pl.log10(resl*1000*qa.convert(freq,"GHz")['value']/672.))/0.0721 01009 confnum=max(1,min(28,confnum)) 01010 conf=str(int(round(confnum))) 01011 if len(conf)<2: conf='0'+conf 01012 antennalist=repodir+"alma.out"+conf+".cfg" 01013 self.msg("converted resolution to antennalist "+antennalist) 01014 01015 if os.path.exists(antennalist): 01016 stnx, stny, stnz, stnd, padnames, nant, telescope = self.readantenna(antennalist) 01017 else: 01018 self.msg("antennalist "+antennalist+" not found",priority="error") 01019 return False 01020 01021 # RI TODO average antenna instead of first? 01022 diam = stnd[0] 01023 antnames=padnames 01024 01025 posobs=me.observatory(telescope) 01026 obs=me.measure(posobs,'WGS84') 01027 #obslat=qa.convert(obs['m1'],'deg')['value'] 01028 #obslon=qa.convert(obs['m0'],'deg')['value'] 01029 #obsalt=qa.convert(obs['m2'],'m')['value'] 01030 01031 01032 if (telescope==None or diam==None): 01033 self.msg("Telescope name and antenna diameter have not been set.",priority="error") 01034 return False 01035 01036 # copied from task_simdata: 01037 01038 self.setcfg(sm, telescope, stnx, stny, stnz, diam, 01039 padnames, posobs) 01040 01041 model_nchan=1 01042 # RI TODO isquantity checks 01043 model_width=qa.quantity(bandwidth) # note: ATM uses band center 01044 01045 # start is center of first channel. for nch=1, that equals center 01046 model_start=qa.quantity(freq) 01047 01048 stokes, feeds = self.polsettings(telescope) 01049 sm.setspwindow(spwname="band1", freq=qa.tos(model_start), 01050 deltafreq=qa.tos(model_width), 01051 freqresolution=qa.tos(model_width), 01052 nchannels=model_nchan, stokes=stokes) 01053 sm.setfeed(mode=feeds, pol=['']) 01054 01055 sm.setlimits(shadowlimit=0.01, elevationlimit='10deg') 01056 sm.setauto(0.0) 01057 01058 obslat=qa.convert(obs['m1'],'deg') 01059 dec=qa.add(obslat, qa.add(qa.quantity("90deg"),qa.mul(elevation,-1))) 01060 01061 sm.setfield(sourcename="src1", 01062 sourcedirection="J2000 00:00:00.00 "+qa.angle(dec)[0], 01063 calcode="OBJ", distance='0m') 01064 reftime = me.epoch('TAI', "2012/01/01/00:00:00") 01065 if integration==None: 01066 integration=qa.mul(etime,0.01) 01067 self.msg("observing for "+qa.tos(etime)+" with integration="+qa.tos(integration)) 01068 sm.settimes(integrationtime=integration, usehourangle=True, 01069 referencetime=reftime) 01070 01071 sm.observe(sourcename="src1", spwname="band1", 01072 starttime=qa.quantity(0, "s"), 01073 stoptime=qa.quantity(etime)); 01074 01075 sm.setdata() 01076 sm.setvp() 01077 01078 eta_p, eta_s, eta_b, eta_t, eta_q, t_rx = self.noisetemp(telescope=telescope,freq=freq) 01079 eta_a = eta_p * eta_s * eta_b * eta_t 01080 if self.verbose: 01081 self.msg('antenna efficiency = ' + str(eta_a),origin="noise") 01082 self.msg('spillover efficiency = ' + str(eta_s),origin="noise") 01083 self.msg('correlator efficiency = ' + str(eta_q),origin="noise") 01084 01085 if pwv==None: 01086 # RI TODO choose based on freq octile 01087 pwv=2.0 01088 01089 # things hardcoded in ALMA etimecalculator, & defaults in simulator.xml 01090 t_ground=270. 01091 t_cmb=2.725 01092 # eta_q = 0.88 01093 # eta_a = 0.95*0.8*eta_s 01094 01095 if telescope=='ALMA' and (qa.convert(freq,"GHz")['value'])>600.: 01096 rxtype=1 # DSB 01097 01098 if method=="tsys-atm": 01099 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q, 01100 antefficiency=eta_a,trx=t_rx, 01101 tground=t_ground,tcmb=t_cmb,pwv=str(pwv)+"mm", 01102 mode="tsys-atm",table=tmpname,rxtype=rxtype) 01103 else: 01104 if method=="tsys-manual": 01105 if not t_sky: 01106 t_sky=200. 01107 self.msg("Warning: Sky brightness temp not set, using 200K",origin="noise",priority="warn") 01108 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q, 01109 antefficiency=eta_a,trx=t_rx,tatmos=t_sky, 01110 tground=t_ground,tcmb=t_cmb,tau=tau0, 01111 mode="tsys-manual",table=tmpname,rxtype=rxtype) 01112 else: 01113 self.msg("Unknown calculation method "+method,priority="error") 01114 return False 01115 01116 if doimnoise: 01117 sm.corrupt() 01118 01119 sm.done() 01120 01121 01122 if doimnoise: 01123 cellsize=qa.quantity(6.e3/250./qa.convert(model_start,"GHz")["value"],"arcsec") # need better cell determination - 250m?! 01124 cellsize=[cellsize,cellsize] 01125 # very light clean - its an empty image! 01126 self.imclean(tmpname+".ms",tmpname, 01127 "csclean",cellsize,[128,128], 01128 "J2000 00:00:00.00 "+qa.angle(dec)[0], 01129 100,"0.01mJy","natural",[],"I") 01130 01131 ia.open(tmpname+".image") 01132 stats= ia.statistics(robust=True, verbose=False,list=False) 01133 ia.done() 01134 imnoise=(stats["rms"][0]) 01135 else: 01136 imnoise=0. 01137 01138 nint = qa.convert(etime,'s')['value'] / qa.convert(integration,'s')['value'] 01139 nbase= 0.5*nant*(nant-1) 01140 01141 if os.path.exists(tmpname+".T.cal"): 01142 tb.open(tmpname+".T.cal") 01143 gain=tb.getcol("CPARAM") 01144 # RI TODO average instead of first? 01145 tb.done() 01146 # gain is per ANT so square for per baseline; 01147 # pick a gain from about the middle of the track 01148 noiseperbase=1./(gain[0][0][0.5*nint*nant].real)**2 01149 else: 01150 noiseperbase=0. 01151 01152 theoreticalnoise=noiseperbase/pl.sqrt(nint*nbase*2) # assume 2-poln 01153 01154 if debug==None: 01155 xx=glob.glob(tmpname+"*") 01156 for k in range(len(xx)): 01157 if os.path.isdir(xx[k]): 01158 cu.removetable(xx[k]) 01159 else: 01160 os.remove(xx[k]) 01161 01162 if doimnoise: 01163 return theoreticalnoise , imnoise 01164 else: 01165 return theoreticalnoise 01166 01167 01168 def setcfg(self, mysm, telescope, x, y, z, diam, 01169 padnames, posobs, mounttype=None): 01170 """ 01171 Sets the antenna positions for the mysm sm instance, which should have 01172 already opened an MS, given 01173 telescope - telescope name 01174 x - array of X positions, i.e. stnx from readantenna 01175 y - array of Y positions, i.e. stny from readantenna 01176 z - array of Z positions, i.e. stnz from readantenna 01177 diam - numpy.array of antenna diameters, i.e. from readantenna 01178 padnames - list of pad names 01179 posobs - The observatory position as a measure. 01180 01181 Optional: 01182 mounttype (Defaults to a guess based on telescope.) 01183 01184 Returns the mounttype that it uses. 01185 01186 Closes mysm and throws a ValueError if it can't set the configuration. 01187 """ 01188 if not mounttype: 01189 mounttype = 'alt-az' 01190 if telescope.upper() in ('ASKAP', # Effectively, if not BIZARRE. 01191 'DRAO', 01192 'WSRT'): 01193 mounttype = 'EQUATORIAL' 01194 elif telescope.upper() in ('LOFAR', 'LWA', 'MOST'): 01195 # And other long wavelength arrays...like that one in WA that 01196 # has 3 orthogonal feeds so they can go to the horizon. 01197 # 01198 # The SKA should go here too once people accept 01199 # that it will have to be correlated as stations. 01200 mounttype = 'BIZARRE' # Ideally it would be 'FIXED' or 01201 # 'GROUND'. 01202 01203 if not mysm.setconfig(telescopename=telescope, x=x, y=y, z=z, 01204 dishdiameter=diam.tolist(), 01205 mount=[mounttype], padname=padnames, 01206 coordsystem='global', referencelocation=posobs): 01207 mysm.close() 01208 raise ValueError("Error setting the configuration") 01209 return mounttype 01210 01211 01212 01213 01214 ########################################################### 01215 #===================== ephemeris ========================== 01216 01217 01218 def ephemeris(self, date, usehourangle=True, direction=None, telescope=None, ms=None): 01219 01220 if direction==None: direction=self.direction 01221 if telescope==None: telescope=self.telescopename 01222 01223 import pdb 01224 01225 # right now, simdata centers at the transit; when that changes, 01226 # or when that becomes optional, then that option needs to be 01227 # stored in the simutil object and used here, to either 01228 # center the plot at transit or not. 01229 # 01230 # direction="J2000 18h00m00.03s -45d59m59.6s" 01231 # refdate="2012/06/21/03:25:00" 01232 01233 # useHourAngle_p means simulate at transit 01234 # TODO: put in reftime parameter, parse 2012/05/21, 2012/05/21/transit, 01235 # and 2012/05/21/22:05:00 separately. 01236 01237 ds=self.direction_splitter(direction) # if list, returns average 01238 src=me.direction(ds[0],ds[1],ds[2]) 01239 01240 me.done() 01241 me.doframe(me.observatory(telescope)) 01242 01243 time=me.epoch('TAI',date) 01244 me.doframe(time) 01245 01246 # what is HA of source at refdate? 01247 offset_ha=qa.convert((me.measure(src,'hadec'))['m0'],'h') 01248 peak=me.epoch("TAI",qa.add(date,qa.mul(-1,offset_ha))) 01249 peaktime_float=peak['m0']['value'] 01250 if usehourangle: 01251 # offset the reftime to be at transit: 01252 time=peak 01253 me.doframe(time) 01254 01255 reftime_float=time['m0']['value'] 01256 reftime_floor=pl.floor(time['m0']['value']) 01257 refdate_str=qa.time(qa.totime(str(reftime_floor)+'d'),form='dmy')[0] 01258 01259 timeinc='15min' # for plotting 01260 timeinc=qa.convert(qa.time(timeinc)[0],'d')['value'] 01261 ntime=int(1./timeinc) 01262 01263 # check for circumpolar: 01264 rset = me.riseset(src) 01265 rise = rset['rise'] 01266 if rise == 'above': 01267 rise = time 01268 rise['m0']['value'] = rise['m0']['value'] - 0.5 01269 settime = time 01270 settime['m0']['value'] = settime['m0']['value'] + 0.5 01271 elif rise == 'below': 01272 raise ValueError(direction + ' is not visible from ' + telescope) 01273 else: 01274 settime = rset['set'] 01275 rise = me.measure(rise['utc'],'tai') 01276 settime = me.measure(settime['utc'],'tai') 01277 01278 # where to start plotting? 01279 offset=-0.5 01280 if settime < time: offset-=0.5 01281 if rise > time: offset+=0.5 01282 time['m0']['value']+=offset 01283 01284 times=[] 01285 az=[] 01286 el=[] 01287 01288 for i in range(ntime): 01289 times.append(time['m0']['value']) 01290 me.doframe(time) 01291 azel=me.measure(src,'azel') 01292 az.append(qa.convert(azel['m0'],'deg')['value']) 01293 el.append(qa.convert(azel['m1'],'deg')['value']) 01294 time['m0']['value']+=timeinc 01295 01296 # self.msg(" ref="+date,origin='ephemeris') 01297 self.msg("rise="+qa.time(rise['m0'],form='dmy')[0],origin='ephemeris') 01298 self.msg(" set="+qa.time(settime['m0'],form='dmy')[0],origin='ephemeris') 01299 01300 pl.plot((pl.array(times)-reftime_floor)*24,el) 01301 # peak=(rise['m0']['value']+settime['m0']['value'])/2 01302 # self.msg("peak="+qa.time('%fd' % peak,form='dmy'),origin='ephemeris') 01303 self.msg("peak="+qa.time('%fd' % reftime_float,form='dmy')[0],origin='ephemeris') 01304 01305 relpeak=peaktime_float-reftime_floor 01306 pl.plot(pl.array([1,1])*24*relpeak,[0,90]) 01307 01308 # if theres an ms, figure out the entire range of observation 01309 if ms: 01310 tb.open(ms+"/OBSERVATION") 01311 timerange=tb.getcol("TIME_RANGE") 01312 tb.done() 01313 obsstart=min(timerange.flat) 01314 obsend=max(timerange.flat) 01315 relstart=me.epoch("UTC",str(obsstart)+"s")['m0']['value']-reftime_floor 01316 relend=me.epoch("UTC",str(obsend)+"s")['m0']['value']-reftime_floor 01317 pl.plot([relstart*24,relend*24],[89,89],'r',linewidth=3) 01318 else: 01319 if self.totaltime>0: 01320 etimeh=qa.convert(self.totaltime,'h')['value'] 01321 pl.plot(pl.array([-0.5,0.5])*etimeh+relpeak*24,[80,80],'r') 01322 01323 pl.xlabel("hours relative to "+refdate_str,fontsize='x-small') 01324 pl.ylabel("elevation",fontsize='x-small') 01325 ax=pl.gca() 01326 l=ax.get_xticklabels() 01327 pl.setp(l,fontsize="x-small") 01328 l=ax.get_yticklabels() 01329 pl.setp(l,fontsize="x-small") 01330 01331 01332 pl.ylim([0,90]) 01333 pl.xlim(pl.array([-12,12])+24*(reftime_float-reftime_floor)) 01334 01335 01336 01337 01338 01339 01340 01341 01342 01343 01344 01345 01346 01347 01348 01349 01350 01351 01352 01353 ########################################################### 01354 #========================================================== 01355 01356 def readantenna(self, antab=None): 01357 ###Helper function to read 4 columns text antenna table X, Y, Z, Diam 01358 f=open(antab) 01359 self.msg("Reading antenna positions from '%s'" % antab,origin="readantenna") 01360 line= ' ' 01361 inx=[] 01362 iny=[] 01363 inz=[] 01364 ind=[] 01365 id=[] # pad id 01366 nant=0 01367 line=' ' 01368 params={} 01369 while (len(line)>0): 01370 try: 01371 line=f.readline() 01372 if line.startswith('#'): 01373 line=line[1:] 01374 paramlist=line.split('=') 01375 ### if the line is a parameter line like coordsys=utms, then it stores that 01376 if (paramlist.__len__() == 2): 01377 if params==None: 01378 params={paramlist[0].strip():paramlist[1].strip()} 01379 else: 01380 params[paramlist[0].strip()]=paramlist[1].strip() 01381 else: 01382 ### ignoring line that has less than 4 elements 01383 ### all coord systems should have x,y,z,diam,id, where xyz varies 01384 #print line.split() 01385 if(len(line.split()) >3): 01386 splitline=line.split() 01387 inx.append(float(splitline[0])) 01388 iny.append(float(splitline[1])) 01389 inz.append(float(splitline[2])) 01390 ind.append(float(splitline[3])) 01391 if len(splitline)>4: 01392 id.append(splitline[4]) 01393 else: 01394 id.append('A%02d'%nant) 01395 nant+=1 01396 except: 01397 break 01398 f.close() 01399 01400 if not params.has_key("observatory"): 01401 self.msg("Must specify observatory in antenna file",origin="readantenna",priority="error") 01402 return -1 01403 else: 01404 self.telescopename=params["observatory"] 01405 if self.verbose: 01406 self.msg("Using observatory= %s" % params["observatory"],origin="readantenna") 01407 01408 if not params.has_key("coordsys"): 01409 self.msg("Must specify coordinate system #coorsys=XYZ|UTM|GEO in antenna file",origin="readantenna",priority="error") 01410 return -1 01411 else: 01412 self.coordsys=params["coordsys"] 01413 01414 if (params["coordsys"].upper()=="XYZ"): 01415 ### earth-centered XYZ i.e. ITRF in casa 01416 stnx=inx 01417 stny=iny 01418 stnz=inz 01419 else: 01420 stnx=[] 01421 stny=[] 01422 stnz=[] 01423 if (params["coordsys"].upper()=="UTM"): 01424 ### expect easting, northing, elevation in m 01425 self.msg("Antenna locations in UTM; will read from file easting, northing, elevation in m",origin="readantenna") 01426 if params.has_key("zone"): 01427 zone=params["zone"] 01428 else: 01429 self.msg("You must specify zone=NN in your antenna file",origin="readantenna",priority="error") 01430 return -1 01431 if params.has_key("datum"): 01432 datum=params["datum"] 01433 else: 01434 self.msg("You must specify datum in your antenna file",origin="readantenna",priority="error") 01435 return -1 01436 if params.has_key("hemisphere"): 01437 nors=params["hemisphere"] 01438 nors=nors[0].upper() 01439 else: 01440 self.msg("You must specify hemisphere=N|S in your antenna file",origin="readantenna",priority="error") 01441 return -1 01442 01443 vsave=self.verbose 01444 for i in range(len(inx)): 01445 x,y,z = self.utm2xyz(inx[i],iny[i],inz[i],int(zone),datum,nors) 01446 if i==1: 01447 self.verbose=False 01448 stnx.append(x) 01449 stny.append(y) 01450 stnz.append(z) 01451 self.verbose=vsave 01452 else: 01453 if (params["coordsys"].upper()[0:3]=="LOC"): 01454 # I'm pretty sure Rob's function only works with lat,lon in degrees; 01455 meobs=me.observatory(params["observatory"]) 01456 if (meobs.__len__()<=1): 01457 self.msg("You need to add "+params["observatory"]+" to the Observatories table in your installation to proceed.",priority="error") 01458 return False,False,False,False,False,params["observatory"] 01459 obs=me.measure(meobs,'WGS84') 01460 obslat=qa.convert(obs['m1'],'deg')['value'] 01461 obslon=qa.convert(obs['m0'],'deg')['value'] 01462 obsalt=qa.convert(obs['m2'],'m')['value'] 01463 if self.verbose: 01464 self.msg("converting local tangent plane coordinates to ITRF using observatory position = %d %d " % (obslat,obslon)) 01465 #foo=self.getdatum(datum,verbose=True) 01466 for i in range(len(inx)): 01467 x,y,z = self.locxyz2itrf(obslat,obslon,inx[i],iny[i],inz[i]+obsalt) 01468 stnx.append(x) 01469 stny.append(y) 01470 stnz.append(z) 01471 else: 01472 if (params["coordsys"].upper()[0:3]=="GEO"): 01473 if params.has_key("datum"): 01474 datum=params["datum"] 01475 else: 01476 self.msg("You must specify zone=NN in your antenna file",origin="readantenna",priority="error") 01477 return -1 01478 if (datum.upper() != "WGS84"): 01479 self.msg("Unfortunately I only can deal with WGS84 right now",origin="readantenna",priority="error") 01480 return -1 01481 self.msg("geodetic coordinates not implemented yet",priority="error") 01482 01483 return (stnx, stny, stnz, pl.array(ind), id, nant, params["observatory"]) 01484 01485 01486 01487 01488 01489 01490 01491 01492 01493 01494 01495 01496 01497 01498 01499 01500 01501 01502 01503 ########################################################### 01504 #==================== geodesy ============================= 01505 01506 01507 def tmgeod(self,n,e,eps,cm,fe,sf,so,r,v0,v2,v4,v6,fn,er,esq): 01508 """ 01509 Transverse Mercator Projection 01510 conversion of grid coords n,e to geodetic coords 01511 revised subroutine of t. vincenty feb. 25, 1985 01512 converted from Fortran R Indebetouw Jan 2009 01513 ********** symbols and definitions *********************** 01514 latitude positive north, longitude positive west. 01515 all angles are in radian measure. 01516 01517 input: 01518 n,e are northing and easting coordinates respectively 01519 er is the semi-major axis of the ellipsoid 01520 esq is the square of the 1st eccentricity 01521 cm is the central meridian of the projection zone 01522 fe is the false easting value at the cm 01523 eps is the square of the 2nd eccentricity 01524 sf is the scale factor at the cm 01525 so is the meridional distance (times the sf) from the 01526 equator to southernmost parallel of lat. for the zone 01527 r is the radius of the rectifying sphere 01528 u0,u2,u4,u6,v0,v2,v4,v6 are precomputed constants for 01529 determination of meridianal dist. from latitude 01530 output: 01531 lat,lon are lat. and long. respectively 01532 conv is convergence 01533 kp point scale factor 01534 01535 the formula used in this subroutine gives geodetic accuracy 01536 within zones of 7 degrees in east-west extent. within state 01537 transverse mercator projection zones, several minor terms of 01538 the equations may be omitted (see a separate ngs publication). 01539 if programmed in full, the subroutine can be used for 01540 computations in surveys extending over two zones. 01541 """ 01542 01543 om=(n-fn+so)/(r*sf) # (northing - flag_north + distance_from_equator) 01544 cosom=pl.cos(om) 01545 foot=om+pl.sin(om)*cosom*(v0+v2*cosom*cosom+v4*cosom**4+v6*cosom**6) 01546 sinf=pl.sin(foot) 01547 cosf=pl.cos(foot) 01548 tn=sinf/cosf 01549 ts=tn*tn 01550 ets=eps*cosf*cosf 01551 rn=er*sf/pl.sqrt(1.-esq*sinf*sinf) 01552 q=(e-fe)/rn 01553 qs=q*q 01554 b2=-tn*(1.+ets)/2. 01555 b4=-(5.+3.*ts+ets*(1.-9.*ts)-4.*ets*ets)/12. 01556 b6=(61.+45.*ts*(2.+ts)+ets*(46.-252.*ts-60.*ts*ts))/360. 01557 b1=1. 01558 b3=-(1.+ts+ts+ets)/6. 01559 b5=(5.+ts*(28.+24.*ts)+ets*(6.+8.*ts))/120. 01560 b7=-(61.+662.*ts+1320.*ts*ts+720.*ts**3)/5040. 01561 lat=foot+b2*qs*(1.+qs*(b4+b6*qs)) 01562 l=b1*q*(1.+qs*(b3+qs*(b5+b7*qs))) 01563 lon=-l/cosf+cm 01564 01565 # compute scale factor 01566 fi=lat 01567 lam = lon 01568 sinfi=pl.sin(fi) 01569 cosfi=pl.cos(fi) 01570 l1=(lam-cm)*cosfi 01571 ls=l1*l1 01572 01573 tn=sinfi/cosfi 01574 ts=tn*tn 01575 01576 # convergence 01577 c1=-tn 01578 c3=(1.+3.*ets+2.*ets**2)/3. 01579 c5=(2.-ts)/15. 01580 conv=c1*l1*(1.+ls*(c3+c5*ls)) 01581 01582 # point scale factor 01583 f2=(1.+ets)/2. 01584 f4=(5.-4.*ts+ets*( 9.-24.*ts))/12. 01585 kp=sf*(1.+f2*ls*(1.+f4*ls)) 01586 01587 return lat,lon,conv,kp 01588 01589 01590 01591 01592 def tconpc(self,sf,orlim,er,esq,rf): 01593 01594 """ 01595 transverse mercator projection *** 01596 conversion of grid coords to geodetic coords 01597 revised subroutine of t. vincenty feb. 25, 1985 01598 converted from fortran r. indebetouw jan 2009 01599 ********** symbols and definitions *********************** 01600 input: 01601 rf is the reciprocal flattening of ellipsoid 01602 esq = e squared (eccentricity?) 01603 er is the semi-major axis for grs-80 01604 sf is the scale factor at the cm 01605 orlim is the southernmost parallel of latitude for which the 01606 northing coord is zero at the cm 01607 01608 output: 01609 eps 01610 so is the meridional distance (times the sf) from the 01611 equator to southernmost parallel of lat. for the zone 01612 r is the radius of the rectifying sphere 01613 u0,u2,u4,u6,v0,v2,v4,v6 are precomputed constants for 01614 determination of meridional dist. from latitude 01615 ********************************************************** 01616 """ 01617 01618 f=1./rf 01619 eps=esq/(1.-esq) 01620 pr=(1.-f)*er 01621 en=(er-pr)/(er+pr) 01622 en2=en*en 01623 en3=en*en*en 01624 en4=en2*en2 01625 01626 c2=-3.*en/2.+9.*en3/16. 01627 c4=15.*en2/16.-15.*en4/32. 01628 c6=-35.*en3/48. 01629 c8=315.*en4/512. 01630 u0=2.*(c2-2.*c4+3.*c6-4.*c8) 01631 u2=8.*(c4-4.*c6+10.*c8) 01632 u4=32.*(c6-6.*c8) 01633 u6=128.*c8 01634 01635 c2=3.*en/2.-27.*en3/32. 01636 c4=21.*en2/16.-55.*en4/32. 01637 c6=151.*en3/96. 01638 c8=1097.*en4/512. 01639 v0=2.*(c2-2.*c4+3.*c6-4.*c8) 01640 v2=8.*(c4-4.*c6+10.*c8) 01641 v4=32.*(c6-6.*c8) 01642 v6=128.*c8 01643 01644 r=er*(1.-en)*(1.-en*en)*(1.+2.25*en*en+(225./64.)*en4) 01645 cosor=pl.cos(orlim) 01646 omo=orlim+pl.sin(orlim)*cosor*(u0+u2*cosor*cosor+u4*cosor**4+u6*cosor**6) 01647 so=sf*r*omo 01648 01649 return eps,r,so,v0,v2,v4,v6 01650 01651 01652 01653 01654 01655 def getdatum(self,datumcode,verbose=False): 01656 """ 01657 local datums and ellipsoids; 01658 take local earth-centered xyz coordinates, add dx,dy,dz to get wgs84 earth-centered 01659 """ 01660 01661 # set equatorial radius and inverse flattening 01662 01663 ellipsoids={'AW':[6377563.396,299.3249647 ,'Airy 1830' ], 01664 'AM':[6377340.189,299.3249647 ,'Modified Airy' ], 01665 'AN':[6378160.0 ,298.25 ,'Australian National' ], 01666 'BR':[6377397.155,299.1528128 ,'Bessel 1841' ], 01667 'CC':[6378206.4 ,294.9786982 ,'Clarke 1866' ], 01668 'CD':[6378249.145,293.465 ,'Clarke 1880' ], 01669 'EA':[6377276.345,300.8017 ,'Everest (India 1830)' ], 01670 'RF':[6378137.0 ,298.257222101,'Geodetic Reference System 1980'], 01671 'HE':[6378200.0 ,298.30 ,'Helmert 1906' ], 01672 'HO':[6378270.0 ,297.00 ,'Hough 1960' ], 01673 'IN':[6378388.0 ,297.00 ,'International 1924' ], 01674 'SA':[6378160.0 ,298.25 ,'South American 1969' ], 01675 'WD':[6378135.0 ,298.26 ,'World Geodetic System 1972' ], 01676 'WE':[6378137.0 ,298.257223563,'World Geodetic System 1984' ]} 01677 01678 datums={ 01679 'AGD66' :[-133, -48, 148,'AN','Australian Geodetic Datum 1966'], 01680 'AGD84' :[-134, -48, 149,'AN','Australian Geodetic Datum 1984'], 01681 'ASTRO' :[-104,-129, 239,'IN','Camp Area Astro (Antarctica)' ], 01682 'CAPE' :[-136,-108,-292,'CD','CAPE (South Africa)' ], 01683 'ED50' :[ -87, -98,-121,'IN','European 1950' ], 01684 'ED79' :[ -86, -98,-119,'IN','European 1979' ], 01685 'GRB36' :[ 375,-111, 431,'AW','Great Britain 1936' ], 01686 'HAWAII':[ 89,-279,-183,'IN','Hawaiian Hawaii (Old)' ], 01687 'KAUAI' :[ 45,-290,-172,'IN','Hawaiian Kauai (Old)' ], 01688 'MAUI' :[ 65,-290,-190,'IN','Hawaiian Maui (Old)' ], 01689 'OAHU' :[ 56,-284,-181,'IN','Hawaiian Oahu (Old)' ], 01690 'INDIA' :[ 289, 734, 257,'EA','Indian' ], 01691 'NAD83' :[ 0, 0, 0,'RF','N. American 1983' ], 01692 'CANADA':[ -10, 158, 187,'CC','N. American Canada 1927' ], 01693 'ALASKA':[ -5, 135, 172,'CC','N. American Alaska 1927' ], 01694 'NAD27' :[ -8, 160, 176,'CC','N. American Conus 1927' ], 01695 'CARIBB':[ -7, 152, 178,'CC','N. American Caribbean' ], 01696 'MEXICO':[ -12, 130, 190,'CC','N. American Mexico' ], 01697 'CAMER' :[ 0, 125, 194,'CC','N. American Central America' ], 01698 'SAM56' :[-288, 175,-376,'IN','South American (Provisional1956)'], 01699 'SAM69' :[ -57, 1 , -41,'SA','South American 1969' ], 01700 'CAMPO' :[-148, 136, 90,'IN','S. American Campo Inchauspe (Argentina)'], 01701 'WGS72' :[ 0, 0 , 4.5,'WD','World Geodetic System - 72' ], 01702 'WGS84' :[ 0, 0 , 0,'WE','World Geodetic System - 84' ]} 01703 01704 if not datums.has_key(datumcode): 01705 self.msg("unknown datum %s" % datumcode,priority="error") 01706 return -1 01707 01708 datum=datums[datumcode] 01709 ellipsoid=datum[3] 01710 01711 if not ellipsoids.has_key(ellipsoid): 01712 self.msg("unknown ellipsoid %s" % ellipsoid,priority="error") 01713 return -1 01714 01715 if verbose: 01716 self.msg("Using %s datum with %s ellipsoid" % (datum[4],ellipsoids[ellipsoid][2]),origin="getdatum") 01717 return datum[0],datum[1],datum[2],ellipsoids[ellipsoid][0],ellipsoids[ellipsoid][1] 01718 01719 01720 01721 01722 01723 def utm2long(self,east,north,zone,datum,nors): 01724 """ 01725 this program converts universal transverse meractor coordinates to gps 01726 converted from fortran by r. indebetouw jan 2009. 01727 ri also added other datums and ellipsoids in a helper function 01728 01729 header from original UTMS fortran program: 01730 * this program was originally written by e. carlson 01731 * subroutines tmgrid, tconst, tmgeod, tconpc, 01732 * were written by t. vincenty, ngs, in july 1984 . 01733 * the orginal program was written in september of 1988. 01734 * 01735 * this program was updated on febuary 16, 1989. the update was 01736 * having the option of saving and *81* record file. 01737 * 01738 * 01739 * this program was updated on april 3, 1990. the following update 01740 * were made: 01741 * 1. change from just a choice of nad27 of nad83 reference 01742 * ellipsoids to; clarke 1866, grs80/wgs84, international, and 01743 * allow for user defined other. 01744 * 2. allow use of latitudes in southern hemisphere and longitudes 01745 * up to 360 degrees west. 01746 * 01747 * this program was updated on december 1, 1993. the following update 01748 * was made: 01749 * 1. the northings would compute the right latitude in the southern 01750 * hemishpere. 01751 * 2. the computed latitude on longidutes would be either in e or w. 01752 * 01753 ***************************************************************** * 01754 * disclaimer * 01755 * * 01756 * this program and supporting information is furnished by the * 01757 * government of the united states of america, and is accepted and * 01758 * used by the recipient with the understanding that the united states * 01759 * government makes no warranties, express or implied, concerning the * 01760 * accuracy, completeness, reliability, or suitability of this * 01761 * program, of its constituent parts, or of any supporting data. * 01762 * * 01763 * the government of the united states of america shall be under no * 01764 * liability whatsoever resulting from any use of this program. this * 01765 * program should not be relied upon as the sole basis for solving a * 01766 * problem whose incorrect solution could result in injury to person * 01767 * or property. * 01768 * * 01769 * this program is property of the government of the united states * 01770 * of america. therefore, the recipient further agrees not to assert * 01771 * proprietary rights therein and not to represent this program to * 01772 * anyone as being other than a government program. * 01773 * * 01774 *********************************************************************** 01775 01776 this is the driver program to compute latitudes and longitudes from 01777 the utms for each zone 01778 01779 input: 01780 northing, easting 01781 zone, datum 01782 nors=N/S 01783 01784 determined according to datum: 01785 er = equatorial radius of the ellipsoid (semi-major axis) 01786 rf = reciprocal of flatting of the ellipsod 01787 f = 01788 esq= e squared 01789 01790 calculated according to longitude and zone: 01791 rad = radian conversion factor 01792 cm = central meridian ( computed using the longitude) 01793 sf = scale factor of central meridian ( always .9996 for utm) 01794 orlim = southernmost parallel of latitude ( always zero for utm) 01795 r, a, b, c, u, v, w = ellipsoid constants used for computing 01796 meridional distance from latitude 01797 so = meridional distance (multiplied by scale factor ) 01798 from the equator to the southernmost parallel of latitude 01799 ( always zero for utm) 01800 01801 """ 01802 01803 rad=180./pl.pi 01804 01805 offx,offy,offz,er,rf = self.getdatum(datum,verbose=self.verbose) 01806 01807 f=1./rf 01808 esq=(2*f-f*f) 01809 01810 # find the central meridian if the zone number is less than 30 01811 01812 if zone < 30 : # ie W long - this code uses W=positive lon 01813 iz=zone 01814 icm=(183-(6*iz)) 01815 cm=float(icm)/rad 01816 ucm=(icm+3)/rad 01817 lcm=(icm-3)/rad 01818 else: 01819 iz=zone 01820 icm=(543 - (6*iz)) 01821 cm= float(icm)/rad 01822 ucm=(icm+3)/rad 01823 lcm=(icm-3)/rad 01824 01825 tol=(5./60.)/rad 01826 01827 if nors == 'S': 01828 fn= 10000000. 01829 else: 01830 fn=0. 01831 01832 fe=500000.0 01833 sf=0.9996 01834 orlim=0.0 01835 01836 found=0 01837 01838 # precompute parameters for this zone: 01839 eps,r,so,v0,v2,v4,v6 = self.tconpc(sf,orlim,er,esq,rf) 01840 01841 # compute the latitudes and longitudes: 01842 lat,lon,conv,kp = self.tmgeod(north,east,eps,cm,fe,sf,so,r,v0,v2,v4,v6,fn,er,esq) 01843 01844 # do the test to see if the longitude is within 5 minutes 01845 # of the boundaries for the zone and if so compute the 01846 # north and easting for the adjacent zone 01847 01848 # if abs(ucm-lam) < tol: 01849 # cm=float(icm+6)/rad 01850 # iz=iz-1 01851 # if iz == 0: 01852 # iz=60 01853 # found=found+1 01854 # lat,lon,conv,kp = tmgeod(n,e,eps,cm,fe,sf,so,r,v0,v2,v4,v6,fn,er,esq) 01855 # 01856 # if abs(lcm-lam) < tol: 01857 # cm=float(icm-6)/rad 01858 # iz=iz+1 01859 # if iz == 61: 01860 # iz=1 01861 # lat,lon,conv,kp = tmgeod(n,e,eps,cm,fe,sf,so,r,v0,v2,v4,v6,fn,er,esq) 01862 # found=found+1 01863 01864 # *** convert to more usual convention of negative lon = W 01865 lon=-lon 01866 01867 if self.verbose: 01868 self.msg(" longitude, latitude = %s %s; conv,kp = %f,%f" % (qa.angle(qa.quantity(lon,"rad"),prec=8)[0],qa.angle(qa.quantity(lat,"rad"),prec=8)[0],conv,kp),origin="utm2long") 01869 01870 return lon,lat 01871 01872 01873 01874 01875 def long2xyz(self,long,lat,elevation,datum): 01876 01877 dx,dy,dz,er,rf = self.getdatum(datum,verbose=False) 01878 01879 f=1./rf 01880 esq=2*f-f**2 01881 nu=er/(1.-esq*(pl.sin(lat))**2) 01882 01883 x=(nu+elevation)*pl.cos(lat)*pl.cos(long) +dx 01884 y=(nu+elevation)*pl.cos(lat)*pl.sin(long) +dy 01885 z=((1.-esq)*nu+elevation)*pl.sin(lat) +dz 01886 01887 return x,y,z 01888 01889 01890 def xyz2long(self,x,y,z,datum): 01891 01892 dx,dy,dz,er,rf = self.getdatum(datum,verbose=False) 01893 01894 f=1./rf 01895 01896 b= ((x-dx)**2 + (y-dy)**2) / er**2 01897 c= (z-dx)**2 / er**2 01898 01899 esq=2*f-f**2 # (a2-b2)/a2 01900 01901 a0=c*(1-esq) 01902 a1=2*a0 01903 efth=esq**2 01904 a2=a0+b-efth 01905 a3=-2.*efth 01906 a4=-efth 01907 01908 b0=4.*a0 01909 b1=3.*a1 01910 b2=2.*a2 01911 b3=a3 01912 01913 # refine/calculate esq 01914 nlqk=esq 01915 for i in range(3): 01916 nlqks = nlqk * nlqk 01917 nlqkc = nlqk * nlqks 01918 nlf = a0*nlqks*nlqks + a1*nlqkc + a2*nlqks + a3*nlqk + a4 01919 nlfprm = b0*nlqkc + b1*nlqks + b2*nlqk + b3 01920 nlqk = nlqk - (nlf / nlfprm) 01921 01922 y0 = (1.+nlqk)*(z-dz) 01923 x0 = pl.sqrt((x-dx)**2 + (y-dy)**2) 01924 lat=pl.arctan2(y0,x0) 01925 lon=pl.arctan2(y-dy,x-dx) 01926 #print x-dx,y-dy,z-dz,x0,y0 01927 01928 return lon,lat 01929 01930 01931 def utm2xyz(self,easting,northing,elevation,zone,datum,nors): 01932 01933 lon,lat = self.utm2long(easting,northing,zone,datum,nors) 01934 x,y,z = self.long2xyz(lon,lat,elevation,datum) 01935 01936 return x,y,z 01937 01938 01939 01940 def locxyz2itrf(self, lat, longitude, locx=0.0, locy=0.0, locz=0.0): 01941 """ 01942 Returns the nominal ITRF (X, Y, Z) coordinates (m) for a point at "local" 01943 (x, y, z) (m) measured at geodetic latitude lat and longitude longitude 01944 (degrees). The ITRF frame used is not the official ITRF, just a right 01945 handed Cartesian system with X going through 0 latitude and 0 longitude, 01946 and Z going through the north pole. The "local" (x, y, z) are measured 01947 relative to the closest point to (lat, longitude) on the WGS84 reference 01948 ellipsoid, with z normal to the ellipsoid and y pointing north. 01949 """ 01950 # from Rob Reid; need to generalize to use any datum... 01951 phi, lmbda = map(pl.radians, (lat, longitude)) 01952 sphi = pl.sin(phi) 01953 a = 6378137.0 # WGS84 equatorial semimajor axis 01954 b = 6356752.3142 # WGS84 polar semimajor axis 01955 ae = pl.arccos(b / a) 01956 N = a / pl.sqrt(1.0 - (pl.sin(ae) * sphi)**2) 01957 01958 # Now you see the connection between the Old Ones and Antarctica... 01959 Nploczcphimlocysphi = (N + locz) * pl.cos(phi) - locy * sphi 01960 01961 clmb = pl.cos(lmbda) 01962 slmb = pl.sin(lmbda) 01963 01964 x = Nploczcphimlocysphi * clmb - locx * slmb 01965 y = Nploczcphimlocysphi * slmb + locx * clmb 01966 z = (N * (b / a)**2 + locz) * sphi + locy * pl.cos(phi) 01967 01968 return x, y, z 01969 01970 01971 01972 01973 def itrf2loc(self, x,y,z, cx,cy,cz): 01974 """ 01975 itrf xyz and COFA cx,cy,cz -> latlon WGS84 01976 """ 01977 clon,clat = self.xyz2long(cx,cy,cz,'WGS84') 01978 ccoslon=pl.cos(clon) 01979 csinlon=pl.sin(clon) 01980 csinlat=pl.sin(clat) 01981 n=x.__len__() 01982 lat=pl.zeros(n) 01983 lon=pl.zeros(n) 01984 01985 # do like MsPlotConvert 01986 for i in range(n): 01987 # translate w/o rotating: 01988 xtrans=x[i]-cx 01989 ytrans=y[i]-cy 01990 ztrans=z[i]-cz 01991 # rotate 01992 lat[i] = (-csinlon*xtrans) + (ccoslon*ytrans) 01993 lon[i] = (-csinlat*ccoslon*xtrans) - (csinlat*csinlon*ytrans) + ztrans 01994 01995 return lat,lon 01996 01997 01998 01999 02000 02001 02002 02003 02004 02005 02006 02007 02008 02009 02010 ########################################################### 02011 02012 def plotants(self,x,y,z,d,name): 02013 # given globals 02014 02015 #stnx, stny, stnz, stnd, nant, telescopename = util.readantenna(antennalist) 02016 cx=pl.mean(x) 02017 cy=pl.mean(y) 02018 cz=pl.mean(z) 02019 lat,lon = self.itrf2loc(x,y,z,cx,cy,cz) 02020 n=lat.__len__() 02021 02022 dolam=0 02023 # TODO convert to klam: (d too) 02024 ### 02025 02026 rg=max(lat)-min(lat) 02027 r2=max(lon)-min(lon) 02028 if r2>rg: 02029 rg=r2 02030 if max(d)>0.01*rg: 02031 pl.plot(lat,lon,',') 02032 #print max(d),ra 02033 for i in range(n): 02034 pl.gca().add_patch(pl.Circle((lat[i],lon[i]),radius=0.5*d[i],fc="#dddd66")) 02035 if n<10: 02036 pl.text(lat[i],lon[i],name[i],horizontalalignment='center',verticalalignment='center') 02037 else: 02038 pl.plot(lat,lon,'o',c="#dddd66") 02039 if n<10: 02040 for i in range(n): 02041 pl.text(lat[i],lon[i],name[i],horizontalalignment='center',fontsize=8) 02042 02043 pl.axis('equal') 02044 #if dolam: 02045 # pl.xlabel("kilolamda") 02046 # pl.ylabel("kilolamda") 02047 02048 02049 02050 02051 02052 02053 02054 02055 02056 02057 02058 02059 02060 02061 02062 02063 02064 02065 02066 ###################################################### 02067 # helper function to get the pixel size from an image (positive increments) 02068 02069 def cellsize(self,image): 02070 ia.open(image) 02071 mycs=ia.coordsys() 02072 ia.done() 02073 increments=mycs.increment(type="direction")['numeric'] 02074 cellx=qa.quantity(abs(increments[0]),mycs.units(type="direction")[0]) 02075 celly=qa.quantity(abs(increments[1]),mycs.units(type="direction")[1]) 02076 xform=mycs.lineartransform(type="direction") 02077 offdiag=max(abs(xform[0,1]),abs(xform[1,0])) 02078 if offdiag > 1e-4: 02079 self.msg("Your image is rotated with respect to Lat/Lon. I can't cope with that yet",priority="error") 02080 cellx=qa.mul(cellx,abs(xform[0,0])) 02081 celly=qa.mul(celly,abs(xform[1,1])) 02082 return [qa.tos(cellx),qa.tos(celly)] 02083 02084 ###################################################### 02085 # helper function to get the spectral size from an image 02086 02087 def spectral(self,image): 02088 ia.open(image) 02089 cs=ia.coordsys() 02090 sh=ia.shape() 02091 ia.done() 02092 spc=cs.findcoordinate("spectral") 02093 if not spc[0]: return (0,0) 02094 model_width=str(cs.increment(type="spectral")['numeric'][0])+cs.units(type="spectral")[0] 02095 model_nchan=sh[spc['pixel']] 02096 return model_nchan,model_width 02097 02098 ###################################################### 02099 02100 def is4d(self, image): 02101 ia.open(image) 02102 s=ia.shape() 02103 if len(s)!=4: return False 02104 cs=ia.coordsys() 02105 ia.done() 02106 dir=cs.findcoordinate("direction") 02107 spc=cs.findcoordinate("spectral") 02108 stk=cs.findcoordinate("stokes") 02109 if not (dir[0] and spc[0] and stk[0]): return False 02110 if dir[1].__len__() != 2: return False 02111 if spc[1].__len__() != 1: return False 02112 if stk[1].__len__() != 1: return False 02113 # they have to be in the correct order too 02114 if stk[1]!=2: return False 02115 if spc[1]!=3: return False 02116 if dir[1][0]!=0: return False 02117 if dir[1][1]!=1: return False 02118 cs.done() 02119 return True 02120 02121 ################################################################## 02122 # fit modelimage into a 4 coordinate image defined by the parameters 02123 02124 # TODO spectral extrapolation and regridding using innchan **** 02125 02126 def modifymodel(self, inimage, outimage, 02127 inbright,indirection,incell,incenter,inwidth,innchan, 02128 flatimage=False): # if nonzero, create mom -1 image 02129 02130 # new ia tool 02131 in_ia=ia.newimagefromfile(inimage) 02132 in_shape=in_ia.shape() 02133 in_csys=in_ia.coordsys() 02134 02135 # pull data first, since ia.stats doesn't work w/o a CS: 02136 if outimage!=inimage: 02137 if self.verbose: self.msg("rearranging input data (may take some time for large cubes)") 02138 arr=in_ia.getchunk() 02139 else: 02140 # TODO move rearrange to inside ia tool, and at least don't do this: 02141 arr=pl.zeros(in_shape) 02142 axmap=[-1,-1,-1,-1] 02143 axassigned=[-1,-1,-1,-1] 02144 02145 02146 02147 # brightness scaling 02148 if (inbright=="unchanged") or (inbright==""): 02149 scalefactor=1. 02150 else: 02151 if self.isquantity(inbright,halt=False): 02152 qinb=qa.quantity(inbright) 02153 if qinb['unit']!='': 02154 # qa doesn't deal universally well with pixels and beams 02155 # so this may fail: 02156 try: 02157 inb=qa.convert(qinb,"Jy/pixel")['value'] 02158 except: 02159 inb=qinb['value'] 02160 self.msg("assuming inbright="+str(inbright)+" means "+str(inb)+" Jy/pixel",priority="warn") 02161 inbright=inb 02162 scalefactor=float(inbright)/pl.nanmax(arr) 02163 02164 # check shape characteristics of the input; 02165 # add degenerate axes as neeed: 02166 02167 in_dir=in_csys.findcoordinate("direction") 02168 in_spc=in_csys.findcoordinate("spectral") 02169 in_stk=in_csys.findcoordinate("stokes") 02170 02171 02172 # first check number of pixel axes and raise to 4 if required 02173 in_nax=arr.shape.__len__() 02174 if in_nax<2: 02175 self.msg("Your input model has fewer than 2 dimensions. Can't proceed",priority="error") 02176 return False 02177 if in_nax==2: 02178 arr=arr.reshape([arr.shape[0],arr.shape[1],1]) 02179 in_shape=arr.shape 02180 in_nax=in_shape.__len__() # which should be 3 02181 if in_nax==3: 02182 arr=arr.reshape([arr.shape[0],arr.shape[1],arr.shape[2],1]) 02183 in_shape=arr.shape 02184 in_nax=in_shape.__len__() # which should be 4 02185 if in_nax>4: 02186 self.msg("model image has more than 4 dimensions. Not sure how to proceed",priority="error") 02187 return False 02188 02189 02190 # make incell a list if not already 02191 if type(incell) == type([]): 02192 incell = map(qa.convert,incell,['arcsec','arcsec']) 02193 else: 02194 incell = qa.abs(qa.convert(incell,'arcsec')) 02195 # incell[0]<0 for RA increasing left 02196 incell = [qa.mul(incell,-1),incell] 02197 # later, we can test validity with qa.compare() 02198 02199 02200 # now parse coordsys: 02201 model_refdir="" 02202 model_cell="" 02203 # look for direction coordinate, with two pixel axes: 02204 if in_dir[0]: 02205 in_ndir = in_dir[1].__len__() 02206 if in_ndir != 2: 02207 self.msg("Mal-formed direction coordinates in modelimage. Discarding and using first two pixel axes for RA and Dec.") 02208 axmap[0]=0 # direction in first two pixel axes 02209 axmap[1]=1 02210 axassigned[0]=0 # coordinate corresponding to first 2 pixel axes 02211 axassigned[1]=0 02212 else: 02213 # we've found direction axes, and may change their increments and direction or not. 02214 dirax=in_dir[1] 02215 axmap[0]=dirax[0] 02216 axmap[1]=dirax[1] 02217 axassigned[dirax[0]]=0 02218 axassigned[dirax[1]]=0 02219 if self.verbose: self.msg("Direction coordinate (%i,%i) parsed" % (axmap[0],axmap[1]),origin="setup model") 02220 # move model_refdir to center of image 02221 model_refpix=[0.5*in_shape[axmap[0]],0.5*in_shape[axmap[1]]] 02222 ra = in_ia.toworld(model_refpix,'q')['quantity']['*'+str(axmap[0]+1)] 02223 dec = in_ia.toworld(model_refpix,'q')['quantity']['*'+str(axmap[1]+1)] 02224 model_refdir= in_csys.referencecode(type="direction",list=False)[0]+" "+qa.formxxx(ra,format='hms',prec=5)+" "+qa.formxxx(dec,format='dms',prec=5) 02225 model_projection=in_csys.projection()['type'] 02226 model_projpars=in_csys.projection()['parameters'] 02227 02228 # cell size from image 02229 increments=in_csys.increment(type="direction")['numeric'] 02230 incellx=qa.quantity(increments[0],in_csys.units(type="direction")[0]) 02231 incelly=qa.quantity(increments[1],in_csys.units(type="direction")[1]) 02232 xform=in_csys.lineartransform(type="direction") 02233 offdiag=max(abs(xform[0,1]),abs(xform[1,0])) 02234 if offdiag > 1e-4: 02235 self.msg("Your image is rotated with respect to Lat/Lon. I can't cope with that yet",priority="error") 02236 incellx=qa.mul(incellx,xform[0,0]) 02237 incelly=qa.mul(incelly,xform[1,1]) 02238 02239 # preserve sign in case input image is backwards (RA increasing right) 02240 model_cell = [qa.convert(incellx,'arcsec'),qa.convert(incelly,'arcsec')] 02241 02242 else: # no valid direction coordinate 02243 axmap[0]=0 # assign direction to first two pixel axes 02244 axmap[1]=1 02245 axassigned[0]=0 # assign coordinate corresponding to first 2 pixel axes 02246 axassigned[1]=0 02247 02248 02249 # try to parse direction using splitter function and override model_refdir 02250 if type(indirection)==type([]): 02251 if len(indirection) > 1: 02252 self.msg("error parsing indirection "+str(indirection)+" -- should be a single direction string") 02253 return False 02254 else: 02255 indirection=indirection[0] 02256 if self.isdirection(indirection,halt=False): 02257 # lon/lat = ra/dec if J2000, =glon/glat if galactic 02258 epoch, lon, lat = self.direction_splitter(indirection) 02259 02260 model_refdir=epoch+qa.formxxx(lon,format='hms',prec=5)+" "+qa.formxxx(lat,format='dms',prec=5) 02261 model_refpix=[0.5*in_shape[axmap[0]],0.5*in_shape[axmap[1]]] 02262 model_projection="SIN" # for indirection we default to SIN. 02263 model_projpars=pl.array([0.,0]) 02264 if self.verbose: self.msg("setting model image direction to indirection = "+model_refdir) 02265 else: 02266 # indirection is not set - is there a direction in the model already? 02267 if not self.isdirection(model_refdir,halt=False): 02268 self.msg("Cannot determine reference direction in model image. Valid 'indirection' parameter must be provided.",priority="error") 02269 return False 02270 02271 02272 # override model_cell? 02273 cell_replaced=False 02274 if self.isquantity(incell[0],halt=False): 02275 if qa.compare(incell[0],"1arcsec"): 02276 model_cell=incell 02277 cell_replaced=True 02278 if self.verbose: self.msg("replacing existing model cell size with incell") 02279 valid_modcell=False 02280 if not cell_replaced: 02281 if self.isquantity(model_cell[0],halt=False): 02282 if qa.compare(model_cell[0],"1arcsec"): 02283 valid_modcell=True 02284 if not valid_modcell: 02285 self.msg("Unable to determine model cell size. Valid 'incell' parameter must be provided.",priority="error") 02286 return False 02287 02288 02289 02290 if self.verbose: 02291 self.msg("model image shape="+str(in_shape),origin="setup model") 02292 self.msg("model pixel = %8.2e x %8.2e arcsec" % (model_cell[0]['value'],model_cell[1]['value']),origin="setup model") 02293 02294 02295 02296 02297 02298 02299 02300 # we've now found or assigned two direction axes, and changed direction and cell if required 02301 # next, work on spectral axis: 02302 02303 model_center="" 02304 model_width="" 02305 # look for a spectral axis: 02306 if in_spc[0]: 02307 #if type(in_spc[1]) == type(1) : 02308 # # should not come here after SWIG switch over 02309 # foo=in_spc[1] 02310 #else: 02311 foo=in_spc[1][0] 02312 if in_spc[1].__len__() > 1: 02313 self.msg("you seem to have two spectral axes",priority="warn") 02314 model_nchan=arr.shape[foo] 02315 axmap[3]=foo 02316 axassigned[foo]=3 02317 model_restfreq=in_csys.restfrequency() 02318 in_startpix=in_csys.referencepixel(type="spectral")['numeric'][0] 02319 model_width=in_csys.increment(type="spectral")['numeric'][0] 02320 model_start=in_csys.referencevalue(type="spectral")['numeric'][0]-in_startpix*model_width 02321 # this maybe can be done more accurately - for nonregular 02322 # grids it may trip things up 02323 # start is center of first channel. for nch=1, that equals center 02324 model_center=model_start+0.5*(model_nchan-1)*model_width 02325 model_width=str(model_width)+in_csys.units(type="spectral")[0] 02326 model_start=str(model_start)+in_csys.units(type="spectral")[0] 02327 model_center=str(model_center)+in_csys.units(type="spectral")[0] 02328 add_spectral_coord=False 02329 if self.verbose: self.msg("Spectral Coordinate %i parsed" % axmap[3],origin="setup model") 02330 else: 02331 # need to add one to the coordsys 02332 add_spectral_coord=True 02333 02334 02335 # override incenter? 02336 center_replaced=False 02337 if self.isquantity(incenter,halt=False): 02338 if qa.compare(incenter,"1Hz"): 02339 if (qa.quantity(incenter))['value']>=0: 02340 model_center=incenter 02341 model_restfreq=model_center 02342 center_replaced=True 02343 if self.verbose: self.msg("setting central frequency to "+incenter) 02344 valid_modcenter=False 02345 if not center_replaced: 02346 if self.isquantity(model_center,halt=False): 02347 if qa.compare(model_center,"1Hz"): 02348 valid_modcenter=True 02349 if not valid_modcenter: 02350 self.msg("Unable to determine model frequency. Valid 'incenter' parameter must be provided.",priority="error") 02351 return False 02352 02353 # override inwidth? 02354 width_replaced=False 02355 if self.isquantity(inwidth,halt=False): 02356 if qa.compare(inwidth,"1Hz"): 02357 if (qa.quantity(inwidth))['value']>=0: 02358 model_width=inwidth 02359 width_replaced=True 02360 if self.verbose: self.msg("setting channel width to "+inwidth) 02361 valid_modwidth=False 02362 if not width_replaced: 02363 if self.isquantity(model_width,halt=False): 02364 if qa.compare(model_width,"1Hz"): 02365 valid_modwidth=True 02366 if not valid_modwidth: 02367 self.msg("Unable to determine model channel or bandwidth. Valid 'inwidth' parameter must be provided.",priority="error") 02368 return False 02369 02370 02371 02372 02373 02374 model_stokes="" 02375 # look for a stokes axis 02376 if in_stk[0]: 02377 model_stokes=in_csys.stokes() 02378 foo=model_stokes[0] 02379 out_nstk=model_stokes.__len__() 02380 for i in range(out_nstk-1): 02381 foo=foo+model_stokes[i+1] 02382 model_stokes=foo 02383 #if type(in_stk[1]) == type(1): 02384 # # should not come here after SWIG switch over 02385 # foo=in_stk[1] 02386 #else: 02387 foo=in_stk[1][0] 02388 if in_stk[1].__len__() > 1: 02389 self.msg("you seem to have two stokes axes",priority="warn") 02390 axmap[2]=foo 02391 axassigned[foo]=2 02392 if in_shape[foo]>4: 02393 self.msg("you appear to have more than 4 Stokes components - please edit your header and/or parameters",priority="error") 02394 return False 02395 add_stokes_coord=False 02396 if self.verbose: self.msg("Stokes Coordinate %i parsed" % axmap[2],origin="setup model") 02397 else: 02398 # need to add one to the coordsys 02399 add_stokes_coord=True 02400 02401 02402 02403 if add_spectral_coord: 02404 # find first unused axis - probably at end, but just in case its not: 02405 i=0 02406 extra_axis=-1 02407 while extra_axis<0 and i<4: 02408 if axassigned[i]<0: extra_axis=i 02409 i+=1 02410 if extra_axis<0: 02411 self.msg("I can't find an unused axis to make Spectral [%i %i %i %i] " % (axassigned[0],axassigned[1],axassigned[2],axassigned[3]),priority="error",origin="setup model") 02412 return False 02413 02414 axmap[3]=extra_axis 02415 axassigned[extra_axis]=3 02416 model_nchan=arr.shape[extra_axis] 02417 02418 02419 if add_stokes_coord: 02420 # find first unused axis - probably at end, but just in case its not: 02421 i=0 02422 extra_axis=-1 02423 while extra_axis<0 and i<4: 02424 if axassigned[i]<0: extra_axis=i 02425 i+=1 02426 if extra_axis<0: 02427 self.msg("I can't find an unused axis to make Stokes [%i %i %i %i] " % (axassigned[0],axassigned[1],axassigned[2],axassigned[3]),priority="error",origin="setup model") 02428 return False 02429 axmap[2]=extra_axis 02430 axassigned[extra_axis]=2 02431 02432 if arr.shape[extra_axis]>4: 02433 self.msg("you have %i Stokes parameters in your potential Stokes axis %i. something is wrong." % (arr.shape[extra_axis],extra_axis),priority="error") 02434 return False 02435 if self.verbose: self.msg("Adding Stokes Coordinate",origin="setup model") 02436 if arr.shape[extra_axis]==4: 02437 model_stokes="IQUV" 02438 if arr.shape[extra_axis]==3: 02439 model_stokes="IQV" 02440 self.msg("setting IQV Stokes parameters from the 4th axis of you model. If that's not what you want, then edit the header",origin="setup model",priority="warn") 02441 if arr.shape[extra_axis]==2: 02442 model_stokes="IQ" 02443 self.msg("setting IQ Stokes parameters from the 4th axis of you model. If that's not what you want, then edit the header",origin="setup model",priority="warn") 02444 if arr.shape[extra_axis]<=1: 02445 model_stokes="I" 02446 02447 02448 02449 02450 # ======================================== 02451 # now we should have 4 assigned pixel axes, and model_cell, model_refdir, model_nchan, 02452 # model_stokes all set 02453 out_nstk=len(model_stokes) 02454 if self.verbose: 02455 self.msg("axis map for model image = %i %i %i %i" % 02456 (axmap[0],axmap[1],axmap[2],axmap[3]),origin="setup model") 02457 02458 modelshape=[in_shape[axmap[0]], in_shape[axmap[1]],out_nstk,model_nchan] 02459 if outimage!=inimage: 02460 ia.fromshape(outimage,modelshape,overwrite=True) 02461 modelcsys=ia.coordsys() 02462 else: 02463 modelcsys=in_ia.coordsys() 02464 modelcsys.setunits(['rad','rad','','Hz']) 02465 modelcsys.setincrement([qa.convert(model_cell[0],modelcsys.units()[0])['value'], # already -1 02466 qa.convert(model_cell[1],modelcsys.units()[1])['value']], 02467 type="direction") 02468 02469 dirm=self.dir_s2m(model_refdir) 02470 lonq=dirm['m0'] 02471 latq=dirm['m1'] 02472 modelcsys.setreferencecode(dirm['refer'],type="direction") 02473 if len(model_projpars)>0: 02474 modelcsys.setprojection(parameters=model_projpars.tolist(),type=model_projection) 02475 else: 02476 modelcsys.setprojection(type=model_projection) 02477 modelcsys.setreferencevalue( 02478 [qa.convert(lonq,modelcsys.units()[0])['value'], 02479 qa.convert(latq,modelcsys.units()[1])['value']], 02480 type="direction") 02481 modelcsys.setreferencepixel(model_refpix,"direction") 02482 if self.verbose: 02483 self.msg("sky model image direction = "+model_refdir) 02484 self.msg("sky model image increment = "+str(model_cell)) 02485 02486 modelcsys.setspectral(refcode="LSRK",restfreq=model_restfreq) 02487 modelcsys.setreferencevalue(qa.convert(model_center,modelcsys.units()[3])['value'],type="spectral") 02488 # modelcsys.setreferencepixel(0.5*model_nchan,type="spectral") # default is middle chan 02489 modelcsys.setreferencepixel(0.5*(model_nchan-1),type="spectral") # but not half-pixel 02490 modelcsys.setincrement(qa.convert(model_width,modelcsys.units()[3])['value'],type="spectral") 02491 02492 02493 # first assure that the csys has the expected order 02494 expected=['Direction', 'Direction', 'Stokes', 'Spectral'] 02495 if modelcsys.axiscoordinatetypes() != expected: 02496 self.msg("internal error with coordinate axis order created by Imager",priority="error") 02497 self.msg(modelcsys.axiscoordinatetypes().__str__(),priority="error") 02498 return False 02499 02500 # more checks: 02501 foo=pl.array(modelshape) 02502 if not (pl.array(arr.shape) == pl.array(foo.take(axmap).tolist())).all(): 02503 self.msg("internal error: I'm confused about the shape if your model data cube",priority="error") 02504 self.msg("have "+foo.take(axmap).__str__()+", want "+in_shape.__str__(),priority="error") 02505 return False 02506 02507 if outimage!=inimage: 02508 ia.setcoordsys(modelcsys.torecord()) 02509 ia.done() 02510 ia.open(outimage) 02511 02512 # now rearrange the pixel axes if ness. 02513 for ax in range(4): 02514 if axmap[ax] != ax: 02515 if self.verbose: self.msg("swapping input axes %i with %i" % (ax,axmap[ax]),origin="setup model") 02516 arr=arr.swapaxes(ax,axmap[ax]) 02517 tmp=axmap[ax] 02518 axmap[ax]=ax 02519 axmap[tmp]=tmp 02520 02521 02522 # there's got to be a better way to remove NaNs: 02523 if outimage!=inimage: 02524 for i0 in range(arr.shape[0]): 02525 for i1 in range(arr.shape[1]): 02526 for i2 in range(arr.shape[2]): 02527 for i3 in range(arr.shape[3]): 02528 foo=arr[i0,i1,i2,i3] 02529 if foo!=foo: arr[i0,i1,i2,i3]=0.0 02530 02531 if self.verbose and outimage!=inimage: 02532 self.msg("model array minmax= %e %e" % (arr.min(),arr.max()),origin="setup model") 02533 self.msg("scaling model brightness by a factor of %f" % scalefactor,origin="setup model") 02534 self.msg("image channel width = %8.2e GHz" % qa.convert(model_width,'GHz')['value'],origin="setup model") 02535 if arr.nbytes > 5e7: 02536 self.msg("your model is large - predicting visibilities may take a while.",priority="warn") 02537 02538 if outimage!=inimage: 02539 ia.putchunk(arr*scalefactor) 02540 ia.setbrightnessunit("Jy/pixel") 02541 ia.close() 02542 in_ia.close() 02543 # outimage should now have correct Coordsys and shape 02544 02545 # make a moment 0 image 02546 if flatimage and outimage!=inimage: 02547 self.flatimage(outimage,verbose=self.verbose) 02548 02549 # returning to the outside world we'll return a positive cell: 02550 model_cell=[qa.abs(model_cell[0]),qa.abs(model_cell[1])] 02551 model_size=[qa.mul(modelshape[0],model_cell[0]), 02552 qa.mul(modelshape[1],model_cell[1])] 02553 02554 02555 return model_refdir,model_cell,model_size,model_nchan,model_center,model_width,model_stokes 02556 02557 02558 02559 02560 02561 ################################################################## 02562 # image/clean subtask 02563 02564 def imclean(self,mstoimage,imagename, 02565 cleanmode,cell,imsize,imcenter,niter,threshold,weighting, 02566 outertaper,stokes,sourcefieldlist="",modelimage="",mask=[]): 02567 from clean import clean 02568 02569 # determine channelization from (first) ms: 02570 if type(mstoimage)==type([]): 02571 ms0=mstoimage[0] 02572 if len(mstoimage)==1: 02573 mstoimage=mstoimage[0] 02574 else: 02575 ms0=mstoimage 02576 02577 tb.open(ms0+"/SPECTRAL_WINDOW") 02578 if tb.nrows() > 1: 02579 self.msg("determining output cube parameters from FIRST of several SPW in MS "+ms0) 02580 freq=tb.getvarcol("CHAN_FREQ")['r1'] 02581 nchan=freq.size 02582 tb.done() 02583 02584 if nchan==1: 02585 chanmode="mfs" 02586 else: 02587 chanmode="channel" 02588 02589 psfmode="clark" 02590 ftmachine="ft" 02591 02592 if cleanmode=="csclean": 02593 imagermode='csclean' 02594 #if cleanmode=="clark": 02595 # imagermode="" 02596 if cleanmode=="mosaic": 02597 imagermode="mosaic" 02598 ftmachine="mosaic" 02599 02600 # in 3.4 clean doesn't accept just any imsize 02601 from cleanhelper import cleanhelper 02602 optsize=[0,0] 02603 optsize[0]=cleanhelper.getOptimumSize(imsize[0]) 02604 nksize=len(imsize) 02605 if nksize==1: # imsize can be a single element or array 02606 optsize[1]=optsize[0] 02607 else: 02608 optsize[1]=cleanhelper.getOptimumSize(imsize[1]) 02609 if((optsize[0] != imsize[0]) or (nksize!=1 and optsize[1] != imsize[1])): 02610 self.msg(str(imsize)+' is not an acceptable imagesize, will use '+str(optsize)+" instead",priority="warn") 02611 imsize=optsize 02612 02613 # print clean inputs no matter what, so user can use them. 02614 # and write a clean.last file 02615 cleanlast=open(imagename+".clean.last","write") 02616 cleanlast.write('taskname = "clean"\n') 02617 02618 self.msg("clean inputs:") 02619 if type(mstoimage)==type([]): 02620 cleanstr="clean(vis="+str(mstoimage)+",imagename='"+imagename+"'" 02621 cleanlast.write('vis = '+str(mstoimage)+'\n') 02622 else: 02623 cleanstr="clean(vis='"+str(mstoimage)+"',imagename='"+imagename+"'" 02624 cleanlast.write('vis = "'+str(mstoimage)+'"\n') 02625 cleanlast.write('imagename = "'+imagename+'"\n') 02626 cleanlast.write('outlierfile = ""\n') 02627 cleanlast.write('field = "'+sourcefieldlist+'"\n') 02628 cleanlast.write('spw = ""\n') 02629 cleanlast.write('selectdata = False\n') 02630 cleanlast.write('timerange = ""\n') 02631 cleanlast.write('uvrange = ""\n') 02632 cleanlast.write('antenna = ""\n') 02633 cleanlast.write('scan = ""\n') 02634 if nchan>1: 02635 cleanstr=cleanstr+",mode='"+chanmode+"',nchan="+str(nchan) 02636 cleanlast.write('mode = "'+chanmode+'"\n') 02637 cleanlast.write('nchan = "'+str(nchan)+'"\n') 02638 else: 02639 cleanlast.write('mode = "mfs"\n') 02640 cleanlast.write('nchan = -1\n') 02641 cleanlast.write('gridmode = ""\n') 02642 cleanlast.write('wprojplanes = 1\n') 02643 cleanlast.write('facets = 1\n') 02644 cleanlast.write('cfcache = "cfcache.dir"\n') 02645 cleanlast.write('painc = 360.0\n') 02646 cleanlast.write('epjtable = ""\n') 02647 #cleanstr=cleanstr+",interpolation='nearest'" # default change 20100518 02648 cleanlast.write('interpolation = "linear"\n') 02649 cleanstr=cleanstr+",niter="+str(niter) 02650 cleanlast.write('niter = '+str(niter)+'\n') 02651 cleanlast.write('gain = 0.1\n') 02652 cleanstr=cleanstr+",threshold='"+str(threshold)+"'" 02653 cleanlast.write('threshold = "'+str(threshold)+'"\n') 02654 cleanstr=cleanstr+",psfmode='"+psfmode+"'" 02655 cleanlast.write('psfmode = "'+psfmode+'"\n') 02656 if imagermode != "": 02657 cleanstr=cleanstr+",imagermode='"+imagermode+"'" 02658 cleanlast.write('imagermode = "'+imagermode+'"\n') 02659 cleanstr=cleanstr+",ftmachine='"+ftmachine+"'" 02660 cleanlast.write('ftmachine = "'+ftmachine+'"\n') 02661 cleanlast.write('mosweight = False\n') 02662 cleanlast.write('scaletype = "SAULT"\n') 02663 cleanlast.write('multiscale = []\n') 02664 cleanlast.write('negcomponent = -1\n') 02665 cleanlast.write('smallscalebias = 0.6\n') 02666 cleanlast.write('interactive = False\n') 02667 if type(mask)==type(" "): 02668 cleanlast.write('mask = "'+mask+'"\n') 02669 cleanstr=cleanstr+",mask='"+mask+"'" 02670 else: 02671 cleanlast.write('mask = '+str(mask)+'\n') 02672 cleanstr=cleanstr+",mask="+str(mask) 02673 cleanlast.write('start = 0\n') 02674 cleanlast.write('width = 1\n') 02675 cleanlast.write('outframe = ""\n') 02676 cleanlast.write('veltype = "radio"\n') 02677 cleanstr=cleanstr+",imsize="+str(imsize)+",cell="+str(map(qa.tos,cell))+",phasecenter='"+str(imcenter)+"'" 02678 cleanlast.write('imsize = '+str(imsize)+'\n'); 02679 cleanlast.write('cell = '+str(map(qa.tos,cell))+'\n'); 02680 cleanlast.write('phasecenter = "'+str(imcenter)+'"\n'); 02681 cleanlast.write('restfreq = ""\n'); 02682 if stokes != "I": 02683 cleanstr=cleanstr+",stokes='"+stokes+"'" 02684 cleanlast.write('stokes = "'+stokes+'"\n'); 02685 cleanlast.write('weighting = "'+weighting+'"\n'); 02686 cleanstr=cleanstr+",weighting='"+weighting+"'" 02687 if weighting == "briggs": 02688 cleanstr=cleanstr+",robust=0.5" 02689 cleanlast.write('robust = 0.5\n'); 02690 robust=0.5 02691 else: 02692 cleanlast.write('robust = 0.0\n'); 02693 robust=0. 02694 02695 taper=False 02696 if len(outertaper) >0: 02697 taper=True 02698 if type(outertaper) == type([]): 02699 if len(outertaper[0])==0: 02700 taper=False 02701 if taper: 02702 uvtaper=True 02703 cleanlast.write('uvtaper = True\n'); 02704 cleanlast.write('outertaper = "'+str(outertaper)+'"\n'); 02705 cleanstr=cleanstr+",uvtaper=True,outertaper="+str(outertaper)+",innertaper=[]" 02706 else: 02707 uvtaper=False 02708 cleanlast.write('uvtaper = False\n'); 02709 cleanlast.write('outertaper = []\n'); 02710 cleanstr=cleanstr+",uvtaper=False" 02711 cleanlast.write('innertaper = []\n'); 02712 if os.path.exists(modelimage): 02713 cleanstr=cleanstr+",modelimage='"+str(modelimage)+"'" 02714 cleanlast.write('modelimage = "'+str(modelimage)+'"\n'); 02715 else: 02716 cleanlast.write('modelimage = ""\n'); 02717 cleanlast.write("restoringbeam = ['']\n"); 02718 cleanlast.write("pbcor = True\n"); 02719 cleanlast.write("minpb = 0.2\n"); 02720 cleanlast.write("calready = True\n"); 02721 cleanlast.write('noise = ""\n'); 02722 cleanlast.write('npixels = 0\n'); 02723 cleanlast.write('npercycle = 100\n'); 02724 cleanlast.write('cyclefactor = 1.5\n'); 02725 cleanlast.write('cyclespeedup = -1\n'); 02726 cleanlast.write('nterms = 1\n'); 02727 cleanlast.write('reffreq = ""\n'); 02728 cleanlast.write('chaniter = False\n'); 02729 cleanstr=cleanstr+")" 02730 if self.verbose: 02731 self.msg(cleanstr,priority="warn") 02732 else: 02733 self.msg(cleanstr) 02734 cleanlast.write("#"+cleanstr+"\n") 02735 cleanlast.close() 02736 02737 clean(vis=mstoimage, imagename=imagename, mode=chanmode, 02738 niter=niter, threshold=threshold, selectdata=False, nchan=nchan, 02739 psfmode=psfmode, imagermode=imagermode, ftmachine=ftmachine, 02740 imsize=imsize, cell=map(qa.tos,cell), phasecenter=imcenter, 02741 stokes=stokes, weighting=weighting, robust=robust, 02742 uvtaper=uvtaper,outertaper=outertaper, pbcor=True, mask=mask, 02743 modelimage=modelimage) 02744 02745 del freq,nchan # something is holding onto the ms in table cache 02746 02747 02748 02749 02750 02751 02752 02753 ################################################### 02754 02755 def flatimage(self,image,complist="",verbose=False): 02756 # flat output 02757 ia.open(image) 02758 imsize=ia.shape() 02759 imcsys=ia.coordsys() 02760 ia.done() 02761 spectax=imcsys.findcoordinate('spectral')[1] 02762 nchan=imsize[spectax] 02763 stokesax=imcsys.findcoordinate('stokes')[1] 02764 nstokes=imsize[stokesax] 02765 02766 flat=image+".flat" 02767 if nchan>1: 02768 if verbose: self.msg("creating moment zero image "+flat,origin="analysis") 02769 ia.open(image) 02770 flat_ia = ia.moments(moments=[-1],outfile=flat,overwrite=True) 02771 ia.done() 02772 flat_ia.close() 02773 del flat_ia 02774 else: 02775 if verbose: self.msg("removing degenerate image axes in "+flat,origin="analysis") 02776 # just remove degenerate axes from image 02777 flat_ia = ia.newimagefromimage(infile=image,outfile=flat,dropdeg=True,overwrite=True) 02778 flat_ia.close() 02779 02780 # seems no way to just drop the spectral and keep the stokes. 02781 if nstokes<=1: 02782 os.rename(flat,flat+".tmp") 02783 ia.open(flat+".tmp") 02784 flat_ia = ia.adddegaxes(outfile=flat,stokes='I',overwrite=True) 02785 ia.done() 02786 flat_ia.close() 02787 shutil.rmtree(flat+".tmp") 02788 del flat_ia 02789 02790 if nstokes>1: 02791 os.rename(flat,flat+".tmp") 02792 po.open(flat+".tmp") 02793 foo=po.stokesi(outfile=flat,stokes='I') 02794 foo.done() 02795 po.done() 02796 shutil.rmtree(flat+".tmp") 02797 02798 imcsys.done() 02799 del imcsys 02800 02801 # add components 02802 if len(complist)>0: 02803 ia.open(flat) 02804 if not os.path.exists(complist): 02805 self.msg("sky component list "+str(complist)+" not found in flatimage",priority="error") 02806 cl.open(complist) 02807 ia.modify(cl.torecord(),subtract=False) 02808 cl.done() 02809 ia.done() 02810 02811 02812 02813 02814 ################################################### 02815 02816 def convimage(self,modelflat,outflat,complist=""): 02817 # regrid flat input to flat output shape and convolve 02818 modelregrid = modelflat+".regrid" 02819 # get outflatcoordsys from outflat 02820 ia.open(outflat) 02821 outflatcs=ia.coordsys() 02822 outflatshape=ia.shape() 02823 # and beam TODO is beam the same in flat as a cube? 02824 beam=ia.restoringbeam() 02825 ia.done() 02826 02827 ia.open(modelflat) 02828 modelflatcs=ia.coordsys() 02829 modelflatshape=ia.shape() 02830 tmpxx=ia.regrid(outfile=modelregrid+'.tmp', overwrite=True, 02831 csys=outflatcs.torecord(),shape=outflatshape) 02832 # im.regrid assumes a surface brightness, or more accurately doesnt 02833 # pay attention to units at all, so we now have to scale 02834 # by the pixel size to have the right values in jy/pixel, 02835 02836 # get pixel size from model image coordsys 02837 tmpxx.done() 02838 increments=outflatcs.increment(type="direction")['numeric'] 02839 incellx=qa.quantity(abs(increments[0]),outflatcs.units(type="direction")[0]) 02840 incelly=qa.quantity(abs(increments[1]),outflatcs.units(type="direction")[1]) 02841 xform=outflatcs.lineartransform(type="direction") 02842 offdiag=max(abs(xform[0,1]),abs(xform[1,0])) 02843 if offdiag > 1e-4: 02844 self.msg("Your image is rotated with respect to Lat/Lon. I can't cope with that yet",priority="error") 02845 incellx=qa.mul(incellx,abs(xform[0,0])) 02846 incelly=qa.mul(incelly,abs(xform[1,1])) 02847 model_cell = [qa.convert(incellx,'arcsec'),qa.convert(incelly,'arcsec')] 02848 02849 # and from outflat (the clean image) 02850 increments=outflatcs.increment(type="direction")['numeric'] 02851 incellx=qa.quantity(abs(increments[0]),outflatcs.units(type="direction")[0]) 02852 incelly=qa.quantity(abs(increments[1]),outflatcs.units(type="direction")[1]) 02853 xform=outflatcs.lineartransform(type="direction") 02854 offdiag=max(abs(xform[0,1]),abs(xform[1,0])) 02855 if offdiag > 1e-4: 02856 self.msg("Your image is rotated with respect to Lat/Lon. I can't cope with that yet",priority="error") 02857 incellx=qa.mul(incellx,abs(xform[0,0])) 02858 incelly=qa.mul(incelly,abs(xform[1,1])) 02859 cell = [qa.convert(incellx,'arcsec'),qa.convert(incelly,'arcsec')] 02860 02861 # image scaling 02862 factor = (qa.convert(cell[0],"arcsec")['value']) 02863 factor *= (qa.convert(cell[1],"arcsec")['value']) 02864 factor /= (qa.convert(model_cell[0],"arcsec")['value']) 02865 factor /= (qa.convert(model_cell[1],"arcsec")['value']) 02866 imrr = ia.imagecalc(modelregrid, 02867 "'%s'*%g" % (modelregrid+'.tmp',factor), 02868 overwrite = True) 02869 shutil.rmtree(modelregrid+".tmp") 02870 if self.verbose: 02871 self.msg("scaling model by pixel area ratio %g" % factor) 02872 02873 # add unresolved components in Jy/pix 02874 # don't do this if you've already done it in flatimage()! 02875 if (os.path.exists(complist)): 02876 cl.open(complist) 02877 imrr.modify(cl.torecord(),subtract=False) 02878 cl.done() 02879 02880 imrr.done() 02881 ia.done() 02882 del imrr 02883 02884 # Convolve model with beam. 02885 convolved = modelregrid + '.conv' 02886 ia.open(modelregrid) 02887 ia.setbrightnessunit("Jy/pixel") 02888 tmpcnv=ia.convolve2d(convolved,major=beam['major'],minor=beam['minor'], 02889 pa=beam['positionangle'],overwrite=True) 02890 ia.done() 02891 02892 #tmpcnv.open(convolved) 02893 tmpcnv.setbrightnessunit("Jy/beam") 02894 tmpcnv.setrestoringbeam(beam=beam) 02895 tmpcnv.done() 02896 02897 02898 def bandname(self, freq): 02899 """ 02900 Given freq in GHz, return the silly and in some cases deliberately 02901 confusing band name that some people insist on using. 02902 """ 02903 # TODO: add a trivia argument to optionally provide the historical 02904 # radar band info from Wikipedia. 02905 band = '' 02906 if freq > 90: # Use the ALMA system. 02907 band = 'band%.0f' % (0.01 * freq) # () are mandatory! 02908 # Now switch to radar band names. Above 40 GHz different groups used 02909 # different names. Wikipedia uses ones from Baytron, a now defunct company 02910 # that made test equipment. 02911 elif freq > 75: # used as a visual sensor for experimental autonomous vehicles 02912 band = 'W' 02913 elif freq > 50: # Very strongly absorbed by atmospheric O2, 02914 band = 'V' # which resonates at 60 GHz. 02915 elif freq >= 40: 02916 band = 'Q' 02917 elif freq >= 26.5: # mapping, short range, airport surveillance; 02918 band = 'Ka' # frequency just above K band (hence 'a') 02919 # Photo radar operates at 34.300 +/- 0.100 GHz. 02920 elif freq >= 18: 02921 band = 'K' # from German kurz, meaning 'short'; limited use due to 02922 # absorption by water vapor, so Ku and Ka were used 02923 # instead for surveillance. Used for detecting clouds 02924 # and at 24.150 +/- 0.100 GHz for speeding motorists. 02925 elif freq >= 12: 02926 band = 'U' # or Ku 02927 elif freq >= 8: # Missile guidance, weather, medium-resolution mapping and ground 02928 band = 'X' # surveillance; in the USA the narrow range 10.525 GHz +/- 25 MHz 02929 # is used for airport radar and short range tracking. Named X 02930 # band because the frequency was a secret during WW2. 02931 elif freq >= 4: 02932 band = 'C' # Satellite transponders; a compromise (hence 'C') between X 02933 # and S bands; weather; long range tracking 02934 elif freq >= 2: 02935 band = 'S' # Moderate range surveillance, air traffic control, 02936 # long-range weather; 'S' for 'short' 02937 elif freq >= 1: 02938 band = 'L' # Long range air traffic control and surveillance; 'L' for 'long' 02939 elif freq >= 0.3: 02940 band = 'UHF' 02941 else: 02942 band = 'P' # for 'previous', applied retrospectively to early radar systems 02943 # Includes HF and VHF. Worse, it leaves no space for me 02944 # to put in rock band easter eggs. 02945 return band 02946 02947 02948 def polsettings(self, telescope, poltype=None, listall=False): 02949 """ 02950 Returns stokes parameters (for use as stokes in sm.setspwindow) 02951 and feed type (for use as mode in sm.setfeed). 02952 02953 If poltype is not specified or recognized, a guess is made using 02954 telescope. Defaults to 'XX YY', 'perfect X Y' 02955 02956 If listall is True, return the options instead. 02957 """ 02958 psets = {'circular': ('RR LL', 'perfect R L'), 02959 'linear': ('XX YY', 'perfect X Y'), 02960 'RR': ('RR', 'perfect R L')} 02961 if listall: 02962 return psets 02963 if poltype not in psets: 02964 poltype = 'linear' 02965 for circ in ('VLA', 'DRAO'): # Done this way to 02966 if circ in telescope.upper(): # catch EVLA. 02967 poltype = 'circular' 02968 return psets[poltype] 02969 02970 ####################################### 02971 # ALMA calcpointings 02972 02973 def applyRotate(self, x, y, tcos, tsin): 02974 return tcos*x-tsin*y, tsin*x+tcos*y 02975 02976 def isEven(self, num): 02977 return (num % 2) == 0 02978 02979 # this was the only algorithm in Cycle 0 - for Cycle 1 it was 02980 # recast as BaseTriangularTiling.getPointings(), the width>height 02981 # case statement was removed, and BaseTriangularTiling was supplemented by 02982 # ShiftTriangularTiling.getPointings() 02983 def getTrianglePoints(self, width, height, angle, spacing): 02984 tcos = pl.cos(angle*pl.pi/180) 02985 tsin = pl.sin(angle*pl.pi/180) 02986 02987 xx=[] 02988 yy=[] 02989 02990 if (width >= height): 02991 wSpacing = spacing 02992 hSpacing = (pl.sqrt(3.) / 2) * spacing 02993 02994 nwEven = int(pl.floor((width / 2) / wSpacing)) 02995 nwOdd = int(pl.floor((width / 2) / wSpacing + 0.5)) 02996 nh = int(pl.floor((height / 2) / hSpacing)) 02997 02998 for ih in pl.array(range(nh*2+1))-nh: 02999 if (self.isEven(ih)): 03000 for iw in pl.array(range(nwEven*2+1))-nwEven: 03001 x,y = self.applyRotate(iw*wSpacing, ih*hSpacing, tcos, tsin) 03002 xx.append(x) 03003 yy.append(y) 03004 else: 03005 for iw in pl.array(range(nwOdd*2+1))-nwOdd: 03006 x,y = self.applyRotate((iw+0.5)*wSpacing, ih*hSpacing, tcos, tsin) 03007 xx.append(x) 03008 yy.append(y) 03009 else: 03010 wSpacing = (pl.sqrt(3) / 2) * spacing 03011 hSpacing = spacing 03012 03013 nw = int(pl.floor((width / 2) / wSpacing)) 03014 nhEven = int(pl.floor((height / 2) / hSpacing)) 03015 nhOdd = int(pl.floor((height / 2) / hSpacing + 0.5)) 03016 03017 for iw in pl.array(range(nw*2+1))-nw: 03018 if (self.isEven(iw)): 03019 for ih in pl.array(range(nhEven*2+1))-nhEven: 03020 x,y = self.applyRotate(iw*wSpacing, ih*hSpacing, tcos, tsin) 03021 xx.append(x) 03022 yy.append(y) 03023 else: 03024 for ih in pl.array(range(nhOdd*2+1))-nhOdd: 03025 x,y = self.applyRotate(iw*wSpacing, (ih+0.5)*hSpacing, tcos, tsin) 03026 xx.append(x) 03027 yy.append(y) 03028 return xx,yy 03029 03030 03031 03032 def getTriangularTiling(self, longlength, latlength, angle, spacing, pb): 03033 03034 # OT if isLandscape, shortside=Latlength 03035 # else isLandscape=false, shortside=Longlength 03036 03037 if longlength > latlength: # OT isLandscape=True (OT uses >= ) 03038 width=longlength # arcsec 03039 height=latlength # arcsec 03040 shortside=height 03041 else: 03042 width=latlength 03043 height=longlength 03044 shortside=height 03045 angle=angle+90 03046 03047 # which tiling? Base or Shifted (OT getTiling) 03048 vSpacing = (pl.sqrt(3) / 2) * spacing 03049 n = pl.ceil(shortside / vSpacing) 03050 03051 if n % 2 == 0: 03052 return self.getShiftTriangularTiling(width, height, angle, spacing, pb) 03053 else: 03054 return self.getBaseTriangularTiling(width, height, angle, spacing, pb) 03055 03056 def needsFiller(self, length, spacing, bs, npoints): 03057 if length > spacing * (npoints - 1) + bs: 03058 return 1 03059 else: 03060 return 0 03061 03062 def getBaseTriangularTiling(self, width, height, angle, spacing, pb): 03063 tcos = pl.cos(angle*pl.pi/180) 03064 tsin = pl.sin(angle*pl.pi/180) 03065 03066 xx=[] 03067 yy=[] 03068 03069 wSpacing = spacing 03070 hSpacing = (pl.sqrt(3.) / 2) * spacing 03071 03072 nwEven = int(pl.floor((width / 2) / wSpacing)) 03073 nwEven += self.needsFiller(width, wSpacing, pb, nwEven*2+1) 03074 03075 nwOdd = int(pl.floor((width / 2) / wSpacing + 0.5)) 03076 nwOdd += self.needsFiller(width, wSpacing, pb, nwOdd*2) 03077 03078 nh = int(pl.floor((height / 2) / hSpacing)) 03079 nh += self.needsFiller(height, hSpacing, pb, nh*2+1) 03080 03081 for ih in pl.array(range(nh*2+1))-nh: 03082 if (self.isEven(ih)): 03083 for iw in pl.array(range(nwEven*2+1))-nwEven: 03084 x,y = self.applyRotate(iw*wSpacing, ih*hSpacing, tcos, tsin) 03085 xx.append(x) 03086 yy.append(-y) # will require additional testing @ angle>0 03087 else: 03088 for iw in pl.array(range(nwOdd*2))-nwOdd: 03089 x,y = self.applyRotate((iw+0.5)*wSpacing, ih*hSpacing, tcos, tsin) 03090 xx.append(x) 03091 yy.append(-y) 03092 03093 return xx,yy 03094 03095 03096 03097 03098 def getShiftTriangularTiling(self, width, height, angle, spacing, pb): 03099 tcos = pl.cos(angle*pl.pi/180) 03100 tsin = pl.sin(angle*pl.pi/180) 03101 03102 xx=[] 03103 yy=[] 03104 03105 wSpacing = spacing 03106 hSpacing = (pl.sqrt(3.) / 2) * spacing 03107 03108 nwEven = int(pl.floor((width / 2) / wSpacing + 0.5)) 03109 nwEven += self.needsFiller(width, wSpacing, pb, nwEven*2) 03110 03111 nwOdd = int(pl.floor((width / 2) / wSpacing )) 03112 nwOdd += self.needsFiller(width, wSpacing, pb, nwOdd*2+1) 03113 03114 nh = int(pl.floor((height - hSpacing) / 2 / hSpacing +1)) 03115 nh += self.needsFiller(height, hSpacing, pb, nh*2) 03116 03117 for ih in pl.array(range(nh*2))-nh: 03118 if (self.isEven(ih)): 03119 for iw in pl.array(range(nwEven*2))-nwEven: 03120 x,y = self.applyRotate((iw+0.5)*wSpacing, (ih+0.5)*hSpacing, tcos, tsin) 03121 xx.append(x) 03122 yy.append(-y) 03123 else: 03124 for iw in pl.array(range(nwOdd*2+1))-nwOdd: 03125 x,y = self.applyRotate(iw*wSpacing, (ih+0.5)*hSpacing, tcos, tsin) 03126 xx.append(x) 03127 yy.append(-y) 03128 03129 return xx,yy
1.8.0
