TsyganenkoModels

Modeling of Earth's Magnetosphere Using Spacecraft Magnetometer Data.

Features

• Magnetic field model (native Julia implementations)
  • [x] T89: A magnetospheric magnetic field model with a warped tail current sheet
  • [x] T96: Effects of the solar wind conditions on the global magnetospheric configuration
  • [x] T01/T02: A model of the near magnetosphere with a dawn-dusk asymmetry
  • [x] TS05/TS04: a dynamical empirical model of the inner storm-time magnetosphere

Installation

using Pkg
Pkg.add("TsyganenkoModels")

Quickstart

The model interface allows you to configure a model once and use it for multiple field calculations:

using TsyganenkoModels
using Dates

# Create model configurations
model_t89 = T89(2)  # Kp level 2
param = (; pdyn=2.0, dst=-87.0, byimf=2.0, bzimf=-5.0)
# pdyn: Solar wind dynamic pressure [nPa]
# dst: Dst index [nT]
# byimf: IMF By [nT]
# bzimf: IMF Bz [nT]

# Calculate fields at position
t = DateTime("1970-01-01T00:01:40")
𝐫 = [1, 2, 3]
ps = -0.533585131  # dipole tilt angle [radians]

# Using dipole tilt angle
B_t89 = T89(2)(𝐫, ps)

3-element GeoCotrans.CoordinateVector{GSM, SpaceDataModel.Cartesian3, Float64, Nothing} with indices SOneTo(3):
   5.367416991592773
  -3.448905444772822
 -17.29126736398316

julia> # Compare with other models
       B_t96 = T96(param)(𝐫, ps)3-element GeoCotrans.CoordinateVector{GSM, SpaceDataModel.Cartesian3, Float64, Nothing} with indices SOneTo(3):
  39.34888430483087
  -9.327382852686146
 -90.39880058123812
julia> B_t01 = T01(param)(𝐫, ps)3-element GeoCotrans.CoordinateVector{GSM, SpaceDataModel.Cartesian3, Float64, Nothing} with indices SOneTo(3):
  25.017949957632542
  -4.606333301128644
 -66.82317110266318
julia> B_ts04 = TS04(param)(𝐫, ps)3-element GeoCotrans.CoordinateVector{GSM, SpaceDataModel.Cartesian3, Float64, Nothing} with indices SOneTo(3):
  13.497153413146576
  -0.13916234543149042
 -24.928503548859755
julia> # Using time (auto-calculates dipole tilt)
       T89(3)(𝐫, t)3-element GeoCotrans.CoordinateVector{GSM, SpaceDataModel.Cartesian3, Float64, Nothing} with indices SOneTo(3):
   7.173912838767684
  -3.7801244076510705
 -19.9277123686887
julia> # Using `TsyIGRF` to combine IGRF14 model and one Tsyganenko model (default T89)
       TsyIGRF()(𝐫, t)3-element GeoCotrans.CoordinateVector{GSM, SpaceDataModel.Cartesian3, Float64, Nothing} with indices SOneTo(3):
 -541.7175541221735
 -576.5916824422504
 -339.77872971738
julia> # This is equivalent to
       IGRF()(GSM(𝐫), t) .+ T89(3)(𝐫, t)3-element GeoCotrans.CoordinateVector{GSM, SpaceDataModel.Cartesian3, Float64, Nothing} with indices SOneTo(3):
 -541.7175541221735
 -576.5916824422504
 -339.77872971738

Comparison with geopack

using TsyganenkoModels
using GeoCotrans
using Test
using Geopack
using Dates

t = DateTime("1970-01-01T00:01:40")
𝐫 = [1, 2, 3]

# using geopack
xgsm, ygsm, zgsm = 𝐫
opt = 2
ut = datetime2unix(t)
ps = Geopack.recalc(ut)
db_py = Geopack.t89(opt, ps, xgsm, ygsm, zgsm)

# using TsyganenkoModels
db_jl = T89(opt)(𝐫, ps)

@test db_py ≈ db_jl

Test Passed

Model Result Comparison

# T01 model comparison
pdyn = 2.0   # Solar wind dynamic pressure [nPa]
dst = -87.0  # Dst index [nT]
byimf = 2.0  # IMF By [nT]
bzimf = -5.0 # IMF Bz [nT]
g1 = 0.0
g2 = 0.0
parmod = [pdyn, dst, byimf, bzimf, g1, g2]

db_t96_py = Geopack.t96(parmod, ps, xgsm, ygsm, zgsm)
db_t96_jl = T96(; pdyn, dst, byimf, bzimf)(𝐫, ps)
@test db_t96_jl ≈ db_t96_py rtol = 1e-6

db_t01_py = Geopack.t01(parmod, ps, xgsm, ygsm, zgsm)
db_t01_jl = T01(; pdyn, dst, byimf, bzimf)(𝐫, ps)

@test db_t01_jl ≈ db_t01_py rtol = 1e-6

Test Passed

The internal magnetic field could be computed using igrf from GeoCotrans.

Here we verify the result with igrf_gsm from Geopack.

julia> using GeoCotrans
julia> b0_py = Geopack.igrf_gsm(xgsm, ygsm, zgsm)3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 -548.892407931754
 -572.8105422196072
 -319.8493868519125
julia> b0_jl = GeoCotrans.igrf(GSM(𝐫), t)3-element GeoCotrans.CoordinateVector{GSM, SpaceDataModel.Cartesian3, Float64, Dates.DateTime} with indices SOneTo(3):
 -548.8914669609412
 -572.8115580345993
 -319.8510173486913
julia> @test b0_jl ≈ b0_py rtol = 1e-5Test Passed
julia> @test GeoCotrans.get_igrf_coeffs(t)[1] ≈ Geopack.load_igrf(ut)[1]Test Passed

Performance Benchmarks

julia> using Chairmarks
julia> @b TsyganenkoModels.T89(2)($𝐫, $ps), Geopack.t89(2, $ps, $𝐫...)(180.000 ns, 51.275 μs (25 allocs: 624 bytes))
julia> @b TsyganenkoModels.T96($pdyn, $dst, $byimf, $bzimf)($𝐫, $ps), Geopack.t96([$pdyn, $dst, $byimf, $bzimf, 0, 0], $ps, $𝐫...)(6.342 μs, 1.193 ms (101 allocs: 2.812 KiB))
julia> @b TsyganenkoModels.T01(; pdyn, dst, byimf, bzimf)($𝐫, $ps), Geopack.t01([pdyn, dst, byimf, bzimf, 0, 0], $ps, $𝐫...)(17.964 μs (5 allocs: 208 bytes), 2.481 ms (148 allocs: 4.078 KiB))
julia> @b GeoCotrans.igrf(GSM($𝐫), t), Geopack.igrf_gsm($𝐫...)(2.795 μs (2 allocs: 80 bytes), 133.779 μs (11 allocs: 192 bytes))

API

• TsyganenkoModels.T01
• TsyganenkoModels.T89
• TsyganenkoModels.T96
• TsyganenkoModels.TS04
• TsyganenkoModels.TsyIGRF
• TsyganenkoModels.TsyganenkoModel
• TsyganenkoModels.dipole_tilt
• TsyganenkoModels.t01
• TsyganenkoModels.t89
• TsyganenkoModels.t96

T01(; pdyn, dst, byimf, bzimf, g1=0.0, g2=0.0)

model = T01(pdyn=2.0, dst=-87.0, byimf=2.0, bzimf=-5.0)
B = model([1, 2, 3], ps)  # Compute field at position

T89(iopt::Integer)

model = T89(2)  # Kp level 2
B = model([1, 2, 3], ps)  # Compute field at position

T96(; pdyn, dst, byimf, bzimf)

model = T96(pdyn=2.0, dst=-87.0, byimf=2.0, bzimf=-5.0)
B = model([1, 2, 3], ps)  # Compute field at position

TS04(; pdyn, dst, byimf, bzimf, w1=0.0, w2=0.0, w3=0.0, w4=0.0, w5=0.0, w6=0.0)

model = TS04(pdyn=2.0, dst=-87.0, byimf=2.0, bzimf=-5.0, w1=0.5)
B = model([1, 2, 3], ps)  # Compute field at position

TsyIGRF(model = T89(iopt=3); name=nothing)

model = TsyIGRF(T89(iopt=3))

# Evaluate at position (in Earth radii, spherical coordinates)
r, θ, φ = 6.6, π/2, 0.0  # 6.6 RE at equator
t = Date(2020, 6, 21)    # Summer solstice 2020
B = model(r, θ, φ, t)    # Returns B in spherical coordinates (Br, Bθ, Bφ) in nT

(m::TsyganenkoModel)(x, y, z, ps, args...; kw...) -> (Bx, By, Bz)

dipole_tilt(time) -> ψ

t01(x, y, z, ps, pdyn, dst, byimf, bzimf, g1 = 0.0, g2 = 0.0) -> (Bx, By, Bz)

t89(𝐫, ps, iopt; cache=nothing) -> GSM(Bx, By, Bz)
t89(x, y, z, ps, iopt)
t89(x, y, z, t::AbstractTime, iopt)

t96(x, y, z, ps, pdyn, dst, byimf, bzimf) -> (Bx, By, Bz)