INSTINCT Code Coverage Report

Directory:	src/
File:	Navigation/GNSS/Satellite/Ephemeris/GalileoEphemeris.cpp
Date:	2025-02-07 16:54:41

	Exec	Total	Coverage
Lines:	150	152	98.7%
Functions:	6	7	85.7%
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) {}
    
      1607
      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,
    
      1607
                                         const double& BGD_E1_E5a, const double& BGD_E1_E5b)
    
          : SatNavData(SatNavData::GalileoEphemeris, toc),
    
      1607
            toc(toc),
    
      1607
            toe(toe),
    
      1607
            IODnav(IODnav),
    
      1607
            a(a),
    
      1607
            sqrt_A(sqrt_A),
    
      1607
            e(e),
    
      1607
            i_0(i_0),
    
      1607
            Omega_0(Omega_0),
    
      1607
            omega(omega),
    
      1607
            M_0(M_0),
    
      1607
            delta_n(delta_n),
    
      1607
            Omega_dot(Omega_dot),
    
      1607
            i_dot(i_dot),
    
      1607
            Cus(Cus),
    
      1607
            Cuc(Cuc),
    
      1607
            Cis(Cis),
    
      1607
            Cic(Cic),
    
      1607
            Crs(Crs),
    
      1607
            Crc(Crc),
    
      1607
            dataSource(dataSource),
    
      1607
            signalAccuracy(signalAccuracy),
    
      1607
            svHealth(svHealth),
    
      1607
            BGD_E1_E5a(BGD_E1_E5a),
    
      1607
            BGD_E1_E5b(BGD_E1_E5b) {}
    
      #ifdef TESTING
    
      7802
      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,
    
      7802
                                         double /* TransmissionTimeOfMessage */, double /* spare2 */, double /* spare3 */, double /* spare4 */)
    
        1/2✓ Branch 1 taken 7802 times.
✗ Branch 2 not taken.

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

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

      15604
            BGD_E1_E5b(BGD_E1_E5b)
    
      7802
      {}
    
      #endif
    
      1987
      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)
    
      1987
          const auto mu = InsConst::GAL::MU;
    
          // Relativistic constant F for clock corrections [s/√m] (-2*√µ/c²)
    
      1987
          const auto F = InsConst::GAL::F;
    
          LOG_DATA("    toe {} (Time of ephemeris)", toe.toGPSweekTow());
    
      1987
          const auto A = sqrt_A * sqrt_A; // Semi-major axis [m]
    
          LOG_DATA("    A {} [m] (Semi-major axis)", A);
    
      1987
          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);
    
      1987
          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 1987 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1987 times.
✗ Branch 6 not taken.

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

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

      3974
              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 3974 times.
✗ Branch 2 not taken.

      3974
              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]
    
      3974
              auto M_k = M_0 + n * t_k;
    
              LOG_DATA("      M_k {} [s] (Mean anomaly)", M_k);
    
              // Eccentric anomaly [rad]
    
      3974
              double E_k = M_k;
    
      3974
              double E_k_old = 0.0;
    
        5/6✓ Branch 1 taken 15550 times.
✓ Branch 2 taken 3974 times.
✓ Branch 3 taken 15550 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 15550 times.
✓ Branch 6 taken 3974 times.

      19524
              for (size_t i = 0; std::abs(E_k - E_k_old) > 1e-13 && i < 10; i++)
    
              {
    
      15550
                  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:
    
      15550
                  E_k = M_k + e * sin(E_k);
    
              }
    
              // Relativistic correction term [s]
    
      3974
              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]
    
      3974
              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 3974 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1586 times.
✓ Branch 6 taken 2388 times.

      3974
              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
    
      3974
              clkDrift = a[1] + a[2] / 2.0 * t_minus_toc;
    
              // Correct transmit time for the satellite clock bias
    
        2/4✓ Branch 2 taken 3974 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3974 times.
✗ Branch 6 not taken.

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

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

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

      3655
          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)
    
      3655
          const auto mu = InsConst::GAL::MU;
    
          // Earth angular velocity [rad/s] (WGS 84 value of the earth's rotation rate)
    
      3655
          const auto Omega_e_dot = InsConst::GAL::omega_ie;
    
          LOG_DATA("    toe {} (Time of ephemeris)", toe.toGPSweekTow());
    
      3655
          const auto A = sqrt_A * sqrt_A; // Semi-major axis [m]
    
          LOG_DATA("    A {} [m] (Semi-major axis)", A);
    
      3655
          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);
    
      3655
          auto n = n_0 + delta_n; // Corrected mean motion [rad/s]
    
          LOG_DATA("    n {} [rad/s] (Corrected mean motion)", n);
    
          // Eccentric anomaly [rad]
    
      3655
          double E_k = 0.0;
    
          // Time difference from ephemeris reference epoch [s]
    
        1/2✓ Branch 1 taken 3655 times.
✗ Branch 2 not taken.

      3655
          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]
    
      3655
          auto M_k = M_0 + n * t_k;
    
          LOG_DATA("    M_k {} [s] (Mean anomaly)", M_k);
    
      3655
          E_k = M_k; // Initial Value [rad]
    
      3655
          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 9714 times.
✓ Branch 2 taken 3655 times.
✓ Branch 3 taken 9714 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 9714 times.
✓ Branch 6 taken 3655 times.

      13369
          for (size_t i = 0; std::abs(E_k - E_k_old) > 1e-13 && i < 10; i++)
    
          {
    
      9714
              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:
    
      9714
              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)
    
      3655
          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);
    
      3655
          auto Phi_k = v_k + omega; // Argument of Latitude [rad]
    
          LOG_DATA("    Phi_k {} [rad] (Argument of Latitude)", Phi_k);
    
          // Second Harmonic Perturbations
    
      3655
          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);
    
      3655
          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);
    
      3655
          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);
    
      3655
          auto u_k = Phi_k + delta_u_k; // Corrected Argument of Latitude [rad]
    
          LOG_DATA("    u_k {} [rad] (Corrected Argument of Latitude)", u_k);
    
      3655
          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);
    
      3655
          auto i_k = i_0 + delta_i_k + i_dot * t_k; // Corrected Inclination [rad]
    
          LOG_DATA("    i_k {} [rad] (Corrected Inclination)", i_k);
    
      3655
          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);
    
      3655
          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 3655 times.
✗ Branch 3 not taken.

      3655
          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]
    
      3655
          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]
    
      3655
          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]
    
      3655
          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 3655 times.
✗ Branch 2 not taken.

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

      3655
          if (calc & Calc_Velocity || calc & Calc_Acceleration)
    
          {
    
              // Eccentric Anomaly Rate [rad/s]
    
      1668
              auto E_k_dot = n / (1 - e * std::cos(E_k));
    
              // True Anomaly Rate [rad/s]
    
      1668
              auto v_k_dot = E_k_dot * std::sqrt(1 - e * e) / (1 - e * std::cos(E_k));
    
              // Corrected Inclination Angle Rate [rad/s]
    
      1668
              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]
    
      1668
              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]
    
      1668
              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]
    
      1668
              auto Omega_k_dot = Omega_dot - Omega_e_dot;
    
              // In-plane x velocity [m/s]
    
      1668
              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]
    
      1668
              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]
    
      1668
              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)
    
      1668
                          - 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]
    
      1668
              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)
    
      1668
                          - 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]
    
      1668
              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 1668 times.
✗ Branch 2 not taken.

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

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

      1668
              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 3655 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3655 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3655 times.
✗ Branch 8 not taken.

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

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

      1982
                 && svHealth.E1B_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
    
        2/2✓ Branch 0 taken 1962 times.
✓ Branch 1 taken 20 times.

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

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

      3982
                 && svHealth.E1B_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK;
    
      }
    
      ✗
      double GalileoEphemeris::calcSatellitePositionVariance() const
    
      {
    
      ✗
          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		1607	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		1607	const double& BGD_E1_E5a, const double& BGD_E1_E5b)
30			: SatNavData(SatNavData::GalileoEphemeris, toc),
31		1607	toc(toc),
32		1607	toe(toe),
33		1607	IODnav(IODnav),
34		1607	a(a),
35		1607	sqrt_A(sqrt_A),
36		1607	e(e),
37		1607	i_0(i_0),
38		1607	Omega_0(Omega_0),
39		1607	omega(omega),
40		1607	M_0(M_0),
41		1607	delta_n(delta_n),
42		1607	Omega_dot(Omega_dot),
43		1607	i_dot(i_dot),
44		1607	Cus(Cus),
45		1607	Cuc(Cuc),
46		1607	Cis(Cis),
47		1607	Cic(Cic),
48		1607	Crs(Crs),
49		1607	Crc(Crc),
50		1607	dataSource(dataSource),
51		1607	signalAccuracy(signalAccuracy),
52		1607	svHealth(svHealth),
53		1607	BGD_E1_E5a(BGD_E1_E5a),
54		1607	BGD_E1_E5b(BGD_E1_E5b) {}
55
56			#ifdef TESTING
57
58		7802	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		7802	double /* TransmissionTimeOfMessage /, double / spare2 /, double / spare3 /, double / spare4 */)
66	1/2 ✓ Branch 1 taken 7802 times. ✗ Branch 2 not taken.	7802	: SatNavData(SatNavData::GalileoEphemeris, InsTime(year, month, day, hour, minute, second, GST)),
67		7802	toc(refTime),
68	1/2 ✓ Branch 1 taken 7802 times. ✗ Branch 2 not taken.	7802	toe(InsTime(0, static_cast<int32_t>(GALWeek), Toe, GST)),
69		7802	IODnav(static_cast<size_t>(IODnav)),
70		7802	a({ svClockBias, svClockDrift, svClockDriftRate }),
71		7802	sqrt_A(sqrt_A),
72		7802	e(e),
73		7802	i_0(i_0),
74		7802	Omega_0(Omega_0),
75		7802	omega(omega),
76		7802	M_0(M_0),
77		7802	delta_n(delta_n),
78		7802	Omega_dot(Omega_dot),
79		7802	i_dot(i_dot),
80		7802	Cus(Cus),
81		7802	Cuc(Cuc),
82		7802	Cis(Cis),
83		7802	Cic(Cic),
84		7802	Crs(Crs),
85		7802	Crc(Crc),
86		7802	dataSource(static_cast<uint16_t>(dataSource)),
87		7802	signalAccuracy(signalAccuracy),
88		7802	svHealth({ .E5a_DataValidityStatus = static_cast<GalileoEphemeris::SvHealth::DataValidityStatus>((static_cast<uint16_t>(svHealth) & 0b000001000) >> 3),
89		7802	.E5b_DataValidityStatus = static_cast<GalileoEphemeris::SvHealth::DataValidityStatus>((static_cast<uint16_t>(svHealth) & 0b001000000) >> 6),
90		7802	.E1B_DataValidityStatus = static_cast<GalileoEphemeris::SvHealth::DataValidityStatus>((static_cast<uint16_t>(svHealth) & 0b000000001) >> 0),
91		7802	.E5a_SignalHealthStatus = static_cast<GalileoEphemeris::SvHealth::SignalHealthStatus>((static_cast<uint16_t>(svHealth) & 0b000110000) >> 4),
92		7802	.E5b_SignalHealthStatus = static_cast<GalileoEphemeris::SvHealth::SignalHealthStatus>((static_cast<uint16_t>(svHealth) & 0b110000000) >> 7),
93		7802	.E1B_SignalHealthStatus = static_cast<GalileoEphemeris::SvHealth::SignalHealthStatus>((static_cast<uint16_t>(svHealth) & 0b000000110) >> 1) }),
94		7802	BGD_E1_E5a(BGD_E1_E5a),
95	1/2 ✓ Branch 2 taken 7802 times. ✗ Branch 3 not taken.	15604	BGD_E1_E5b(BGD_E1_E5b)
96		7802	{}
97
98			#endif
99
100		1987	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		1987	const auto mu = InsConst::GAL::MU;
105			// Relativistic constant F for clock corrections [s/√m] (-2*√µ/c²)
106		1987	const auto F = InsConst::GAL::F;
107
108			LOG_DATA(" toe {} (Time of ephemeris)", toe.toGPSweekTow());
109
110		1987	const auto A = sqrt_A * sqrt_A; // Semi-major axis [m]
111			LOG_DATA(" A {} [m] (Semi-major axis)", A);
112		1987	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		1987	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 1987 times. ✗ Branch 3 not taken. ✓ Branch 5 taken 1987 times. ✗ Branch 6 not taken.	1987	InsTime transTime0 = recvTime - std::chrono::duration<double>(dist / InsConst::C);
119
120		1987	InsTime transTime = transTime0;
121			LOG_DATA(" Iterating Time at transmission");
122		1987	double dt_sv = 0.0;
123		1987	double clkDrift = 0.0;
124
125	2/2 ✓ Branch 0 taken 3974 times. ✓ Branch 1 taken 1987 times.	5961	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 3974 times. ✗ Branch 2 not taken.	3974	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 3974 times. ✗ Branch 2 not taken.	3974	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		3974	auto M_k = M_0 + n * t_k;
139			LOG_DATA(" M_k {} [s] (Mean anomaly)", M_k);
140
141			// Eccentric anomaly [rad]
142		3974	double E_k = M_k;
143		3974	double E_k_old = 0.0;
144
145	5/6 ✓ Branch 1 taken 15550 times. ✓ Branch 2 taken 3974 times. ✓ Branch 3 taken 15550 times. ✗ Branch 4 not taken. ✓ Branch 5 taken 15550 times. ✓ Branch 6 taken 3974 times.	19524	for (size_t i = 0; std::abs(E_k - E_k_old) > 1e-13 && i < 10; i++)
146			{
147		15550	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		15550	E_k = M_k + e * sin(E_k);
149			}
150
151			// Relativistic correction term [s]
152		3974	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		3974	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 3974 times. ✗ Branch 3 not taken. ✓ Branch 5 taken 1586 times. ✓ Branch 6 taken 2388 times.	3974	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		3974	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 3974 times. ✗ Branch 3 not taken. ✓ Branch 5 taken 3974 times. ✗ Branch 6 not taken.	3974	transTime = transTime0 - std::chrono::duration<double>(dt_sv);
168			}
169			LOG_DATA(" transTime {} (Time at transmission)", transTime.toGPSweekTow());
170
171		1987	return { .transmitTime = transTime, .bias = dt_sv, .drift = clkDrift };
172			}
173
174		3655	Orbit::PosVelAccel GalileoEphemeris::calcSatelliteData(const InsTime& transTime, Orbit::Calc calc) const
175			{
176	2/4 ✓ Branch 1 taken 3655 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 3655 times. ✗ Branch 5 not taken.	3655	Eigen::Vector3d e_pos = Eigen::Vector3d::Zero();
177	2/4 ✓ Branch 1 taken 3655 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 3655 times. ✗ Branch 5 not taken.	3655	Eigen::Vector3d e_vel = Eigen::Vector3d::Zero();
178	2/4 ✓ Branch 1 taken 3655 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 3655 times. ✗ Branch 5 not taken.	3655	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		3655	const auto mu = InsConst::GAL::MU;
183			// Earth angular velocity [rad/s] (WGS 84 value of the earth's rotation rate)
184		3655	const auto Omega_e_dot = InsConst::GAL::omega_ie;
185
186			LOG_DATA(" toe {} (Time of ephemeris)", toe.toGPSweekTow());
187
188		3655	const auto A = sqrt_A * sqrt_A; // Semi-major axis [m]
189			LOG_DATA(" A {} [m] (Semi-major axis)", A);
190		3655	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		3655	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		3655	double E_k = 0.0;
197
198			// Time difference from ephemeris reference epoch [s]
199	1/2 ✓ Branch 1 taken 3655 times. ✗ Branch 2 not taken.	3655	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		3655	auto M_k = M_0 + n * t_k;
204			LOG_DATA(" M_k {} [s] (Mean anomaly)", M_k);
205
206		3655	E_k = M_k; // Initial Value [rad]
207		3655	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 9714 times. ✓ Branch 2 taken 3655 times. ✓ Branch 3 taken 9714 times. ✗ Branch 4 not taken. ✓ Branch 5 taken 9714 times. ✓ Branch 6 taken 3655 times.	13369	for (size_t i = 0; std::abs(E_k - E_k_old) > 1e-13 && i < 10; i++)
211			{
212		9714	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		9714	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		3655	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		3655	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		3655	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		3655	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		3655	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		3655	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		3655	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		3655	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		3655	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		3655	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 3655 times. ✗ Branch 3 not taken.	3655	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		3655	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		3655	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		3655	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 3655 times. ✗ Branch 2 not taken.	3655	e_pos = Eigen::Vector3d{ x_k, y_k, z_k };
259
260	5/6 ✓ Branch 1 taken 1987 times. ✓ Branch 2 taken 1668 times. ✗ Branch 4 not taken. ✓ Branch 5 taken 1987 times. ✓ Branch 6 taken 1668 times. ✓ Branch 7 taken 1987 times.	3655	if (calc & Calc_Velocity \|\| calc & Calc_Acceleration)
261			{
262			// Eccentric Anomaly Rate [rad/s]
263		1668	auto E_k_dot = n / (1 - e * std::cos(E_k));
264			// True Anomaly Rate [rad/s]
265		1668	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		1668	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		1668	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		1668	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		1668	auto Omega_k_dot = Omega_dot - Omega_e_dot;
274			// In-plane x velocity [m/s]
275		1668	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		1668	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		1668	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		1668	- 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		1668	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		1668	- 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		1668	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 1668 times. ✗ Branch 2 not taken.	1668	if (calc & Calc_Velocity)
288			{
289	1/2 ✓ Branch 1 taken 1668 times. ✗ Branch 2 not taken.	1668	e_vel = Eigen::Vector3d{ vx_k, vy_k, vz_k };
290			}
291
292	2/2 ✓ Branch 1 taken 834 times. ✓ Branch 2 taken 834 times.	1668	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 3655 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 3655 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 3655 times. ✗ Branch 8 not taken.	7310	.e_accel = e_accel };
312			}
313
314		2000	bool GalileoEphemeris::isHealthy() const
315			{
316		2000	return svHealth.E5a_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
317	1/2 ✓ Branch 0 taken 1982 times. ✗ Branch 1 not taken.	1982	&& svHealth.E5b_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
318	1/2 ✓ Branch 0 taken 1982 times. ✗ Branch 1 not taken.	1982	&& svHealth.E1B_DataValidityStatus == SvHealth::DataValidityStatus::NavigationDataValid
319	2/2 ✓ Branch 0 taken 1962 times. ✓ Branch 1 taken 20 times.	1982	&& svHealth.E5a_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK
320	1/2 ✓ Branch 0 taken 1962 times. ✗ Branch 1 not taken.	1962	&& svHealth.E5b_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK
321	3/4 ✓ Branch 0 taken 1982 times. ✓ Branch 1 taken 18 times. ✓ Branch 2 taken 1962 times. ✗ Branch 3 not taken.	3982	&& svHealth.E1B_SignalHealthStatus == SvHealth::SignalHealthStatus::SignalOK;
322			}
323
324		✗	double GalileoEphemeris::calcSatellitePositionVariance() const
325			{
326		✗	return std::pow(signalAccuracy, 2);
327			}
328
329			} // namespace NAV
330