GitHub

Field Line Tracing

Basic Usage

using GeoCotrans, OrdinaryDiffEqTsit5, Dates

# Define time for field evaluation
t = DateTime(2020, 1, 1)

# Trace a field line starting from [3, 0, 0] Earth radii (GEO coordinates)
sol = trace(GEO(3.0, 0.0, 0.0), t, Tsit5())
retcode: Terminated
Interpolation: specialized 4th order "free" interpolation
t: 12-element Vector{Float64}:
 0.0
 0.0009999885482375654
 0.029020707803271383
 0.14198902734225044
 0.3538815232377085
 0.6456004378522566
 1.0275240680500026
 1.5121468182029143
 2.1038792660290904
 2.7979652750977886
 3.205958564065051
 3.205958564065051
u: 12-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [3.0, 0.0, 0.0]
 [3.000024874434897, -0.00014500231793512463, 0.0009891069741289906]
 [3.0003221759221472, -0.004199525924009845, 0.028712445065702575]
 [2.993710455916287, -0.020325429467939554, 0.14026976756264384]
 [2.948291901630456, -0.049090843155660445, 0.34484880212466446]
 [2.8211600666510206, -0.08386276728637002, 0.6041330674995906]
 [2.563770907700326, -0.11821263063834815, 0.8820714803865584]
 [2.1403030907166247, -0.14156972113591745, 1.1116464607902048]
 [1.558300338510598, -0.14084784845796552, 1.1991042585184868]
 [0.8863678177637029, -0.10566984174806507, 1.0478584134473157]
 [0.5381603392062427, -0.07300839837305084, 0.8396744744675901]
 [0.5381603392062427, -0.07300839837305084, 0.8396744744675901]

Plotting with Makie

Here's a complete example showing how to trace and visualize magnetic field lines using GLMakie:

using GeoCotrans, OrdinaryDiffEqTsit5, Dates
using CairoMakie

# Time for field evaluation
t = DateTime(2020, 1, 1)

# Create figure
fig = Figure(size = (800, 800))
ax = Axis3(fig[1, 1],
    xlabel = "X (Re)",
    ylabel = "Y (Re)",
    zlabel = "Z (Re)",
    title = "IGRF Field Lines at $(Date(t))",
    aspect = :data
)

# Draw Earth (unit sphere)
θ_sphere = range(0, π, length=30)
φ_sphere = range(0, 2π, length=60)
x_sphere = [sin(θ) * cos(φ) for θ in θ_sphere, φ in φ_sphere]
y_sphere = [sin(θ) * sin(φ) for θ in θ_sphere, φ in φ_sphere]
z_sphere = [cos(θ) for θ in θ_sphere, φ in φ_sphere]
surface!(ax, x_sphere, y_sphere, z_sphere, color = :lightblue, alpha = 0.8)

# Trace field lines from different starting positions
L_values = [2.0, 3.0, 4.0, 5.0]  # L-shell values
colors = [:red, :orange, :green, :blue]

for (L, color) in zip(L_values, colors)
    # Trace in both directions from the equator
    for dir in [1, -1]
        sol = trace(GEO(L, 0.0, 0.0), t, Tsit5(); dir = dir)
        # Plot the field line
        plot!(ax, sol; color = color, linewidth = 2, idxs = (1, 2, 3),
            label = dir == 1 ? "L = $L" : nothing)
    end
end

axislegend(ax, position = :rt) # Add legend
limits!(ax, -6, 6, -6, 6, -6, 6) # Set view limits

fig
Example block output

2D Visualization

using GeoCotrans, OrdinaryDiffEqTsit5, Dates
using CairoMakie

t = DateTime(2020, 1, 1)

fig = Figure(size = (1000, 500))

# Left panel: Meridional plane (X-Z)
ax1 = Axis(fig[1, 1],
    xlabel = "X (Re)", ylabel = "Z (Re)",
    title = "Field Lines in Meridional Plane",
    aspect = DataAspect()
)

# Draw Earth cross-section
θ_circle = range(0, 2π, length=100)
lines!(ax1, cos.(θ_circle), sin.(θ_circle), color = :lightblue, linewidth = 3)

# Trace field lines
for L in [2.0, 3.0, 4.0, 5.0, 6.0]
    for dir in [1, -1]
        sol = trace(GEO(L, 0.0, 0.0), t, Tsit5(); dir = dir)
        lines!(ax1, sol, idxs = (1, 3), color = :darkblue, linewidth = 1.5)
    end
end

limits!(ax1, -7, 7, -7, 7)

# Right panel: Equatorial plane (X-Y)
ax2 = Axis(fig[1, 2],
    xlabel = "X (Re)", ylabel = "Y (Re)",
    title = "Field Lines in Equatorial Plane",
    aspect = DataAspect()
)

# Draw Earth cross-section
lines!(ax2, cos.(θ_circle), sin.(θ_circle), color = :lightblue, linewidth = 3)

# Starting points slightly off equator to see field line structure
for L in [3.0, 4.0, 5.0]
    for φ in range(0, 2π, length=13)[1:12]
        y0, x0 = L .* sincos(φ)
        # Start slightly above equator
        sol = trace(GEO(x0, y0, 0.1), t, Tsit5(); dir = 1)
        lines!(ax2, sol, idxs = (1, 2), color = :darkgreen, linewidth = 1)
    end
end

limits!(ax2, -7, 7, -7, 7)

fig
Example block output

API Reference

GeoCotrans.FieldLineTracing.trace — Function
trace(pos, t, solver; kwargs...) :: ODESolution

Trace a magnetic field line through pos (Cartesian, Re) at time t with a SciML ODE solver, parameterized by arc length.

Keyword Arguments

  • model = IGRF(): Magnetic field model to use
  • dir = 1: Tracing direction (+1 for parallel to B, -1 for anti-parallel)
  • r0 = 1.0: Inner radial boundary (Earth radii)
  • rlim = 10.0: Outer radial boundary (Earth radii)
  • maxs = 100.0: Maximum arc length for integration
  • from = getcsys(pos): Input coordinate system (Reference frame and coordinate representation)
  • callback = nothing: extra SciML callback(s), combined with the boundary callbacks; u is in the model's native frame
  • Additional keyword arguments are passed to solve()

Example

using GeoCotrans, OrdinaryDiffEqTsit5, Dates
sol = trace(GEO(3.0, 0.0, 0.0), DateTime(2020, 1, 1), Tsit5())
sol.u  # positions [Re] along the line; sol.t is arc length
source
GeoCotrans.FieldLineTracing.find_magequator — Function
find_magequator(pos, t, solver; kwargs...) -> (; pos, Bmin, s) or nothing

Find the magnetic equator, the |B| minimum along the field line through pos. Traces downhill in |B| and stops at the root of d|B|/ds (a ContinuousCallback), so the result is the nearest local minimum; s is the signed arc length from pos (positive along B). Returns nothing if a boundary (r0, rlim, maxs) is reached first (open field line).

Keyword arguments are those of trace; raise rlim for nightside lines in external models, whose equator can lie beyond 10 Re.

Example

using GeoCotrans, OrdinaryDiffEqTsit5, Dates
eq = find_magequator(GEO(3.0, 0.0, 0.5), DateTime(2020, 1, 1), Tsit5())
eq.pos, eq.Bmin
source