""" Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009. Added to signal Chain by Joab Apaza, ROJ, Jun 2023. """ import numpy import time import os from scipy.special import lpmn from schainpy.model.utils import Astro_Coords 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. Added to signal Chain by Joab Apaza, ROJ, Jun 2023. """ 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([-37,-37,37,37])): """ 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] heta = (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) elif self.site==2: #AMISR title_site = 'AMISR 14' coord_site = [-76.874913, -11.953371, 0.52984] theta = (0.0977)*numpy.pi/180. # 0.0977 delta = 0.110*numpy.pi/180 # 0.11 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 = 2.5 nlon = int(int(maglimits[2] - maglimits[0])/grid_res + 1) nlat = int(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 = int(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]= 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,26): 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. """ base_path = os.path.dirname(os.path.abspath(__file__)) filename = os.path.join(base_path, "igrf13coeffs.txt") period_v, g_v, h_v = self.__readIGRFfile(filename) g2 = numpy.zeros((14,14,26)) h2 = numpy.zeros((14,14,26)) g2[1:14,:,:] = g_v h2[1:14,:,:] = h_v g = numpy.transpose(g2, (2,0,1)) h = numpy.transpose(h2, (2,0,1)) periods = period_v.copy() return periods, g, h def rotvector(self,vector,axis=1,ang=0): """ rotvector function returns the new vector generated rotating the rectagular coords. Parameters ---------- vector = A lineal 3-elements array (x,y,z). axis = A integer to specify the axis used to rotate the coord systems. The default value is 1. axis = 1 -> 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