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