Module ephemeris_library.rotation_libraries.lib_epop_geopack
Expand source code
import os
import sys
import datetime
import numpy as np
import pandas as pd
# For consistancy with legacy code & values, all interfaces & time-values use UT/UTC
# init_igrf() is no longer called on import, as to allow for pathing options.
# The function is called with it's default path when recalc is called.
# Please note that all function now accept NP-columns, rather than singleton values.
# This allows for far faster processing.
########
## CONFIG:
DEFAULT_IGRF_PATH = os.path.join( os.path.dirname( __file__ ), "igrf13coeffs.txt" ) # Default IGRF file deployed with the Frame_Rotation library.
########
### GLOBALS
# Pre-instantiate all global variables as to be able to check if instantion has occurred.
# Variables loaded by load_igrf
igrf = None # A dictionary ['g','h'] of dims 105x26 each, containing the 'g'/'h' entries from the IGRF constants file.
nmn = None # The number of arguments to each g/h entry.
mns = None # A incrementing 105x2 of triangluar-pair vales of the n/m columns for the 'g/h' entries in the IGRF constants file.
nyear = None # Integer. Number of years that the IGRF constants file covers
years = None # List of Integers, the years that are covered.
yruts = None # The Unix-time entries for the midgnight start of each year in 'years'
# The coefficients used for the actual frame translation. Set up by recalc().
# A small amount of documentation for these can be found in the recalc() function.
st0,ct0,sl0,cl0,ctcl,stcl, \
ctsl,stsl,sfi,cfi,sps,cps, \
shi,chi,hi,psi,xmut,ds3,cgst,sgst,ba, \
a11,a21,a31,a12,a22,a32,a13,a23,a33, \
e11,e21,e31,e12,e22,e32,e13,e23,e33 = [None] * 39
# The IGRF coefficients.
# Loaded by load_igrf(), used by any magnetic-field derivation (IGRF, MAG frame)
g, h, rec = [None] * 3
###############
## Supplementary Functions:
def _Sanitize_Input( recalc_dependant, *Value_Columns ):
# Given input columns, try to sanitize them into 1d-arrays.
# Throw exceptions when needed.
# Always returns a list, so will require indexing for single columns.
if len(Value_Columns) == 0:
raise ValueError("_Sanitize_Input given an empty list of data to sanitize!")
parsed_Columns = []
# First make sure everything is in a 1d.
# Happily recasts lists, np-arrays, and singleton values.
for i in range(len( Value_Columns )):
parsed_Columns.append( np.atleast_1d( Value_Columns[i] ) )
# And check that all given lists are of 1d, and of the same length.
First_Column_Shape = parsed_Columns[0].shape
for Column in parsed_Columns:
if Column.ndim != 1:
raise ValueError("epopGeopack function inputs must be either singleton or 1-d columns!")
if Column.shape != First_Column_Shape:
raise ValueError("epopGeopack function inputs must be of equal length!")
if recalc_dependant:
if st0 is None:
raise ValueError("This function requires that recalc() be called first!")
elif st0.shape != Column.shape:
raise ValueError("This function requires that recalc() be called with matching column-length!")
return parsed_Columns
###############
## Primary functions
def init_igrf( IGRF_Filepath = DEFAULT_IGRF_PATH ):
# Initialize the IGRF coefficients and related coefs.
global igrf, nmn,mns, nyear,years,yruts
# Read the IGRF file into a pandas dataframe
IGRF_df = pd.read_csv( IGRF_Filepath, sep='\s+', skiprows=3)
IGRF_np = IGRF_df.to_numpy()
Year_Str = IGRF_df.columns[3:] # Extract the year list from the column name, then parse to integers.
Year_Int = np.array([Year[:4] for Year in Year_Str], dtype=int) # The [:4] scrubs the float, and the final entry's year range
Year_Int[-1] += 5 # And then add the '+5' range s.t. the "2020-5" == 2025
# And finally write that to the global var.
years = Year_Int
nyear = len(Year_Int)
# From the final line of the IGRF file, extract the order of the polynomial, and the number of arguments used to generate it.
k = IGRF_np[-1,1] + 1 # k = 14 at time of writing.
nmn = (k+1) * k // 2 # nmn= 105
igrf = np.zeros(( nmn, nyear, 2), dtype=float)
mns = np.zeros(( nmn, 2), dtype=float)
l = 0
for i in range(k):
for j in range(i + 1):
mns[l, :] = [i, j]
l += 1
# For each row in the file:
for Row in IGRF_np:
if Row[0] == 'g': i = 0
else: i = 1
n, m = Row[1:3].astype(int)
mn = n * (n + 1) // 2 + m # The row-index of the related pair in the 'mns' array.
igrf[mn, :, i] = Row[3:]
# I have no idea what this is supposed to accomplish.
igrf[:, -1, :] = igrf[:, -2, :] + igrf[:, -1, :] * 5
# And then build the 'year-ut-second' variable, the number of seconds since the unix epoch for each year-column.
yruts = np.zeros(nyear, dtype=float)
Unix_Epoch = datetime.datetime(1970,1,1)
for i in range(nyear):
yruts[i] = (datetime.datetime(years[i],1,1) - Unix_Epoch).total_seconds()
# Translate the 'igrf' var into a dictionary of 'g' & 'h', of 2d arrays of 105x26 (default).
igrf = {'g': igrf[:,:,0], 'h': igrf[:,:,1]}
def load_igrf( ut ):
if igrf is None: init_igrf()
ut = _Sanitize_Input( False, ut )[0] # Make sure the array is at least 1d.
# Seek the list of indecies s.t. can use the values on either side to extrapolate between.
Indecies = np.searchsorted( yruts, ut, side='right' ) - 1
# In the event that values are inserted -after- the end of the
# given constants, extrapolate using the [-2] to [-1] constants.
Indecies[Indecies < 0 ] = 0
Indecies[Indecies >= len( yruts )-1 ] = len( yruts ) - 2
g0 = igrf['g'][:, Indecies ]
h0 = igrf['h'][:, Indecies ]
g1 = igrf['g'][:,Indecies+1]
h1 = igrf['h'][:,Indecies+1]
# Calculate the time-fractions to use.
Interp_Start = yruts[Indecies]
Interp_End = yruts[Indecies+1]
Weight_Right = (ut - Interp_Start) / (Interp_End - Interp_Start) # The time-fraction of the right-side components
Weight_Left = 1 - Weight_Right # The time-fraction of the left-side components
g_Estimate = g0 * Weight_Left + g1 * Weight_Right
h_Estimate = h0 * Weight_Left + h1 * Weight_Right
# Transpose s.t. they're row-indexed (Eg: g[0] rather than g[:,0])
g_Estimate, h_Estimate = g_Estimate.T, h_Estimate.T
return g_Estimate, h_Estimate
def recalc( ut, vxgse=-400,vygse=0,vzgse=0 ):
# Make sure the 'ut' timestamps are given as a parsable singleton or list.
ut = _Sanitize_Input( False, ut )[0]
if igrf is None: init_igrf()
# Documentation from Fortran's code.
# 1. PREPARES ELEMENTS OF ROTATION MATRICES FOR TRANSFORMATIONS OF VECTORS BETWEEN
# SEVERAL COORDINATE SYSTEMS, MOST FREQUENTLY USED IN SPACE PHYSICS.
#
# 2. PREPARES COEFFICIENTS USED IN THE CALCULATION OF THE MAIN GEOMAGNETIC FIELD
# (IGRF MODEL)
#
#-----INPUT PARAMETERS:
# VGSEX,VGSEY,VGSEZ - GSE (GEOCENTRIC SOLAR-ECLIPTIC) COMPONENTS OF THE OBSERVED
# SOLAR WIND FLOW VELOCITY (IN KM/S)
# st0/ct0 - sin/cos of teta0 (colat in geo?). sl0/cl0 - sin/cos of lambda0 (longitude in geo?).
# ctcl/stcl - . ctsl/stsl - . geo x and z?
# sfi/cfi/xmut - rotate angle between mag and sm and its sin/cos.
# sps/cps/psi - tilt angle and its sin/cos.
# shi/chi/hi - rotate angle between gse to gsm and its sin/cos.
# a[11,...,33] - matrix converts geo to gsm.
# cgst/sgst - cos/sin of gst.
# ds3.
# ba(6).
global st0,ct0,sl0,cl0,ctcl,stcl,ctsl,stsl,sfi,cfi,sps,cps, \
shi,chi,hi,psi,xmut,ds3,cgst,sgst,ba, \
a11,a21,a31,a12,a22,a32,a13,a23,a33, \
e11,e21,e31,e12,e22,e32,e13,e23,e33
# The common block /geopack2/ contains coefficients of the IGRF field model, calculated
# for a given year and day from their standard epoch values. the array rec contains
# coefficients used in the recursion relations for legendre associate polynomials.
# common /geopack2/ g(105),h(105),rec(105)
# Compute the m,n related coefficients (following the notation in Davis 2004):
# 1. The Schmidt quasi-normalization: S_n,m, which normalizes the associated Legendre polynomials P_n^m
# to the Guassian normalized associated Legendre polynomials P^n,m.
#
# Since g_n^m * P_n^m should be constant, mutiplying S_n,m to g_n^m is equivalently converting P_n^m to P^n,m.
# The benefit of doing so is that P^n,m follows simple recursive form, c.f. Eq (16 a-c) and Eq(17 a-b).
# 2. rec[mn], which used in the recursion relation.
global g,h,rec # Used in the igrf_geo() & dip()
# Build the 14-degree order Legendre polynomial coefficients
rec = np.array([])
for n in range( 14 ):
m = np.arange( n+1 )
rec = np.append(rec, (n-m)*(n+m)/((2*n+1)*(2*n-1)) )
# WH- The 'rec' array is a 105-len array of floats with increasing subarrays of [Val*,0], with the Val-subarray starting with a zero-length entry
# Eg: [ 0, 1/3,0, 4/15,1/5,0,....]
# Fetch the IGRF coefficients for the given timestamps
g, h = load_igrf( ut ) # Returns two lists of shape Nx105, where N is the input 'ut'-len.
# Multiply the g/h IGRF values by the 'Schmidt Normalization Factors'
s = 1
_Factor = [0.]
for n in range(1, 14):
s *= (2*n - 1) / n
_Factor.append( s )
p = s
for m in range(1, n+1):
if m == 1: aa = 2
else: aa = 1
p *= np.sqrt( aa*(n-m+1) / (n+m))
_Factor.append( p )
# Then multiply the g/h coefficients by the array of '_Factor'
g *= _Factor
h *= _Factor
g10=-g[:,1]
g11=-g[:,2]
h11=-h[:,2]
# Now calculate the components of the unit vector EzMAG in the GEO system, 'parallel to geodipole axis':
# sin(teta0)*cos(lambda0), sin(teta0)*sin(lambda0), and cos(teta0)
# st0 * cl0 st0 * sl0 ct0
sq = g11**2 + h11**2
sqq = np.sqrt(sq)
sqr = np.sqrt(g10**2+sq)
sl0 = h11/sqq
cl0 = g11/sqq
st0 = sqq/sqr
ct0 = g10/sqr
stcl = st0*cl0
stsl = st0*sl0
ctsl = ct0*sl0
ctcl = ct0*cl0
# Fetch the solar data at each saught time.
# slong never used.
gst, slong, srasn, sdec, obliq = sun(ut)
# Each component returned is an array of floats, len(ut)
# All vectors are expressed in GEI.
# xgse_[xyz] (s[123]) are the components of the unit vector exgsm=exgse in GEI,
# pointing from the earth's center to the sun:
xgse_x = np.cos(srasn) * np.cos(sdec)
xgse_y = np.sin(srasn) * np.cos(sdec)
xgse_z = np.sin(sdec)
# zgse_[xyz] in GEI has the components (0,-sin(delta),cos(delta)) = (0.,-0.397823,0.917462);
# Here delta = 23.44214 deg for the epoch 1978 (see the book by Gurevich or other astronomical handbooks).
# Here the most accurate time-dependent formula is used:
zgse_x = 0.
zgse_y = -np.sin(obliq) #HERE
zgse_z = np.cos(obliq)
# ygse_[xyz] (dy[123]) = zgse_[xyz] x xgsm_[xyz] in GEI:
ygse_x = zgse_y * xgse_z - zgse_z * xgse_y
ygse_y = zgse_z * xgse_x - zgse_x * xgse_z
ygse_z = zgse_x * xgse_y - zgse_y * xgse_x
# zsm_[xyz] (dip[123]) are the components of the unit vector zsm=zmag in GEI:
cgst = np.cos(gst)
sgst = np.sin(gst)
zsm_x = stcl * cgst - stsl * sgst
zsm_y = stcl * sgst + stsl * cgst
zsm_z = ct0
# xgsw_[xyz] (x[123]) in GEI.
v1 = -1 / np.sqrt(vxgse*vxgse+vygse*vygse+vzgse*vzgse)
xgsw_x = (vxgse*xgse_x + vygse*ygse_x + vzgse*zgse_x)*v1
xgsw_y = (vxgse*xgse_y + vygse*ygse_y + vzgse*zgse_y)*v1
xgsw_z = (vxgse*xgse_z + vygse*ygse_z + vzgse*zgse_z)*v1
# ygsw (y[123]) = zsm x xgsw in GEI.
ygsw_x = zsm_y*xgsw_z - zsm_z*xgsw_y
ygsw_y = zsm_z*xgsw_x - zsm_x*xgsw_z
ygsw_z = zsm_x*xgsw_y - zsm_y*xgsw_x
y = np.sqrt(ygsw_x*ygsw_x+ygsw_y*ygsw_y+ygsw_z*ygsw_z)
ygsw_x = ygsw_x/y
ygsw_y = ygsw_y/y
ygsw_z = ygsw_z/y
# zgsw (z[123]) = xgsw x ygsw in GEI.
zgsw_x = xgsw_y*ygsw_z-xgsw_z*ygsw_y
zgsw_y = xgsw_z*ygsw_x-xgsw_x*ygsw_z
zgsw_z = xgsw_x*ygsw_y-xgsw_y*ygsw_x
# xgsm = xgse in GEI.
xgsm_x,xgsm_y,xgsm_z = xgse_x,xgse_y,xgse_z
# ygsm = zsm x xgsm in GEI.
ygsm_x = zsm_y*xgsm_z - zsm_z*xgsm_y
ygsm_y = zsm_z*xgsm_x - zsm_x*xgsm_z
ygsm_z = zsm_x*xgsm_y - zsm_y*xgsm_x
y=np.sqrt(ygsm_x*ygsm_x+ygsm_y*ygsm_y+ygsm_z*ygsm_z)
ygsm_x = ygsm_x/y
ygsm_y = ygsm_y/y
ygsm_z = ygsm_z/y
# ezgsm = exgsm x eygsm in GEI.
zgsm_x = xgse_y*ygsm_z-xgse_z*ygsm_y
zgsm_y = xgse_z*ygsm_x-xgse_x*ygsm_z
zgsm_z = xgse_x*ygsm_y-xgse_y*ygsm_x
# The elements of the matrix gse to gsm are the scalar products:
# chi=em22=(eygsm,eygse), shi=em23=(eygsm,ezgse), em32=(ezgsm,eygse)=-em23, and em33=(ezgsm,ezgse)=em22
chi = ygsm_x*ygse_x + ygsm_y*ygse_y + ygsm_z*ygse_z
shi = ygsm_x*zgse_x + ygsm_y*zgse_y + ygsm_z*zgse_z
hi = np.arcsin(shi)
# elements of the matrix gsm to gsw are the scalar products:
# e11 = (exgsm, exgsw) e12 = (exgsm, eygsw) e13 = (exgsm, ezgsw)
# e21 = (eygsm, exgsw) e22 = (eygsm, eygsw) e23 = (eygsm, ezgsw)
# e31 = (ezgsm, exgsw) e32 = (ezgsm, eygsw) e33 = (ezgsm, ezgsw)
e11 = xgsm_x*xgsw_x + xgsm_y*xgsw_y + xgsm_z*xgsw_z
e12 = xgsm_x*ygsw_x + xgsm_y*ygsw_y + xgsm_z*ygsw_z
e13 = xgsm_x*zgsw_x + xgsm_y*zgsw_y + xgsm_z*zgsw_z
e21 = ygsm_x*xgsw_x + ygsm_y*xgsw_y + ygsm_z*xgsw_z
e22 = ygsm_x*ygsw_x + ygsm_y*ygsw_y + ygsm_z*ygsw_z
e23 = ygsm_x*zgsw_x + ygsm_y*zgsw_y + ygsm_z*zgsw_z
e31 = zgsm_x*xgsw_x + zgsm_y*xgsw_y + zgsm_z*xgsw_z
e32 = zgsm_x*ygsw_x + zgsm_y*ygsw_y + zgsm_z*ygsw_z
e33 = zgsm_x*zgsw_x + zgsm_y*zgsw_y + zgsm_z*zgsw_z
# Solar magnetic (sm) to Geocentric Solar Magnetospheric (gsm) rotation.
# Tilt angle: psi=arcsin(ezsm . exgsm)
sps = zsm_x*xgse_x + zsm_y*xgse_y + zsm_z*xgse_z
cps = np.sqrt(1 - sps**2)
psi = np.arcsin(sps) # Not used elsewhere, but kept for consistency.
# The elements of the matrix mag to sm are the scalar products:
# cfi=gm22=(eysm,eymag), sfi=gm23=(eysm,exmag); They can be derived as follows:
# In geo the vectors exmag and eymag have the components (ct0*cl0,ct0*sl0,-st0) and (-sl0,cl0,0), respectively.
# Hence, in gei the components are:
# exmag: ct0*cl0*cos(gst)-ct0*sl0*sin(gst)
# ct0*cl0*sin(gst)+ct0*sl0*cos(gst)
# -st0
# eymag: -sl0*cos(gst)-cl0*sin(gst)
# -sl0*sin(gst)+cl0*cos(gst)
# 0
# The components of eysm in gei were found above as ysm_in_geix, ysm_in_geiy, and ysm_in_geiz;
# Now we only have to combine the quantities into scalar products:
xmag_x = ct0*(cl0*cgst - sl0*sgst)
xmag_y = ct0*(cl0*sgst + sl0*cgst)
xmag_z = -st0
ymag_x = -(sl0*cgst + cl0*sgst)
ymag_y = -(sl0*sgst - cl0*cgst)
cfi = ygsm_x*ymag_x + ygsm_y*ymag_y
sfi = ygsm_x*xmag_x + ygsm_y*xmag_y + ygsm_z*xmag_z
xmut=(np.arctan2(sfi,cfi)+3.1415926536)*3.8197186342
# The elements of the matrix geo to gsm are the scalar products:
# a11=(exgeo,exgsm), a12=(eygeo,exgsm), a13=(ezgeo,exgsm),
# a21=(exgeo,eygsm), a22=(eygeo,eygsm), a23=(ezgeo,eygsm),
# a31=(exgeo,ezgsm), a32=(eygeo,ezgsm), a33=(ezgeo,ezgsm),
# All the unit vectors in brackets are already defined in gei:
# xgeo=(cgst,sgst,0), ygeo=(-sgst,cgst,0), zgeo=(0,0,1)
# and therefore:
a11 = xgsm_x*cgst + xgse_y*sgst
a12 = -xgsm_x*sgst + xgse_y*cgst
a13 = xgsm_z
a21 = ygsm_x*cgst + ygsm_y *sgst
a22 = -ygsm_x*sgst + ygsm_y *cgst
a23 = ygsm_z
a31 = zgsm_x*cgst + zgsm_y*sgst
a32 = -zgsm_x*sgst + zgsm_y*cgst
a33 = zgsm_z
return psi
# Returns gst,slong,srasn,sdec,obliq used in recalc.
def sun( ut ):
"""
Calculates four quantities necessary for coordinate transformations
which depend on sun position (and, hence, on universal time and season)
Based on http://aa.usno.navy.mil/faq/docs/SunApprox.php and http://aa.usno.navy.mil/faq/docs/GAST.php
:param ut: ut sec, can be array.
:return: gst,slong,srasn,sdec. gst - greenwich mean sidereal time, slong - longitude along ecliptic
srasn - right ascension, sdec - declination of the sun (radians)
obliq - mean oblique of the ecliptic (radian)
Python version by Sheng Tian
"""
ut = np.atleast_1d(ut)
# Various constants.
twopi = 2*np.pi
jd2000 = 2451545.0 # Reduced Julian Date of the J2K epoch.
# convert to Julian date.
t0_jd = 2440587.5 # in day, 0 of Julian day.
secofdaygsw_x = 1./86400 # 1/sec of day.
t_jd = ut*secofdaygsw_x+t0_jd
# d = mjd - mj2000.
d = t_jd-jd2000
d = np.squeeze(d)
# mean obliquity of the ecliptic, e. (Rad)
# e = 23.439 - 0.00000036*d ; in degrees.
e = 0.4090877233749509 - 6.2831853e-9*d
# mean anomaly of the Sun, g.
# g = 357.529 + 0.98560028*d ; in degree.
g = 6.2400582213628066 + 0.0172019699945780*d
g = np.mod(g, twopi)
# mean longitude of the Sun, q.
# q = 280.459 + 0.98564736*d ; in degree.
q = 4.8949329668507771 + 0.0172027916955899*d
q = np.mod(q, twopi)
# geocentric apparent ecliptic longitude, l.
# l = q + 1.915 sin g + 0.020 sin 2g ; in degree.
l = q + 0.0334230551756914*np.sin(g) + 0.0003490658503989*np.sin(2*g)
# vl - q, mean longitude of the sun.
# vl = np.mod(279.696678+0.9856473354*dj,360.)/rad
# q = np.mod(4.881627937990388+0.01720279126623886*dj, twopi)
# g, mean anomaly of the sun.
# g = np.mod(358.475845+0.985600267*dj,360.)/rad
# g = np.mod(6.256583784118852+0.017201969767685215*dj, twopi)
# slong - l, geocentric apparent ecliptic longitude.
# slong = (vl + (1.91946-0.004789*t)*sin(g) + 0.020094*sin(2*g))/rad
# l = q+(0.03350089686033036-2.2884002156881157e-09*dj)*np.sin(g)+0.0003507064598957406*np.sin(2*g)
# l = np.mod(l, twopi)
# obliq - e, mean obliquity of the ecliptic.
# obliq = (23.45229-0.0130125*t)/rad
# e = 0.40931967763254096-6.217959450123535e-09*dj
# sin(d) = sin(e) * sin(L)
sind = np.sin(e)*np.sin(l)
sdec = np.arcsin(sind)
# tan(RA) = cos(e)*sin(L)/cos(L)
srasn = np.arctan2(np.cos(e)*np.sin(l), np.cos(l))
srasn = np.mod(srasn, twopi)
# http://aa.usno.navy.mil/faq/docs/GAST.php
# gst - gmst, greenwich mean sidereal time.
# gst = np.mod(279.690983+.9856473354*dj+360.*fday+180.,360.)/rad
# gst = np.mod(4.881528541489487+0.01720279126623886*dj+twopi*fday+np.pi, twopi)
# gmst = 18.697374558 + 24.06570982441908*d # in hour
gmst = 4.894961212735792 + 6.30038809898489*d # in rad
gmst = np.mod(gmst, twopi)
return gmst,l,srasn,sdec,e
#########################
# Simple frame tranlsations.
# Loads from the globally assigned variables.
def gswgsm( p1, p2, p3, j ):
"""
Converts gsm to gsw coordinates or vice versa.
j>0 j<0
input: xgsm,ygsm,zgsm xgsw,ygsw,zgsw
output: xgsw,ygsw,zgsw xgsm,ygsm,zgsm
:param p1,p2,p3: input position
:param j: flag
:return: output position
"""
if j > 0:
xgsw, ygsw, zgsw = p1, p2, p3
xgsm = xgsw*e11 + ygsw*e12 + zgsw*e13
ygsm = xgsw*e21 + ygsw*e22 + zgsw*e23
zgsm = xgsw*e31 + ygsw*e32 + zgsw*e33
return xgsm,ygsm,zgsm
else:
xgsm, ygsm, zgsm = p1, p2, p3
xgsw = xgsm*e11 + ygsm*e21 + zgsm*e31
ygsw = xgsm*e12 + ygsm*e22 + zgsm*e32
zgsw = xgsm*e13 + ygsm*e23 + zgsm*e33
return xgsw,ygsw,zgsw
def geomag( p1, p2, p3, j ):
"""
Converts geographic (geo) to dipole (mag) coordinates or vice versa.
j>0 j<0
input: j,xgeo,ygeo,zgeo j,xmag,ymag,zmag
output: xmag,ymag,zmag xgeo,ygeo,zgeo
:param p1,p2,p3: input position
:param j: flag
:return: output position
"""
if j > 0:
xgeo, ygeo, zgeo = p1, p2, p3
xmag = xgeo*ctcl + ygeo*ctsl - zgeo*st0
ymag = -xgeo*sl0 + ygeo*cl0
zmag = xgeo*stcl + ygeo*stsl + zgeo*ct0
return xmag, ymag, zmag
else:
xmag, ymag, zmag = p1, p2, p3
xgeo = xmag*ctcl - ymag*sl0 + zmag*stcl
ygeo = xmag*ctsl + ymag*cl0 + zmag*stsl
zgeo = -xmag*st0 + zmag*ct0
return xgeo, ygeo, zgeo
def geigeo( p1, p2, p3, j ):
"""
Converts equatorial inertial (gei) to geographical (geo) coords or vice versa.
j>0 j<0
input: xgei,ygei,zgei xgeo,ygeo,zgeo
output: xgeo,ygeo,zgeo xgei,ygei,zgei
:param p1,p2,p3: input position
:param j: flag
:return: output position
"""
# Uses the 'cgst' & 'sgst' globals set by recalc.
if j > 0:
xgei,ygei,zgei = p1, p2, p3
xgeo = xgei*cgst + ygei*sgst
ygeo = ygei*cgst - xgei*sgst
zgeo = zgei
return xgeo,ygeo,zgeo
else:
xgeo,ygeo,zgeo = p1, p2, p3
xgei = xgeo*cgst - ygeo*sgst
ygei = ygeo*cgst + xgeo*sgst
zgei = zgeo
return xgei,ygei,zgei
def magsm( p1, p2, p3, j ):
"""
Converts dipole (mag) to solar magnetic (sm) coordinates or vice versa
j>0 j<0
input: xmag,ymag,zmag xsm, ysm, zsm
output: xsm, ysm, zsm xmag,ymag,zmag
:param p1,p2,p3: input position
:param j: flag
:return: output position
"""
# Uses the 'sfi' & 'cfi' globals generated by recalc.
if j > 0:
xmag,ymag,zmag = p1, p2, p3
xsm = xmag*cfi - ymag*sfi
ysm = xmag*sfi + ymag*cfi
zsm = zmag
return xsm,ysm,zsm
else:
xsm,ysm,zsm = p1, p2, p3
xmag = xsm*cfi + ysm*sfi
ymag = ysm*cfi - xsm*sfi
zmag = zsm
return xmag,ymag,zmag
def gsmgse( p1, p2, p3, j ):
"""
converts geocentric solar magnetospheric (gsm) coords to solar ecliptic (gse) ones or vice versa.
j>0 j<0
input: xgsm,ygsm,zgsm xgse,ygse,zgse
output: xgse,ygse,zgse xgsm,ygsm,zgsm
:param p1,p2,p3: input position
:param j: flag
:return: output position
"""
# Uses the 'shi' & 'chi' global vars set by recalc.
if j > 0:
xgsm,ygsm,zgsm = p1, p2, p3
xgse = xgsm
ygse = ygsm*chi - zgsm*shi
zgse = ygsm*shi + zgsm*chi
return xgse,ygse,zgse
else:
xgse,ygse,zgse = p1, p2, p3
xgsm = xgse
ygsm = ygse*chi + zgse*shi
zgsm = zgse*chi - ygse*shi
return xgsm,ygsm,zgsm
def smgsm( p1, p2, p3, j ):
"""
Converts solar magnetic (sm) to geocentric solar magnetospheric (gsm) coordinates or vice versa.
j>0 j<0
input: xsm, ysm, zsm xgsm,ygsm,zgsm
output: xgsm,ygsm,zgsm xsm, ysm, zsm
:param p1,p2,p3: input position
:param j: flag
:return: output position
"""
# Uses the 'cps' & 'sps' values set by recalc.
if j > 0:
xsm,ysm,zsm = p1, p2, p3
xgsm = xsm*cps + zsm*sps
ygsm = ysm
zgsm = zsm*cps - xsm*sps
return xgsm,ygsm,zgsm
else:
xgsm,ygsm,zgsm = p1, p2, p3
xsm = xgsm*cps - zgsm*sps
ysm = ygsm
zsm = xgsm*sps + zgsm*cps
return xsm,ysm,zsm
def geogsm( p1, p2, p3, j ):
"""
Converts geographic (geo) to geocentric solar magnetospheric (gsm) coordinates or vice versa.
j>0 j<0
input: j,xgeo,ygeo,zgeo j,xgsm,ygsm,zgsm
output: xgsm,ygsm,zgsm xgeo,ygeo,zgeo
:param p1,p2,p3: input position
:param j: flag
:return: output position
"""
# Uses the a11-a33 rotation matrix columns given by recalc.
if j > 0:
xgeo,ygeo,zgeo = p1, p2, p3
xgsm = a11*xgeo + a12*ygeo + a13*zgeo
ygsm = a21*xgeo + a22*ygeo + a23*zgeo
zgsm = a31*xgeo + a32*ygeo + a33*zgeo
return xgsm,ygsm,zgsm
else:
xgsm,ygsm,zgsm = p1, p2, p3
xgeo = a11*xgsm + a21*ygsm + a31*zgsm
ygeo = a12*xgsm + a22*ygsm + a32*zgsm
zgeo = a13*xgsm + a23*ygsm + a33*zgsm
return xgeo,ygeo,zgeo
def geodgeo( p1, p2, j):
"""
This subroutine (1) converts vertical local height (altitude) h and geodetic
latitude xmu into geocentric coordinates r and theta (geocentric radial
distance and colatitude, respectively; also known as ecef coordinates),
as well as (2) performs the inverse transformation from {r,theta} to {h,xmu}.
The subroutine uses world geodetic system wgs84 parameters for the earth's
ellipsoid. the angular quantities (geo colatitude theta and geodetic latitude
xmu) are in radians, and the distances (geocentric radius r and altitude h
above the earth's ellipsoid) are in kilometers.
if j>0, the transformation is made from geodetic to geocentric coordinates using simple direct equations.
if j<0, the inverse transformation from geocentric to geodetic coordinates is made by means of a fast iterative algorithm.
j>0 j<0
input: h,xmu r,theta
output: r,theta h,xmu
Author: N.A. Tsyganenko
Date: Dec 5, 2007
:param h: Geodetic Altitude in km.
:param xmu: Geodetic latitude in radian.
:param r: Geocentric distance in km.
:param theta: Spherical co-latitude in radian.
"""
# WGS84 constants:
a = 6378.137 # Kilometers, semi-major axis
b = 6356.752314140 # Kilometers, semi-minor axis
Earth_Radius = 6378.137 # Kilometers. Differentiated from 'a' to follow legacy code.
if j > 0: # GEOD->GEO
# Rename constants
h, phi = p1, p2 # h is kilometers above earth's surface. Phi is geodetic latitude (to earth's surface).
# Convert using the GEOD->Cartesian->GEO, assuming a longitude (tau) of zero.
# As according to https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
# Pre-compute N(phi)
N_Phi = a**2 / ( a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2 )**0.5
X = ( N_Phi + h ) * np.cos(phi)
Z = (b**2 / a**2 * N_Phi + h ) * np.sin(phi)
#Y=0, as using latitude zero. sin(0)=0
# And then compute the cartesian radius & latitude from X,Z.
Radius = (X**2 + Z**2)**0.5 # Kilometers cartesian.
Theta = np.arccos( Z / Radius ) # Radians co latitude (North Pole=0, South Pole=pi)
return Radius, Theta
else: # GEO->GEOD
# WH - I was unable to find a clearer explaination of the algorithm used, so this is just a readable port of the legacy code.
r,theta = p1,p2 # Kilometres from center of earth, geographic co-latitude.
phi = np.pi/2 - theta # Input 'theta' constant is a co-latitude, so 90-Theta=Latitutde.
phi1 = phi.copy() # Rad. Current itteration's estimate of geographic angle to the surface-level of the geodetic sphere.
dphi = 0 # Rad. Current itteration's error between the geographic latitude found by the itteration, and the input saught geographic latitude.
tol = 1e-8 # Rad. Accepted tolerance of dphi before breaking early.
h = 0 # Kilometers. Return value, the height above ground.
xmu = 0 # Radians. Return value, the geodetic latitude.
# Pre-calculate the 'Beta' value.
beta = a**2 / b**2 - 1
# Iterate through the function 100x, or if the iteration step, for all values, is LT the tolerable error.
for n in range(100):
# There's a known error (1e-11) in the 'Beta' value used by geopack.
# Find the geodetic latitude of the earth's surface of the current itteration's estimate (phi1)
arg = np.sin(phi1) * ( 1 + beta ) / np.sqrt( 1 + beta * ( 2 + beta ) * np.sin(phi1)**2 )
xmu = np.arcsin( arg )
# This block calculates this estimate's 'h' value above ground-level.
rs = Earth_Radius / np.sqrt( 1 + beta * np.sin(phi1)**2)
cosfims = np.cos( phi1 - xmu )
h = np.sqrt( ( rs*cosfims )**2 + r**2 - rs**2 ) - rs*cosfims
# From this itteration's XMU+h geodetic coordinates, calulate the current cartesian coordinates.
# Uses a variant provided in GEOD->GEO, basically uses the ground-projection angle to find the cartesian position of the surface,
# Then adds on the surface-to-space coordinates.
x = rs * np.cos( phi1 ) + h * np.cos( xmu )
z = rs * np.sin( phi1 ) + h * np.sin( xmu )
# These is no Y-component, as we're assuming longitude=0. Only need to consider the latitudes.
# Calculate the current angle-error between the geographic latitude found in reverse by this itteration,
# And then apply this error to the next itteration's 'phi1' seed.
dphi = np.arctan2(z,x) - phi
phi1 -= dphi
# And if all value's iterative adjustments are less than the tolerable error, can break early.
if np.all( np.abs(dphi) < tol ): break
# After processing is completed, print the run-count, & max, average dphi error in arcseconds.
ArcSecErr = np.abs( dphi * 180/np.pi * 3600 )
#print("Loop Count: {:d} - Max Error {:.6f} arcsec, Avg error {:.6f} arcsec".format( n, np.max(ArcSecErr), np.mean(ArcSecErr) ))
#print( (datetime.datetime.now() - Start_Time).total_seconds() )
# for i in range(len( dphi )):
# if max(np.abs(dphi)) == np.abs(dphi[i]):
# print( p1[i], phi[i]*180/np.pi )
return h,xmu
def sphcar( p1, p2, p3, j ):
"""
Converts spherical coords into cartesian ones and vice versa (theta and phi in radians).
j>0 j<0
input: r,theta,phi x,y,z
output: x,y,z r,theta,phi
:param r,theta,phi:
:param x,y,z:
:param j:
:return:
Note: at the poles (x=0 and y=0) we assume phi=0 (when converting from cartesian to spherical coords, i.e., for j<0)
Last mofification: April 1, 2003 (only some notation changes and more comments added)
Author: N.A. Tsyganenko
"""
# For some blasted reason, Geopack has decied that Theta and Phi map to Lat/Lon rather than Lon/Lat, as damn near every other
# translation system uses. Output matches Geopack.
#
if j > 0:
r,theta,phi = p1, p2, p3
x = r * np.sin(theta) * np.cos(phi)
y = r * np.sin(theta) * np.sin(phi)
z = r * np.cos(theta)
return x, y, z
else:
x, y, z = p1, p2, p3
r = np.sqrt( x**2 + y**2 + z**2 )
theta = np.arccos ( z / r )
phi = np.arctan2( y , x )
return r, theta, phi
##################################
# Magnetic Field-Dependant translations & functions.
# Due to the primary function igrf_geo() not yet being ported, all functions below are prototyped, but not implemented.
# Each will raise a 'NotImplementedError' exception
def bspcar( theta, phi, br, btheta, bphi ):
raise NotImplementedError()
"""
Calculates cartesian field components from spherical ones.
:param theta,phi: spherical angles of the point in radians
:param br,btheta,bphi: spherical components of the field
:return: bx,by,bz. cartesian components of the field
# Last mofification: April 1, 2003 (only some notation changes and more comments added)
# Author: N.A. Tsyganenko
"""
be = br*np.sin(theta) + btheta*np.cos(theta)
bx = br*np.sin(theta) + btheta*np.cos(theta)*np.cos(phi) - bphi*np.sin(phi)
by = be*np.sin(phi) + bphi*np.cos(phi)
bz = br*np.cos(theta) - btheta*np.sin(theta)
return bx,by,bz
def bcarsp( x, y, z, bx, by, bz ):
raise NotImplementedError()
"""
Calculates spherical field components from those in cartesian system
:param x,y,z: cartesian components of the position vector
:param bx,by,bz: cartesian components of the field vector
:return: br,btheta,bphi. Local spherical components of the field vector
Last mofification: April 1, 2003 (only some notation changes and more comments added)
Author: N.A. Tsyganenko
"""
rho = np.sqrt(x**2 + y**2)
r = np.sqrt(x**2 + y**2 + z**2)
_Nonzero_Mask = rho != 0
cphi = np.ones ( len(rho), dtype=float )
sphi = np.zeros( len(rho), dtype=float )
# And only calculate for entries that will not divide-by-zero.
cphi [_Nonzero_Mask] = x[_Nonzero_Mask] / rho[_Nonzero_Mask]
sphi [_Nonzero_Mask] = y[_Nonzero_Mask] / rho[_Nonzero_Mask]
ct = z / r
st = rho / r
br = ( x*bx + y*by + z*bz ) / r
btheta = ( bx*cphi + by*sphi ) * ct - bz*st
bphi = by*cphi - bx*sphi
return br, btheta, bphi
def igrf_gsw(xgsw,ygsw,zgsw):
raise NotImplementedError()
"""
Calculates components of the main (internal) geomagnetic field in the geocentric solar
magnetospheric coordinate system, using IAGA international geomagnetic reference model
coefficients (e.g., http://www.ngdc.noaa.gov/iaga/vmod/igrf.html revised: 22 march, 2005)
Before the first call of this subroutine, or if the date/time
was changed, the model coefficients and GEO-GSW rotation matrix elements should be updated
by calling the subroutine recalc
Python version by Sheng Tian
:param xgsw,ygsw,zgsw: cartesian GSW coordinates (in units Re=6371.2 km)
:return: hxgsw,hygsw,hzgsw. Cartesian GSW components of the main geomagnetic field in nanotesla
"""
xgsm, ygsm, zgsm = gswgsm( xgsw, ygsw, zgsw, 1)
bxgsm,bygsm,bzgsm = igrf_gsm( xgsm, ygsm, zgsm)
bxgsw,bygsw,gzgsw = gswgsm( bxgsm,bygsm,bzgsm, -1)
return bxgsw,bygsw,gzgsw
def igrf_gsm(xgsm,ygsm,zgsm):
raise NotImplementedError()
"""
Calculates components of the main (internal) geomagnetic field in the geocentric solar
magnetospheric coordinate system, using IAGA international geomagnetic reference model
coefficients (e.g., http://www.ngdc.noaa.gov/iaga/vmod/igrf.html revised: 22 march, 2005)
Before the first call of this subroutine, or if the date/time
was changed, the model coefficients and GEO-GSM rotation matrix elements should be updated
by calling the subroutine recalc
Python version by Sheng Tian
:param xgsm,ygsm,zgsm: cartesian GSM coordinates (in units Re=6371.2 km)
:return: hxgsm,hygsm,hzgsm. Cartesian GSM components of the main geomagnetic field in nanotesla
"""
xgeo, ygeo, zgeo = geogsm ( xgsm, ygsm, zgsm, -1 )
r, theta, phi = sphcar ( xgeo, ygeo, zgeo, -1 )
br, btheta, bphi = igrf_geo( r, theta, phi)
bxgeo,bygeo,bzgeo = bspcar ( theta, phi, br, btheta, bphi)
bxgsm,bygsm,gzgsm = geogsm ( bxgeo, bygeo,bzgeo, 1)
return bxgsm, bygsm, gzgsm
def igrf_geo(r,theta,phi):
raise NotImplementedError()
"""
Calculates components of the main (internal) geomagnetic field in the spherical geographic
(geocentric) coordinate system, using IAGA international geomagnetic reference model
coefficients (e.g., http://www.ngdc.noaa.gov/iaga/vmod/igrf.html, revised: 22 march, 2005)
Before the first call of this subroutine, or if the time was changed,
the model coefficients should be updated by calling the subroutine recalc
Python version by Sheng Tian
:param r: spherical geographic (geocentric) coordinates: radial distance r in units Re=6371.2 km
:param theta: colatitude theta in radians
:param phi: longitude phi in radians
:return: br, btheta, bphi. Spherical components of the main geomagnetic field in nanotesla
(positive br outward, btheta southward, bphi eastward)
"""
# Uses the global 'g', 'h', & 'rec' variables.
# g & h are per-recalc value.
# rec is a constant linear array of 105. (As of IGRF_13 in 2021-03)
# No simple way of doing this in BLAS, as varying Earth-Radii can cause the 'k' array to be of uneven size.
# Thus, for sanity, seperating this into a loop.
# TODO. Get this readable & fast.
Return_Mag_Rad_Theta_Phi = []
for i in range( len(r) ):
This_R, This_Theta, This_Phi = r[i] / 6.7e6, theta[i], phi[i]
This_g, This_h = g[i], h[i]
ct = np.cos( This_Theta )
st = np.sin( This_Theta )
# Minimum Theta value.
minst = 1e-5
if np.abs(st) < minst: smlst = True
else: smlst = False
# In this new version, the optimal value of the parameter nm (maximal order of the spherical
# harmonic expansion) is not user-prescribed, but calculated inside the subroutine, based
# on the value of the radial distance r:
irp3 = np.int64(This_R+2)
nm = np.int64(3+30/irp3)
if nm > 13: nm = 13
k = nm+1
# r dependence is encapsulated here.
a = np.empty(k)
b = np.empty(k)
ar = 1/This_R # a/r
a[0] = ar*ar # a[n] = (a/r)^(n+2).
b[0] = a[0] # b[n] = (n+1)(a/r)^(n+2)
for n in range(1,k):
a[n] = a[n-1]*ar
b[n] = a[n]*(n+1)
# t - short for theta, f - short for phi.
br,bt,bf = [0.]*3
d,p = [0.,1]
# m = 0. P^n,0
m = 0
smf,cmf = [0.,1]
p1,d1,p2,d2 = [p,d,0.,0]
l0 = 0
mn = l0
for n in range(m,k):
w = This_g[mn]*cmf + This_h[mn]*smf
br += b[n]*w*p1 # p1 is P^n,m.
bt -= a[n]*w*d1 # d1 is dP^n,m/dt.
xk = rec[mn]
# Eq 16c and its derivative on theta.
d0 = ct*d1-st*p1-xk*d2 # dP^n,m/dt = ct*dP^n-1,m/dt - st*P_n-1,m - K^n,m*dP^n-2,m/dt
p0 = ct*p1-xk*p2 # P^n,m = ct*P^n-1,m - K^n,m*P^n-2,m
d2,p2,d1 = [d1,p1,d0]
p1 = p0
mn += n+1
# Eq 16b and its derivative on theta.
d = st*d+ct*p # dP^m,m/dt = st*dP^m-1,m-1/dt + ct*P^m-1,m-1
p = st*p # P^m,m = st*P^m-1,m-1
# Similarly for P^n,m
l0 = 0
for m in range(1,k): # sum over m
smf = np.sin(m*This_Phi) # sin(m*phi)
cmf = np.cos(m*This_Phi) # cos(m*phi)
p1,d1,p2,d2 = [p,d,0.,0]
tbf = 0.
l0 += m+1
mn = l0
for n in range(m,k): # sum over n
w=This_g[mn]*cmf+This_h[mn]*smf # [g^n,m*cos(m*phi)+h^n,m*sin(m*phi)]
br += b[n]*w*p1
bt -= a[n]*w*d1
tp = p1
if smlst: tp = d1
tbf += a[n]*(This_g[mn]*smf-This_h[mn]*cmf)*tp
xk = rec[mn]
d0 = ct*d1-st*p1-xk*d2 # dP^n,m/dt = ct*dP^n-1,m/dt - st*P_n-1,m - K^n,m*dP^n-2,m/dt
p0 = ct*p1-xk*p2 # P^n,m = ct*P^n-1,m - K^n,m*P^n-2,m
d2,p2,d1 = [d1,p1,d0]
p1=p0
mn += n+1
d = st*d+ct*p
p = st*p
# update B_phi.
tbf *= m
bf += tbf
if smlst:
if ct < 0.: bf = -bf
else: bf /= st
print( br, br, bf )
Return_Mag_Rad_Theta_Phi.append( [br, bt, bf] )
#return Return_Mag_Rad_Theta_Phi
# Translate the return list into a NP array and split into bRad bTheta bPhi
bRad, bTheta, bPhi = np.squeeze( np.hsplit( Return_Mag_Rad_Theta_Phi, 3 ) )
return bRad, bTheta, bPhi
def dip(xgsm,ygsm,zgsm):
raise NotImplementedError()
"""
Calculates gsm components of a geodipole field with the dipole moment
corresponding to the epoch, specified by calling subroutine recalc (should be
invoked before the first use of this one and in case the date/time was changed).
:param xgsm,ygsm,zgsm: GSM coordinates in Re (1 Re = 6371.2 km)
:return: bxgsm,bygsm,gzgsm. Field components in gsm system, in nanotesla.
Last modification: May 4, 2005.
Author: N. A. Tsyganenko
"""
# Uses the following globals from recalc():
# sps,cps, g, h
dipmom = np.sqrt( g[:,1]**2 + g[:,2]**2 + h[:,2]**2 )
p = xgsm**2
u = zgsm**2
v = 3*zgsm*xgsm
t = ygsm**2
q = dipmom / np.sqrt( p + t + u )**5
bxgsm = q * (( t + u - 2*p )*sps - v*cps )
bygsm = q * ( xgsm*sps + zgsm*cps ) * -3*ygsm
bzgsm = q * (( p + t - 2*u )*cps - v*sps )
return bxgsm,bygsm,bzgsm
def dip_gsw(xgsw,ygsw,zgsw):
raise NotImplementedError()
"""
Calculates gsm components of a geodipole field with the dipole moment
corresponding to the epoch, specified by calling subroutine recalc (should be
invoked before the first use of this one and in case the date/time was changed).
:param xgsw,ygsw,zgsw: GSW coordinates in Re (1 Re = 6371.2 km)
:return: bxgsm,bygsm,gzgsm. Field components in gsm system, in nanotesla.
Author: Sheng Tian
"""
xgsm, ygsm, zgsm = gswgsm( xgsw, ygsw, zgsw, 1 )
bxgsm,bygsm,bzgsm = dip ( xgsm, ygsm, zgsm )
bxgsw,bygsw,bzgsw = gswgsm( bxgsm, bygsm, bzgsm,-1)
return bxgsw,bygsw,bzgsw
Functions
def bcarsp(x, y, z, bx, by, bz)
-
Expand source code
def bcarsp( x, y, z, bx, by, bz ): raise NotImplementedError() """ Calculates spherical field components from those in cartesian system :param x,y,z: cartesian components of the position vector :param bx,by,bz: cartesian components of the field vector :return: br,btheta,bphi. Local spherical components of the field vector Last mofification: April 1, 2003 (only some notation changes and more comments added) Author: N.A. Tsyganenko """ rho = np.sqrt(x**2 + y**2) r = np.sqrt(x**2 + y**2 + z**2) _Nonzero_Mask = rho != 0 cphi = np.ones ( len(rho), dtype=float ) sphi = np.zeros( len(rho), dtype=float ) # And only calculate for entries that will not divide-by-zero. cphi [_Nonzero_Mask] = x[_Nonzero_Mask] / rho[_Nonzero_Mask] sphi [_Nonzero_Mask] = y[_Nonzero_Mask] / rho[_Nonzero_Mask] ct = z / r st = rho / r br = ( x*bx + y*by + z*bz ) / r btheta = ( bx*cphi + by*sphi ) * ct - bz*st bphi = by*cphi - bx*sphi return br, btheta, bphi
def bspcar(theta, phi, br, btheta, bphi)
-
Expand source code
def bspcar( theta, phi, br, btheta, bphi ): raise NotImplementedError() """ Calculates cartesian field components from spherical ones. :param theta,phi: spherical angles of the point in radians :param br,btheta,bphi: spherical components of the field :return: bx,by,bz. cartesian components of the field # Last mofification: April 1, 2003 (only some notation changes and more comments added) # Author: N.A. Tsyganenko """ be = br*np.sin(theta) + btheta*np.cos(theta) bx = br*np.sin(theta) + btheta*np.cos(theta)*np.cos(phi) - bphi*np.sin(phi) by = be*np.sin(phi) + bphi*np.cos(phi) bz = br*np.cos(theta) - btheta*np.sin(theta) return bx,by,bz
def dip(xgsm, ygsm, zgsm)
-
Expand source code
def dip(xgsm,ygsm,zgsm): raise NotImplementedError() """ Calculates gsm components of a geodipole field with the dipole moment corresponding to the epoch, specified by calling subroutine recalc (should be invoked before the first use of this one and in case the date/time was changed). :param xgsm,ygsm,zgsm: GSM coordinates in Re (1 Re = 6371.2 km) :return: bxgsm,bygsm,gzgsm. Field components in gsm system, in nanotesla. Last modification: May 4, 2005. Author: N. A. Tsyganenko """ # Uses the following globals from recalc(): # sps,cps, g, h dipmom = np.sqrt( g[:,1]**2 + g[:,2]**2 + h[:,2]**2 ) p = xgsm**2 u = zgsm**2 v = 3*zgsm*xgsm t = ygsm**2 q = dipmom / np.sqrt( p + t + u )**5 bxgsm = q * (( t + u - 2*p )*sps - v*cps ) bygsm = q * ( xgsm*sps + zgsm*cps ) * -3*ygsm bzgsm = q * (( p + t - 2*u )*cps - v*sps ) return bxgsm,bygsm,bzgsm
def dip_gsw(xgsw, ygsw, zgsw)
-
Expand source code
def dip_gsw(xgsw,ygsw,zgsw): raise NotImplementedError() """ Calculates gsm components of a geodipole field with the dipole moment corresponding to the epoch, specified by calling subroutine recalc (should be invoked before the first use of this one and in case the date/time was changed). :param xgsw,ygsw,zgsw: GSW coordinates in Re (1 Re = 6371.2 km) :return: bxgsm,bygsm,gzgsm. Field components in gsm system, in nanotesla. Author: Sheng Tian """ xgsm, ygsm, zgsm = gswgsm( xgsw, ygsw, zgsw, 1 ) bxgsm,bygsm,bzgsm = dip ( xgsm, ygsm, zgsm ) bxgsw,bygsw,bzgsw = gswgsm( bxgsm, bygsm, bzgsm,-1) return bxgsw,bygsw,bzgsw
def geigeo(p1, p2, p3, j)
-
Converts equatorial inertial (gei) to geographical (geo) coords or vice versa. j>0 j<0 input: xgei,ygei,zgei xgeo,ygeo,zgeo output: xgeo,ygeo,zgeo xgei,ygei,zgei
:param p1,p2,p3: input position :param j: flag :return: output position
Expand source code
def geigeo( p1, p2, p3, j ): """ Converts equatorial inertial (gei) to geographical (geo) coords or vice versa. j>0 j<0 input: xgei,ygei,zgei xgeo,ygeo,zgeo output: xgeo,ygeo,zgeo xgei,ygei,zgei :param p1,p2,p3: input position :param j: flag :return: output position """ # Uses the 'cgst' & 'sgst' globals set by recalc. if j > 0: xgei,ygei,zgei = p1, p2, p3 xgeo = xgei*cgst + ygei*sgst ygeo = ygei*cgst - xgei*sgst zgeo = zgei return xgeo,ygeo,zgeo else: xgeo,ygeo,zgeo = p1, p2, p3 xgei = xgeo*cgst - ygeo*sgst ygei = ygeo*cgst + xgeo*sgst zgei = zgeo return xgei,ygei,zgei
def geodgeo(p1, p2, j)
-
This subroutine (1) converts vertical local height (altitude) h and geodetic latitude xmu into geocentric coordinates r and theta (geocentric radial distance and colatitude, respectively; also known as ecef coordinates), as well as (2) performs the inverse transformation from {r,theta} to {h,xmu}.
The subroutine uses world geodetic system wgs84 parameters for the earth's ellipsoid. the angular quantities (geo colatitude theta and geodetic latitude xmu) are in radians, and the distances (geocentric radius r and altitude h above the earth's ellipsoid) are in kilometers.
if j>0, the transformation is made from geodetic to geocentric coordinates using simple direct equations. if j<0, the inverse transformation from geocentric to geodetic coordinates is made by means of a fast iterative algorithm.
j>0 j<0
input: h,xmu r,theta output: r,theta h,xmu
Author: N.A. Tsyganenko Date: Dec 5, 2007
:param h: Geodetic Altitude in km. :param xmu: Geodetic latitude in radian. :param r: Geocentric distance in km. :param theta: Spherical co-latitude in radian.
Expand source code
def geodgeo( p1, p2, j): """ This subroutine (1) converts vertical local height (altitude) h and geodetic latitude xmu into geocentric coordinates r and theta (geocentric radial distance and colatitude, respectively; also known as ecef coordinates), as well as (2) performs the inverse transformation from {r,theta} to {h,xmu}. The subroutine uses world geodetic system wgs84 parameters for the earth's ellipsoid. the angular quantities (geo colatitude theta and geodetic latitude xmu) are in radians, and the distances (geocentric radius r and altitude h above the earth's ellipsoid) are in kilometers. if j>0, the transformation is made from geodetic to geocentric coordinates using simple direct equations. if j<0, the inverse transformation from geocentric to geodetic coordinates is made by means of a fast iterative algorithm. j>0 j<0 input: h,xmu r,theta output: r,theta h,xmu Author: N.A. Tsyganenko Date: Dec 5, 2007 :param h: Geodetic Altitude in km. :param xmu: Geodetic latitude in radian. :param r: Geocentric distance in km. :param theta: Spherical co-latitude in radian. """ # WGS84 constants: a = 6378.137 # Kilometers, semi-major axis b = 6356.752314140 # Kilometers, semi-minor axis Earth_Radius = 6378.137 # Kilometers. Differentiated from 'a' to follow legacy code. if j > 0: # GEOD->GEO # Rename constants h, phi = p1, p2 # h is kilometers above earth's surface. Phi is geodetic latitude (to earth's surface). # Convert using the GEOD->Cartesian->GEO, assuming a longitude (tau) of zero. # As according to https://en.wikipedia.org/wiki/Geographic_coordinate_conversion # Pre-compute N(phi) N_Phi = a**2 / ( a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2 )**0.5 X = ( N_Phi + h ) * np.cos(phi) Z = (b**2 / a**2 * N_Phi + h ) * np.sin(phi) #Y=0, as using latitude zero. sin(0)=0 # And then compute the cartesian radius & latitude from X,Z. Radius = (X**2 + Z**2)**0.5 # Kilometers cartesian. Theta = np.arccos( Z / Radius ) # Radians co latitude (North Pole=0, South Pole=pi) return Radius, Theta else: # GEO->GEOD # WH - I was unable to find a clearer explaination of the algorithm used, so this is just a readable port of the legacy code. r,theta = p1,p2 # Kilometres from center of earth, geographic co-latitude. phi = np.pi/2 - theta # Input 'theta' constant is a co-latitude, so 90-Theta=Latitutde. phi1 = phi.copy() # Rad. Current itteration's estimate of geographic angle to the surface-level of the geodetic sphere. dphi = 0 # Rad. Current itteration's error between the geographic latitude found by the itteration, and the input saught geographic latitude. tol = 1e-8 # Rad. Accepted tolerance of dphi before breaking early. h = 0 # Kilometers. Return value, the height above ground. xmu = 0 # Radians. Return value, the geodetic latitude. # Pre-calculate the 'Beta' value. beta = a**2 / b**2 - 1 # Iterate through the function 100x, or if the iteration step, for all values, is LT the tolerable error. for n in range(100): # There's a known error (1e-11) in the 'Beta' value used by geopack. # Find the geodetic latitude of the earth's surface of the current itteration's estimate (phi1) arg = np.sin(phi1) * ( 1 + beta ) / np.sqrt( 1 + beta * ( 2 + beta ) * np.sin(phi1)**2 ) xmu = np.arcsin( arg ) # This block calculates this estimate's 'h' value above ground-level. rs = Earth_Radius / np.sqrt( 1 + beta * np.sin(phi1)**2) cosfims = np.cos( phi1 - xmu ) h = np.sqrt( ( rs*cosfims )**2 + r**2 - rs**2 ) - rs*cosfims # From this itteration's XMU+h geodetic coordinates, calulate the current cartesian coordinates. # Uses a variant provided in GEOD->GEO, basically uses the ground-projection angle to find the cartesian position of the surface, # Then adds on the surface-to-space coordinates. x = rs * np.cos( phi1 ) + h * np.cos( xmu ) z = rs * np.sin( phi1 ) + h * np.sin( xmu ) # These is no Y-component, as we're assuming longitude=0. Only need to consider the latitudes. # Calculate the current angle-error between the geographic latitude found in reverse by this itteration, # And then apply this error to the next itteration's 'phi1' seed. dphi = np.arctan2(z,x) - phi phi1 -= dphi # And if all value's iterative adjustments are less than the tolerable error, can break early. if np.all( np.abs(dphi) < tol ): break # After processing is completed, print the run-count, & max, average dphi error in arcseconds. ArcSecErr = np.abs( dphi * 180/np.pi * 3600 ) #print("Loop Count: {:d} - Max Error {:.6f} arcsec, Avg error {:.6f} arcsec".format( n, np.max(ArcSecErr), np.mean(ArcSecErr) )) #print( (datetime.datetime.now() - Start_Time).total_seconds() ) # for i in range(len( dphi )): # if max(np.abs(dphi)) == np.abs(dphi[i]): # print( p1[i], phi[i]*180/np.pi ) return h,xmu
def geogsm(p1, p2, p3, j)
-
Converts geographic (geo) to geocentric solar magnetospheric (gsm) coordinates or vice versa. j>0 j<0 input: j,xgeo,ygeo,zgeo j,xgsm,ygsm,zgsm output: xgsm,ygsm,zgsm xgeo,ygeo,zgeo
:param p1,p2,p3: input position :param j: flag :return: output position
Expand source code
def geogsm( p1, p2, p3, j ): """ Converts geographic (geo) to geocentric solar magnetospheric (gsm) coordinates or vice versa. j>0 j<0 input: j,xgeo,ygeo,zgeo j,xgsm,ygsm,zgsm output: xgsm,ygsm,zgsm xgeo,ygeo,zgeo :param p1,p2,p3: input position :param j: flag :return: output position """ # Uses the a11-a33 rotation matrix columns given by recalc. if j > 0: xgeo,ygeo,zgeo = p1, p2, p3 xgsm = a11*xgeo + a12*ygeo + a13*zgeo ygsm = a21*xgeo + a22*ygeo + a23*zgeo zgsm = a31*xgeo + a32*ygeo + a33*zgeo return xgsm,ygsm,zgsm else: xgsm,ygsm,zgsm = p1, p2, p3 xgeo = a11*xgsm + a21*ygsm + a31*zgsm ygeo = a12*xgsm + a22*ygsm + a32*zgsm zgeo = a13*xgsm + a23*ygsm + a33*zgsm return xgeo,ygeo,zgeo
def geomag(p1, p2, p3, j)
-
Converts geographic (geo) to dipole (mag) coordinates or vice versa. j>0 j<0 input: j,xgeo,ygeo,zgeo j,xmag,ymag,zmag output: xmag,ymag,zmag xgeo,ygeo,zgeo
:param p1,p2,p3: input position :param j: flag :return: output position
Expand source code
def geomag( p1, p2, p3, j ): """ Converts geographic (geo) to dipole (mag) coordinates or vice versa. j>0 j<0 input: j,xgeo,ygeo,zgeo j,xmag,ymag,zmag output: xmag,ymag,zmag xgeo,ygeo,zgeo :param p1,p2,p3: input position :param j: flag :return: output position """ if j > 0: xgeo, ygeo, zgeo = p1, p2, p3 xmag = xgeo*ctcl + ygeo*ctsl - zgeo*st0 ymag = -xgeo*sl0 + ygeo*cl0 zmag = xgeo*stcl + ygeo*stsl + zgeo*ct0 return xmag, ymag, zmag else: xmag, ymag, zmag = p1, p2, p3 xgeo = xmag*ctcl - ymag*sl0 + zmag*stcl ygeo = xmag*ctsl + ymag*cl0 + zmag*stsl zgeo = -xmag*st0 + zmag*ct0 return xgeo, ygeo, zgeo
def gsmgse(p1, p2, p3, j)
-
converts geocentric solar magnetospheric (gsm) coords to solar ecliptic (gse) ones or vice versa. j>0 j<0 input: xgsm,ygsm,zgsm xgse,ygse,zgse output: xgse,ygse,zgse xgsm,ygsm,zgsm
:param p1,p2,p3: input position :param j: flag :return: output position
Expand source code
def gsmgse( p1, p2, p3, j ): """ converts geocentric solar magnetospheric (gsm) coords to solar ecliptic (gse) ones or vice versa. j>0 j<0 input: xgsm,ygsm,zgsm xgse,ygse,zgse output: xgse,ygse,zgse xgsm,ygsm,zgsm :param p1,p2,p3: input position :param j: flag :return: output position """ # Uses the 'shi' & 'chi' global vars set by recalc. if j > 0: xgsm,ygsm,zgsm = p1, p2, p3 xgse = xgsm ygse = ygsm*chi - zgsm*shi zgse = ygsm*shi + zgsm*chi return xgse,ygse,zgse else: xgse,ygse,zgse = p1, p2, p3 xgsm = xgse ygsm = ygse*chi + zgse*shi zgsm = zgse*chi - ygse*shi return xgsm,ygsm,zgsm
def gswgsm(p1, p2, p3, j)
-
Converts gsm to gsw coordinates or vice versa. j>0 j<0 input: xgsm,ygsm,zgsm xgsw,ygsw,zgsw output: xgsw,ygsw,zgsw xgsm,ygsm,zgsm
:param p1,p2,p3: input position :param j: flag :return: output position
Expand source code
def gswgsm( p1, p2, p3, j ): """ Converts gsm to gsw coordinates or vice versa. j>0 j<0 input: xgsm,ygsm,zgsm xgsw,ygsw,zgsw output: xgsw,ygsw,zgsw xgsm,ygsm,zgsm :param p1,p2,p3: input position :param j: flag :return: output position """ if j > 0: xgsw, ygsw, zgsw = p1, p2, p3 xgsm = xgsw*e11 + ygsw*e12 + zgsw*e13 ygsm = xgsw*e21 + ygsw*e22 + zgsw*e23 zgsm = xgsw*e31 + ygsw*e32 + zgsw*e33 return xgsm,ygsm,zgsm else: xgsm, ygsm, zgsm = p1, p2, p3 xgsw = xgsm*e11 + ygsm*e21 + zgsm*e31 ygsw = xgsm*e12 + ygsm*e22 + zgsm*e32 zgsw = xgsm*e13 + ygsm*e23 + zgsm*e33 return xgsw,ygsw,zgsw
def igrf_geo(r, theta, phi)
-
Expand source code
def igrf_geo(r,theta,phi): raise NotImplementedError() """ Calculates components of the main (internal) geomagnetic field in the spherical geographic (geocentric) coordinate system, using IAGA international geomagnetic reference model coefficients (e.g., http://www.ngdc.noaa.gov/iaga/vmod/igrf.html, revised: 22 march, 2005) Before the first call of this subroutine, or if the time was changed, the model coefficients should be updated by calling the subroutine recalc Python version by Sheng Tian :param r: spherical geographic (geocentric) coordinates: radial distance r in units Re=6371.2 km :param theta: colatitude theta in radians :param phi: longitude phi in radians :return: br, btheta, bphi. Spherical components of the main geomagnetic field in nanotesla (positive br outward, btheta southward, bphi eastward) """ # Uses the global 'g', 'h', & 'rec' variables. # g & h are per-recalc value. # rec is a constant linear array of 105. (As of IGRF_13 in 2021-03) # No simple way of doing this in BLAS, as varying Earth-Radii can cause the 'k' array to be of uneven size. # Thus, for sanity, seperating this into a loop. # TODO. Get this readable & fast. Return_Mag_Rad_Theta_Phi = [] for i in range( len(r) ): This_R, This_Theta, This_Phi = r[i] / 6.7e6, theta[i], phi[i] This_g, This_h = g[i], h[i] ct = np.cos( This_Theta ) st = np.sin( This_Theta ) # Minimum Theta value. minst = 1e-5 if np.abs(st) < minst: smlst = True else: smlst = False # In this new version, the optimal value of the parameter nm (maximal order of the spherical # harmonic expansion) is not user-prescribed, but calculated inside the subroutine, based # on the value of the radial distance r: irp3 = np.int64(This_R+2) nm = np.int64(3+30/irp3) if nm > 13: nm = 13 k = nm+1 # r dependence is encapsulated here. a = np.empty(k) b = np.empty(k) ar = 1/This_R # a/r a[0] = ar*ar # a[n] = (a/r)^(n+2). b[0] = a[0] # b[n] = (n+1)(a/r)^(n+2) for n in range(1,k): a[n] = a[n-1]*ar b[n] = a[n]*(n+1) # t - short for theta, f - short for phi. br,bt,bf = [0.]*3 d,p = [0.,1] # m = 0. P^n,0 m = 0 smf,cmf = [0.,1] p1,d1,p2,d2 = [p,d,0.,0] l0 = 0 mn = l0 for n in range(m,k): w = This_g[mn]*cmf + This_h[mn]*smf br += b[n]*w*p1 # p1 is P^n,m. bt -= a[n]*w*d1 # d1 is dP^n,m/dt. xk = rec[mn] # Eq 16c and its derivative on theta. d0 = ct*d1-st*p1-xk*d2 # dP^n,m/dt = ct*dP^n-1,m/dt - st*P_n-1,m - K^n,m*dP^n-2,m/dt p0 = ct*p1-xk*p2 # P^n,m = ct*P^n-1,m - K^n,m*P^n-2,m d2,p2,d1 = [d1,p1,d0] p1 = p0 mn += n+1 # Eq 16b and its derivative on theta. d = st*d+ct*p # dP^m,m/dt = st*dP^m-1,m-1/dt + ct*P^m-1,m-1 p = st*p # P^m,m = st*P^m-1,m-1 # Similarly for P^n,m l0 = 0 for m in range(1,k): # sum over m smf = np.sin(m*This_Phi) # sin(m*phi) cmf = np.cos(m*This_Phi) # cos(m*phi) p1,d1,p2,d2 = [p,d,0.,0] tbf = 0. l0 += m+1 mn = l0 for n in range(m,k): # sum over n w=This_g[mn]*cmf+This_h[mn]*smf # [g^n,m*cos(m*phi)+h^n,m*sin(m*phi)] br += b[n]*w*p1 bt -= a[n]*w*d1 tp = p1 if smlst: tp = d1 tbf += a[n]*(This_g[mn]*smf-This_h[mn]*cmf)*tp xk = rec[mn] d0 = ct*d1-st*p1-xk*d2 # dP^n,m/dt = ct*dP^n-1,m/dt - st*P_n-1,m - K^n,m*dP^n-2,m/dt p0 = ct*p1-xk*p2 # P^n,m = ct*P^n-1,m - K^n,m*P^n-2,m d2,p2,d1 = [d1,p1,d0] p1=p0 mn += n+1 d = st*d+ct*p p = st*p # update B_phi. tbf *= m bf += tbf if smlst: if ct < 0.: bf = -bf else: bf /= st print( br, br, bf ) Return_Mag_Rad_Theta_Phi.append( [br, bt, bf] ) #return Return_Mag_Rad_Theta_Phi # Translate the return list into a NP array and split into bRad bTheta bPhi bRad, bTheta, bPhi = np.squeeze( np.hsplit( Return_Mag_Rad_Theta_Phi, 3 ) ) return bRad, bTheta, bPhi
def igrf_gsm(xgsm, ygsm, zgsm)
-
Expand source code
def igrf_gsm(xgsm,ygsm,zgsm): raise NotImplementedError() """ Calculates components of the main (internal) geomagnetic field in the geocentric solar magnetospheric coordinate system, using IAGA international geomagnetic reference model coefficients (e.g., http://www.ngdc.noaa.gov/iaga/vmod/igrf.html revised: 22 march, 2005) Before the first call of this subroutine, or if the date/time was changed, the model coefficients and GEO-GSM rotation matrix elements should be updated by calling the subroutine recalc Python version by Sheng Tian :param xgsm,ygsm,zgsm: cartesian GSM coordinates (in units Re=6371.2 km) :return: hxgsm,hygsm,hzgsm. Cartesian GSM components of the main geomagnetic field in nanotesla """ xgeo, ygeo, zgeo = geogsm ( xgsm, ygsm, zgsm, -1 ) r, theta, phi = sphcar ( xgeo, ygeo, zgeo, -1 ) br, btheta, bphi = igrf_geo( r, theta, phi) bxgeo,bygeo,bzgeo = bspcar ( theta, phi, br, btheta, bphi) bxgsm,bygsm,gzgsm = geogsm ( bxgeo, bygeo,bzgeo, 1) return bxgsm, bygsm, gzgsm
def igrf_gsw(xgsw, ygsw, zgsw)
-
Expand source code
def igrf_gsw(xgsw,ygsw,zgsw): raise NotImplementedError() """ Calculates components of the main (internal) geomagnetic field in the geocentric solar magnetospheric coordinate system, using IAGA international geomagnetic reference model coefficients (e.g., http://www.ngdc.noaa.gov/iaga/vmod/igrf.html revised: 22 march, 2005) Before the first call of this subroutine, or if the date/time was changed, the model coefficients and GEO-GSW rotation matrix elements should be updated by calling the subroutine recalc Python version by Sheng Tian :param xgsw,ygsw,zgsw: cartesian GSW coordinates (in units Re=6371.2 km) :return: hxgsw,hygsw,hzgsw. Cartesian GSW components of the main geomagnetic field in nanotesla """ xgsm, ygsm, zgsm = gswgsm( xgsw, ygsw, zgsw, 1) bxgsm,bygsm,bzgsm = igrf_gsm( xgsm, ygsm, zgsm) bxgsw,bygsw,gzgsw = gswgsm( bxgsm,bygsm,bzgsm, -1) return bxgsw,bygsw,gzgsw
def init_igrf(IGRF_Filepath='/mnt/c/Users/Warren Holley/Desktop/ephemeris_library_3/ephemeris_library/rotation_libraries/igrf13coeffs.txt')
-
Expand source code
def init_igrf( IGRF_Filepath = DEFAULT_IGRF_PATH ): # Initialize the IGRF coefficients and related coefs. global igrf, nmn,mns, nyear,years,yruts # Read the IGRF file into a pandas dataframe IGRF_df = pd.read_csv( IGRF_Filepath, sep='\s+', skiprows=3) IGRF_np = IGRF_df.to_numpy() Year_Str = IGRF_df.columns[3:] # Extract the year list from the column name, then parse to integers. Year_Int = np.array([Year[:4] for Year in Year_Str], dtype=int) # The [:4] scrubs the float, and the final entry's year range Year_Int[-1] += 5 # And then add the '+5' range s.t. the "2020-5" == 2025 # And finally write that to the global var. years = Year_Int nyear = len(Year_Int) # From the final line of the IGRF file, extract the order of the polynomial, and the number of arguments used to generate it. k = IGRF_np[-1,1] + 1 # k = 14 at time of writing. nmn = (k+1) * k // 2 # nmn= 105 igrf = np.zeros(( nmn, nyear, 2), dtype=float) mns = np.zeros(( nmn, 2), dtype=float) l = 0 for i in range(k): for j in range(i + 1): mns[l, :] = [i, j] l += 1 # For each row in the file: for Row in IGRF_np: if Row[0] == 'g': i = 0 else: i = 1 n, m = Row[1:3].astype(int) mn = n * (n + 1) // 2 + m # The row-index of the related pair in the 'mns' array. igrf[mn, :, i] = Row[3:] # I have no idea what this is supposed to accomplish. igrf[:, -1, :] = igrf[:, -2, :] + igrf[:, -1, :] * 5 # And then build the 'year-ut-second' variable, the number of seconds since the unix epoch for each year-column. yruts = np.zeros(nyear, dtype=float) Unix_Epoch = datetime.datetime(1970,1,1) for i in range(nyear): yruts[i] = (datetime.datetime(years[i],1,1) - Unix_Epoch).total_seconds() # Translate the 'igrf' var into a dictionary of 'g' & 'h', of 2d arrays of 105x26 (default). igrf = {'g': igrf[:,:,0], 'h': igrf[:,:,1]}
def load_igrf(ut)
-
Expand source code
def load_igrf( ut ): if igrf is None: init_igrf() ut = _Sanitize_Input( False, ut )[0] # Make sure the array is at least 1d. # Seek the list of indecies s.t. can use the values on either side to extrapolate between. Indecies = np.searchsorted( yruts, ut, side='right' ) - 1 # In the event that values are inserted -after- the end of the # given constants, extrapolate using the [-2] to [-1] constants. Indecies[Indecies < 0 ] = 0 Indecies[Indecies >= len( yruts )-1 ] = len( yruts ) - 2 g0 = igrf['g'][:, Indecies ] h0 = igrf['h'][:, Indecies ] g1 = igrf['g'][:,Indecies+1] h1 = igrf['h'][:,Indecies+1] # Calculate the time-fractions to use. Interp_Start = yruts[Indecies] Interp_End = yruts[Indecies+1] Weight_Right = (ut - Interp_Start) / (Interp_End - Interp_Start) # The time-fraction of the right-side components Weight_Left = 1 - Weight_Right # The time-fraction of the left-side components g_Estimate = g0 * Weight_Left + g1 * Weight_Right h_Estimate = h0 * Weight_Left + h1 * Weight_Right # Transpose s.t. they're row-indexed (Eg: g[0] rather than g[:,0]) g_Estimate, h_Estimate = g_Estimate.T, h_Estimate.T return g_Estimate, h_Estimate
def magsm(p1, p2, p3, j)
-
Converts dipole (mag) to solar magnetic (sm) coordinates or vice versa j>0 j<0 input: xmag,ymag,zmag xsm, ysm, zsm output: xsm, ysm, zsm xmag,ymag,zmag
:param p1,p2,p3: input position :param j: flag :return: output position
Expand source code
def magsm( p1, p2, p3, j ): """ Converts dipole (mag) to solar magnetic (sm) coordinates or vice versa j>0 j<0 input: xmag,ymag,zmag xsm, ysm, zsm output: xsm, ysm, zsm xmag,ymag,zmag :param p1,p2,p3: input position :param j: flag :return: output position """ # Uses the 'sfi' & 'cfi' globals generated by recalc. if j > 0: xmag,ymag,zmag = p1, p2, p3 xsm = xmag*cfi - ymag*sfi ysm = xmag*sfi + ymag*cfi zsm = zmag return xsm,ysm,zsm else: xsm,ysm,zsm = p1, p2, p3 xmag = xsm*cfi + ysm*sfi ymag = ysm*cfi - xsm*sfi zmag = zsm return xmag,ymag,zmag
def recalc(ut, vxgse=-400, vygse=0, vzgse=0)
-
Expand source code
def recalc( ut, vxgse=-400,vygse=0,vzgse=0 ): # Make sure the 'ut' timestamps are given as a parsable singleton or list. ut = _Sanitize_Input( False, ut )[0] if igrf is None: init_igrf() # Documentation from Fortran's code. # 1. PREPARES ELEMENTS OF ROTATION MATRICES FOR TRANSFORMATIONS OF VECTORS BETWEEN # SEVERAL COORDINATE SYSTEMS, MOST FREQUENTLY USED IN SPACE PHYSICS. # # 2. PREPARES COEFFICIENTS USED IN THE CALCULATION OF THE MAIN GEOMAGNETIC FIELD # (IGRF MODEL) # #-----INPUT PARAMETERS: # VGSEX,VGSEY,VGSEZ - GSE (GEOCENTRIC SOLAR-ECLIPTIC) COMPONENTS OF THE OBSERVED # SOLAR WIND FLOW VELOCITY (IN KM/S) # st0/ct0 - sin/cos of teta0 (colat in geo?). sl0/cl0 - sin/cos of lambda0 (longitude in geo?). # ctcl/stcl - . ctsl/stsl - . geo x and z? # sfi/cfi/xmut - rotate angle between mag and sm and its sin/cos. # sps/cps/psi - tilt angle and its sin/cos. # shi/chi/hi - rotate angle between gse to gsm and its sin/cos. # a[11,...,33] - matrix converts geo to gsm. # cgst/sgst - cos/sin of gst. # ds3. # ba(6). global st0,ct0,sl0,cl0,ctcl,stcl,ctsl,stsl,sfi,cfi,sps,cps, \ shi,chi,hi,psi,xmut,ds3,cgst,sgst,ba, \ a11,a21,a31,a12,a22,a32,a13,a23,a33, \ e11,e21,e31,e12,e22,e32,e13,e23,e33 # The common block /geopack2/ contains coefficients of the IGRF field model, calculated # for a given year and day from their standard epoch values. the array rec contains # coefficients used in the recursion relations for legendre associate polynomials. # common /geopack2/ g(105),h(105),rec(105) # Compute the m,n related coefficients (following the notation in Davis 2004): # 1. The Schmidt quasi-normalization: S_n,m, which normalizes the associated Legendre polynomials P_n^m # to the Guassian normalized associated Legendre polynomials P^n,m. # # Since g_n^m * P_n^m should be constant, mutiplying S_n,m to g_n^m is equivalently converting P_n^m to P^n,m. # The benefit of doing so is that P^n,m follows simple recursive form, c.f. Eq (16 a-c) and Eq(17 a-b). # 2. rec[mn], which used in the recursion relation. global g,h,rec # Used in the igrf_geo() & dip() # Build the 14-degree order Legendre polynomial coefficients rec = np.array([]) for n in range( 14 ): m = np.arange( n+1 ) rec = np.append(rec, (n-m)*(n+m)/((2*n+1)*(2*n-1)) ) # WH- The 'rec' array is a 105-len array of floats with increasing subarrays of [Val*,0], with the Val-subarray starting with a zero-length entry # Eg: [ 0, 1/3,0, 4/15,1/5,0,....] # Fetch the IGRF coefficients for the given timestamps g, h = load_igrf( ut ) # Returns two lists of shape Nx105, where N is the input 'ut'-len. # Multiply the g/h IGRF values by the 'Schmidt Normalization Factors' s = 1 _Factor = [0.] for n in range(1, 14): s *= (2*n - 1) / n _Factor.append( s ) p = s for m in range(1, n+1): if m == 1: aa = 2 else: aa = 1 p *= np.sqrt( aa*(n-m+1) / (n+m)) _Factor.append( p ) # Then multiply the g/h coefficients by the array of '_Factor' g *= _Factor h *= _Factor g10=-g[:,1] g11=-g[:,2] h11=-h[:,2] # Now calculate the components of the unit vector EzMAG in the GEO system, 'parallel to geodipole axis': # sin(teta0)*cos(lambda0), sin(teta0)*sin(lambda0), and cos(teta0) # st0 * cl0 st0 * sl0 ct0 sq = g11**2 + h11**2 sqq = np.sqrt(sq) sqr = np.sqrt(g10**2+sq) sl0 = h11/sqq cl0 = g11/sqq st0 = sqq/sqr ct0 = g10/sqr stcl = st0*cl0 stsl = st0*sl0 ctsl = ct0*sl0 ctcl = ct0*cl0 # Fetch the solar data at each saught time. # slong never used. gst, slong, srasn, sdec, obliq = sun(ut) # Each component returned is an array of floats, len(ut) # All vectors are expressed in GEI. # xgse_[xyz] (s[123]) are the components of the unit vector exgsm=exgse in GEI, # pointing from the earth's center to the sun: xgse_x = np.cos(srasn) * np.cos(sdec) xgse_y = np.sin(srasn) * np.cos(sdec) xgse_z = np.sin(sdec) # zgse_[xyz] in GEI has the components (0,-sin(delta),cos(delta)) = (0.,-0.397823,0.917462); # Here delta = 23.44214 deg for the epoch 1978 (see the book by Gurevich or other astronomical handbooks). # Here the most accurate time-dependent formula is used: zgse_x = 0. zgse_y = -np.sin(obliq) #HERE zgse_z = np.cos(obliq) # ygse_[xyz] (dy[123]) = zgse_[xyz] x xgsm_[xyz] in GEI: ygse_x = zgse_y * xgse_z - zgse_z * xgse_y ygse_y = zgse_z * xgse_x - zgse_x * xgse_z ygse_z = zgse_x * xgse_y - zgse_y * xgse_x # zsm_[xyz] (dip[123]) are the components of the unit vector zsm=zmag in GEI: cgst = np.cos(gst) sgst = np.sin(gst) zsm_x = stcl * cgst - stsl * sgst zsm_y = stcl * sgst + stsl * cgst zsm_z = ct0 # xgsw_[xyz] (x[123]) in GEI. v1 = -1 / np.sqrt(vxgse*vxgse+vygse*vygse+vzgse*vzgse) xgsw_x = (vxgse*xgse_x + vygse*ygse_x + vzgse*zgse_x)*v1 xgsw_y = (vxgse*xgse_y + vygse*ygse_y + vzgse*zgse_y)*v1 xgsw_z = (vxgse*xgse_z + vygse*ygse_z + vzgse*zgse_z)*v1 # ygsw (y[123]) = zsm x xgsw in GEI. ygsw_x = zsm_y*xgsw_z - zsm_z*xgsw_y ygsw_y = zsm_z*xgsw_x - zsm_x*xgsw_z ygsw_z = zsm_x*xgsw_y - zsm_y*xgsw_x y = np.sqrt(ygsw_x*ygsw_x+ygsw_y*ygsw_y+ygsw_z*ygsw_z) ygsw_x = ygsw_x/y ygsw_y = ygsw_y/y ygsw_z = ygsw_z/y # zgsw (z[123]) = xgsw x ygsw in GEI. zgsw_x = xgsw_y*ygsw_z-xgsw_z*ygsw_y zgsw_y = xgsw_z*ygsw_x-xgsw_x*ygsw_z zgsw_z = xgsw_x*ygsw_y-xgsw_y*ygsw_x # xgsm = xgse in GEI. xgsm_x,xgsm_y,xgsm_z = xgse_x,xgse_y,xgse_z # ygsm = zsm x xgsm in GEI. ygsm_x = zsm_y*xgsm_z - zsm_z*xgsm_y ygsm_y = zsm_z*xgsm_x - zsm_x*xgsm_z ygsm_z = zsm_x*xgsm_y - zsm_y*xgsm_x y=np.sqrt(ygsm_x*ygsm_x+ygsm_y*ygsm_y+ygsm_z*ygsm_z) ygsm_x = ygsm_x/y ygsm_y = ygsm_y/y ygsm_z = ygsm_z/y # ezgsm = exgsm x eygsm in GEI. zgsm_x = xgse_y*ygsm_z-xgse_z*ygsm_y zgsm_y = xgse_z*ygsm_x-xgse_x*ygsm_z zgsm_z = xgse_x*ygsm_y-xgse_y*ygsm_x # The elements of the matrix gse to gsm are the scalar products: # chi=em22=(eygsm,eygse), shi=em23=(eygsm,ezgse), em32=(ezgsm,eygse)=-em23, and em33=(ezgsm,ezgse)=em22 chi = ygsm_x*ygse_x + ygsm_y*ygse_y + ygsm_z*ygse_z shi = ygsm_x*zgse_x + ygsm_y*zgse_y + ygsm_z*zgse_z hi = np.arcsin(shi) # elements of the matrix gsm to gsw are the scalar products: # e11 = (exgsm, exgsw) e12 = (exgsm, eygsw) e13 = (exgsm, ezgsw) # e21 = (eygsm, exgsw) e22 = (eygsm, eygsw) e23 = (eygsm, ezgsw) # e31 = (ezgsm, exgsw) e32 = (ezgsm, eygsw) e33 = (ezgsm, ezgsw) e11 = xgsm_x*xgsw_x + xgsm_y*xgsw_y + xgsm_z*xgsw_z e12 = xgsm_x*ygsw_x + xgsm_y*ygsw_y + xgsm_z*ygsw_z e13 = xgsm_x*zgsw_x + xgsm_y*zgsw_y + xgsm_z*zgsw_z e21 = ygsm_x*xgsw_x + ygsm_y*xgsw_y + ygsm_z*xgsw_z e22 = ygsm_x*ygsw_x + ygsm_y*ygsw_y + ygsm_z*ygsw_z e23 = ygsm_x*zgsw_x + ygsm_y*zgsw_y + ygsm_z*zgsw_z e31 = zgsm_x*xgsw_x + zgsm_y*xgsw_y + zgsm_z*xgsw_z e32 = zgsm_x*ygsw_x + zgsm_y*ygsw_y + zgsm_z*ygsw_z e33 = zgsm_x*zgsw_x + zgsm_y*zgsw_y + zgsm_z*zgsw_z # Solar magnetic (sm) to Geocentric Solar Magnetospheric (gsm) rotation. # Tilt angle: psi=arcsin(ezsm . exgsm) sps = zsm_x*xgse_x + zsm_y*xgse_y + zsm_z*xgse_z cps = np.sqrt(1 - sps**2) psi = np.arcsin(sps) # Not used elsewhere, but kept for consistency. # The elements of the matrix mag to sm are the scalar products: # cfi=gm22=(eysm,eymag), sfi=gm23=(eysm,exmag); They can be derived as follows: # In geo the vectors exmag and eymag have the components (ct0*cl0,ct0*sl0,-st0) and (-sl0,cl0,0), respectively. # Hence, in gei the components are: # exmag: ct0*cl0*cos(gst)-ct0*sl0*sin(gst) # ct0*cl0*sin(gst)+ct0*sl0*cos(gst) # -st0 # eymag: -sl0*cos(gst)-cl0*sin(gst) # -sl0*sin(gst)+cl0*cos(gst) # 0 # The components of eysm in gei were found above as ysm_in_geix, ysm_in_geiy, and ysm_in_geiz; # Now we only have to combine the quantities into scalar products: xmag_x = ct0*(cl0*cgst - sl0*sgst) xmag_y = ct0*(cl0*sgst + sl0*cgst) xmag_z = -st0 ymag_x = -(sl0*cgst + cl0*sgst) ymag_y = -(sl0*sgst - cl0*cgst) cfi = ygsm_x*ymag_x + ygsm_y*ymag_y sfi = ygsm_x*xmag_x + ygsm_y*xmag_y + ygsm_z*xmag_z xmut=(np.arctan2(sfi,cfi)+3.1415926536)*3.8197186342 # The elements of the matrix geo to gsm are the scalar products: # a11=(exgeo,exgsm), a12=(eygeo,exgsm), a13=(ezgeo,exgsm), # a21=(exgeo,eygsm), a22=(eygeo,eygsm), a23=(ezgeo,eygsm), # a31=(exgeo,ezgsm), a32=(eygeo,ezgsm), a33=(ezgeo,ezgsm), # All the unit vectors in brackets are already defined in gei: # xgeo=(cgst,sgst,0), ygeo=(-sgst,cgst,0), zgeo=(0,0,1) # and therefore: a11 = xgsm_x*cgst + xgse_y*sgst a12 = -xgsm_x*sgst + xgse_y*cgst a13 = xgsm_z a21 = ygsm_x*cgst + ygsm_y *sgst a22 = -ygsm_x*sgst + ygsm_y *cgst a23 = ygsm_z a31 = zgsm_x*cgst + zgsm_y*sgst a32 = -zgsm_x*sgst + zgsm_y*cgst a33 = zgsm_z return psi
def smgsm(p1, p2, p3, j)
-
Converts solar magnetic (sm) to geocentric solar magnetospheric (gsm) coordinates or vice versa. j>0 j<0 input: xsm, ysm, zsm xgsm,ygsm,zgsm output: xgsm,ygsm,zgsm xsm, ysm, zsm
:param p1,p2,p3: input position :param j: flag :return: output position
Expand source code
def smgsm( p1, p2, p3, j ): """ Converts solar magnetic (sm) to geocentric solar magnetospheric (gsm) coordinates or vice versa. j>0 j<0 input: xsm, ysm, zsm xgsm,ygsm,zgsm output: xgsm,ygsm,zgsm xsm, ysm, zsm :param p1,p2,p3: input position :param j: flag :return: output position """ # Uses the 'cps' & 'sps' values set by recalc. if j > 0: xsm,ysm,zsm = p1, p2, p3 xgsm = xsm*cps + zsm*sps ygsm = ysm zgsm = zsm*cps - xsm*sps return xgsm,ygsm,zgsm else: xgsm,ygsm,zgsm = p1, p2, p3 xsm = xgsm*cps - zgsm*sps ysm = ygsm zsm = xgsm*sps + zgsm*cps return xsm,ysm,zsm
def sphcar(p1, p2, p3, j)
-
Converts spherical coords into cartesian ones and vice versa (theta and phi in radians). j>0 j<0 input: r,theta,phi x,y,z output: x,y,z r,theta,phi
:param r,theta,phi: :param x,y,z: :param j: :return:
Note: at the poles (x=0 and y=0) we assume phi=0 (when converting from cartesian to spherical coords, i.e., for j<0) Last mofification: April 1, 2003 (only some notation changes and more comments added) Author: N.A. Tsyganenko
Expand source code
def sphcar( p1, p2, p3, j ): """ Converts spherical coords into cartesian ones and vice versa (theta and phi in radians). j>0 j<0 input: r,theta,phi x,y,z output: x,y,z r,theta,phi :param r,theta,phi: :param x,y,z: :param j: :return: Note: at the poles (x=0 and y=0) we assume phi=0 (when converting from cartesian to spherical coords, i.e., for j<0) Last mofification: April 1, 2003 (only some notation changes and more comments added) Author: N.A. Tsyganenko """ # For some blasted reason, Geopack has decied that Theta and Phi map to Lat/Lon rather than Lon/Lat, as damn near every other # translation system uses. Output matches Geopack. # if j > 0: r,theta,phi = p1, p2, p3 x = r * np.sin(theta) * np.cos(phi) y = r * np.sin(theta) * np.sin(phi) z = r * np.cos(theta) return x, y, z else: x, y, z = p1, p2, p3 r = np.sqrt( x**2 + y**2 + z**2 ) theta = np.arccos ( z / r ) phi = np.arctan2( y , x ) return r, theta, phi
def sun(ut)
-
Calculates four quantities necessary for coordinate transformations which depend on sun position (and, hence, on universal time and season) Based on http://aa.usno.navy.mil/faq/docs/SunApprox.php and http://aa.usno.navy.mil/faq/docs/GAST.php
:param ut: ut sec, can be array. :return: gst,slong,srasn,sdec. gst - greenwich mean sidereal time, slong - longitude along ecliptic srasn - right ascension, sdec - declination of the sun (radians) obliq - mean oblique of the ecliptic (radian)
Python version by Sheng Tian
Expand source code
def sun( ut ): """ Calculates four quantities necessary for coordinate transformations which depend on sun position (and, hence, on universal time and season) Based on http://aa.usno.navy.mil/faq/docs/SunApprox.php and http://aa.usno.navy.mil/faq/docs/GAST.php :param ut: ut sec, can be array. :return: gst,slong,srasn,sdec. gst - greenwich mean sidereal time, slong - longitude along ecliptic srasn - right ascension, sdec - declination of the sun (radians) obliq - mean oblique of the ecliptic (radian) Python version by Sheng Tian """ ut = np.atleast_1d(ut) # Various constants. twopi = 2*np.pi jd2000 = 2451545.0 # Reduced Julian Date of the J2K epoch. # convert to Julian date. t0_jd = 2440587.5 # in day, 0 of Julian day. secofdaygsw_x = 1./86400 # 1/sec of day. t_jd = ut*secofdaygsw_x+t0_jd # d = mjd - mj2000. d = t_jd-jd2000 d = np.squeeze(d) # mean obliquity of the ecliptic, e. (Rad) # e = 23.439 - 0.00000036*d ; in degrees. e = 0.4090877233749509 - 6.2831853e-9*d # mean anomaly of the Sun, g. # g = 357.529 + 0.98560028*d ; in degree. g = 6.2400582213628066 + 0.0172019699945780*d g = np.mod(g, twopi) # mean longitude of the Sun, q. # q = 280.459 + 0.98564736*d ; in degree. q = 4.8949329668507771 + 0.0172027916955899*d q = np.mod(q, twopi) # geocentric apparent ecliptic longitude, l. # l = q + 1.915 sin g + 0.020 sin 2g ; in degree. l = q + 0.0334230551756914*np.sin(g) + 0.0003490658503989*np.sin(2*g) # vl - q, mean longitude of the sun. # vl = np.mod(279.696678+0.9856473354*dj,360.)/rad # q = np.mod(4.881627937990388+0.01720279126623886*dj, twopi) # g, mean anomaly of the sun. # g = np.mod(358.475845+0.985600267*dj,360.)/rad # g = np.mod(6.256583784118852+0.017201969767685215*dj, twopi) # slong - l, geocentric apparent ecliptic longitude. # slong = (vl + (1.91946-0.004789*t)*sin(g) + 0.020094*sin(2*g))/rad # l = q+(0.03350089686033036-2.2884002156881157e-09*dj)*np.sin(g)+0.0003507064598957406*np.sin(2*g) # l = np.mod(l, twopi) # obliq - e, mean obliquity of the ecliptic. # obliq = (23.45229-0.0130125*t)/rad # e = 0.40931967763254096-6.217959450123535e-09*dj # sin(d) = sin(e) * sin(L) sind = np.sin(e)*np.sin(l) sdec = np.arcsin(sind) # tan(RA) = cos(e)*sin(L)/cos(L) srasn = np.arctan2(np.cos(e)*np.sin(l), np.cos(l)) srasn = np.mod(srasn, twopi) # http://aa.usno.navy.mil/faq/docs/GAST.php # gst - gmst, greenwich mean sidereal time. # gst = np.mod(279.690983+.9856473354*dj+360.*fday+180.,360.)/rad # gst = np.mod(4.881528541489487+0.01720279126623886*dj+twopi*fday+np.pi, twopi) # gmst = 18.697374558 + 24.06570982441908*d # in hour gmst = 4.894961212735792 + 6.30038809898489*d # in rad gmst = np.mod(gmst, twopi) return gmst,l,srasn,sdec,e