geopack

2024-01-29 22:15| 来源: 网络整理| 查看: 265

Python 中的 geopack 和 Tsyganenko 模型

作者：盛天，大学。明尼苏达州，[email protected]

该 pythongeopack集成了最初用 Fortran 编写的两个模块：thegeopack和 Tsyganenko 模型（T89、T96、T01 和 T04）。Fortrangeopack05可在https://ccmc.gsfc.nasa.gov/modelweb/magnetos/data-based/Geopack_2005.html和http://geo.phys.spbu.ru/~tsyganenko/Geopack-2008获得.html。他们在 IDL 中的 DLM 可在http://ampere.jhuapl.edu/code/idl_geopack.html获得。作为和的重要补充，Tsyganenko 模型在 Fortran 中可用，网址为https://ccmc.gsfc.nasa.gov/models/modelinfo.php?model=Tsyganenko%20Magnetic%20Field。geopack08geopack05geopack08

附上测试结果./test_geopack1.md以证明 Pythongeopack返回的输出与 Fortran 和 IDL 对应的输出相同。然而，在用户级别是不可见的，内部实现了一些改进：

使用了最新的 IGRF 系数，涵盖了从 1900 年到 2025 年的时间范围。超出此范围的年份是有效输入，相应的 IGRF 系数将被外推，而 Fortran 和 IDL 版本不能很好地外推。

Pythongeopack中的 IGRF 系数是毫秒节奏时间的平滑函数，而 Fortran 中的系数是每日的geopack。

igrf_gsm更改为igrf_geo加上适当的坐标变换的包装。Fortran 版本中有很多地方是复制粘贴代码页的。虽然不美观，但我还是让它们住在 python 版本中，因为修复它们需要付出巨大的努力。但是，igrf_geo 是一个明显且易于修复的地方，所以我做到了。

gotoFortrangeopack和 Tsyganenko 模型中的所有语句都被删除。

添加了Agswgsm以支持中引入的新 GSW 坐标geopack08。

安装

该包需要预先安装 Python 并依赖于numpy和scipy包。我只geopack在 Python 3.6 的 Mac OS 上测试过 Python。在其他平台和其他 Python 版本上的性能尚不清楚。

geopack要通过安装 Python ，请在终端中pip键入。> pip3 install geopack

要安装最新版本，请在 Mac（希望是 Linux）上手动安装：

在https://github.com/tsssss/geopack/下载最新包。解压，打开终端，cd进入解压目录 python3 setup.py install通过在终端中键入将包安装到 Python 关于geopack08和的注释T07d

Python 版本与 Fortrangeopack05 和geopack08. 的主要变化geopack08是一个名为的新坐标GSW，它与广泛使用GSM但更适合研究尾部物理的坐标相似。为了向后兼容geopack05，Python 版本仍然使用GSM作为向量的主坐标。不过，为了保持更新geopack08，Python版本提供了新的坐标变换功能GSWGSM，让用户可以轻松切换到自己喜欢的坐标。一个新的 TsyganenkoT07d模型已经发布了一个新的算法。对 T07d 的支持正在开发中。

G 和 W 参数注意事项

有两个 G 参数用作 T01 模型的可选输入。Tsyganenko (2001) 中有定义。同样，有六个 W 参数用作 T04 模型的可选输入，在 Tsyganenko (2005) 中定义。python版本不支持G和W参数的计算。对于感兴趣的用户，这里是https://rbsp-ect.newmexicoconsortium.org/data_pub/QinDenton/上的 Qin-Denton W 和 G 参数的链接。感谢 Shawn Young 博士提供参考资料和相关信息。

回想一下，有一些潜在的方法可以实现 G 和 W 参数。但请理解该配套没有任何资金支持。我通常在暑假或寒假期间进行重大更新，这样更容易找到空闲时间。对于熟悉 G 和 W 参数的用户，如果您对在包中实现它们的解决方案有任何建议或想法，请告诉我！

获取时间标签的示例

每个新的时间步都需要更新模型。使用的时间是 unix 时间戳，即从 1970-01-01/00:00 开始的秒数。以下是 Python 中的一些示例，用于从直观的输入中获取时间。

# 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 用法

这是一个关于如何导入包和调用函数的简短示例。下一节列出了所有功能的详细说明。

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

这是导入包和引用函数的另一种方法。

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

导入包的另一种方法。

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)

要使用中的功能geopack08，用户可以在 GSE 中提供太阳风磁场并在 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) 封装接口

Pythongeopack遵循 Python 的方式：函数参数都是输入参数并返回输出。（这与 Fortran 和 IDL 非常不同geopack。）

更改为新的感兴趣时间时

recalc. 重新计算给定时间的偶极倾斜角（和其他内部参数）。

Example ps = recalc(ut) ps = recalc(ut, vxgse,vygse,vzgse) Input ut: The given time in the universal time in second. vxgse,vygse,vzgse: The solar wind velocity in GSE. If they are omitted, a default value of [-400,0,0] is used so that GSM and GSW are the same. Return ps: Dipole tilt angle in radian (defined in GSM, not GSW).

获取内模磁场

dip. 在 GSM 坐标中计算给定位置和时间的偶极子模型的内部磁场（时间相关性由处理recalc）。

Example bxgsm,bygsm,bzgsm = dip(xgsm,ygsm,zgsm) Input xgsm,ygsm,zgsm: The given position in cartesian GSM coordinate in Re (earth radii, 1 Re = 6371.2 km). Return bxgsm,bygsm,bzgsm: Cartesian GSM components of the internal magnetic field in nT.

dip_gsw. 在 GSW 坐标中计算给定位置和时间的偶极子模型的内部磁场（时间相关性由处理recalc）。

Example bxgsw,bygsw,bzgsw = dip_gsw(xgsw,ygsw,zgsw) Input xgsw,ygsw,zgsw: The given position in cartesian GSW coordinate in Re (earth radii, 1 Re = 6371.2 km). Return bxgsw,bygsw,bzgsw: Cartesian GSW components of the internal magnetic field in nT.

igrf_gsm. 在 GSM 坐标中计算给定位置和时间的 IGRF 模型 ( http://www.ngdc.noaa.gov/iaga/vmod/igrf.html )的内部磁场。

Example bxgsm,bygsm,bzgsm = igrf_gsm(xgsm,ygsm,zgsm) Input xgsm,ygsm,zgsm: The given position in cartesian GSM coordinate in Re (earth radii, 1 Re = 6371.2 km). Return bxgsm,bygsm,bzgsm: Cartesian GSM components of the internal magnetic field in nT.

igrf_gsw. 在 GSW 坐标中计算给定位置和时间的 IGRF 模型 ( http://www.ngdc.noaa.gov/iaga/vmod/igrf.html )的内部磁场。

Example bxgsw,bygsw,bzgsw = igrf_gsw(xgsw,ygsw,zgsw) Input xgsw,ygsw,zgsw: The given position in cartesian GSW coordinate in Re (earth radii, 1 Re = 6371.2 km). Return bxgsw,bygsw,bzgsw: Cartesian GSW components of the internal magnetic field in nT.

igrf_geo. 在 GEO 坐标中计算给定位置和时间的 IGRF 模型 ( http://www.ngdc.noaa.gov/iaga/vmod/igrf.html )的内部磁场。

Example br,btheta,bphi = igrf_gsm(r,theta,phi) Input r,theta,phi: The given position in spherical GEO coordinate. r is the radia distance in Re; theta is the co-latitude in radian; phi is the longitude in radian. Return br,btheta,bphi: Spherical GSM components of the internal magnetic field in nT. br is outward; btheta is southward; bphi is eastward.

获取外部模型磁场

由 Tsyganenko 博士开发的四个模型（T89、T96、T01 和 T04）在该包中实现。

t89. 在 GSM 坐标中计算给定位置和时间的 T89 模型的外部磁场。

Example bxgsm,bygsm,bzgsm = t89(par, ps, xgsm,ygsm,zgsm) Input par: A model parameter. It is an integer (1-7) maps to the Kp index | par | 1 | 2 | 3 | 4 | 5 | 6 | 7 | | Kp | 0,0+ | 1-,1,1+ | 2-,2,2+ | 3-,3,3+ | 4-,4,4+ | 5-,5,5+ | > 6- | ps: Dipole tilt angle in radian. xgsm,ygsm,zgsm: The given position in cartesian GSM coordinate in Re (earth radii, 1 Re = 6371.2 km).

t96. 在 GSM 坐标中计算给定位置和时间的 T96 模型的外部磁场。

Example bxgsm,bygsm,bzgsm = t96(par, ps, xgsm,ygsm,zgsm) Input ps: Dipole tilt angle in radian. xgsm,ygsm,zgsm: The given position in cartesian GSM coordinate in Re (earth radii, 1 Re = 6371.2 km). par: A model paramter. It is a 10-element array, whose elements are (1-10) | par | 1 | 2 | 3-4 | 5-10 | | Var | Pdyn | Dst | ByIMF,BzIMF | not used | where Pdyn is the solar wind dynamic pressure in nPa; Dst is the Dst index in nT; ByImf,BzImf are the y and z components of the IMF (interplanetary magnetif field) in GSM.

t01. 在 GSM 坐标中计算给定位置和时间的 T01 模型的外部磁场。

Example bxgsm,bygsm,bzgsm = t01(par, ps, xgsm,ygsm,zgsm) Input ps: Dipole tilt angle in radian. xgsm,ygsm,zgsm: The given position in cartesian GSM coordinate in Re (earth radii, 1 Re = 6371.2 km). par: A model paramter. It is a 10-element array, whose elements are (1-10) | par | 1 | 2 | 3-4 | 5-6 | 7-10 | | Var | Pdyn | Dst | ByIMF,BzIMF | G1,G2 | not used | where Pdyn is the solar wind dynamic pressure in nPa; Dst is the Dst index in nT; ByImf,BzImf are the y and z components of the IMF (interplanetary magnetif field) in GSM; G1,G2 are two indices defined in Tsyganenko (2001). 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)

t04. 在 GSM 坐标中计算给定位置和时间的 T04 模型的外部磁场。

Example bxgsm,bygsm,bzgsm = t04(par, ps, xgsm,ygsm,zgsm) Input ps: Dipole tilt angle in radian. xgsm,ygsm,zgsm: The given position in cartesian GSM coordinate in Re (earth radii, 1 Re = 6371.2 km). par: A model paramter. It is a 10-element array, whose elements are (1-10) | par | 1 | 2 | 3-4 | 5-10 | | Var | Pdyn | Dst | ByIMF,BzIMF | W1 to W6 | where Pdyn is the solar wind dynamic pressure in nPa; Dst is the Dst index in nT; ByImf,BzImf are the y and z components of the IMF (interplanetary magnetif field) in GSM; W1,W2,...,W6 are six indices defined in Tsyganenko (2005). 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.

注意：所有 4 个模型共享相同的接口，但含义par有很大不同。

在坐标之间转换笛卡尔向量

支持的坐标有：GEO、GEI、MAG、GSM、GSE 和 SM。它们在 Hapgood (1992) 中定义。并且在 Hones+(1986) 中定义的 GSW 被添加到geopack_08. 坐标变换的函数是： geomag, geigeo, magsm, gsmgse, smgsm, geogsm, gswgsm. 它们共享相同的界面，因此将它们一起解释。

Usage b1,b2,b3 = geomag(h1,h2,h3, flag) Example xmag,ymag,zmag = geomag(xgeo,ygeo,zgeo, 1) xgeo,ygeo,zgeo = geomag(xmag,ymag,zmag, -1) ... Input and Return h1,h2,h3: Cartesian components of a vector in "coord1" b1,b2,b3: Cartesian components of the vector in "coord2" flag: flag > 0 -- coord1 to coord2; flag < 0 -- coord2 to coord1

此外geodgeo，在海拔（千米）/大地纬度（弧度）和地心距离（千米）/纬度（弧度）之间转换位置。

Usage b1,b2 = geodgeo(h1,h2, flag) Example rgeo,thetageo = geodgeo(hgeod,xmugeod, 1) hgeod,xmugeod = geodgeo(rgeo,thetageo, -1) Input and Return h1,h2: Components of a vector in "coord1" b1,b2: Components of a vector in "coord2" flag: flag > 0 -- coord1 to coord2; flag < 0 -- coord2 to coord1

沿模型磁场追踪：trace

Example x1gsm,y1gsm,z1gsm = trace(x0gsm,y0gsm,z0gsm, dir, rlim, r0, par, exname, inname) Input x0gsm,y0gsm,z0gsm: The given position in cartesian GSM coordinate in Re (earth radii, 1 Re = 6371.2 km). dir: Direction of tracing. dir = -1 for parallel; dir = 1 for anti-parallel. rlim: Maximum tracing radius in Re. Default value is 10 Re. r0: Minimum tracing radius in Re. Default value is 1 Re. inname: A string specifies the internal model, one of 'dipole','igrf'. The default value is 'igrf'. exname: A string specifies the external model, one of 't89','t96','t01','t04'. The default value is 't89' and its par is default to be 2. par: The model parameter. Its dimension and the meaning depend on the external model. Please check the interface of the models for details.

未出现在上述列表中的函数被视为内部函数。对于它们的用法，高级用户可以查看 Python 的源代码geopack。

参考

马萨诸塞州哈普古德 (1992)。空间物理坐标变换：用户指南。行星与空间科学，40（5），711-717。http://doi.org/10.1016/0032-0633(92)90012-D

NA Tsyganenko，基于数据的新近磁层磁场模型：1. 数学结构。2. 参数化和拟合观测（提交给 JGR，2001 年 7 月）

NA Tsyganenko 和 MI Sitnov，模拟强地磁暴期间内磁层的动力学，J. Geophys。Res., v. 110 (A3), A03208, doi: 10.1029/2004JA010798, 2005。

【本文地址】

公司简介

联系我们