GNSS Measurement Simulation

This page outlines the receiver-side GNSS simulation flow in C++. GnssConstellation stores GNSS satellite ephemerides and transmitter metadata. GNSSMeasurements owns the receiver-dependent work: light-time iteration, visibility, relativistic range corrections, plasma delays, channel selection, and observable generation.

Time and State Conventions

Epochs passed to GNSSMeasurements::Compute and GNSSMeasurements::Precompute are receiver signal-reception epochs. The scale is explicit. A cislunar simulation will often use TDB for the receiver coordinate timeline while keeping GNSS ephemeris tables in TAI:

GnssMeasurementOptions options;
options.receive_time_scale = Time::TDB;     // receive epochs you pass in
options.ephemeris_time_scale = Time::TAI;   // GnssConstellation tables
options.frame = Frame::MOON_CI;             // common inertial state frame

The receiver state must contain position, velocity, clock bias, and clock drift at the configured indices:

options.indices.position = 0;
options.indices.velocity = 3;
options.indices.clock_bias = 6;
options.indices.clock_drift = 7;
options.clock_bias_unit = ClockBiasUnit::SECONDS;

ClockBiasUnit::METERS and ClockBiasUnit::KILOMETERS are also supported when a filter state should use range-like clock units for numerical scaling.

Constellation Setup

The constellation is a satellite catalog. It does not decide which satellites are visible to a receiver.

auto gps = MakePtr<GnssConstellation>(GnssConst::GPS);
gps->SetSatelliteStates(prns, t_tai, rv_eci_by_prn);
gps->SetupTransmitters();

SetSatelliteStates stores the tabulated ephemerides in TAI seconds and fits a piecewise-Chebyshev state model. Runtime satellite interpolation uses the Chebyshev model.

Online Flow

Online computation is for filters and closed-loop simulations where the receiver state is known one epoch at a time.

GNSSMeasurements gnss(gps);
gnss.SetFrequency(GnssFreq::L1);
gnss.SetOptions(options);

GnssOccludingBody earth;
earth.radius_m = R_EARTH;
earth.position_m = Vec3::Zero();
gnss.SetOccludingBodies({earth});

gnss.SetSunPositionProvider([](Real t_rx) {
  return Vec3(0.0, AU, 0.0);
});
gnss.SetBoresightTargetProvider([](Real t_rx) {
  return Vec3::Zero();
});

MatXd H;
GNSSMeasurementsEpoch epoch = gnss.Compute(t_rx, receiver_state, &H);

For each PRN, Compute:

  1. Converts the receiver receive epoch to the ephemeris scale.

  2. Iterates transmit time with geometric light time.

  3. Builds a GnssChannel containing transmit state, receive/transmit epochs, ephemeris, C/N0, and correction terms.

  4. Applies receiver-dependent visibility and C/N0 thresholding.

  5. Computes pseudorange, Doppler, carrier phase, covariance, and optional Jacobian.

Relativistic and Plasma Corrections

The measurement model includes transmitter clock relativity, Shapiro delay, and optional ionosphere/plasma delays.

options.apply_transmitter_relativity = true;

options.apply_shapiro_delay = true;
options.shapiro_mu = GM_EARTH;
options.shapiro_body_position = Vec3::Zero();

options.apply_ionosphere_plasma_delay = true;
options.default_ionosphere_plasma_delay_m = 0.0;

GnssIonospherePlasmaRayTraceOptions plasma;
plasma.config.kp = kp_index;
plasma.config.step_size = step_size_km;
plasma.config.cutoff_r = cutoff_radius_km;
plasma.config.gradn_dx = gradient_step_km;
plasma.config.integ_method = integration_method;
plasma.config.correction_method = correction_method;
plasma.config.correction = apply_ray_correction;
plasma.config.fine_correction = apply_fine_correction;
plasma.config.corr_tol = correction_tolerance_m;
plasma.config.use_fortran_gcpm = use_fortran_gcpm;
plasma.config.compute_higher_order = compute_higher_order;
plasma.config.use_adaptive_step = use_adaptive_step;
plasma.config.straight_ray = use_straight_ray;
plasma.raytrace_frame = Frame::ECEF;
plasma.raytrace_epoch_scale = Time::UTC;
plasma.iri_model = pecsim::IRIModel::IRI_2020;
gnss.SetIonospherePlasmaRayTraceOptions(plasma);

Shapiro delay is non-dispersive and is added to both pseudorange and carrier phase range. Ionosphere/plasma delay is represented as a positive code delay in meters: pseudorange uses +delay and carrier phase uses -delay. The GNSS raytrace hook calls pecsim::trace_ray and uses the returned TEC delay. The caller owns the plasma simulation settings; the measurement class does not pick Kp, step size, cutoff radius, integrator, or correction policy.

For cislunar GNSS filtering runs, the recommended high-throughput workflow is to precompute all receiver-specific links first, then run the expensive GCPM/IRI delay calculation as a separate Python multiprocessing step. See GNSS Filtering Pipeline for the full staged pipeline.

Precompute Flow

Precompute is for generating a measurement sequence before filtering. It uses the same channel-building logic as online computation, but it can call a batch ionosphere/plasma model after all epochs and links are known.

std::vector<Real> receive_times = {t0, t1, t2};
std::vector<State> receiver_states = {x0, x1, x2};

gnss.SetBatchCustomIonospherePlasmaDelayModel(
    [](const std::vector<Real>& times,
       const std::vector<State>& states,
       const std::vector<std::vector<GnssChannel>>& channels) {
      std::vector<std::vector<Real>> delays(times.size());
      for (size_t i = 0; i < times.size(); i++) {
        delays[i].resize(channels[i].size());
        for (size_t j = 0; j < channels[i].size(); j++) {
          delays[i][j] = 0.0;  // batch model result in meters
        }
      }
      return delays;
    });

std::vector<GNSSMeasurementsEpoch> epochs =
    gnss.Precompute(receive_times, receiver_states, true);

The batch provider receives:

  • receiver receive epochs in options.receive_time_scale

  • receiver states in options.frame

  • receiver-specific channels with light-time transmit epochs already solved

This is the right hook for ionosphere/plasma models that need the full simulation interval, endpoint geometry, cached environmental grids, or shared line-of-sight integrations. When the batch provider is installed, Precompute skips the online plasma hook during channel construction and fills each channel delay from the batch result.

Using Precomputed Data in a Filter

Each GNSSMeasurementsEpoch contains the channels and the stacked observables:

for (const GNSSMeasurementsEpoch& epoch : epochs) {
  const VecXd& y = epoch.values;
  const MatXd& R = epoch.covariance;
  const MatXd& H = epoch.jacobian;

  for (const GnssChannel& channel : epoch.channels) {
    int prn = channel.prn;
    Real t_tx = channel.transmit_time;
    Real shapiro_m = channel.shapiro_delay_m;
    Real plasma_m = channel.ionosphere_plasma_delay_m;
  }

End-to-End Project Script

The release-facing GNSS filtering example is under projects/GNSS_Filtering and is implemented by the core LunarGnssODTSSimulation under cpp/lupnt/simulations/LunarGnssODTS. It builds the LuPNT simulation executable, precomputes links, optionally precomputes GCPM delays, runs the Monte Carlo filter, and generates plots:

pixi run run-gnss-pipeline

Skip the delay batch only when the configured delay table already exists, or when using a config with plasma.simulate_truth: false:

pixi run run-gnss-pipeline --skip-delays
}