##// END OF EJS Templates
schain
RSS
merge with amisr tmp
merge with amisr tmp
joabAM - - Load All Authors

File last commit:

r1738:e8767bd42c24
r1738:e8767bd42c24
Show More
BField.py
507 lines | 18.0 KiB | text/x-python | PythonLexer
/ schainpy / model / utils / BField.py
History | Source | Raw |Copy content |Copy permalink
joabAM
merge with amisr tmp
 r1738 """
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 - 1900) % 5
ntime = int((time0 - 1900 - dtime)/5)
else:
# Estimating coefficients for times > 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