search

tsssss / geopack

Python version of geopack and Tsyganenko models
MIT License
30 stars 12 forks source link
geopack physics python python3 space-physics tsyganenko-model

readme

The geopack and Tsyganenko models in Python

Author: Sheng Tian, UCLA, ts0110@atmos.ucla.edu

This python geopack has integrated two modules originally written in Fortran: the geopack and the Tsyganenko models (T89, T96, T01, and T04). The Fortran geopack05 is available at https://ccmc.gsfc.nasa.gov/modelweb/magnetos/data-based/Geopack_2005.html and geopack08 is available at http://geo.phys.spbu.ru/~tsyganenko/Geopack-2008.html. Their DLM in IDL is available at http://ampere.jhuapl.edu/code/idl_geopack.html. As a crucial complement to geopack05 and geopack08, the Tsyganenko models are available in Fortran at https://ccmc.gsfc.nasa.gov/models/modelinfo.php?model=Tsyganenko%20Magnetic%20Field.

Test results are attached in ./test_geopack1.md to demonstrate that the Python geopack returns the same outputs as the Fortran and IDL counterparts. However, invisible at the user level, several improvements have been internally implemented:

  1. The latest IGRF coefficients are used, which cover the time range from 1900 to 2025. Years beyond this range are valid inputs and the corresponding IGRF coefficients will be extrapolated, whereas the Fortran and IDL versions do not extrapolate well if at all.

  2. The IGRF coefficients in the Python geopack are time series at a milli-second cadence, whereas the coefficients are daily in the Fortran geopack.

  3. igrf_gsm is changed to a wrapper of igrf_geo plus the proper coordinate transforms. There are many places in the Fortran version where pages of codes are copied and pasted. Though not aesthetically pleasing, I let them live in the Python version, because it requires tremendous effort to fix them all. However, the igrf_geo is the one place that is obvious and easy to fix, so I did it.

  4. All goto statements in the Fortran geopack and Tsyganenko models are eliminated.

  5. A gswgsm is added to support the new GSW coordinate introduced in geopack08.

Installation

The package requires Python pre-installed and depends on the numpy and scipy packages. I've only tested the Python geopack on Mac OS in Python 3.6. Performance on other platform and other versions of Python is unclear.

To install the Python geopack through pip, type > pip3 install geopack in the terminal.

To install the latest version, manually install it on a Mac (and hopefully Linux):

  1. Download the latest package at https://github.com/tsssss/geopack/.
  2. Unzip it, open a terminal, and cd to the unzipped directory
  3. Install the package to Python by typing python3 setup.py install in the terminal

Notes on geopack08 and T07d

The Python version of geopack tries to be compatible with both Fortran geopack05 and geopack08. The major change of geopack08 is a new coordinate called GSW, which is similar to the widely used GSM but more suitable for studying tail physics. To be backward compatible with geopack05, the Python version still uses GSM as the major coordinate for vectors. However, to keep updated with geopack08, the Python version provides a new coordinate transform function GSWGSM, so that users can easily switch to their favorite coordinate. A new Tsyganenko T07d model has been released with a new algorithm. Support for T07d is under development.

Notes on the G and W parameters

There are two G parameters used as optional inputs to the T01 model. Their definitions are in Tsyganenko (2001). Similarly, there are six W parameters used as optional inputs to the T04 model, as defined in Tsyganenko (2005). The Python version does not support the calculations of the G and W parameters. For users interested, here is the link for the Qin-Denton W and G parameters at https://rbsp-ect.newmexicoconsortium.org/data_pub/QinDenton/. Thanks to Dr Shawn Young for providing the references and relevant information.

Back in my mind, there are some potential ways to implement the G and W parameters. But please do understand that the package does not have any funding support. I usually do major updates during summer or winter break, when it's easier to find spare time. For users who are familiar with the G and W parameters, let me know if you have any suggestions or ideas on solutions to implement them in the package!

Example of getting the time tag

The model needs to be updated for each new time step. The time used is the unix timestamp, which is the seconds from 1970-01-01/00:00. Here are some examples in Python to get the time from intuitive inputs.

# Test for 2001-01-02/03:04:05 UT
import datetime
from dateutil import parser

# From date and time
t1 = datetime.datetime(2001,1,2,3,4,5)
t0 = datetime.datetime(1970,1,1)
ut = (t1-t0).total_seconds()
print(ut)
978404645.0

# From string, need the package dateutil
t1 = parser.parse('2001-01-02/03:04:05')
ut = (t1-t0).total_seconds()
print(ut)
978404645.0

Usage

Here is a short example of importing the package and call functions. A detailed explanation of all functions is listed in the next section.

from geopack import geopack, t89

ut = 100    # 1970-01-01/00:01:40 UT.
xgsm,ygsm,zgsm = [1,2,3]
ps = geopack.recalc(ut)
b0xgsm,b0ygsm,b0zgsm = geopack.dip(xgsm,ygsm,zgsm)          # calc dipole B in GSM.
dbxgsm,dbygsm,dbzgsm = t89.t89(2, ps, xgsm,ygsm,zgsm)       # calc T89 dB in GSM.
bxgsm,bygsm,bzgsm = [b0xgsm+dbxgsm,b0ygsm+dbygsm,b0zgsm+dbzgsm]
print(bxgsm,bygsm,bzgsm)
-539.5083883330017 -569.5906371610358 -338.8680547453352

And here is another way to import the package and refer to the functions.

import geopack

ut = 100    # 1970-01-01/00:01:40 UT.
xgsm,ygsm,zgsm = [1,2,3]
ps = geopack.geopack.recalc(ut)
b0xgsm,b0ygsm,b0zgsm = geopack.geopack.dip(xgsm,ygsm,zgsm)
dbxgsm,dbygsm,dbzgsm = geopack.t89.t89(2, ps, xgsm,ygsm,zgsm)
print(b0xgsm,b0ygsm,b0zgsm)
-544.425907831383 -565.7731166717405 -321.43413443108597

Another way to import the package.

import geopack.geopack as gp

ut = 100    # 1970-01-01/00:01:40 UT.
xgsm,ygsm,zgsm = [2,1,100]
ps = gp.recalc(ut)
xgsm,ygsm,zgsm = gp.geogsm(2,1,100, 1)
print(xgsm,ygsm,zgsm)
(-41.00700906453125, -19.962123759781406, 89.0221254665413)

To use the feature in geopack08, users can supply the solar wind magnetic field in GSE and express vectors in GSW

from geopack import geopack

ut = 100    # 1970-01-01/00:01:40 UT.
xgsm,ygsm,zgsm = [1,2,3]
vgse = [-400,0,10]                                       # solar wind velocity in GSE.
ps = geopack.recalc(ut, *vgse)                           # init with time & SW velocity.
# or use ps = geopack.recalc(ut, vgse[0],vgse[1],vgse[2])
xgsw,ygsw,zgsw = gswgsm(xgsm,ygsm,zgsm, -1)              # convert position to GSW.
b0xgsw,b0ygsw,b0zgsw = geopack.dip_gsw(xgsw,ygsw,zgsw)   # calc dipole B in GSW.
print(b0xgsw,b0ygsw,b0zgsw)
-540.5392569443875 -560.7296994901754 -336.47913346240205
print((geopack.gswgsm(b0xgsw,b0ygsw,b0zgsw, 1)))         # dipole B in GSM.
(-544.4259078313833, -565.7731166717405, -321.4341344310859)

Package Interface

The Python geopack follows the Python way: function parameters are all input parameters and the outputs are returned. (This is very different from the Fortran and IDL geopack.)

Functions do not appear in the above list are considered as internal functions. For usages of them, advanced users can check the source code of the Python geopack.

References

Hapgood, M. A. (1992). Space physics coordinate transformations: A user guide. Planetary and Space Science, 40(5), 711–717. http://doi.org/10.1016/0032-0633(92)90012-D

N. A. Tsyganenko, A new data-based model of the near magnetosphere magnetic field: 1. Mathematical structure. 2. Parameterization and fitting to observations (submitted to JGR, July 2001)

N. A. Tsyganenko and M. I. Sitnov, Modeling the dynamics of the inner magnetosphere during strong geomagnetic storms, J. Geophys. Res., v. 110 (A3), A03208, doi: 10.1029/2004JA010798, 2005.