INSTINCT Code Coverage Report

Directory:	src/
File:	Navigation/GNSS/Satellite/Ephemeris/GalileoEphemeris.cpp
Date:	2025-07-19 10:51:51

	Exec	Total	Coverage
Lines:	152	152	100.0%
Functions:	7	7	100.0%
Branches:	54	86	62.8%

  
      Line
      Branch
      Exec
      Source
    
      // This file is part of INSTINCT, the INS Toolkit for Integrated
    
      // Navigation Concepts and Training by the Institute of Navigation of
    
      // the University of Stuttgart, Germany.
    
      //
    
      // This Source Code Form is subject to the terms of the Mozilla Public
    
      // License, v. 2.0. If a copy of the MPL was not distributed with this
    
      // file, You can obtain one at https://mozilla.org/MPL/2.0/.
    
      #include "GalileoEphemeris.hpp"
    
      #include "Navigation/Constants.hpp"
    
      #include "Navigation/GNSS/Functions.hpp"
    
      #include "util/Logger.hpp"
    
      namespace NAV
    
      {
    
      28
      GalileoEphemeris::GalileoEphemeris(const InsTime& toc)
    
      28
          : SatNavData(SatNavData::GalileoEphemeris, toc) {}
    
      1633
      GalileoEphemeris::GalileoEphemeris(const InsTime& toc, const InsTime& toe,
    
                                         const size_t& IODnav,
    
                                         const std::array<double, 3>& a,
    
                                         const double& sqrt_A, const double& e, const double& i_0, const double& Omega_0, const double& omega, const double& M_0,
    
                                         const double& delta_n, const double& Omega_dot, const double& i_dot, const double& Cus, const double& Cuc,
    
                                         const double& Cis, const double& Cic, const double& Crs, const double& Crc,
    
                                         const std::bitset<10>& dataSource, const double& signalAccuracy, const SvHealth& svHealth,
    
      1633
                                         const double& BGD_E1_E5a, const double& BGD_E1_E5b)
    
          : SatNavData(SatNavData::GalileoEphemeris, toc),
    
      1633
            toc(toc),
    
      1633
            toe(toe),
    
      1633
            IODnav(IODnav),
    
      1633
            a(a),
    
      1633
            sqrt_A(sqrt_A),
    
      1633
            e(e),
    
      1633
            i_0(i_0),
    
      1633
            Omega_0(Omega_0),
    
      1633
            omega(omega),
    
      1633
            M_0(M_0),
    
      1633
            delta_n(delta_n),
    
      1633
            Omega_dot(Omega_dot),
    
      1633
            i_dot(i_dot),
    
      1633
            Cus(Cus),
    
      1633
            Cuc(Cuc),
    
      1633
            Cis(Cis),
    
      1633
            Cic(Cic),
    
      1633
            Crs(Crs),
    
      1633
            Crc(Crc),
    
      1633
            dataSource(dataSource),
    
      1633
            signalAccuracy(signalAccuracy),
    
      1633
            svHealth(svHealth),
    
      1633
            BGD_E1_E5a(BGD_E1_E5a),
    
      1633
            BGD_E1_E5b(BGD_E1_E5b) {}
    
      #ifdef TESTING
    
      8556
      GalileoEphemeris::GalileoEphemeris(int32_t year, int32_t month, int32_t day, int32_t hour, int32_t minute, double second, double svClockBias, double svClockDrift, double svClockDriftRate,
    
                                         double IODnav, double Crs, double delta_n, double M_0,
    
                                         double Cuc, double e, double Cus, double sqrt_A,
    
                                         double Toe, double Cic, double Omega_0, double Cis,
    
                                         double i_0, double Crc, double omega, double Omega_dot,
    
                                         double i_dot, double dataSource, double GALWeek, double /* spare1 */,
    
                                         double signalAccuracy, double svHealth, double BGD_E1_E5a, double BGD_E1_E5b,
    
      8556
                                         double /* TransmissionTimeOfMessage */, double /* spare2 */, double /* spare3 */, double /* spare4 */)
    
        1/2✓ Branch 1 taken 8556 times.
✗ Branch 2 not taken.

      8556
          : SatNavData(SatNavData::GalileoEphemeris, InsTime(year, month, day, hour, minute, second, GST)),
    
      8556
            toc(refTime),
    
        1/2✓ Branch 1 taken 8556 times.
✗ Branch 2 not taken.

      8556
            toe(InsTime(0, static_cast<int32_t>(GALWeek), Toe, GST)),
    
      8556
            IODnav(static_cast<size_t>(IODnav)),
    
      8556
            a({ svClockBias, svClockDrift, svClockDriftRate }),
    
      8556
            sqrt_A(sqrt_A),
    
      8556
            e(e),
    
      8556
            i_0(i_0),
    
      8556
            Omega_0(Omega_0),
    
      8556
            omega(omega),
    
      8556
            M_0(M_0),
    
      8556
            delta_n(delta_n),
    
      8556
            Omega_dot(Omega_dot),
    
      8556
            i_dot(i_dot),
    
      8556
            Cus(Cus),
    
      8556
            Cuc(Cuc),
    
      8556
            Cis(Cis),
    
      8556
            Cic(Cic),
    
      8556
            Crs(Crs),
    
      8556
            Crc(Crc),
    
      8556
            dataSource(static_cast<uint16_t>(dataSource)),
    
      8556
            signalAccuracy(signalAccuracy),
    
      8556
            svHealth({ .E5a_DataValidityStatus = static_cast<GalileoEphemeris::SvHealth::DataValidityStatus>((static_cast<uint16_t>(svHealth) & 0b000001000) >> 3),
    
      8556
                       .E5b_DataValidityStatus = static_cast<GalileoEphemeris::SvHealth::DataValidityStatus>((static_cast<uint16_t>(svHealth) & 0b001000000) >> 6),
    
      8556
                       .E1B_DataValidityStatus = static_cast<GalileoEphemeris::SvHealth::DataValidityStatus>((static_cast<uint16_t>(svHealth) & 0b000000001) >> 0),
    
      8556
                       .E5a_SignalHealthStatus = static_cast<GalileoEphemeris::SvHealth::SignalHealthStatus>((static_cast<uint16_t>(svHealth) & 0b000110000) >> 4),
    
      8556
                       .E5b_SignalHealthStatus = static_cast<GalileoEphemeris::SvHealth::SignalHealthStatus>((static_cast<uint16_t>(svHealth) & 0b110000000) >> 7),
    
      8556
                       .E1B_SignalHealthStatus = static_cast<GalileoEphemeris::SvHealth::SignalHealthStatus>((static_cast<uint16_t>(svHealth) & 0b000000110) >> 1) }),
    
      8556
            BGD_E1_E5a(BGD_E1_E5a),
    
        1/2✓ Branch 2 taken 8556 times.
✗ Branch 3 not taken.

      17112
            BGD_E1_E5b(BGD_E1_E5b)
    
      8556
      {}
    
      #endif
    
      122361
      Clock::Corrections GalileoEphemeris::calcClockCorrections(const InsTime& recvTime, double dist, const Frequency& freq) const
    
      {
    
          LOG_DATA("Calc Sat Clock corrections at receiver time {}", recvTime.toGPSweekTow());
    
          // Earth gravitational constant [m³/s²] (WGS 84 value of the earth's gravitational constant for GPS user)
    
      122361
          const auto mu = InsConst::GAL::MU;
    
          // Relativistic constant F for clock corrections [s/√m] (-2*√µ/c²)
    
      122361
          const auto F = InsConst::GAL::F;
    
          LOG_DATA("    toe {} (Time of ephemeris)", toe.toGPSweekTow());
    
      122361
          const auto A = sqrt_A * sqrt_A; // Semi-major axis [m]
    
          LOG_DATA("    A {} [m] (Semi-major axis)", A);
    
      122361
          auto n_0 = std::sqrt(mu / std::pow(A, 3)); // Computed mean motion [rad/s]
    
          LOG_DATA("    n_0 {} [rad/s] (Computed mean motion)", n_0);
    
      122362
          auto n = n_0 + delta_n; // Corrected mean motion [rad/s]
    
          LOG_DATA("    n {} [rad/s] (Corrected mean motion)", n);
    
          // Time at transmission
    
        2/4✓ Branch 2 taken 122362 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 122363 times.
✗ Branch 6 not taken.

      122362
          InsTime transTime0 = recvTime - std::chrono::duration<double>(dist / InsConst::C);
    
      122363
          InsTime transTime = transTime0;
    
          LOG_DATA("    Iterating Time at transmission");
    
      122363
          double dt_sv = 0.0;
    
      122363
          double clkDrift = 0.0;
    
        2/2✓ Branch 0 taken 244726 times.
✓ Branch 1 taken 122360 times.

      367086
          for (size_t i = 0; i < 2; i++)
    
          {
    
              LOG_DATA("      transTime {} (Time at transmission)", transTime.toGPSweekTow());
    
              // [s]
    
        1/2✓ Branch 1 taken 244722 times.
✗ Branch 2 not taken.

      244726
              auto t_minus_toc = static_cast<double>((transTime - toc).count());
    
              LOG_DATA("      transTime - toc {} [s]", t_minus_toc);
    
              // Time difference from ephemeris reference epoch [s]
    
        1/2✓ Branch 1 taken 244723 times.
✗ Branch 2 not taken.

      244722
              double t_k = static_cast<double>((transTime - toe).count());
    
              LOG_DATA("      transTime - toe {} [s] (t_k = Time difference from ephemeris reference epoch)", t_k);
    
              // Mean anomaly [rad]
    
      244723
              auto M_k = M_0 + n * t_k;
    
              LOG_DATA("      M_k {} [s] (Mean anomaly)", M_k);
    
              // Eccentric anomaly [rad]
    
      244723
              double E_k = M_k;
    
      244723
              double E_k_old = 0.0;
    
        5/6✓ Branch 1 taken 972905 times.
✓ Branch 2 taken 244712 times.
✓ Branch 3 taken 972905 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 972902 times.
✓ Branch 6 taken 244715 times.

      1217625
              for (size_t i = 0; std::abs(E_k - E_k_old) > 1e-13 && i < 10; i++)
    
              {
    
      972902
                  E_k_old = E_k; // Kepler’s equation ( Mk = E_k − e sin E_k ) may be solved for Eccentric anomaly (E_k) by iteration:
    
      972902
                  E_k = M_k + e * sin(E_k);
    
              }
    
              // Relativistic correction term [s]
    
      244715
              double dt_r = F * e * sqrt_A * std::sin(E_k);
    
              LOG_DATA("      dt_r {} [s] (Relativistic correction term)", dt_r);
    
              // SV PRN code phase time offset [s]
    
      244715
              dt_sv = a[0] + a[1] * t_minus_toc + a[2] * std::pow(t_minus_toc, 2) + dt_r;
    
              // See GAL-ICD-2.0 GAL ICD, ch. 5.1.5, p.47
    
        3/4✓ Branch 2 taken 244720 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 61773 times.
✓ Branch 6 taken 182950 times.

      244722
              dt_sv -= ratioFreqSquared(E01, freq, -128, -128) * (freq == E05 ? BGD_E1_E5a : BGD_E1_E5b); // TODO: Check again
    
              LOG_DATA("      dt_sv {} [s] (SV PRN code phase time offset)", dt_sv);
    
              // Groves ch. 9.3.1, eq. 9.78, p. 391
    
      244723
              clkDrift = a[1] + a[2] / 2.0 * t_minus_toc;
    
              // Correct transmit time for the satellite clock bias
    
        2/4✓ Branch 2 taken 244721 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 244723 times.
✗ Branch 6 not taken.

      244724
              transTime = transTime0 - std::chrono::duration<double>(dt_sv);
    
          }
    
          LOG_DATA("      transTime {} (Time at transmission)", transTime.toGPSweekTow());
    
      122360
          return { .transmitTime = transTime, .bias = dt_sv, .drift = clkDrift };
    
      }
    
      124029
      Orbit::PosVelAccel GalileoEphemeris::calcSatelliteData(const InsTime& transTime, Orbit::Calc calc) const
    
      {
    
        2/4✓ Branch 1 taken 124030 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124031 times.
✗ Branch 5 not taken.

      124029
          Eigen::Vector3d e_pos = Eigen::Vector3d::Zero();
    
        2/4✓ Branch 1 taken 124031 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124031 times.
✗ Branch 5 not taken.

      124031
          Eigen::Vector3d e_vel = Eigen::Vector3d::Zero();
    
        2/4✓ Branch 1 taken 124031 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124031 times.
✗ Branch 5 not taken.

      124031
          Eigen::Vector3d e_accel = Eigen::Vector3d::Zero();
    
          LOG_DATA("Calc Sat Position at transmit time {}", transTime.toGPSweekTow());
    
          // Earth gravitational constant [m³/s²] (WGS 84 value of the earth's gravitational constant for GPS user)
    
      124031
          const auto mu = InsConst::GAL::MU;
    
          // Earth angular velocity [rad/s] (WGS 84 value of the earth's rotation rate)
    
      124031
          const auto Omega_e_dot = InsConst::GAL::omega_ie;
    
          LOG_DATA("    toe {} (Time of ephemeris)", toe.toGPSweekTow());
    
      124031
          const auto A = sqrt_A * sqrt_A; // Semi-major axis [m]
    
          LOG_DATA("    A {} [m] (Semi-major axis)", A);
    
      124031
          auto n_0 = std::sqrt(mu / std::pow(A, 3)); // Computed mean motion [rad/s]
    
          LOG_DATA("    n_0 {} [rad/s] (Computed mean motion)", n_0);
    
      124031
          auto n = n_0 + delta_n; // Corrected mean motion [rad/s]
    
          LOG_DATA("    n {} [rad/s] (Corrected mean motion)", n);
    
          // Eccentric anomaly [rad]
    
      124031
          double E_k = 0.0;
    
          // Time difference from ephemeris reference epoch [s]
    
        1/2✓ Branch 1 taken 124031 times.
✗ Branch 2 not taken.

      124031
          double t_k = static_cast<double>((transTime - toe).count());
    
          LOG_DATA("    t_k {} [s] (Time difference from ephemeris reference epoch)", t_k);
    
          // Mean anomaly [rad]
    
      124031
          auto M_k = M_0 + n * t_k;
    
          LOG_DATA("    M_k {} [s] (Mean anomaly)", M_k);
    
      124031
          E_k = M_k; // Initial Value [rad]
    
      124031
          double E_k_old = 0.0;
    
          LOG_DATA("    Iterating E_k");
    
          LOG_DATA("      E_k {} [rad] (Eccentric anomaly)", E_k);
    
        5/6✓ Branch 1 taken 306105 times.
✓ Branch 2 taken 124028 times.
✓ Branch 3 taken 306105 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 306104 times.
✓ Branch 6 taken 124029 times.

      430135
          for (size_t i = 0; std::abs(E_k - E_k_old) > 1e-13 && i < 10; i++)
    
          {
    
      306104
              E_k_old = E_k;                                                         // Kepler’s equation ( Mk = E_k − e sin E_k ) may be solved for Eccentric anomaly (E_k) by iteration:
    
      306104
              E_k = E_k + (M_k - E_k + e * std::sin(E_k)) / (1 - e * std::cos(E_k)); // - Refined Value, minimum of three iterations, (j=1,2,3)
    
              LOG_DATA("      E_k {} [rad] (Eccentric anomaly)", E_k);               // - Final Value (radians)
    
          }
    
          // auto v_k = 2.0 * std::atan(std::sqrt((1.0 + e) / (1.0 - e)) * std::tan(E_k / 2.0)); // True Anomaly (unambiguous quadrant) [rad] (GPS ICD algorithm)
    
          // auto v_k = std::atan2(std::sqrt(1 - e * e) * std::sin(E_k) / (1 - e * std::cos(E_k)), (std::cos(E_k) - e) / (1 - e * std::cos(E_k))); // True Anomaly [rad] (GALILEO ICD algorithm)
    
      124029
          auto v_k = std::atan2(std::sqrt(1 - e * e) * std::sin(E_k), (std::cos(E_k) - e)); // True Anomaly [rad] // simplified, since the denominators cancel out
    
          LOG_DATA("    v_k {} [rad] (True Anomaly (unambiguous quadrant))", v_k);
    
      124029
          auto Phi_k = v_k + omega; // Argument of Latitude [rad]
    
          LOG_DATA("    Phi_k {} [rad] (Argument of Latitude)", Phi_k);
    
          // Second Harmonic Perturbations
    
      124029
          auto delta_u_k = Cus * std::sin(2 * Phi_k) + Cuc * std::cos(2 * Phi_k); // Argument of Latitude Correction [rad]
    
          LOG_DATA("    delta_u_k {} [rad] (Argument of Latitude Correction)", delta_u_k);
    
      124029
          auto delta_r_k = Crs * std::sin(2 * Phi_k) + Crc * std::cos(2 * Phi_k); // Radius Correction [m]
    
          LOG_DATA("    delta_r_k {} [m] (Radius Correction)", delta_r_k);
    
      124029
          auto delta_i_k = Cis * std::sin(2 * Phi_k) + Cic * std::cos(2 * Phi_k); // Inclination Correction [rad]
    
          LOG_DATA("    delta_i_k {} [rad] (Inclination Correction)", delta_i_k);
    
      124029
          auto u_k = Phi_k + delta_u_k; // Corrected Argument of Latitude [rad]
    
          LOG_DATA("    u_k {} [rad] (Corrected Argument of Latitude)", u_k);
    
      124029
          auto r_k = A * (1 - e * std::cos(E_k)) + delta_r_k; // Corrected Radius [m]
    
          LOG_DATA("    r_k {} [m] (Corrected Radius)", r_k);
    
      124029
          auto i_k = i_0 + delta_i_k + i_dot * t_k; // Corrected Inclination [rad]
    
          LOG_DATA("    i_k {} [rad] (Corrected Inclination)", i_k);
    
      124029
          auto x_k_op = r_k * std::cos(u_k); // Position in orbital plane [m]
    
          LOG_DATA("    x_k_op {} [m] (Position in orbital plane)", x_k_op);
    
      124029
          auto y_k_op = r_k * std::sin(u_k); // Position in orbital plane [m]
    
          LOG_DATA("    y_k_op {} [m] (Position in orbital plane)", y_k_op);
    
          // Corrected longitude of ascending node [rad]
    
        1/2✓ Branch 2 taken 124031 times.
✗ Branch 3 not taken.

      124029
          auto Omega_k = Omega_0 + (Omega_dot - Omega_e_dot) * t_k - Omega_e_dot * static_cast<double>(toe.toGPSweekTow().tow);
    
          LOG_DATA("    Omega_k {} [rad] (Corrected longitude of ascending node)", Omega_k);
    
          // Earth-fixed x coordinates [m]
    
      124031
          auto x_k = x_k_op * std::cos(Omega_k) - y_k_op * std::cos(i_k) * std::sin(Omega_k);
    
          LOG_DATA("    x_k {} [m] (Earth-fixed x coordinates)", x_k);
    
          // Earth-fixed y coordinates [m]
    
      124031
          auto y_k = x_k_op * std::sin(Omega_k) + y_k_op * std::cos(i_k) * std::cos(Omega_k);
    
          LOG_DATA("    y_k {} [m] (Earth-fixed y coordinates)", y_k);
    
          // Earth-fixed z coordinates [m]
    
      124031
          auto z_k = y_k_op * std::sin(i_k);
    
          LOG_DATA("    z_k {} [m] (Earth-fixed z coordinates)", z_k);
    
        1/2✓ Branch 1 taken 124031 times.
✗ Branch 2 not taken.

      124031
          e_pos = Eigen::Vector3d{ x_k, y_k, z_k };
    
        5/6✓ Branch 1 taken 1987 times.
✓ Branch 2 taken 122044 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1987 times.
✓ Branch 6 taken 122044 times.
✓ Branch 7 taken 1987 times.

      124031
          if (calc & Calc_Velocity || calc & Calc_Acceleration)
    
          {
    
              // Eccentric Anomaly Rate [rad/s]
    
      122044
              auto E_k_dot = n / (1 - e * std::cos(E_k));
    
              // True Anomaly Rate [rad/s]
    
      122044
              auto v_k_dot = E_k_dot * std::sqrt(1 - e * e) / (1 - e * std::cos(E_k));
    
              // Corrected Inclination Angle Rate [rad/s]
    
      122044
              auto i_k_dot = i_dot + 2 * v_k_dot * (Cis * std::cos(2 * Phi_k) - Cic * std::sin(2 * Phi_k));
    
              // Corrected Argument of Latitude Rate [rad/s]
    
      122044
              auto u_k_dot = v_k_dot + 2 * v_k_dot * (Cus * std::cos(2 * Phi_k) - Cuc * std::sin(2 * Phi_k));
    
              // Corrected Radius Rate [m/s]
    
      122044
              auto r_k_dot = e * A * E_k_dot * std::sin(E_k) + 2 * v_k_dot * (Crs * std::cos(2 * Phi_k) - Crc * std::sin(2 * Phi_k));
    
              // Longitude of Ascending Node Rate [rad/s]
    
      122044
              auto Omega_k_dot = Omega_dot - Omega_e_dot;
    
              // In-plane x velocity [m/s]
    
      122044
              auto vx_k_op = r_k_dot * std::cos(u_k) - r_k * u_k_dot * std::sin(u_k);
    
              // In-plane y velocity [m/s]
    
      122044
              auto vy_k_op = r_k_dot * std::sin(u_k) + r_k * u_k_dot * std::cos(u_k);
    
              // Earth-Fixed x velocity [m/s]
    
      122044
              auto vx_k = -x_k_op * Omega_k_dot * std::sin(Omega_k) + vx_k_op * std::cos(Omega_k) - vy_k_op * std::sin(Omega_k) * std::cos(i_k)
    
      122044
                          - y_k_op * (Omega_k_dot * std::cos(Omega_k) * std::cos(i_k) - i_k_dot * std::sin(Omega_k) * std::sin(i_k));
    
              // Earth-Fixed y velocity [m/s]
    
      122044
              auto vy_k = x_k_op * Omega_k_dot * std::cos(Omega_k) + vx_k_op * std::sin(Omega_k) + vy_k_op * std::cos(Omega_k) * std::cos(i_k)
    
      122044
                          - y_k_op * (Omega_k_dot * std::sin(Omega_k) * std::cos(i_k) + i_k_dot * std::cos(Omega_k) * std::sin(i_k));
    
              // Earth-Fixed z velocity [m/s]
    
      122044
              auto vz_k = vy_k_op * std::sin(i_k) + y_k_op * i_k_dot * std::cos(i_k);
    
        1/2✓ Branch 1 taken 122044 times.
✗ Branch 2 not taken.

      122044
              if (calc & Calc_Velocity)
    
              {
    
        1/2✓ Branch 1 taken 122044 times.
✗ Branch 2 not taken.

      122044
                  e_vel = Eigen::Vector3d{ vx_k, vy_k, vz_k };
    
              }
    
        2/2✓ Branch 1 taken 834 times.
✓ Branch 2 taken 121210 times.

      122044
              if (calc & Calc_Acceleration)
    
              {
    
                  // Oblate Earth acceleration Factor [m/s^2]
    
      834
                  auto F = -(3.0 / 2.0) * InsConst::GPS::J2 * (mu / std::pow(r_k, 2)) * std::pow(InsConst::GPS::R_E / r_k, 2);
    
                  // Earth-Fixed x acceleration [m/s^2]
    
      834
                  auto ax_k = -mu * (x_k / std::pow(r_k, 3)) + F * ((1.0 - 5.0 * std::pow(z_k / r_k, 2)) * (x_k / r_k))
    
      834
                              + 2 * vy_k * Omega_e_dot + x_k * std::pow(Omega_e_dot, 2);
    
                  // Earth-Fixed y acceleration [m/s^2]
    
      834
                  auto ay_k = -mu * (y_k / std::pow(r_k, 3)) + F * ((1.0 - 5.0 * std::pow(z_k / r_k, 2)) * (y_k / r_k))
    
      834
                              + 2 * vx_k * Omega_e_dot + y_k * std::pow(Omega_e_dot, 2);
    
                  // Earth-Fixed z acceleration [m/s^2]
    
      834
                  auto az_k = -mu * (z_k / std::pow(r_k, 3)) + F * ((3.0 - 5.0 * std::pow(z_k / r_k, 2)) * (z_k / r_k));
    
        1/2✓ Branch 1 taken 834 times.
✗ Branch 2 not taken.

      834
                  e_accel = Eigen::Vector3d{ ax_k, ay_k, az_k };
    
              }
    
          }
    
          return { .e_pos = e_pos,
    
                   .e_vel = e_vel,
    
        3/6✓ Branch 1 taken 124030 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124030 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 124029 times.
✗ Branch 8 not taken.

      248060
                   .e_accel = e_accel };
    
      }
    
      122372
      bool GalileoEphemeris::isHealthy() const
    
      {
    
      122372
          return svHealth.E5a_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
    
        1/2✓ Branch 0 taken 122354 times.
✗ Branch 1 not taken.

      122354
                 && svHealth.E5b_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
    
        1/2✓ Branch 0 taken 122354 times.
✗ Branch 1 not taken.

      122354
                 && svHealth.E1B_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
    
        2/2✓ Branch 0 taken 122333 times.
✓ Branch 1 taken 21 times.

      122354
                 && svHealth.E5a_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK
    
        1/2✓ Branch 0 taken 122333 times.
✗ Branch 1 not taken.

      122333
                 && svHealth.E5b_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK
    
        3/4✓ Branch 0 taken 122354 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 122333 times.
✗ Branch 3 not taken.

      244726
                 && svHealth.E1B_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK;
    
      }
    
      36304
      double GalileoEphemeris::calcSatellitePositionVariance() const
    
      {
    
      36304
          return std::pow(signalAccuracy, 2);
    
      }
    
      } // namespace NAV

Line	Branch	Exec	Source
1			// This file is part of INSTINCT, the INS Toolkit for Integrated
2			// Navigation Concepts and Training by the Institute of Navigation of
3			// the University of Stuttgart, Germany.
4			//
5			// This Source Code Form is subject to the terms of the Mozilla Public
6			// License, v. 2.0. If a copy of the MPL was not distributed with this
7			// file, You can obtain one at https://mozilla.org/MPL/2.0/.
8
9			#include "GalileoEphemeris.hpp"
10
11			#include "Navigation/Constants.hpp"
12			#include "Navigation/GNSS/Functions.hpp"
13
14			#include "util/Logger.hpp"
15
16			namespace NAV
17			{
18
19		28	GalileoEphemeris::GalileoEphemeris(const InsTime& toc)
20		28	: SatNavData(SatNavData::GalileoEphemeris, toc) {}
21
22		1633	GalileoEphemeris::GalileoEphemeris(const InsTime& toc, const InsTime& toe,
23			const size_t& IODnav,
24			const std::array<double, 3>& a,
25			const double& sqrt_A, const double& e, const double& i_0, const double& Omega_0, const double& omega, const double& M_0,
26			const double& delta_n, const double& Omega_dot, const double& i_dot, const double& Cus, const double& Cuc,
27			const double& Cis, const double& Cic, const double& Crs, const double& Crc,
28			const std::bitset<10>& dataSource, const double& signalAccuracy, const SvHealth& svHealth,
29		1633	const double& BGD_E1_E5a, const double& BGD_E1_E5b)
30			: SatNavData(SatNavData::GalileoEphemeris, toc),
31		1633	toc(toc),
32		1633	toe(toe),
33		1633	IODnav(IODnav),
34		1633	a(a),
35		1633	sqrt_A(sqrt_A),
36		1633	e(e),
37		1633	i_0(i_0),
38		1633	Omega_0(Omega_0),
39		1633	omega(omega),
40		1633	M_0(M_0),
41		1633	delta_n(delta_n),
42		1633	Omega_dot(Omega_dot),
43		1633	i_dot(i_dot),
44		1633	Cus(Cus),
45		1633	Cuc(Cuc),
46		1633	Cis(Cis),
47		1633	Cic(Cic),
48		1633	Crs(Crs),
49		1633	Crc(Crc),
50		1633	dataSource(dataSource),
51		1633	signalAccuracy(signalAccuracy),
52		1633	svHealth(svHealth),
53		1633	BGD_E1_E5a(BGD_E1_E5a),
54		1633	BGD_E1_E5b(BGD_E1_E5b) {}
55
56			#ifdef TESTING
57
58		8556	GalileoEphemeris::GalileoEphemeris(int32_t year, int32_t month, int32_t day, int32_t hour, int32_t minute, double second, double svClockBias, double svClockDrift, double svClockDriftRate,
59			double IODnav, double Crs, double delta_n, double M_0,
60			double Cuc, double e, double Cus, double sqrt_A,
61			double Toe, double Cic, double Omega_0, double Cis,
62			double i_0, double Crc, double omega, double Omega_dot,
63			double i_dot, double dataSource, double GALWeek, double /* spare1 */,
64			double signalAccuracy, double svHealth, double BGD_E1_E5a, double BGD_E1_E5b,
65		8556	double /* TransmissionTimeOfMessage /, double / spare2 /, double / spare3 /, double / spare4 */)
66	1/2 ✓ Branch 1 taken 8556 times. ✗ Branch 2 not taken.	8556	: SatNavData(SatNavData::GalileoEphemeris, InsTime(year, month, day, hour, minute, second, GST)),
67		8556	toc(refTime),
68	1/2 ✓ Branch 1 taken 8556 times. ✗ Branch 2 not taken.	8556	toe(InsTime(0, static_cast<int32_t>(GALWeek), Toe, GST)),
69		8556	IODnav(static_cast<size_t>(IODnav)),
70		8556	a({ svClockBias, svClockDrift, svClockDriftRate }),
71		8556	sqrt_A(sqrt_A),
72		8556	e(e),
73		8556	i_0(i_0),
74		8556	Omega_0(Omega_0),
75		8556	omega(omega),
76		8556	M_0(M_0),
77		8556	delta_n(delta_n),
78		8556	Omega_dot(Omega_dot),
79		8556	i_dot(i_dot),
80		8556	Cus(Cus),
81		8556	Cuc(Cuc),
82		8556	Cis(Cis),
83		8556	Cic(Cic),
84		8556	Crs(Crs),
85		8556	Crc(Crc),
86		8556	dataSource(static_cast<uint16_t>(dataSource)),
87		8556	signalAccuracy(signalAccuracy),
88		8556	svHealth({ .E5a_DataValidityStatus = static_cast<GalileoEphemeris::SvHealth::DataValidityStatus>((static_cast<uint16_t>(svHealth) & 0b000001000) >> 3),
89		8556	.E5b_DataValidityStatus = static_cast<GalileoEphemeris::SvHealth::DataValidityStatus>((static_cast<uint16_t>(svHealth) & 0b001000000) >> 6),
90		8556	.E1B_DataValidityStatus = static_cast<GalileoEphemeris::SvHealth::DataValidityStatus>((static_cast<uint16_t>(svHealth) & 0b000000001) >> 0),
91		8556	.E5a_SignalHealthStatus = static_cast<GalileoEphemeris::SvHealth::SignalHealthStatus>((static_cast<uint16_t>(svHealth) & 0b000110000) >> 4),
92		8556	.E5b_SignalHealthStatus = static_cast<GalileoEphemeris::SvHealth::SignalHealthStatus>((static_cast<uint16_t>(svHealth) & 0b110000000) >> 7),
93		8556	.E1B_SignalHealthStatus = static_cast<GalileoEphemeris::SvHealth::SignalHealthStatus>((static_cast<uint16_t>(svHealth) & 0b000000110) >> 1) }),
94		8556	BGD_E1_E5a(BGD_E1_E5a),
95	1/2 ✓ Branch 2 taken 8556 times. ✗ Branch 3 not taken.	17112	BGD_E1_E5b(BGD_E1_E5b)
96		8556	{}
97
98			#endif
99
100		122361	Clock::Corrections GalileoEphemeris::calcClockCorrections(const InsTime& recvTime, double dist, const Frequency& freq) const
101			{
102			LOG_DATA("Calc Sat Clock corrections at receiver time {}", recvTime.toGPSweekTow());
103			// Earth gravitational constant [m³/s²] (WGS 84 value of the earth's gravitational constant for GPS user)
104		122361	const auto mu = InsConst::GAL::MU;
105			// Relativistic constant F for clock corrections [s/√m] (-2*√µ/c²)
106		122361	const auto F = InsConst::GAL::F;
107
108			LOG_DATA(" toe {} (Time of ephemeris)", toe.toGPSweekTow());
109
110		122361	const auto A = sqrt_A * sqrt_A; // Semi-major axis [m]
111			LOG_DATA(" A {} [m] (Semi-major axis)", A);
112		122361	auto n_0 = std::sqrt(mu / std::pow(A, 3)); // Computed mean motion [rad/s]
113			LOG_DATA(" n_0 {} [rad/s] (Computed mean motion)", n_0);
114		122362	auto n = n_0 + delta_n; // Corrected mean motion [rad/s]
115			LOG_DATA(" n {} [rad/s] (Corrected mean motion)", n);
116
117			// Time at transmission
118	2/4 ✓ Branch 2 taken 122362 times. ✗ Branch 3 not taken. ✓ Branch 5 taken 122363 times. ✗ Branch 6 not taken.	122362	InsTime transTime0 = recvTime - std::chrono::duration<double>(dist / InsConst::C);
119
120		122363	InsTime transTime = transTime0;
121			LOG_DATA(" Iterating Time at transmission");
122		122363	double dt_sv = 0.0;
123		122363	double clkDrift = 0.0;
124
125	2/2 ✓ Branch 0 taken 244726 times. ✓ Branch 1 taken 122360 times.	367086	for (size_t i = 0; i < 2; i++)
126			{
127			LOG_DATA(" transTime {} (Time at transmission)", transTime.toGPSweekTow());
128
129			// [s]
130	1/2 ✓ Branch 1 taken 244722 times. ✗ Branch 2 not taken.	244726	auto t_minus_toc = static_cast<double>((transTime - toc).count());
131			LOG_DATA(" transTime - toc {} [s]", t_minus_toc);
132
133			// Time difference from ephemeris reference epoch [s]
134	1/2 ✓ Branch 1 taken 244723 times. ✗ Branch 2 not taken.	244722	double t_k = static_cast<double>((transTime - toe).count());
135			LOG_DATA(" transTime - toe {} [s] (t_k = Time difference from ephemeris reference epoch)", t_k);
136
137			// Mean anomaly [rad]
138		244723	auto M_k = M_0 + n * t_k;
139			LOG_DATA(" M_k {} [s] (Mean anomaly)", M_k);
140
141			// Eccentric anomaly [rad]
142		244723	double E_k = M_k;
143		244723	double E_k_old = 0.0;
144
145	5/6 ✓ Branch 1 taken 972905 times. ✓ Branch 2 taken 244712 times. ✓ Branch 3 taken 972905 times. ✗ Branch 4 not taken. ✓ Branch 5 taken 972902 times. ✓ Branch 6 taken 244715 times.	1217625	for (size_t i = 0; std::abs(E_k - E_k_old) > 1e-13 && i < 10; i++)
146			{
147		972902	E_k_old = E_k; // Kepler’s equation ( Mk = E_k − e sin E_k ) may be solved for Eccentric anomaly (E_k) by iteration:
148		972902	E_k = M_k + e * sin(E_k);
149			}
150
151			// Relativistic correction term [s]
152		244715	double dt_r = F * e * sqrt_A * std::sin(E_k);
153			LOG_DATA(" dt_r {} [s] (Relativistic correction term)", dt_r);
154
155			// SV PRN code phase time offset [s]
156		244715	dt_sv = a[0] + a[1] * t_minus_toc + a[2] * std::pow(t_minus_toc, 2) + dt_r;
157
158			// See GAL-ICD-2.0 GAL ICD, ch. 5.1.5, p.47
159	3/4 ✓ Branch 2 taken 244720 times. ✗ Branch 3 not taken. ✓ Branch 5 taken 61773 times. ✓ Branch 6 taken 182950 times.	244722	dt_sv -= ratioFreqSquared(E01, freq, -128, -128) * (freq == E05 ? BGD_E1_E5a : BGD_E1_E5b); // TODO: Check again
160
161			LOG_DATA(" dt_sv {} [s] (SV PRN code phase time offset)", dt_sv);
162
163			// Groves ch. 9.3.1, eq. 9.78, p. 391
164		244723	clkDrift = a[1] + a[2] / 2.0 * t_minus_toc;
165
166			// Correct transmit time for the satellite clock bias
167	2/4 ✓ Branch 2 taken 244721 times. ✗ Branch 3 not taken. ✓ Branch 5 taken 244723 times. ✗ Branch 6 not taken.	244724	transTime = transTime0 - std::chrono::duration<double>(dt_sv);
168			}
169			LOG_DATA(" transTime {} (Time at transmission)", transTime.toGPSweekTow());
170
171		122360	return { .transmitTime = transTime, .bias = dt_sv, .drift = clkDrift };
172			}
173
174		124029	Orbit::PosVelAccel GalileoEphemeris::calcSatelliteData(const InsTime& transTime, Orbit::Calc calc) const
175			{
176	2/4 ✓ Branch 1 taken 124030 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 124031 times. ✗ Branch 5 not taken.	124029	Eigen::Vector3d e_pos = Eigen::Vector3d::Zero();
177	2/4 ✓ Branch 1 taken 124031 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 124031 times. ✗ Branch 5 not taken.	124031	Eigen::Vector3d e_vel = Eigen::Vector3d::Zero();
178	2/4 ✓ Branch 1 taken 124031 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 124031 times. ✗ Branch 5 not taken.	124031	Eigen::Vector3d e_accel = Eigen::Vector3d::Zero();
179
180			LOG_DATA("Calc Sat Position at transmit time {}", transTime.toGPSweekTow());
181			// Earth gravitational constant [m³/s²] (WGS 84 value of the earth's gravitational constant for GPS user)
182		124031	const auto mu = InsConst::GAL::MU;
183			// Earth angular velocity [rad/s] (WGS 84 value of the earth's rotation rate)
184		124031	const auto Omega_e_dot = InsConst::GAL::omega_ie;
185
186			LOG_DATA(" toe {} (Time of ephemeris)", toe.toGPSweekTow());
187
188		124031	const auto A = sqrt_A * sqrt_A; // Semi-major axis [m]
189			LOG_DATA(" A {} [m] (Semi-major axis)", A);
190		124031	auto n_0 = std::sqrt(mu / std::pow(A, 3)); // Computed mean motion [rad/s]
191			LOG_DATA(" n_0 {} [rad/s] (Computed mean motion)", n_0);
192		124031	auto n = n_0 + delta_n; // Corrected mean motion [rad/s]
193			LOG_DATA(" n {} [rad/s] (Corrected mean motion)", n);
194
195			// Eccentric anomaly [rad]
196		124031	double E_k = 0.0;
197
198			// Time difference from ephemeris reference epoch [s]
199	1/2 ✓ Branch 1 taken 124031 times. ✗ Branch 2 not taken.	124031	double t_k = static_cast<double>((transTime - toe).count());
200			LOG_DATA(" t_k {} [s] (Time difference from ephemeris reference epoch)", t_k);
201
202			// Mean anomaly [rad]
203		124031	auto M_k = M_0 + n * t_k;
204			LOG_DATA(" M_k {} [s] (Mean anomaly)", M_k);
205
206		124031	E_k = M_k; // Initial Value [rad]
207		124031	double E_k_old = 0.0;
208			LOG_DATA(" Iterating E_k");
209			LOG_DATA(" E_k {} [rad] (Eccentric anomaly)", E_k);
210	5/6 ✓ Branch 1 taken 306105 times. ✓ Branch 2 taken 124028 times. ✓ Branch 3 taken 306105 times. ✗ Branch 4 not taken. ✓ Branch 5 taken 306104 times. ✓ Branch 6 taken 124029 times.	430135	for (size_t i = 0; std::abs(E_k - E_k_old) > 1e-13 && i < 10; i++)
211			{
212		306104	E_k_old = E_k; // Kepler’s equation ( Mk = E_k − e sin E_k ) may be solved for Eccentric anomaly (E_k) by iteration:
213		306104	E_k = E_k + (M_k - E_k + e * std::sin(E_k)) / (1 - e * std::cos(E_k)); // - Refined Value, minimum of three iterations, (j=1,2,3)
214			LOG_DATA(" E_k {} [rad] (Eccentric anomaly)", E_k); // - Final Value (radians)
215			}
216
217			// auto v_k = 2.0 * std::atan(std::sqrt((1.0 + e) / (1.0 - e)) * std::tan(E_k / 2.0)); // True Anomaly (unambiguous quadrant) [rad] (GPS ICD algorithm)
218			// auto v_k = std::atan2(std::sqrt(1 - e * e) * std::sin(E_k) / (1 - e * std::cos(E_k)), (std::cos(E_k) - e) / (1 - e * std::cos(E_k))); // True Anomaly [rad] (GALILEO ICD algorithm)
219		124029	auto v_k = std::atan2(std::sqrt(1 - e * e) * std::sin(E_k), (std::cos(E_k) - e)); // True Anomaly [rad] // simplified, since the denominators cancel out
220			LOG_DATA(" v_k {} [rad] (True Anomaly (unambiguous quadrant))", v_k);
221		124029	auto Phi_k = v_k + omega; // Argument of Latitude [rad]
222			LOG_DATA(" Phi_k {} [rad] (Argument of Latitude)", Phi_k);
223
224			// Second Harmonic Perturbations
225		124029	auto delta_u_k = Cus * std::sin(2 * Phi_k) + Cuc * std::cos(2 * Phi_k); // Argument of Latitude Correction [rad]
226			LOG_DATA(" delta_u_k {} [rad] (Argument of Latitude Correction)", delta_u_k);
227		124029	auto delta_r_k = Crs * std::sin(2 * Phi_k) + Crc * std::cos(2 * Phi_k); // Radius Correction [m]
228			LOG_DATA(" delta_r_k {} [m] (Radius Correction)", delta_r_k);
229		124029	auto delta_i_k = Cis * std::sin(2 * Phi_k) + Cic * std::cos(2 * Phi_k); // Inclination Correction [rad]
230			LOG_DATA(" delta_i_k {} [rad] (Inclination Correction)", delta_i_k);
231
232		124029	auto u_k = Phi_k + delta_u_k; // Corrected Argument of Latitude [rad]
233			LOG_DATA(" u_k {} [rad] (Corrected Argument of Latitude)", u_k);
234		124029	auto r_k = A * (1 - e * std::cos(E_k)) + delta_r_k; // Corrected Radius [m]
235			LOG_DATA(" r_k {} [m] (Corrected Radius)", r_k);
236		124029	auto i_k = i_0 + delta_i_k + i_dot * t_k; // Corrected Inclination [rad]
237			LOG_DATA(" i_k {} [rad] (Corrected Inclination)", i_k);
238
239		124029	auto x_k_op = r_k * std::cos(u_k); // Position in orbital plane [m]
240			LOG_DATA(" x_k_op {} [m] (Position in orbital plane)", x_k_op);
241		124029	auto y_k_op = r_k * std::sin(u_k); // Position in orbital plane [m]
242			LOG_DATA(" y_k_op {} [m] (Position in orbital plane)", y_k_op);
243
244			// Corrected longitude of ascending node [rad]
245	1/2 ✓ Branch 2 taken 124031 times. ✗ Branch 3 not taken.	124029	auto Omega_k = Omega_0 + (Omega_dot - Omega_e_dot) * t_k - Omega_e_dot * static_cast<double>(toe.toGPSweekTow().tow);
246			LOG_DATA(" Omega_k {} [rad] (Corrected longitude of ascending node)", Omega_k);
247
248			// Earth-fixed x coordinates [m]
249		124031	auto x_k = x_k_op * std::cos(Omega_k) - y_k_op * std::cos(i_k) * std::sin(Omega_k);
250			LOG_DATA(" x_k {} [m] (Earth-fixed x coordinates)", x_k);
251			// Earth-fixed y coordinates [m]
252		124031	auto y_k = x_k_op * std::sin(Omega_k) + y_k_op * std::cos(i_k) * std::cos(Omega_k);
253			LOG_DATA(" y_k {} [m] (Earth-fixed y coordinates)", y_k);
254			// Earth-fixed z coordinates [m]
255		124031	auto z_k = y_k_op * std::sin(i_k);
256			LOG_DATA(" z_k {} [m] (Earth-fixed z coordinates)", z_k);
257
258	1/2 ✓ Branch 1 taken 124031 times. ✗ Branch 2 not taken.	124031	e_pos = Eigen::Vector3d{ x_k, y_k, z_k };
259
260	5/6 ✓ Branch 1 taken 1987 times. ✓ Branch 2 taken 122044 times. ✗ Branch 4 not taken. ✓ Branch 5 taken 1987 times. ✓ Branch 6 taken 122044 times. ✓ Branch 7 taken 1987 times.	124031	if (calc & Calc_Velocity \|\| calc & Calc_Acceleration)
261			{
262			// Eccentric Anomaly Rate [rad/s]
263		122044	auto E_k_dot = n / (1 - e * std::cos(E_k));
264			// True Anomaly Rate [rad/s]
265		122044	auto v_k_dot = E_k_dot * std::sqrt(1 - e * e) / (1 - e * std::cos(E_k));
266			// Corrected Inclination Angle Rate [rad/s]
267		122044	auto i_k_dot = i_dot + 2 * v_k_dot * (Cis * std::cos(2 * Phi_k) - Cic * std::sin(2 * Phi_k));
268			// Corrected Argument of Latitude Rate [rad/s]
269		122044	auto u_k_dot = v_k_dot + 2 * v_k_dot * (Cus * std::cos(2 * Phi_k) - Cuc * std::sin(2 * Phi_k));
270			// Corrected Radius Rate [m/s]
271		122044	auto r_k_dot = e * A * E_k_dot * std::sin(E_k) + 2 * v_k_dot * (Crs * std::cos(2 * Phi_k) - Crc * std::sin(2 * Phi_k));
272			// Longitude of Ascending Node Rate [rad/s]
273		122044	auto Omega_k_dot = Omega_dot - Omega_e_dot;
274			// In-plane x velocity [m/s]
275		122044	auto vx_k_op = r_k_dot * std::cos(u_k) - r_k * u_k_dot * std::sin(u_k);
276			// In-plane y velocity [m/s]
277		122044	auto vy_k_op = r_k_dot * std::sin(u_k) + r_k * u_k_dot * std::cos(u_k);
278			// Earth-Fixed x velocity [m/s]
279		122044	auto vx_k = -x_k_op * Omega_k_dot * std::sin(Omega_k) + vx_k_op * std::cos(Omega_k) - vy_k_op * std::sin(Omega_k) * std::cos(i_k)
280		122044	- y_k_op * (Omega_k_dot * std::cos(Omega_k) * std::cos(i_k) - i_k_dot * std::sin(Omega_k) * std::sin(i_k));
281			// Earth-Fixed y velocity [m/s]
282		122044	auto vy_k = x_k_op * Omega_k_dot * std::cos(Omega_k) + vx_k_op * std::sin(Omega_k) + vy_k_op * std::cos(Omega_k) * std::cos(i_k)
283		122044	- y_k_op * (Omega_k_dot * std::sin(Omega_k) * std::cos(i_k) + i_k_dot * std::cos(Omega_k) * std::sin(i_k));
284			// Earth-Fixed z velocity [m/s]
285		122044	auto vz_k = vy_k_op * std::sin(i_k) + y_k_op * i_k_dot * std::cos(i_k);
286
287	1/2 ✓ Branch 1 taken 122044 times. ✗ Branch 2 not taken.	122044	if (calc & Calc_Velocity)
288			{
289	1/2 ✓ Branch 1 taken 122044 times. ✗ Branch 2 not taken.	122044	e_vel = Eigen::Vector3d{ vx_k, vy_k, vz_k };
290			}
291
292	2/2 ✓ Branch 1 taken 834 times. ✓ Branch 2 taken 121210 times.	122044	if (calc & Calc_Acceleration)
293			{
294			// Oblate Earth acceleration Factor [m/s^2]
295		834	auto F = -(3.0 / 2.0) * InsConst::GPS::J2 * (mu / std::pow(r_k, 2)) * std::pow(InsConst::GPS::R_E / r_k, 2);
296			// Earth-Fixed x acceleration [m/s^2]
297		834	auto ax_k = -mu * (x_k / std::pow(r_k, 3)) + F * ((1.0 - 5.0 * std::pow(z_k / r_k, 2)) * (x_k / r_k))
298		834	+ 2 * vy_k * Omega_e_dot + x_k * std::pow(Omega_e_dot, 2);
299			// Earth-Fixed y acceleration [m/s^2]
300		834	auto ay_k = -mu * (y_k / std::pow(r_k, 3)) + F * ((1.0 - 5.0 * std::pow(z_k / r_k, 2)) * (y_k / r_k))
301		834	+ 2 * vx_k * Omega_e_dot + y_k * std::pow(Omega_e_dot, 2);
302			// Earth-Fixed z acceleration [m/s^2]
303		834	auto az_k = -mu * (z_k / std::pow(r_k, 3)) + F * ((3.0 - 5.0 * std::pow(z_k / r_k, 2)) * (z_k / r_k));
304
305	1/2 ✓ Branch 1 taken 834 times. ✗ Branch 2 not taken.	834	e_accel = Eigen::Vector3d{ ax_k, ay_k, az_k };
306			}
307			}
308
309			return { .e_pos = e_pos,
310			.e_vel = e_vel,
311	3/6 ✓ Branch 1 taken 124030 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 124030 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 124029 times. ✗ Branch 8 not taken.	248060	.e_accel = e_accel };
312			}
313
314		122372	bool GalileoEphemeris::isHealthy() const
315			{
316		122372	return svHealth.E5a_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
317	1/2 ✓ Branch 0 taken 122354 times. ✗ Branch 1 not taken.	122354	&& svHealth.E5b_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
318	1/2 ✓ Branch 0 taken 122354 times. ✗ Branch 1 not taken.	122354	&& svHealth.E1B_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
319	2/2 ✓ Branch 0 taken 122333 times. ✓ Branch 1 taken 21 times.	122354	&& svHealth.E5a_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK
320	1/2 ✓ Branch 0 taken 122333 times. ✗ Branch 1 not taken.	122333	&& svHealth.E5b_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK
321	3/4 ✓ Branch 0 taken 122354 times. ✓ Branch 1 taken 18 times. ✓ Branch 2 taken 122333 times. ✗ Branch 3 not taken.	244726	&& svHealth.E1B_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK;
322			}
323
324		36304	double GalileoEphemeris::calcSatellitePositionVariance() const
325			{
326		36304	return std::pow(signalAccuracy, 2);
327			}
328
329			} // namespace NAV
330