Coordinate transformations and field line tracing

GeoCotrans.jl compared with IRBEM and PySPEDAS

GeoCotrans.jl is a native Julia implementation of geocentric coordinate transformations and IGRF field line tracing. Here we check it against IRBEM.jl (Fortran) and PySPEDAS (Python), and benchmark the three.

NoteTakeaway

GeoCotrans agrees closely with IRBEM and PySPEDAS and is an order of magnitude faster. It interpolates IGRF coefficients and the Sun’s direction more finely, which accounts for the small differences.

Code 
using GeoCotrans, IRBEM
using PySPEDAS, PySPEDAS.PythonCall
using DimensionalData, Dates, Chairmarks, Test
using OrdinaryDiffEqTsit5

irbem_cotrans(pos, in, out) = coord_trans(in, out, collect(dims(pos, Ti)), parent(pos)')'

THEMIS-A positions from the PySPEDAS test cases:

Code 
@py import pyspedas.cotrans_tools.tests.test_cotrans: CotransTestCases
CotransTestCases().test_cotrans()

pyspedas.projects.themis.state(["2010-02-25/00:00:00", "2010-02-25/23:59:59"], probe = "a", time_clip = true)
tha_pos = set(DimArray(PySPEDAS.get_data("tha_pos")), Dim{:time} => Ti)
tha_pos_gse = set(DimArray(PySPEDAS.get_data("tha_pos_gse")), Dim{:time} => Ti)
pyimport("logging").getLogger().setLevel("WARNING")

Agreement

Code 
jl = gei2geo(tha_pos)
py = PySPEDAS.get_data("tha_pos_new_geo")
@test jl ≈ parent(py)
@test jl ≈ irbem_cotrans(tha_pos, "GEI", "GEO")
Test Passed
Code 
jl = gei2gsm(tha_pos)
pyspedas.cotrans("tha_pos", "tha_pos_new_gsm", coord_in = "GEI", coord_out = "GSM")
@test jl ≈ parent(PySPEDAS.get_data("tha_pos_new_gsm")) rtol = 1e-5
@test jl ≈ irbem_cotrans(tha_pos, "GEI", "GSM") rtol = 1e-3
Test Passed
Code 
jl = gse2gsm(tha_pos_gse)
pyspedas.cotrans("tha_pos_gse", "tha_pos_new_gsm", coord_in = "GSE", coord_out = "GSM")
@test jl ≈ parent(PySPEDAS.get_data("tha_pos_new_gsm")) rtol = 1e-5
@test jl ≈ irbem_cotrans(tha_pos_gse, "GSE", "GSM") rtol = 1e-3
Test Passed

Benchmark

GeoCotrans, IRBEM, PySPEDAS:

Code 
@b gei2geo($tha_pos), irbem_cotrans($tha_pos, "GEI", "GEO"), pyspedas.cotrans("tha_pos", "tha_pos_new_geo", coord_in = "GEI", coord_out = "GEO")
(17.583 μs (7 allocs: 48.797 KiB), 964.042 μs (32 allocs: 132.875 KiB), 1.133 ms (16 allocs: 336 bytes))
Code 
@b gei2gsm($tha_pos), irbem_cotrans($tha_pos, "GEI", "GSM"), pyspedas.cotrans("tha_pos", "tha_pos_new_gsm", coord_in = "GEI", coord_out = "GSM")
(116.083 μs (7 allocs: 48.797 KiB), 1.870 ms (32 allocs: 132.875 KiB), 3.395 ms (16 allocs: 336 bytes))
Code 
@b gse2gsm($tha_pos_gse), irbem_cotrans($tha_pos_gse, "GSE", "GSM"), pyspedas.cotrans("tha_pos_gse", "tha_pos_new_gsm", coord_in = "GSE", coord_out = "GSM")
(105.708 μs (7 allocs: 48.797 KiB), 1.905 ms (32 allocs: 132.875 KiB), 2.388 ms (16 allocs: 336 bytes))

Field line tracing and magnetic equator

GeoCotrans.find_magequator traces the field line through a point and locates the |B| minimum, like IRBEM’s find_magequator with IGRF only (kext = 0).

Code 
t = DateTime(2020, 1, 1)
pos = GEO(3.0, 0.5, 1.2)
jl_eq = GeoCotrans.find_magequator(pos, t, Tsit5())
ir_eq = IRBEM.find_magequator(t, pos; kext = 0)
@test jl_eq.pos ≈ ir_eq.XGEO rtol = 1e-3
@test jl_eq.Bmin ≈ ir_eq.Bmin rtol = 1e-3
Test Passed
Code 
@b GeoCotrans.find_magequator($pos, $t, Tsit5()), IRBEM.find_magequator($t, $pos; kext = 0)
(17.333 μs (72 allocs: 4.906 KiB), 85.542 μs (54 allocs: 2.250 KiB))

Tracing the whole field line:

Code 
@b trace($pos, $t, Tsit5()), IRBEM.trace_field_line($t, $pos; kext = 0)
(11.958 μs (68 allocs: 7.094 KiB), 1.181 ms (71 allocs: 117.672 KiB))
Back to top