import os import time import numpy import scipy import base64 import datetime from io import BytesIO import scipy.interpolate from matplotlib.figure import Figure from utils import TimeTools from utils.patterns import select_pattern from utils import Misc_Routines from utils import Astro_Coords from utils import gaussfit attenuation = numpy.array([[[-21.25,-15.25,-9.25,-3.25,3.25,9.25,15.25,21.25], [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25]], [[21.25,21.25,21.25,21.25,21.25,21.25,21.25,21.25], [15.25,15.25,15.25,15.25,15.25,15.25,15.25,15.25], [09.25,09.25,09.25,09.25,09.25,09.25,09.25,09.25], [03.25,03.25,03.25,03.25,03.25,03.25,03.25,03.25], [-03.25,-03.25,-03.25,-03.25,-03.25,-03.25,-03.25,-03.25], [-09.25,-09.25,-09.25,-09.25,-09.25,-09.25,-09.25,-09.25], [-15.25,-15.25,-15.25,-15.25,-15.25,-15.25,-15.25,-15.25], [-21.25,-21.25,-21.25,-21.25,-21.25,-21.25,-21.25,-21.25]]]) class BField(): def __init__(self,year=None,doy=None,site=1,heights=None,alpha_i=90): """ BField class creates an object to get the Magnetic field for a specific date and height(s). Parameters ---------- year = A scalar giving the desired year. If the value is None (default value) then the current year will be used. doy = A scalar giving the desired day of the year. If the value is None (default va- lue) then the current doy will be used. site = An integer to choose the geographic coordinates of the place where the magne- tic field will be computed. The default value is over Jicamarca (site=1) heights = An array giving the heights (km) where the magnetic field will be modeled By default the magnetic field will be computed at 100, 500 and 1000km. alpha_i = Angle to interpolate the magnetic field. Modification History -------------------- Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009. """ tmp = time.localtime() if year==None: year = tmp[0] if doy==None: doy = tmp[7] self.year = year self.doy = doy self.site = site if heights is None: heights = numpy.array([100,500,1000]) else: heights = numpy.array(heights) self.heights = heights self.alpha_i = alpha_i def getBField(self,maglimits=numpy.array([-7,-7,7,7])): """ getBField models the magnetic field for a different heights in a specific date. Parameters ---------- maglimits = An 4-elements array giving ..... The default value is [-7,-7,7,7]. Return ------ dcos = An 4-dimensional array giving the directional cosines of the magnetic field over the desired place. alpha = An 3-dimensional array giving the angle of the magnetic field over the desi- red place. Modification History -------------------- Converted to Python by Freddy R. Galindo, ROJ, 07 October 2009. """ x_ant = numpy.array([1,0,0]) y_ant = numpy.array([0,1,0]) z_ant = numpy.array([0,0,1]) if self.site==0: title_site = "Magnetic equator" coord_site = numpy.array([-76+52./60.,-11+57/60.,0.5]) elif self.site==1: title_site = 'Jicamarca' coord_site = [-76-52./60.,-11-57/60.,0.5] theta = (45+5.35)*numpy.pi/180. # (50.35 and 1.46 from Fleish Thesis) delta = -1.46*numpy.pi/180 x_ant1 = numpy.roll(self.rotvector(self.rotvector(x_ant,1,delta),3,theta),1) y_ant1 = numpy.roll(self.rotvector(self.rotvector(y_ant,1,delta),3,theta),1) z_ant1 = numpy.roll(self.rotvector(self.rotvector(z_ant,1,delta),3,theta),1) ang0 = -1*coord_site[0]*numpy.pi/180. ang1 = coord_site[1]*numpy.pi/180. x_ant = self.rotvector(self.rotvector(x_ant1,2,ang1),3,ang0) y_ant = self.rotvector(self.rotvector(y_ant1,2,ang1),3,ang0) z_ant = self.rotvector(self.rotvector(z_ant1,2,ang1),3,ang0) else: # print "No defined Site. Skip..." return None nhei = self.heights.size pt_intercep = numpy.zeros((nhei,2)) nfields = 1 grid_res = 0.5 nlon = int(numpy.int(maglimits[2] - maglimits[0])/grid_res + 1) nlat = int(numpy.int(maglimits[3] - maglimits[1])/grid_res + 1) location = numpy.zeros((nlon,nlat,2)) mlon = numpy.atleast_2d(numpy.arange(nlon)*grid_res + maglimits[0]) mrep = numpy.atleast_2d(numpy.zeros(nlat) + 1) location0 = numpy.dot(mlon.transpose(),mrep) mlat = numpy.atleast_2d(numpy.arange(nlat)*grid_res + maglimits[1]) mrep = numpy.atleast_2d(numpy.zeros(nlon) + 1) location1 = numpy.dot(mrep.transpose(),mlat) location[:,:,0] = location0 location[:,:,1] = location1 alpha = numpy.zeros((nlon,nlat,nhei)) rr = numpy.zeros((nlon,nlat,nhei,3)) dcos = numpy.zeros((nlon,nlat,nhei,2)) global first_time first_time = None for ilon in numpy.arange(nlon): for ilat in numpy.arange(nlat): outs = self.__bdotk(self.heights, self.year + self.doy/366., coord_site[1], coord_site[0], coord_site[2], coord_site[1]+location[ilon,ilat,1], location[ilon,ilat,0]*720./180.) alpha[ilon, ilat,:] = outs[1] rr[ilon, ilat,:,:] = outs[3] mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose() tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(x_ant)) tmp = tmp.sum(axis=1) dcos[ilon,ilat,:,0] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1)) mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose() tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(y_ant)) tmp = tmp.sum(axis=1) dcos[ilon,ilat,:,1] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1)) return dcos, alpha, nlon, nlat def __bdotk(self,heights,tm,gdlat=-11.95,gdlon=-76.8667,gdalt=0.0,decd=-12.88, ham=-4.61666667): global first_time # Mean Earth radius in Km WGS 84 a_igrf = 6371.2 bk = numpy.zeros(heights.size) alpha = numpy.zeros(heights.size) bfm = numpy.zeros(heights.size) rr = numpy.zeros((heights.size,3)) rgc = numpy.zeros((heights.size,3)) ObjGeodetic = Astro_Coords.Geodetic(gdlat,gdalt) [gclat,gcalt] = ObjGeodetic.change2geocentric() gclat = gclat*numpy.pi/180. gclon = gdlon*numpy.pi/180. # Antenna position from center of Earth ca_vector = [numpy.cos(gclat)*numpy.cos(gclon),numpy.cos(gclat)*numpy.sin(gclon),numpy.sin(gclat)] ca_vector = gcalt*numpy.array(ca_vector) dec = decd*numpy.pi/180. # K vector respect to the center of earth. klon = gclon + ham*numpy.pi/720. k_vector = [numpy.cos(dec)*numpy.cos(klon),numpy.cos(dec)*numpy.sin(klon),numpy.sin(dec)] k_vector = numpy.array(k_vector) for ih in numpy.arange(heights.size): # Vector from Earth's center to volume of interest rr[ih,:] = k_vector*heights[ih] cv_vector = numpy.squeeze(ca_vector) + rr[ih,:] cv_gcalt = numpy.sqrt(numpy.sum(cv_vector**2.)) cvxy = numpy.sqrt(numpy.sum(cv_vector[0:2]**2.)) radial = cv_vector/cv_gcalt east = numpy.array([-1*cv_vector[1],cv_vector[0],0])/cvxy comp1 = east[1]*radial[2] - radial[1]*east[2] comp2 = east[2]*radial[0] - radial[2]*east[0] comp3 = east[0]*radial[1] - radial[0]*east[1] north = -1*numpy.array([comp1, comp2, comp3]) rr_k = cv_vector - numpy.squeeze(ca_vector) u_rr = rr_k/numpy.sqrt(numpy.sum(rr_k**2.)) cv_gclat = numpy.arctan2(cv_vector[2],cvxy) cv_gclon = numpy.arctan2(cv_vector[1],cv_vector[0]) bhei = cv_gcalt-a_igrf blat = cv_gclat*180./numpy.pi blon = cv_gclon*180./numpy.pi bfield = self.__igrfkudeki(bhei,tm,blat,blon) B = (bfield[0]*north + bfield[1]*east - bfield[2]*radial)*1.0e-5 bfm[ih] = numpy.sqrt(numpy.sum(B**2.)) #module bk[ih] = numpy.sum(u_rr*B) alpha[ih] = numpy.arccos(bk[ih]/bfm[ih])*180/numpy.pi rgc[ih,:] = numpy.array([cv_gclon, cv_gclat, cv_gcalt]) return bk, alpha, bfm, rr, rgc def __igrfkudeki(self,heights,time,latitude,longitude,ae=6371.2): """ __igrfkudeki calculates the International Geomagnetic Reference Field for given in- put conditions based on IGRF2005 coefficients. Parameters ---------- heights = Scalar or vector giving the height above the Earth of the point in ques- tion in kilometers. time = Scalar or vector giving the decimal year of time in question (e.g. 1991.2). latitude = Latitude of point in question in decimal degrees. Scalar or vector. longitude = Longitude of point in question in decimal degrees. Scalar or vector. ae = first_time = Return ------ bn = be = bd = bmod = balpha = first_time = Modification History -------------------- Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009. """ global first_time global gs, hs, nvec, mvec, maxcoef heights = numpy.atleast_1d(heights) time = numpy.atleast_1d(time) latitude = numpy.atleast_1d(latitude) longitude = numpy.atleast_1d(longitude) if numpy.max(latitude)==90: # print "Field calculations are not supported at geographic poles" pass # output arrays bn = numpy.zeros(heights.size) be = numpy.zeros(heights.size) bd = numpy.zeros(heights.size) if first_time==None:first_time=0 time0 = time[0] if time!=first_time: #print "Getting coefficients for", time0 [periods,g,h ] = self.__readIGRFcoeff() top_year = numpy.max(periods) nperiod = (top_year - 1900)/5 + 1 maxcoef = 10 if time0>=2000:maxcoef = 12 # Normalization array for Schmidt fucntions multer = numpy.zeros((2+maxcoef,1+maxcoef)) + 1 for cn in (numpy.arange(maxcoef)+1): for rm in (numpy.arange(cn)+1): tmp = numpy.arange(2*rm) + cn - rm + 1. multer[rm+1,cn] = ((-1.)**rm)*numpy.sqrt(2./tmp.prod()) schmidt = multer[1:,1:].transpose() # n and m arrays nvec = numpy.atleast_2d(numpy.arange(maxcoef)+2) mvec = numpy.atleast_2d(numpy.arange(maxcoef+1)).transpose() # Time adjusted igrf g and h with Schmidt normalization # IGRF coefficient arrays: g0(n,m), n=1, maxcoeff,m=0, maxcoeff, ... if time0 top_year dtime = (time0 - top_year) + 5 ntime = g[:,0,0].size - 2 g0 = g[ntime,1:maxcoef+1,:maxcoef+1] h0 = h[ntime,1:maxcoef+1,:maxcoef+1] gdot = g[ntime+1,1:maxcoef+1,:maxcoef+1]-g[ntime,1:maxcoef+1,:maxcoef+1] hdot = h[ntime+1,1:maxcoef+1,:maxcoef+1]-h[ntime,1:maxcoef+1,:maxcoef+1] gs = (g0 + dtime*(gdot/5.))*schmidt[:maxcoef,0:maxcoef+1] hs = (h0 + dtime*(hdot/5.))*schmidt[:maxcoef,0:maxcoef+1] first_time = time0 for ii in numpy.arange(heights.size): # Height dependence array rad = (ae/(ae+height))**(n+3) rad = numpy.atleast_2d((ae/(ae + heights[ii]))**(nvec+1)) # Sin and Cos of m times longitude phi arrays mphi = mvec*longitude[ii]*numpy.pi/180. cosmphi = numpy.atleast_2d(numpy.cos(mphi)) sinmphi = numpy.atleast_2d(numpy.sin(mphi)) # Cos of colatitude theta c = numpy.cos((90 - latitude[ii])*numpy.pi/180.) # Legendre functions p(n,m|c) [p,dp]= scipy.special.lpmn(maxcoef+1,maxcoef+1,c) p = p[:,:-1].transpose() s = numpy.sqrt((1. - c)*(1 + c)) # Generate derivative array dpdtheta = -s*dpdc dpdtheta = c*p/s for m in numpy.arange(maxcoef+2): dpdtheta[:,m] = m*dpdtheta[:,m] dpdtheta = dpdtheta + numpy.roll(p,-1,axis=1) # Extracting arrays required for field calculations p = p[1:maxcoef+1,:maxcoef+1] dpdtheta = dpdtheta[1:maxcoef+1,:maxcoef+1] # Weigh p and dpdtheta with gs and hs coefficients. gp = gs*p hp = hs*p gdpdtheta = gs*dpdtheta hdpdtheta = hs*dpdtheta # Calcultate field components matrix0 = numpy.dot(gdpdtheta,cosmphi) matrix1 = numpy.dot(hdpdtheta,sinmphi) bn[ii] = numpy.dot(rad,(matrix0 + matrix1)) matrix0 = numpy.dot(hp,(mvec*cosmphi)) matrix1 = numpy.dot(gp,(mvec*sinmphi)) be[ii] = numpy.dot((-1*rad),((matrix0 - matrix1)/s)) matrix0 = numpy.dot(gp,cosmphi) matrix1 = numpy.dot(hp,sinmphi) bd[ii] = numpy.dot((-1*nvec*rad),(matrix0 + matrix1)) bmod = numpy.sqrt(bn**2. + be**2. + bd**2.) btheta = numpy.arctan(bd/numpy.sqrt(be**2. + bn**2.))*180/numpy.pi balpha = numpy.arctan(be/bn)*180./numpy.pi #bn : north #be : east #bn : radial #bmod : module return bn, be, bd, bmod, btheta, balpha def str2num(self, datum): try: return int(datum) except: try: return float(datum) except: return datum def __readIGRFfile(self, filename): list_years=[] for i in range(1,24): list_years.append(1895.0 + i*5) epochs=list_years epochs.append(epochs[-1]+5) nepochs = numpy.shape(epochs) gg = numpy.zeros((13,14,nepochs[0]),dtype=float) hh = numpy.zeros((13,14,nepochs[0]),dtype=float) coeffs_file=open(filename) lines=coeffs_file.readlines() coeffs_file.close() for line in lines: items = line.split() g_h = items[0] n = self.str2num(items[1]) m = self.str2num(items[2]) coeffs = items[3:] for i in range(len(coeffs)-1): coeffs[i] = self.str2num(coeffs[i]) #coeffs = numpy.array(coeffs) ncoeffs = numpy.shape(coeffs)[0] if g_h == 'g': # print n," g ",m gg[n-1,m,:]=coeffs elif g_h=='h': # print n," h ",m hh[n-1,m,:]=coeffs # else : # continue # Ultimo Reordenamiento para almacenar . gg[:,:,nepochs[0]-1] = gg[:,:,nepochs[0]-2] + 5*gg[:,:,nepochs[0]-1] hh[:,:,nepochs[0]-1] = hh[:,:,nepochs[0]-2] + 5*hh[:,:,nepochs[0]-1] # return numpy.array([gg,hh]) periods = numpy.array(epochs) g = gg h = hh return periods, g, h def __readIGRFcoeff(self,filename="igrf10coeffs.dat"): """ __readIGRFcoeff reads the coefficients from a binary file which is located in the folder "resource." Parameter --------- filename = A string to specify the name of the file which contains thec coeffs. The default value is "igrf10coeffs.dat" Return ------ periods = A lineal array giving... g1 = h1 = Modification History -------------------- Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009. """ # # igrfile = sys.path[-1] + os.sep + "resource" + os.sep + filename # igrfile = os.path.join('./resource',filename) # f = open(igrfile,'rb') # #f = open(os.getcwd() + os.sep + "resource" + os.sep + filename,'rb') # # # Reading SkyNoise Power (lineal scale) # periods = numpy.fromfile(f,numpy.dtype([('var',' Around "x" axis = 2 -> Around "y" axis = 3 -> Around "z" ang = Angle of rotation (in radians). The default value is zero. Return ------ rotvector = A lineal array of 3 elements giving the new coordinates. Modification History -------------------- Converted to Python by Freddy R. Galindo, ROJ, 01 October 2009. """ if axis==1: t = [[1,0,0],[0,numpy.cos(ang),numpy.sin(ang)],[0,-numpy.sin(ang),numpy.cos(ang)]] elif axis==2: t = [[numpy.cos(ang),0,-numpy.sin(ang)],[0,1,0],[numpy.sin(ang),0,numpy.cos(ang)]] elif axis==3: t = [[numpy.cos(ang),numpy.sin(ang),0],[-numpy.sin(ang),numpy.cos(ang),0],[0,0,1]] rotvector = numpy.array(numpy.dot(numpy.array(t),numpy.array(vector))) return rotvector class AntPatternPlot: def __init__(self): """ AntPatternPlot creates an object to call methods to plot the antenna pattern. Modification History -------------------- Created by Freddy Galindo, ROJ, 06 October 2009. """ self.fig = Figure(figsize=(8,8), facecolor='white') self.ax = self.fig.add_subplot(111) def contPattern(self,iplot=0,gpath='',filename='',mesg='',amp=None ,x=None ,y=None ,getCut=None,title='', save=False): """ contPattern plots a contour map of the antenna pattern. Parameters ---------- iplot = A integer to specify if the plot is the first, second, ... The default va- lue is 0. Examples -------- >> Over_Jro.JroPattern(pattern=2).contPattern() Modification history -------------------- Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. """ if getCut == 1: return xmax = numpy.max(x) xmin = numpy.min(x) ymax = numpy.max(y) ymin = numpy.min(y) levels = numpy.array([1e-3,1e-2,1e-1,0.5,1.0]) tmp = numpy.round(10*numpy.log10(levels),decimals=1) labels = list(range(5)) for i in numpy.arange(5): labels[i] = str(numpy.int(tmp[i])) colors = ((0,0,1.),(0,170/255.,0),(127/255.,1.,0),(1.,109/255.,0),(128/255.,0,0)) CS = self.ax.contour(x,y,amp.transpose(),levels,colors=colors) fmt = {} for l,s in zip(CS.levels,labels): fmt[l] = s self.ax.annotate('Ng',xy=(-0.05,1.04),xytext=(0.01,0.962),xycoords='axes fraction',arrowprops=dict(facecolor='black', width=1.,shrink=0.2),fontsize=15.) self.ax.annotate(mesg,xy=(0,0),xytext=(0.01,0.01),xycoords='figure fraction') self.ax.clabel(CS,CS.levels,inline=True,fmt=fmt,fontsize=10) self.ax.set_xlim(xmin,xmax) self.ax.set_ylim(ymin,ymax) self.ax.set_title("Total Pattern: " + title) self.ax.set_xlabel("West to South") self.ax.set_ylabel("West to North") self.ax.grid(True) def plotRaDec(self,gpath=None,filename=None,jd=2452640.5,ra_obs=None,xg=None,yg=None,x=None,y=None, save=True): """ plotRaDec draws right ascension and declination lines on a JRO plane. This function must call after conPattern. Parameters ---------- jd = A scalar giving the Julian date. ra_obs = Scalar giving the right ascension of the observatory. xg = A 3-element array to specify .. yg = A 3-element array to specify .. Examples -------- >> Over_Jro.JroPattern(pattern=2).contPattern() >> Over_Jro.JroPattern(pattern=2).plotRaDec(jd=jd,ra_obs=ra_obs,xg=xg,yg=yg) Modification history -------------------- Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. """ # Finding RA of observatory for a specific date if ra_obs is None:ra_obs = numpy.array([23.37060849]) if xg is None:xg = numpy.array([0.62918474,-0.77725579,0.]) if yg is None:yg = numpy.array([0.77700346,0.62898048,0.02547905]) # Getting HA and DEC axes mindec = -28; maxdec = 4; incdec = 2. ndec = numpy.int((maxdec - mindec)/incdec) + 1 minha = -20; maxha = 20; incha = 2. nha = numpy.int((maxha - minha)/incha) + 1 #mcosx = numpy.zeros((nha,ndec)) #mcosy = numpy.zeros((nha,ndec)) ha_axes = numpy.reshape(numpy.arange(nha)*incha + minha,(nha,1)) ones_dec = numpy.reshape(numpy.zeros(ndec) + 1,(ndec,1)) ha_axes = numpy.dot(ha_axes,ones_dec.transpose()) ha_axes2 = numpy.array(ra_obs - ha_axes) dec_axes = numpy.reshape(numpy.arange(ndec)*incdec + mindec,(ndec,1)) ones_ra = numpy.reshape(numpy.zeros(nha) + 1,(nha,1)) dec_axes = numpy.dot(ones_ra,dec_axes.transpose()) dec_axes2 = numpy.array(dec_axes) ObjHor = Astro_Coords.Equatorial(ha_axes2,dec_axes2,jd) [alt,az,ha] = ObjHor.change2AltAz() z = numpy.transpose(alt)*Misc_Routines.CoFactors.d2r ; z = z.flatten() az = numpy.transpose(az)*Misc_Routines.CoFactors.d2r ; az = az.flatten() vect = numpy.array([numpy.cos(z)*numpy.sin(az),numpy.cos(z)*numpy.cos(az),numpy.sin(z)]) xg = numpy.atleast_2d(xg) dcosx = numpy.array(numpy.dot(xg,vect)) yg = numpy.atleast_2d(yg) dcosy = numpy.array(numpy.dot(yg,vect)) mcosx = dcosx.reshape(ndec,nha) mcosy = dcosy.reshape(ndec,nha) # Defining NAN for points outof limits. xmax = numpy.max(x) xmin = numpy.min(x) ymax = numpy.max(y) ymin = numpy.min(y) factor = 1.3 noval = numpy.where((mcosx>(xmax*factor)) | (mcosx<(xmin*factor))) if noval[0].size>0:mcosx[noval] = numpy.nan noval = numpy.where((mcosy>(ymax*factor)) | (mcosy<(ymin*factor))) if noval[0].size>0:mcosy[noval] = numpy.nan # Plotting HA and declination grid. iha0 = numpy.int((0 - minha)/incha) idec0 = numpy.int((-14 - mindec)/incdec) colorgrid = (1.,109/255.,0) self.ax.plot(mcosx.transpose(),mcosy.transpose(),color=colorgrid,linestyle='--', lw=0.5) for idec in numpy.arange(ndec): if idec != idec0: valx = (mcosx[idec,iha0]<=xmax) & (mcosx[idec,iha0]>=xmin) valy = (mcosy[idec,iha0]<=ymax) & (mcosy[idec,iha0]>=ymin) if valx & valy: text = str(numpy.int(mindec + incdec*idec))+'$^o$' self.ax.text(mcosx[idec,iha0],mcosy[idec,iha0],text) self.ax.plot(mcosx,mcosy,color=colorgrid,linestyle='--',lw=0.5) for iha in numpy.arange(nha): if iha != iha0: valx = (mcosx[idec0,iha]<=xmax) & (mcosx[idec0,iha]>=xmin) valy = (mcosy[idec0,iha]<=ymax) & (mcosy[idec0,iha]>=ymin) if valx & valy: text = str(4*numpy.int(minha + incha*iha))+"'" self.ax.text(mcosx[idec0,iha],mcosy[idec0,iha],text) if save: save_fig = os.path.join(gpath,filename) self.fig.savefig(save_fig,format='png') def plotBField(self,gpath,filename,dcos,alpha, nlon, nlat, dcosxrange, dcosyrange, heights, alpha_i, save=False): """ plotBField draws the magnetic field in a directional cosines plot. Parameters ---------- dcos = An 4-dimensional array giving the directional cosines of the magnetic field over the desired place. alpha = An 3-dimensional array giving the angle of the magnetic field over the desi- red place. nlon = An integer to specify the number of elements per longitude. nlat = An integer to specify the number of elements per latitude. dcosxrange = A 2-element array giving the range of the directional cosines in the "x" axis. dcosyrange = A 2-element array giving the range of the directional cosines in the "y" axis. heights = An array giving the heights (km) where the magnetic field will be modeled By default the magnetic field will be computed at 100, 500 and 1000km. alpha_i = Angle to interpolate the magnetic field. Modification History -------------------- Converted to Python by Freddy R. Galindo, ROJ, 07 October 2009. """ handles = [] objects = [] colors = ['k','m','c','b','g','r','y'] marker = ['-+','-*','-D','-x','-s','->','-o','-^'] alpha_location = numpy.zeros((nlon,2,heights.size)) for ih in numpy.arange(heights.size): alpha_location[:,0,ih] = dcos[:,0,ih,0] for ilon in numpy.arange(nlon): myx = (alpha[ilon,:,ih])[::-1] myy = (dcos[ilon,:,ih,0])[::-1] tck = scipy.interpolate.splrep(myx,myy,s=0) mydcosx = scipy.interpolate.splev(alpha_i,tck,der=0) myx = (alpha[ilon,:,ih])[::-1] myy = (dcos[ilon,:,ih,1])[::-1] tck = scipy.interpolate.splrep(myx,myy,s=0) mydcosy = scipy.interpolate.splev(alpha_i,tck,der=0) alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy]) ObjFig, = self.ax.plot(alpha_location[:,0,ih],alpha_location[:,1,ih], marker[ih % 8],color=colors[numpy.int(ih/8)],ms=4.5,lw=0.5) handles.append(ObjFig) objects.append(numpy.str(heights[ih]) + ' km') legend = self.ax.legend(handles, objects,loc="lower right", numpoints=1, handlelength=0.3, handletextpad=0.02, borderpad=0.3, labelspacing=0.1) self.ax.add_artist(legend) def plotCelestial(self, jd, main_dec, tod, maxha_min, objects, glat, glon, xg, yg, dcosxrange, dcosyrange): self.tod = tod self.dcosx_sun = 1 self.dcosy_sun = 1 self.ha_sun = 1 self.time_sun = 1 self.dcosx_moon = 1 self.dcosy_moon = 1 self.ha_moon = 1 self.time_moon = 1 self.dcosx_hydra = 1 self.dcosy_hydra = 1 self.ha_hydra = 1 self.time_hydra = 1 self.dcosx_galaxy = 1 self.dcosy_galaxy = 1 self.ha_galaxy = 1 self.time_galaxy = 1 tod = self.tod maxlev = 24; minlev = 0; maxcol = 39; mincol = 10 handles = [] titles = ['$Sun$','$Moon$','$Hydra$','$Galaxy$'] marker = ['--^','--s','--*','--o'] # Getting RGB table to plot celestial object over Jicamarca colortable = ['olive', 'indigo', 'cyan', 'red'] labels = [] for io in objects: ObjBodies = Astro_Coords.CelestialBodies() if io==0: [ra,dec,sunlon,sunobliq] = ObjBodies.sunpos(jd) elif io==1: [ra,dec,dist,moonlon,moonlat] = ObjBodies.moonpos(jd) elif io==2: [ra,dec] = ObjBodies.hydrapos() elif io==3: [maxra,ra] = ObjBodies.skynoise_jro(dec_cut=main_dec) ra = maxra*15. dec = main_dec ObjEq = Astro_Coords.Equatorial(ra, dec, jd, lat=glat, lon=glon) [alt, az, ha] = ObjEq.change2AltAz() vect = numpy.array([az,alt]).transpose() vect = Misc_Routines.Vector(vect,direction=0).Polar2Rect() dcosx = numpy.array(numpy.dot(vect,xg)) dcosy = numpy.array(numpy.dot(vect,yg)) wrap = numpy.where(ha>=180.) if wrap[0].size>0: ha[wrap] = ha[wrap] - 360. val = numpy.where((numpy.abs(ha))<=(maxha_min*0.25)) if val[0].size>2: tod_1 = tod*1. shift_1 = numpy.where(tod>12.) tod_1[shift_1] = tod_1[shift_1] - 24. tod_2 = tod*1. shift_2 = numpy.where(tod<12.) tod_2[shift_2] = tod_2[shift_2] + 24. diff0 = numpy.nanmax(tod[val]) - numpy.nanmin(tod[val]) diff1 = numpy.nanmax(tod_1[val]) - numpy.nanmin(tod_1[val]) diff2 = numpy.nanmax(tod_2[val]) - numpy.nanmin(tod_2[val]) if ((diff0<=diff1) & (diff0<=diff2)): tod_0 = tod elif ((diff1> [pattern, mean_pos] = JroPattern(pattern=2).getPattern() >> print meanpos [ 8.08728085e-14 -4.78193873e-14] Modification history -------------------- Developed by Jorge L. Chau. Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. """ if (self.fftopt>0) and (self.getcut>0): #print "Conflict bewteen fftopt and getcut" #print "To get a cut of the antenna pattern uses ffopt=0" return None, None if (self.fftopt==0): # Getting antenna pattern using the array method self.pattern = self.__usingArray(rx=1) if (self.justrx==0):self.pattern = self.pattern*self.__usingArray(rx=0) elif (self.fftopt>0): # Getting antenna pattern using FFT method self.pattern = self.__usingFFT(rx=1) if (self.justrx==0):self.pattern = self.pattern*self.__usingFFT(rx=0) self.maxpattern = numpy.nanmax(self.pattern) self.norpattern = self.pattern/self.maxpattern if self.getcut==0:self.__getBeamPars() def __usingArray(self,rx): """ __usingArray method returns the Jicamarca antenna pattern computed using array model pattern = dipolepattern x modulepattern Parameters ---------- rx = Set to 1 to use the Rx information. Otherwise set to 0 for Tx. Return ------ pattern = An array giving the modelled antenna pattern using the array model. Modification history -------------------- Developed by Jorge L. Chau. Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. """ if rx==1: ues = self.uesrx phase = self.phaserx gain = self.gainrx elif rx==0: ues = self.uestx phase = self.phasetx gain = self.gaintx ues = ues*360./self.airwl phase = phase*360./self.airwl for ii in range(4): if ii==0:dim = numpy.array([4,0,8,4]) # WEST elif ii==1:dim = numpy.array([0,0,4,4]) # NORTH elif ii==2:dim = numpy.array([0,4,4,8]) # EAST elif ii==3:dim = numpy.array([4,4,8,8]) # SOUTH xi = dim[0]; xf = dim[2]; yi = dim[1]; yf = dim[3] phase[xi:xf,yi:yf] = phase[xi:xf,yi:yf] + ues[ii] phase = -phase ar = self.eomwl*numpy.array([[0.5,6., 24.5],[0.5,6.,24.5]]) nr = numpy.array([[12.,4.,2.],[12.,4.,2.]]) lr = 0.25*self.eomwl*numpy.array([[0,0.,0],[0.,0,0]]) # Computing module and dipole patterns. pattern = (numpy.abs(self.__dipPattern(ar,nr,lr)*self.__modPattern(phase,gain)))**2 return pattern def __usingFFT(self,rx): """ __usingFFT method returns the Jicamarca antenna pattern computed using The Fast Fou- rier Transform. pattern = iFFT(FFT(gain*EXP(j*phase))) Parameters ---------- rx = Set to 1 to use the Rx information. Otherwise set to 0 for Tx. Return ------ pattern = An array giving the modelled antenna pattern using the array model. Modification history -------------------- Developed by Jorge L. Chau. Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. """ if rx==1: ues = self.uesrx phase = self.phaserx gain = self.gainrx elif rx==0: ues = self.uestx phase = self.phasetx gain = self.gaintx ues = ues*360./self.airwl phase = phase*360./self.airwl for ii in range(4): if ii==0:dim = numpy.array([4,0,8,4]) # WEST elif ii==1:dim = numpy.array([0,0,4,4]) # NORTH elif ii==2:dim = numpy.array([0,4,4,8]) # EAST elif ii==3:dim = numpy.array([4,4,8,8]) # SOUTH xi = dim[0]; xf = dim[2]; yi = dim[1]; yf = dim[3] phase[xi:xf,yi:yf] = phase[xi:xf,yi:yf] + ues[ii] phase = -phase delta_x = self.eomwl/2. delta_y = self.eomwl/2. nxfft = 2048 nyfft = 2048 dcosx = (numpy.arange(nxfft) - (0.5*nxfft))/(nxfft*delta_x)*self.eomwl dcosy = (numpy.arange(nyfft) - (0.5*nyfft))/(nyfft*delta_y)*self.eomwl fft_gain = numpy.zeros((nxfft,nyfft)) fft_phase = numpy.zeros((nxfft,nyfft)) nx = 8 ny = 8 ndx =12 ndy =12 for iy in numpy.arange(ny): for ix in numpy.arange(nx): ix1 = nxfft/2-self.nx/2*ndx+ix*ndx if ix<(nx/2):ix1 = ix1 - 1 if ix>=(nx/2):ix1 = ix1 + 1 iy1 = nyfft/2-ny/2*ndx+iy*ndy if iy<(ny/2):iy1 = iy1 - 1 if iy>=(ny/2):iy1 = iy1 + 1 fft_gain[ix1:ix1+ndx-1,iy1:iy1+ndy-1] = gain[ix,ny-1-iy] fft_phase[ix1:ix1+ndx-1,iy1:iy1+ndy-1] = phase[ix,ny-1-iy] fft_phase = fft_phase*Misc_Routines.CoFactors.d2r pattern = numpy.abs(numpy.fft.fft2(fft_gain*numpy.exp(numpy.complex(0,1)*fft_phase)))**2 pattern = numpy.fft.fftshift(pattern) xvals = numpy.where((dcosx>=(numpy.min(self.dcosx))) & (dcosx<=(numpy.max(self.dcosx)))) yvals = numpy.where((dcosy>=(numpy.min(self.dcosy))) & (dcosy<=(numpy.max(self.dcosy)))) pattern = pattern[xvals[0][0]:xvals[0][-1],yvals[0][0]:yvals[0][-1]] return pattern def __dipPattern(self,ar,nr,lr): """ _dipPattern function computes the dipole's pattern to the Jicamarca radar. The next equation defines the pattern as a function of the mainlobe direction: sincx = SIN(k/2*n0x*(a0x*SIN(phi)*COS(alpha)))/SIN(k/2*(a0x*SIN(phi)*COS(alpha))) sincy = SIN(k/2*n0y*(a0y*SIN(phi)*SIN(alpha)))/SIN(k/2*(a0y*SIN(phi)*SIN(alpha))) A0(phi,alpha) = sincx*sincy Parameters ---------- ar = ? nr = ? lr = ? Return ------ dipole = An array giving antenna pattern from the dipole point of view.. Modification history -------------------- Developed by Jorge L. Chau. Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. """ dipole = numpy.zeros((self.nx,self.ny),dtype=complex) for iy in range(self.ny): for ix in range(self.nx): yindex = iy*(self.getcut==0) + ix*(self.getcut==1) argx = ar[0,0]*self.dcosx[ix] - lr[0,0] if argx == 0.0: junkx = nr[0,0] else: junkx = numpy.sin(0.5*self.kk*nr[0,0]*argx)/numpy.sin(0.5*self.kk*argx) argy = ar[1,0]*self.dcosy[yindex] - lr[1,0] if argy == 0.0: junky = nr[1,0] else: junky = numpy.sin(0.5*self.kk*nr[1,0]*argy)/numpy.sin(0.5*self.kk*argy) dipole[ix,iy] = junkx*junky return dipole def __modPattern(self,phase,gain): """ ModPattern computes the module's pattern to the Jicamarca radar. The next equation defines the pattern as a function mainlobe direction: phasex = pos(x)*SIN(phi)*COS(alpha) phasey = pos(y)*SIN(phi)*SIN(alpha) A1(phi,alpha) = TOTAL(gain*EXP(COMPLEX(0,k*(phasex+phasey)+phase))) Parameters ---------- phase = Bidimensional array (8x8) giving the phase (in meters) of each module. gain = Bidimensional array (8x8) giving to define modules will be active (ones) and which will not (zeros). Return ------ module = An array giving antenna pattern from the module point of view.. Modification history -------------------- Developed by Jorge L. Chau. Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. """ pos = self.eomwl*attenuation posx = pos[0,:,:] posy = pos[1,:,:] phase = phase*Misc_Routines.CoFactors.d2r module = numpy.zeros((self.nx,self.ny),dtype=complex) for iy in range(self.ny): for ix in range(self.nx): yindex = iy*(self.getcut==0) + ix*(self.getcut==1) phasex = posx*self.dcosx[ix] phasey = posy*self.dcosy[yindex] tmp = gain*numpy.exp(numpy.complex(0,1.)*(self.kk*(phasex+phasey)+phase)) module[ix,iy] = tmp.sum() return module def __getBeamPars(self): """ _getBeamPars computes the main-beam parameters of the antenna. Modification history -------------------- Developed by Jorge L. Chau. Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. """ dx = self.dcosx[1] - self.dcosx[0] dy = self.dcosy[1] - self.dcosy[0] amp = self.norpattern xx = numpy.resize(self.dcosx,(self.nx,self.nx)).transpose() yy = numpy.resize(self.dcosy,(self.ny,self.ny)) mm0 = amp[numpy.where(amp > 0.5)] xx0 = xx[numpy.where(amp > 0.5)] yy0 = yy[numpy.where(amp > 0.5)] xc = numpy.sum(mm0*xx0)/numpy.sum(mm0) yc = numpy.sum(mm0*yy0)/numpy.sum(mm0) rc = numpy.sqrt(mm0.size*dx*dy/numpy.pi) nnx = numpy.where(numpy.abs(self.dcosx - xc) < rc) nny = numpy.where(numpy.abs(self.dcosy - yc) < rc) mm1 = amp[numpy.min(nnx):numpy.max(nnx)+1,numpy.min(nny):numpy.max(nny)+1] xx1 = self.dcosx[numpy.min(nnx):numpy.max(nnx)+1] yy1 = self.dcosy[numpy.min(nny):numpy.max(nny)+1] # fitting data into the main beam. params = gaussfit.fitgaussian(mm1) # Tranforming from indexes to axis' values xcenter = xx1[0] + (((xx1[xx1.size-1] - xx1[0])/(xx1.size -1))*(params[1])) ycenter = yy1[0] + (((yy1[yy1.size-1] - yy1[0])/(yy1.size -1))*(params[2])) xwidth = ((xx1[xx1.size-1] - xx1[0])/(xx1.size-1))*(params[3])*(1/Misc_Routines.CoFactors.d2r) ywidth = ((yy1[yy1.size-1] - yy1[0])/(yy1.size-1))*(params[4])*(1/Misc_Routines.CoFactors.d2r) meanwx = (xwidth*ywidth) meanpos = numpy.array([xcenter,ycenter]) #print 'Position: %f %f' %(xcenter,ycenter) #print 'Widths: %f %f' %(xwidth, ywidth) #print 'BWHP: %f' %(2*numpy.sqrt(2*meanwx)*numpy.sqrt(-numpy.log(0.5))) self.meanpos = meanpos class overJroShow: __serverdocspath = '' __tmpDir = '' def __init__(self, title='', heights=None): self.year = None self.month = None self.dom = None self.pattern = None self.maxphi = None self.heights = None self.filename = None self.showType = None self.path = None self.objects = None self.nptsx = 101 self.nptsy = 101 self.fftopt = 0 self.site = 1 self.dcosx = 1 self.dcosy = 1 self.dcosxrange = None self.dcosyrange = None self.maxha_min= 0. self.show_object = None self.dcosx_mag = None self.dcosy_mag = None self.ha_mag = None self.time_mag = None self.main_dec = None self.ObjC = None self.ptitle = title self.path4plotname = None self.plotname0 = None self.plotname1 = None self.plotname2 = None self.scriptHeaders = 0 self.glat = -11.95 self.glon = -76.8667 self.UT = 5 #timezone self.glat = -11.951481 self.glon = -76.874383 if heights is None: self.heights = numpy.array([100.,500.,1000.]) else: self.heights = numpy.array(heights) def initParameters1(self): gui=1 if self.pattern==None: if gui==1: self.filename = self.filename.split(',') pattern = numpy.atleast_1d(self.pattern) filename = numpy.atleast_1d(self.filename) npatterns = numpy.max(numpy.array([pattern.size,filename.size])) self.pattern = numpy.resize(pattern,npatterns) self.filename = numpy.resize(filename,npatterns) self.doy = datetime.datetime(self.year,self.month,self.dom).timetuple().tm_yday if self.objects==None: self.objects=numpy.zeros(5) else: tmp = numpy.atleast_1d(self.objects) self.objects = numpy.zeros(5) self.objects[0:tmp.size] = tmp self.show_object = self.objects self.maxha_min = 4*self.maxphi*numpy.sqrt(2)*1.25 #ROJ geographic coordinates and time zone self.UT = 5 #timezone self.glat = -11.951481 self.glon = -76.874383 self.junkjd = TimeTools.Time(self.year,self.month,self.dom).change2julday() self.junklst = TimeTools.Julian(self.junkjd).change2lst(longitude=self.glon) # Finding RA of observatory for a specific date self.ra_obs = self.junklst*Misc_Routines.CoFactors.h2d def initParameters(self, dt): self.year = dt.year self.month = dt.month self.dom = dt.day # Defining plot filenames self.path4plotname = os.path.join(self.__serverdocspath,self.__tmpDir) self.plotname0 = 'over_jro_0_%i.png'% (time.time()) #plot pattern & objects self.plotname1 = 'over_jro_1_%i.png'% (time.time()) #plot antenna cuts self.plotname2 = 'over_jro_2_%i.png'% (time.time()) #plot sky noise # Defining antenna axes respect to geographic coordinates (See Ochs report). # alfa = 1.46*Misc_Routines.CoFactors.d2r # theta = 51.01*Misc_Routines.CoFactors.d2r alfa = 1.488312*Misc_Routines.CoFactors.d2r th = 6.166710 + 45.0 theta = th*Misc_Routines.CoFactors.d2r self.maxha_min = 4*self.maxphi*numpy.sqrt(2)*1.25 sina = numpy.sin(alfa) cosa = numpy.cos(alfa) MT1 = numpy.array([[1,0,0],[0,cosa,-sina],[0,sina,cosa]]) sinb = numpy.sin(theta) cosb = numpy.cos(theta) MT2 = numpy.array([[cosb,sinb,0],[-sinb,cosb,0],[0,0,1]]) self.MT3 = numpy.array(numpy.dot(MT2, MT1)).transpose() self.xg = numpy.dot(self.MT3.transpose(),numpy.array([1,0,0])) self.yg = numpy.dot(self.MT3.transpose(),numpy.array([0,1,0])) self.zg = numpy.dot(self.MT3.transpose(),numpy.array([0,0,1])) def plotPattern2(self, date, phases, gain_tx, gain_rx, ues, just_rx, angle): # Plotting Antenna patterns. self.maxphi = angle self.initParameters(date) self.doy = datetime.datetime(date.year,date.month,date.day).timetuple().tm_yday self.junkjd = TimeTools.Time(self.year,self.month,self.dom).change2julday() self.junklst = TimeTools.Julian(self.junkjd).change2lst(longitude=self.glon) self.ra_obs = self.junklst*Misc_Routines.CoFactors.h2d date = TimeTools.Time(date.year,date.month,date.day).change2strdate(mode=2) mesg = 'Over Jicamarca: ' + date[0] ObjAnt = JroPattern(pattern=0, filename=None, path=None, nptsx=self.nptsx, nptsy=self.nptsy, maxphi=angle, fftopt=self.fftopt, phases=phases, gain_tx=gain_tx, gain_rx=gain_rx, ues=ues, just_rx=just_rx ) self.pattern_plot = AntPatternPlot() self.pattern_plot.contPattern(iplot=0, gpath=self.path4plotname, filename=self.plotname0, mesg=mesg, amp=ObjAnt.norpattern, x=ObjAnt.dcosx, y=ObjAnt.dcosy, getCut=ObjAnt.getcut, title=self.ptitle, save=False) self.pattern_plot.plotRaDec(gpath=self.path4plotname, filename=self.plotname0, jd=self.junkjd, ra_obs=self.ra_obs, xg=self.xg, yg=self.yg, x=ObjAnt.dcosx, y=ObjAnt.dcosy, save=False) vect_ant = numpy.array([ObjAnt.meanpos[0],ObjAnt.meanpos[1],numpy.sqrt(1-numpy.sum(ObjAnt.meanpos**2.))]) vect_geo = numpy.dot(scipy.linalg.inv(self.MT3),vect_ant) vect_polar = Misc_Routines.Vector(numpy.array(vect_geo),direction=1).Polar2Rect() [ra,dec,ha] = Astro_Coords.AltAz(vect_polar[1],vect_polar[0],self.junkjd).change2equatorial() self.main_dec = dec def plotBfield(self): ObjB = BField(self.year,self.doy,1,self.heights) [dcos, alpha, nlon, nlat] = ObjB.getBField() self.pattern_plot.plotBField('', '',dcos,alpha,nlon,nlat, self.dcosxrange, self.dcosyrange, ObjB.heights, ObjB.alpha_i, save=False) def plotCelestial(self, objects): ntod = 24.*16. tod = numpy.arange(ntod)/ntod*24. [month,dom] = TimeTools.Doy2Date(self.year, self.doy).change2date() jd = TimeTools.Time(self.year, month, dom, tod+self.UT).change2julday() self.pattern_plot.plotCelestial( jd, self.main_dec, tod, self.maxha_min, objects, self.glat, self.glon, self.xg, self.yg, self.dcosxrange, self.dcosyrange ) def plotAntennaCuts(self): # print "Drawing antenna cuts" incha = 0.05 # min nha = numpy.int32(2*self.maxha_min/incha) + 1. newha = numpy.arange(nha)/nha*2.*self.maxha_min - self.maxha_min nha_star = numpy.int32(200./incha) star_ha = (numpy.arange(nha_star) - (nha_star/2))*nha_star #Init ObjCut for PatternCutPlot() view_objects = numpy.where(self.show_object>0) subplots = len(view_objects[0]) ObjCut = Graphics_OverJro.PatternCutPlot(subplots) for io in (numpy.arange(5)): if self.show_object[io]==2: if io==0: if self.dcosx_mag.size!=0: dcosx = self.dcosx_mag dcosy = self.dcosy_mag dcosz = 1 - numpy.sqrt(dcosx**2. + dcosy**2.) # Finding rotation of B respec to antenna coords. [mm,bb] = scipy.polyfit(dcosx,dcosy,1) alfa = 0.0 theta = -1.*numpy.arctan(mm) sina = numpy.sin(alfa); cosa = numpy.cos(alfa) MT1 = [[1,0,0],[0,cosa,-sina],[0,sina,cosa]] MT1 = numpy.array(MT1) sinb = numpy.sin(theta); cosb = numpy.cos(theta) MT2 = [[cosb,sinb,0],[-sinb,cosb,0],[0,0,1]] MT2 = numpy.array(MT2) MT3_mag = numpy.dot(MT2, MT1) MT3_mag = numpy.array(MT3_mag).transpose() # Getting dcos respec to B coords vector = numpy.array([dcosx,dcosy,dcosz]) nvector = numpy.dot(MT3_mag,vector) nvector = numpy.array(nvector).transpose() ## print 'Rotation (deg) %f'%(theta/Misc_Routines.CoFactors.d2r) yoffset = numpy.sum(nvector[:,1])/nvector[:,1].size # print 'Dcosyoffset %f'%(yoffset) ha = self.ha_mag*4. time = self.time_mag width_star = 0.1 # half width in minutes otitle = 'B Perp. cut' # else: # print "No B perp. over Observatory" # # elif io==1: if self.ObjC.dcosx_sun.size!=0: dcosx = self.ObjC.dcosx_sun dcosy = self.ObjC.dcosy_sun ha = self.ObjC.ha_sun*4.0 time = self.ObjC.time_sun width_star = 2. # half width in minutes otitle = 'Sun cut' # else: # print "Sun is not passing over Observatory" elif io==2: if self.ObjC.dcosx_moon.size!=0: dcosx = self.ObjC.dcosx_moon dcosy = self.ObjC.dcosy_moon ha = self.ObjC.ha_moon*4 time = self.ObjC.time_moon m_distance = 404114.6 # distance to the Earth in km m_diameter = 1734.4 # diameter in km. width_star = numpy.arctan(m_distance/m_diameter) width_star = width_star/2./Misc_Routines.CoFactors.d2r*4. otitle = 'Moon cut' # else: # print "Moon is not passing over Observatory" elif io==3: if self.ObjC.dcosx_hydra.size!=0: dcosx = self.ObjC.dcosx_hydra dcosy = self.ObjC.dcosy_hydra ha = self.ObjC.ha_hydra*4. time = self.ObjC.time_hydra width_star = 0.25 # half width in minutes otitle = 'Hydra cut' # else: # print "Hydra is not passing over Observatory" elif io==4: if self.ObjC.dcosx_galaxy.size!=0: dcosx = self.ObjC.dcosx_galaxy dcosy = self.ObjC.dcosy_galaxy ha = self.ObjC.ha_galaxy*4. time = self.ObjC.time_galaxy width_star = 25. # half width in minutes otitle = 'Galaxy cut' # else: # print "Galaxy center is not passing over Jicamarca" # # hour = numpy.int32(time) mins = numpy.int32((time - hour)*60.) secs = numpy.int32(((time - hour)*60. - mins)*60.) ObjT = TimeTools.Time(self.year,self.month,self.dom,hour,mins,secs) subtitle = ObjT.change2strdate() star_cut = numpy.exp(-(star_ha/width_star)**2./2.) pol = scipy.polyfit(ha,dcosx,3.) polx = numpy.poly1d(pol); newdcosx = polx(newha) pol = scipy.polyfit(ha,dcosy,3.) poly = numpy.poly1d(pol);newdcosy = poly(newha) patterns = [] for icut in numpy.arange(self.pattern.size): # Getting Antenna cut. Obj = JroPattern(dcosx=newdcosx, dcosy=newdcosy, getcut=1, pattern=self.pattern[icut], path=self.path, filename=self.filename[icut]) Obj.getPattern() patterns.append(Obj.pattern) ObjCut.drawCut(io, patterns, self.pattern.size, newha, otitle, subtitle, self.ptitle) ObjCut.saveFig(self.path4plotname,self.plotname1) def skyNoise(jd, ut=-5.0, longitude=-76.87, filename='/app/utils/galaxy.txt'): """ hydrapos returns RA and Dec provided by Bill Coles (Oct 2003). Parameters ---------- jd = The julian date of the day (and time), scalar or vector. dec_cut = A scalar giving the declination to get a cut of the skynoise over Jica- marca. The default value is -11.95. filename = A string to specify name the skynoise file. The default value is skynoi- se_jro.dat Return ------ maxra = The maximum right ascension to the declination used to get a cut. ra = The right ascension. Examples -------- >> [maxra,ra] = skynoise_jro() >> print maxra, ra 139.45 -12.0833333333 Modification history -------------------- Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. """ # Defining date to compute SkyNoise. [year, month, dom, hour, mis, secs] = TimeTools.Julian(jd).change2time() is_dom = (month==9) & (dom==21) if is_dom: tmp = jd jd = TimeTools.Time(year,9,22).change2julian() dom = 22 # Reading SkyNoise g = os.getcwd() f = open(filename) lines = f.read() f.close() nlines = 99 lines = lines.split('\n') data = numpy.zeros((2,nlines))*numpy.float32(0.) for ii in numpy.arange(nlines): line = numpy.array([lines[ii][0:6],lines[ii][6:]]) data[:,ii] = numpy.float32(line) # Getting SkyNoise to the date desired. otime = data[0,:]*60.0 opowr = data[1,:] hour = numpy.array([0,23]) mins = numpy.array([0,59]) secs = numpy.array([0,59]) LTrange = TimeTools.Time(year,month,dom,hour,mins,secs).change2julday() LTtime = LTrange[0] + numpy.arange(1440)*((LTrange[1] - LTrange[0])/(1440.-1)) lst = TimeTools.Julian(LTtime + (-3600.*ut/86400.)).change2lst() ipowr = lst*0.0 # Interpolating using scipy (inside max and min "x") otime = otime/3600. val = numpy.where((lst>numpy.min(otime)) & (lstnumpy.max(otime)) if uval[0].size>0: ii = numpy.min(uval[0]) m = (ipowr[ii-1] - ipowr[ii-2])/(lst[ii-1] - lst[ii-2]) b = ipowr[ii-1] - m*lst[ii-1] ipowr[uval] = m*lst[uval] + b if is_dom: lst = numpy.roll(lst,4) ipowr = numpy.roll(ipowr,4) new_lst = numpy.int32(lst*3600.) new_pow = ipowr return ipowr, LTtime, lst def skynoise_plot(year, month, day): """ getPlot method creates a skynoise map over Jicamarca for a desired day. Additionally save a PNG file of this plot. Examples -------- >> SkyNoisePlot(skypower,skytime,skytime_lst).getPlot() Modification history -------------------- Created by Freddy R. Galindo, ROJ, 18 October 2009. """ # Working with the time before to plot the SkyNoise julian = TimeTools.Time(year, month, day).change2julday() power, tm, time_lst = skyNoise(julian) secs = TimeTools.Julian(tm).change2secs() secs_lst = time_lst*1. if time_lst.argmin()>0: secs_lst[time_lst.argmin():] = secs_lst[time_lst.argmin():] + 24. secs_lst = secs_lst*3600. label_secs = time.localtime(secs[power.argmax()] + time.timezone) label_secs_lst = time.localtime(secs_lst[power.argmax()] + time.timezone) # Creating labels for main x-labelticks (Local Time): snow = TimeTools.Time(year, month, day, 0, 0, 0).change2secs() today = secs - snow xticks_dpos = [] xticks_dval = [] for ix in [0,120,240,360,480,600,720,840,960,1080,1200,1320,1439]: xticks_dpos.append(today[ix]) time_tuple = time.localtime(today[ix] + time.timezone) xticks_dval.append(time.strftime('%H:%M',time_tuple)) if ix==1439:xticks_dval[12] = '' # Creating labels for secondary x-labelticks (Sidereal Time): xticks_upos = [secs_lst[0],secs_lst[-1]] xticks_uval = ['',''] for ix in [0,120,240,360,480,600,720,840,960,1080,1200,1320]: index = numpy.argmin(numpy.abs(time_lst - (ix/60. + 1.))) xticks_upos.append(secs_lst[index]) time_tuple = time.localtime(secs_lst[index] + time.timezone) if time_tuple[4]==59: tmp_list = list(time_tuple) tmp_list[4] = 0 tmp_list[3] = tmp_list[3] + 1 time_tuple = tuple(tmp_list) xticks_uval.append(time.strftime('%H:%M',time_tuple)) # Creating SkyNoise Map. fig = Figure() fig.subplots_adjust(top=0.80) main_mesg = "SKY BRIGHTNESS AT 50Mhz - Date: " main_mesg = main_mesg + time.strftime("%d-%b-%Y (%j)",label_secs) main_mesg = main_mesg + "\nGalaxy Pass at " + time.strftime("%H:%M:%S LT",label_secs) main_mesg = main_mesg + time.strftime(" (%H:%M:%S LST)",label_secs_lst) fig.suptitle(main_mesg, fontsize=12) axl = fig.add_subplot(1,1,1) axl.plot(today,power,'b-') axr = axl.twinx() axl.set_xlabel('Local Time (Hours)') axl.set_ylabel('Signal Strengh (mA)') axr.set_ylabel('Temperature [K]x10^3') axl.set_xlim(0,24) axl.set_ylim(0,115) axr.set_ylim(0,50) axl.set_xticks(xticks_dpos) axl.set_xticklabels(xticks_dval, fontsize=8) axl.grid() axt = axl.twiny() axt.set_xlabel('Local Sidereal Time (Hours)') axt.set_xlim(secs_lst[0],secs_lst[-1]) axt.set_xticks(xticks_upos) axt.set_xticklabels(xticks_uval, fontsize=8) buf = BytesIO() fig.savefig(buf, format="png") return buf def overjro_plot(pattern_id, date, angle, height, bodys): pattern = select_pattern(pattern = pattern_id) ob = overJroShow(pattern['title'], heights=height) ob.plotPattern2( date, pattern['phase'], pattern['gaintx'], pattern['gainrx'], pattern['ues'], pattern['justrx'], angle ) if 'bfield' in bodys: ob.plotBfield() objects = [] for body in bodys: if body=='sun': objects.append(0) elif body=='moon': objects.append(1) elif body=='hydra': objects.append(2) elif body=='galaxy': objects.append(3) if objects: ob.plotCelestial(objects) buf = BytesIO() ob.pattern_plot.fig.savefig(buf, format="png") return buf