Ionosphere, Plasmasphere, and Ray Tracing

Purpose

This specification defines the mathematical contract for the ionospheric and plasmaspheric group-delay model applied to terrestrial GNSS signals received by lunar (cislunar) users. Unlike near-Earth links, lunar GNSS paths traverse long, near-horizontal trajectories through the ionosphere and plasmasphere; the delay is therefore computed by numerically ray-tracing the bent signal path and integrating the electron density along it.

The primary modeling reference is Keidai Iiyama and Grace Gao, “Ionospheric and Plasmaspheric Delay Characterization for Lunar Terrestrial GNSS Receivers with Global Core Plasma Model,” submitted to NAVIGATION: Journal of the Institute of Navigation, arXiv:2510.10059 (2025), and the author’s PhD thesis Chapter 5. Equation numbers (5.x) below refer to that chapter.

The ray tracer lives in the pecsim namespace under cpp/lupnt/environment/plasma/tec/raytrace.{h,cc} (pecsim::trace_ray, RayTraceConfig, PathProfile). Electron density comes from GCPM v2.4 (cpp/lupnt/environment/plasma/gcpm/gcpm_interface.{h,cc}) wrapping the IRI ionosphere (gcpm/iri_interface.{h,cc}); the magnetic field comes from IGRF-14 (cpp/lupnt/environment/plasma/igrf/igrf_interface.{h,cc}). The delay is consumed by the GNSS measurement model through GNSSMeasurements::ComputeIonospherePlasmaRayTraceDelay (cpp/lupnt/measurements/gnss_measurement.cc), specified as \(\Delta\rho_\mathrm{plasma}\) in GNSS Measurement Model.

Note

The entire plasma module is km-based. Positions x1, x2 passed to trace_ray are in km, the step size and cutoff radius are in km, the Earth radius is RE = 6371.0 km, and the speed of light is C = 299792.458 km/s (cpp/lupnt/environment/plasma/core/constants.h). Electron density is evaluated internally in cm-3 (GCPM) and converted to m-3 before entering the delay integrals; one TEC unit is TECU = 1e16 electrons/m2. Final delays are returned in meters. Callers such as GNSSMeasurements divide their SI (meter) positions by 1000 before calling trace_ray.

Coordinate, Unit, and Time Contract

The ray tracer works in a geocentric Earth-fixed Cartesian frame. The GNSS coupling converts receiver/transmitter states from options.frame into raytrace_frame (default Frame::ECEF) and passes an epoch in raytrace_epoch_scale (default Time::UTC), expressed as seconds from J2000:

\[x_1 = r_S^\mathrm{ecef}\,[\mathrm{km}], \qquad x_2 = r_R^\mathrm{ecef}\,[\mathrm{km}], \qquad t = t_\mathrm{rx}^\mathrm{UTC}\ (\text{s since J2000}).\]

Electron-density evaluation converts the geocentric position into solar magnetic (SM) coordinates required by GCPM (cpp/lupnt/environment/plasma/tec/raytrace.cc :: compute_ne):

std::array<double,3> pos_gei_arr = {pos_geo[0]/RE, pos_geo[1]/RE, pos_geo[2]/RE};
geo_to_sm(itime, pos_gei_arr, pos_sm);       // geographic -> solar magnetic
cart_to_pol(pos_sm, alatr, along, r);        // r in RE, alatr in rad
double amlt = 12.0 + along / AMLTRAD;         // magnetic local time [hours]
if (amlt > 24.0) amlt -= 24.0;
out = gcpm_v24_fortran(datetime, r, amlt, alatr, kp);   // out[0] = ne [cm^-3]

The SM frame (thesis, following the plasma-physics convention) has \(X_\mathrm{SM}\) in the Sun-Earth plane, \(Z_\mathrm{SM}\) aligned with the geomagnetic pole, and \(Y_\mathrm{SM}\) completing the right-hand system.

Electron Density Model

The electron density is supplied by the Global Core Plasma Model (GCPM v2.4), an empirical model of the inner magnetosphere that stitches together region-specific submodels for the ionosphere (based on the International Reference Ionosphere, IRI), the plasmasphere, the plasmapause, the magnetospheric trough, and the polar cap, with hyperbolic-tangent transitions. As a mathematical contract, the density along the path is

\[n_e = n_e\big(r,\ \phi_\mathrm{SM},\ \mathrm{MLT};\ K_p,\ R_{12},\ t\big),\]

where \(r\) is geocentric radius, \(\phi_\mathrm{SM}\) is SM latitude, and \(\mathrm{MLT}\) is magnetic local time. The inputs are:

  • \(K_p\) – planetary geomagnetic index (config.kp). It sets the plasmapause radius: higher \(K_p\) erodes the plasmasphere, contracting the plasmapause from ~4-6 \(R_E\) (quiet) to ~2-3 \(R_E\) (active).

  • \(R_{12}\) / \(F_{10.7}\) – solar activity (config.rz12). The IRI submodel derives \(F_{10.7}\) and the ionospheric index \(IG_{12}\) from \(R_{12}\) via the PHaRLAP relations (thesis Eqs. 5.19-5.20):

    \[F_{10.7} = 63.75 + 0.728\,R_{12} + 0.00089\,R_{12}^2,\]
    \[IG_{12} = -12.349154 + 1.4683266\,R_{12} - 0.00267690893\,R_{12}^2.\]

  • \(t\) – epoch (date/time), which sets IRI solar-activity parameters and the SM frame orientation.

GCPM returns total electron density in cm-3; compute_ne converts to m-3:

\[n_e\,[\mathrm{m}^{-3}] = 10^{6}\, n_e\,[\mathrm{cm}^{-3}].\]

Note

config.kp must be set to a non-negative value; compute_ne throws if kp < 0, and an IRI model must be selected (set_iri_model) or it throws. config.rz12 follows IRI’s convention: > 0 uses the value directly (\(0 < R_{12} \le 200\)), -1 uses IRI’s historical/projected \(R_{12}\) with the storm model, and -2 the same without the storm model. trace_ray applies config.rz12 to the active IRI option (IRI2007Option/IRI2020Option) before tracing. GCPM/IRI model only the inner-magnetospheric cold-plasma environment; contributions outside the cutoff radius (default \(4\,R_E\)) are neglected, where densities are much lower.

Note

GCPM cannot be evaluated for arbitrary future epochs (it requires historical solar/geomagnetic parameters). The thesis baseline uses 2025-01-01 12:00 UTC as a proxy epoch for a 2027 scenario; the characterization spans \(K_p \in \{1,3,5,6.667,9\}\) and \(R_{12} \in \{10,50,100,150,200\}\).

Magnetic Field Model

The magnetic field, needed only for the second- and third-order terms, is evaluated with IGRF-14 in geographic spherical coordinates (raytrace.cc :: compute_B):

cart_to_pol(pos_geo_arr, lat_rad, lon_rad, r);
std::vector<double> B = igrf14(lat_rad, lon_rad, r, decimal_year);   // nT

The projection of the field on the ray direction enters the delay through \(\cos\theta\), the angle between the ray path and the field line:

\[B = \lVert \mathbf{B} \rVert, \qquad \cos\theta = \frac{\mathbf{B}\cdot\hat{s}}{B}.\]

Note

The delay coefficients in the code assume \(B\) in nanotesla. This is the source of the numeric difference from the thesis symbolic coefficients, which are written for \(B\) in tesla (see the note under Refractive Index and Group Delay).

Refractive Index and Group Delay

For a right-hand circularly polarized GNSS signal the phase refractive index \(n\) and group refractive index \(n_\mathrm{gr}\) are (thesis Eqs. 5.1-5.2)

\[n = 1 - \frac{f_p^2}{2f^2} - \frac{f_p^2 f_g \cos\theta}{2f^3} - \frac{f_p^2}{4f^4}\!\left[\frac{f_p^2}{2} + f_g^2(1+\cos^2\theta)\right],\]
\[n_\mathrm{gr} = 1 + \frac{f_p^2}{2f^2} + \frac{f_p^2 f_g \cos\theta}{f^3} + \frac{3f_p^2}{4f^4}\!\left[\frac{f_p^2}{2} + f_g^2(1+\cos^2\theta)\right],\]

with plasma frequency \(f_p = n_e e^2/(4\pi\epsilon_0 m)\) and gyro frequency \(f_g = eB/(2\pi m)\). The delay is dispersive (\(\propto 1/f^2\) at leading order). The ionospheric/plasmaspheric group delay along the path (thesis Eqs. 5.4-5.7) is

\[d_{I_\mathrm{gr}} = \int_\Gamma (n_\mathrm{gr}-1)\,ds = \frac{p}{f^2} + \frac{q}{f^3} + \frac{u}{f^4},\]
\[p = 40.3 \int_\Gamma n_e\, ds = 40.3\,\mathrm{TEC} = 40.3\,(\mathrm{TEC}_\mathrm{LOS} + \Delta\mathrm{TEC}_\mathrm{bend}),\]
\[q = -2.2566\times10^{12} \int_\Gamma n_e\, B \cos\theta\, ds,\]
\[u = 2437 \int_\Gamma n_e^2\, ds + 4.74\times10^{22} \int_\Gamma n_e\, B^2 (1+\cos^2\theta)\, ds .\]

The integral of \(n_e\) along the path is the TEC. Because the path \(\Gamma\) is curved, the TEC splits into a line-of-sight part \(\mathrm{TEC}_\mathrm{LOS}\) and a bending increment \(\Delta\mathrm{TEC}_\mathrm{bend}\) (the difference between the bent-path and straight-line TEC). The first-order term dominates; the second and third orders are below ~1% of the first for terrestrial and LEO users.

The phase index used to bend the ray is raytrace.cc :: refractive_index_neB, and the accumulators for the three group-delay orders are set in raytrace.cc :: ray_derivative:

// refractive_index_neB (phase index that bends the ray)
double p = p_coeff * ne_m3;                          // 40.3 * ne
double q = q_coeff * ne_m3 * B * cos_theta;          // 2.2566e3 * ne*B*cos (B in nT)
double u = u1_coeff*ne_m3*ne_m3
         + u2_coeff*ne_m3*B*B*(1 + cos_theta*cos_theta);
double n = 1 - (p/f2 + q/(2*f3) + u/(3*f4));

// ray_derivative: per-arc group-delay accumulators (s in km -> *1000)
dz(7) = ne_m3 * 1000;                                // d(TEC)/ds
dz(8) = q_coeff * ne_m3 * B * cos_theta * 1000;      // 2nd-order integrand
dz(9) = 1000*(u1_coeff*ne_m3*ne_m3
           + u2_coeff*ne_m3*B*B*(1 + cos_theta*cos_theta));   // 3rd-order integrand

where the coefficients are p_coeff = 40.3, q_coeff = 2.2566e3, u1_coeff = 2437, u2_coeff = 4.74e4.

Note

The code coefficients are scaled for \(B\) in nanotesla, so q_coeff = 2.2566e3 and u2_coeff = 4.74e4 correspond to the thesis SI-tesla values \(2.2566\times10^{12}\) and \(4.74\times10^{22}\) (factors \(10^{-9}\) and \(10^{-18}\)). The code stores q with a positive coefficient and applies the sign of the delay through the phase index (\(-q/(2f^3)\) in refractive_index_neB) and through the sign of \(\cos\theta\); the thesis absorbs the sign into the symbolic definition of \(q\) in Eq. (5.6). Both the phase index (\(q/(2f^3)\), \(u/(3f^4)\)) and the group-delay accumulators (\(q/f^3\), \(u/f^4\)) match Eqs. (5.1)-(5.4) – the group form omits the \(\tfrac12\) and \(\tfrac13\) phase factors.

Ray-Tracing Algorithm

Signal propagation is governed by the vector eikonal equation (thesis Eqs. 5.11-5.13):

\[\frac{dt}{ds} = \frac{n}{c}, \qquad \frac{d\mathbf{r}}{ds} = \hat{s}, \qquad \frac{d\hat{s}}{ds} = \frac{\nabla n - \hat{s}(\hat{s}\cdot\nabla n)}{n},\]

subject to the transmitter/receiver boundary conditions \(\mathbf{r}(t_\mathrm{tx})=\mathbf{r}_\mathrm{tx}\), \(\mathbf{r}(t_\mathrm{rx})=\mathbf{r}_\mathrm{rx}\). Given the reception epoch/position, the transmit epoch is solved by the standard light-time iteration (thesis Eq. 5.16), consistent with GNSS Measurement Model.

The state vector integrated along arc length \(s\) is \(z = [\,\mathbf{r},\ \hat{s},\ t,\ \mathrm{TEC},\ I_2,\ I_3\,]\) (10 components), and the refractive-index gradient is a central finite difference with step config.gradn_dx (default 1 km). The eikonal right-hand side is raytrace.cc :: ray_derivative:

dz.segment<3>(0) = s_hat;                                   // dr/ds = s_hat
dz.segment<3>(3) = (grad_n - s_hat*grad_n.dot(s_hat)) / n;  // ds_hat/ds  (Eq. 5.13)
dz(6) = 1.0 / C;                                            // dt/ds ~ 1/c

Note

The geometric propagation time is accumulated as dz(6) = 1/C (vacuum) rather than n/c; the ionospheric delay is then added separately as (tec_delay + second_delay + third_delay)/c. This is a deliberate machine-precision choice (the code comment notes that the per-step TEC time increment is \(\sim n_e\cdot 10^{-23}\)s, too small to co-accumulate with the epoch), and is numerically equivalent to integrating \(dt/ds = n/c\).

Shooting method with cutoff split. The path is divided into a plasmasphere section and a vacuum section at the cutoff radius config.cutoff_r (default \(4\,R_E \approx 25484\) km): the ray is assumed straight once it exits the cutoff. Integration marches from the transmitter with the fourth-order Runge-Kutta or Euler integrator (config.integ_method) until \(r \ge\) cutoff, then shoots straight to the receiver (thesis Eqs. 5.17-5.18):

\[s_f = s_\mathrm{exit} + \Delta s, \qquad \Delta s = (\mathbf{x}_\mathrm{rx} - \mathbf{x}_\mathrm{exit})\cdot \hat{s}_\mathrm{exit},\]
// propagate_ray: exit at cutoff, project remaining distance onto s_hat_exit
Vec3d delta_x2 = x2 - pos_now;
double delta_s = delta_x2.dot(dir_now);
sf = s + delta_s;
Vec3d final_pos = pos_now + dir_now * (sf - s);

An adaptive step size is used inside the plasmasphere (raytrace.cc :: adjust_stepsize), keyed on altitude \(\mathrm{alt}=r-R_E\):

if (alt < 1000)      step_size = 10.0;   // km
else if (alt < 4000) step_size = 20.0;
else                 step_size = 100.0;

Initial-direction correction (the iterative bent-path solve). The unknown initial ray direction \(\hat{s}_0\) is adjusted so the shot lands on the receiver. trace_ray runs up to 20 coarse iterations that (a) refresh the propagation time from the accumulated delay (correct_proptime) and (b) re-optimize \(\hat{s}_0\) over the 2-D azimuth/elevation of the initial direction with a Nelder-Mead simplex minimizing the terminal miss \(\lVert\bar{\mathbf{x}}_\mathrm{rx} - \mathbf{x}_\mathrm{rx}\rVert\) (correct_dir_neldermead). To keep cost low, the coarse pass reuses the stored per-step refractive index and its gradient (store_mat) instead of re-querying GCPM at every point. An optional finer pass (config.fine_correction) refines with Nelder-Mead or a finite-difference Newton solve (correct_dir_newton) to the tolerance config.corr_tol (default 1 m). Without correction the terminal error is 10-100 km due to bending; with correction it is driven to sub-meter-100 m.

// trace_ray: coarse loop -- proptime update then direction re-solve
correct_proptime(prop_time, dir, sf, x1, x2, epoch_utc_rx, config, zf);
correct_dir_neldermead(prop_time, dir, sf, x1, x2, epoch_utc_rx, config, zf,
                       store_mat, /*coarse_mode=*/true, debug_corr);

Note

config.correction_method must be "neldermead" or "newton"; trace_ray throws otherwise. The struct default is "grid" (which would throw), so a caller enabling correction must set this field. The thrown message text ("Supported methods are 'naive' and 'newton'") is stale relative to the actual check.

Delay Decomposition (PathProfile)

After the corrected path is integrated, compute_path_profile assembles the delay terms in meters from the accumulators (z(7..9)) and the geometry:

double tec_delay_m    = p_coeff * z(7) / freq2;          // 40.3 * TEC_bent / f^2
double second_delay_m = z(8) / freq3;
double third_delay_m  = z(9) / freq4;
pp.dist_bend_m        = (sf - dist) * 1000;              // excess path length
pp.tec_delay_bend_m   = tec_delay_m - tec_delay_straight_m;

These map onto the thesis decomposition (Eqs. 5.8-5.10). The excess path length from bending is

\[d_\mathrm{len} = \int_\Gamma ds - \rho, \qquad \rho = \lVert \mathbf{r}_\mathrm{tx} - \mathbf{r}_\mathrm{rx} \rVert,\]

and the total group delay relative to the LOS distance is

\[d_\mathrm{total} = d_{I_\mathrm{gr}} + d_\mathrm{len} = \underbrace{40.3\,\mathrm{TEC}_\mathrm{LOS}}_{d_{I1,\mathrm{LOS}}} + \underbrace{q/f^3}_{d_{I2}} + \underbrace{u/f^4}_{d_{I3}} + \underbrace{40.3\,\Delta\mathrm{TEC}_\mathrm{bend}}_{d_{I1,B}} + d_\mathrm{len}.\]

The PathProfile fields carry each term: tec_delay_m (\(40.3\,\mathrm{TEC}\) along the bent path, i.e. \(d_{I1,\mathrm{LOS}}+d_{I1,B}\)), tec_delay_bend_m (\(d_{I1,B}\)), second_delay_m (\(d_{I2}\)), third_delay_m (\(d_{I3}\)), dist_bend_m (\(d_\mathrm{len}\)), and total_delay_m (\(d_\mathrm{total}\)).

Configuration Knobs (RayTraceConfig)

The complete set of knobs (raytrace.h), with defaults:

  • freq_Hz – carrier frequency [Hz]; L1/L2/L5 constants are provided (freq_L1 = 1575.42e6 etc.). All delays scale as \(1/f^2\) to leading order, so L5 (1176.45 MHz) sees ~1.8x the L1 delay.

  • step_size = 10.0 – nominal integrator step [km].

  • use_adaptive_step = true – enable the 10/20/100 km altitude schedule.

  • integ_method = "Euler" – "Euler" or "RK4".

  • correction = true – run the bent-path direction correction.

  • fine_correction = false – run the additional fine pass.

  • correction_method = "grid" – must be set to "neldermead" or "newton" when correction is on (see note above).

  • corr_tol = 1.0 – terminal miss tolerance [m].

  • cutoff_r = 4*RE – plasmasphere/vacuum boundary [km].

  • gradn_dx = 1.0 – central-difference step for \(\nabla n\) [km].

  • kp = -1 – geomagnetic \(K_p\) (must be set \(\ge 0\)).

  • rz12 = -1.0 – solar \(R_{12}\) (see IRI convention above).

  • compute_higher_order = true – include second/third orders (needs IGRF).

  • straight_ray = false – assume a straight path (no bending; skips coarse correction, zeros \(\nabla n\)), for the straight-line TEC baseline.

  • use_fortran_gcpm = true – use the Fortran GCPM v2.4 backend.

Coupling into the GNSS Measurement Model

When options.apply_ionosphere_plasma_delay is true and ray-trace options have been registered via SetIonospherePlasmaRayTraceOptions, the GNSS model calls the tracer per link (gnss_measurement.cc :: ComputeIonospherePlasmaRayTraceDelay):

pecsim::RayTraceConfig config = ray_opts.config;
if (ray_opts.use_channel_frequency) config.freq_Hz = GnssFrequencyHz(freq);
if (ray_opts.iri_model != pecsim::IRIModel::NONE) set_iri_model(ray_opts.iri_model);
Real t_raytrace = ConvertGnssEpoch(receive_time, options_.receive_time_scale,
                                   ray_opts.raytrace_epoch_scale);
// r_rx, r_tx converted options_.frame -> ray_opts.raytrace_frame, then to km
pecsim::PathProfile profile = pecsim::trace_ray(
    t_raytrace.val(), ToVec3d(r_tx_ray)/1000.0, ToVec3d(r_rx_ray)/1000.0, config);
Real delay_m = profile.tec_delay_m;
if (ray_opts.delay_mode == GnssIonospherePlasmaRayTraceDelayMode::TEC_PLUS_HIGHER_ORDER)
  delay_m += profile.second_delay_m + profile.third_delay_m;

The returned \(\Delta\rho_\mathrm{plasma}\) is the positive code delay of GNSS Measurement Model, entering the pseudorange as \(+\Delta\rho_\mathrm{plasma}\) and the carrier phase as \(-\Delta\rho_\mathrm{plasma}\). The GnssIonospherePlasmaRayTraceOptions fields are: config (the tracer knobs; GNSSMeasurements never fills \(K_p\), integrator, or step sizes), raytrace_frame (default ECEF), raytrace_epoch_scale (default UTC), use_channel_frequency (default true – overwrite config.freq_Hz with the channel frequency), iri_model, and delay_mode.

Note

delay_mode defaults to TEC_ONLY, so the returned delay is profile.tec_delay_m (the first-order bent-path TEC delay, \(d_{I1,\mathrm{LOS}}+d_{I1,B}\)) alone; TEC_PLUS_HIGHER_ORDER adds \(d_{I2}+d_{I3}\). The geometric excess-path term \(d_\mathrm{len}\) (profile.dist_bend_m) is computed but not added to the pseudorange delay; its contribution to timing is already the \(\Delta\mathrm{TEC}_\mathrm{bend}\) part of tec_delay_m, and the sub-meter geometric term is left out of the range equation.

Batch vs Online Precompute

Ray tracing is expensive (a GCPM call per step per iteration), so two paths exist. The online path calls trace_ray link-by-link inside Compute. The batch/precompute path decouples the tracer from filtering: in the staged Lunar GNSS ODTS scenario, Precompute writes all light-time-corrected links with zero plasma delay, a Python driver (precompute_delays.py) fills each epoch/channel delay with GCPM/IRI ray tracing, and Run merges the delay file back before forming observables (see the Online and Precompute Equivalence section of GNSS Measurement Model). The batch tracer is exposed through cpp/lupnt/environment/plasma/tec/raytrace_wrapper.{h,cc} (pecsim::run_raytrace_batch, load_raytrace_data / RaytraceData), which stores per-ray total_delays, tec_delays, second_delays, third_delays, dist_bend_m, tec_delay_bend_m, tecus, and min_alts. The \(K_p\) time series consumed by these runs is fetched by python/pylupnt/plasma/kp_loader.py from GFZ Potsdam and cached under the plasma data directory.

Model Boundaries and Notes

  • Km-unit convention. All tracer inputs/outputs are km / km-s based; SI callers scale by 1000 at the boundary.

  • Magnetic-field units. Delay coefficients assume \(B\) in nT; the thesis symbolic coefficients are the SI-tesla equivalents.

  • Second/third order are small. For terrestrial and LEO users they are below ~1% of the first-order delay; the thesis lunar characterization finds them several orders of magnitude below \(d_{I1}\) (Table 5.2), and the GNSS default (TEC_ONLY) omits them.

  • Cutoff at \(4\,R_E\). Plasma outside the inner magnetosphere is neglected and the path is straight beyond the cutoff.

  • Time integration. dt/ds uses the vacuum \(1/c\) with the delay added separately, equivalent to \(n/c\) but numerically robust.

  • Excess path length \(d_\mathrm{len}\) is characterized but not added to the measurement delay.

  • Validity ranges. GCPM/IRI require historical solar/geomagnetic inputs; the characterization covers \(K_p \in [1,9]\) and \(R_{12}\in[10,200]\), and delays are strongly altitude-dependent (mean first-order delay ~50 m for tangential altitudes below 500 km, falling below ~0.2 m above 10,000 km).