Add relative orbit dynamics with Carter's STM

settings/sample_satellite/satellite_sub.ini

@@ -104,7 +104,7 @@ relative_orbit_update_method = 1
 relative_dynamics_model_type = 0
 // STM Relative Dynamics model type (only valid for STM update)
 // 0: HCW, 1: Melton, 2: SS, 3: Sabatini, 4: Carter, 5: Yamanaka-Ankersen
-stm_model_type = 5
+stm_model_type = 4
 // Initial satellite position relative to the reference satellite in LVLH frame[m]
 // * The coordinate system is defined at [PLANET_SELECTION] in SampleSimBase.ini
 initial_relative_position_lvlh_m(0) = 0.0

src/dynamics/orbit/relative_orbit.cpp

@@ -99,6 +99,9 @@ void RelativeOrbit::InitializeStmMatrix(orbit::StmModel stm_model_type, const Or
   double f_ref_rad = phi_rad - arg_perigee_rad;
 
   switch (stm_model_type) {
+    case orbit::StmModel::kCarter:
+      relative_orbit_carter_.CalculateInitialInverseMatrix(gravity_constant_m3_s2, f_ref_rad, &reference_oe);
+      break;
     case orbit::StmModel::kYamakawaAnkersen:
       relative_orbit_yamanaka_ankersen_.CalculateInitialInverseMatrix(f_ref_rad, &reference_oe);
       break;
@@ -143,7 +146,7 @@ void RelativeOrbit::CalculateStm(orbit::StmModel stm_model_type, const Orbit* re
       break;
     }
     case orbit::StmModel::kCarter: {
-      stm_ = orbit::CalcCarterStm(reference_sat_orbit_radius, gravity_constant_m3_s2, f_ref_rad, &reference_oe);
+      stm_ = relative_orbit_carter_.CalculateSTM(gravity_constant_m3_s2, elapsed_sec, f_ref_rad, &reference_oe);
       break;
     }
     case orbit::StmModel::kYamakawaAnkersen: {

src/dynamics/orbit/relative_orbit.hpp

@@ -9,6 +9,7 @@
 #include <math_physics/math/ordinary_differential_equation.hpp>
 #include <math_physics/orbit/relative_orbit_models.hpp>
 #include <math_physics/orbit/relative_orbit_yamanaka_ankersen.hpp>
+#include <math_physics/orbit/relative_orbit_carter.hpp>
 #include <simulation/multiple_spacecraft/relative_information.hpp>
 #include <string>
 
@@ -87,6 +88,7 @@ class RelativeOrbit : public Orbit, public math::OrdinaryDifferentialEquation<6>
   orbit::StmModel stm_model_type_;                                          //!< State Transition Matrix model type
   RelativeInformation* relative_information_;                               //!< Relative information
   orbit::RelativeOrbitYamanakaAnkersen relative_orbit_yamanaka_ankersen_;  //!< Relative Orbit Calcilater with Yamanaka-Ankersen's STM
+  orbit::RelativeOrbitCarter relative_orbit_carter_;                       //!< Relative Orbit Calcilater with Carter's STM
 
   /**
    * @fn InitializeState

src/math_physics/CMakeLists.txt

@@ -33,6 +33,7 @@ add_library(${PROJECT_NAME} STATIC
   orbit/kepler_orbit.cpp
   orbit/relative_orbit_models.cpp
   orbit/relative_orbit_yamanaka_ankersen.cpp
+  orbit/relative_orbit_carter.cpp
   orbit/interpolation_orbit.cpp
   orbit/sgp4/sgp4ext.cpp
   orbit/sgp4/sgp4io.cpp

src/math_physics/orbit/relative_orbit_carter.cpp

@@ -0,0 +1,143 @@
+/**
+ * @file relative_orbit_carter.cpp
+ * @brief Functions to calculate Carter's STM for relative orbit
+ */
+#include "relative_orbit_carter.hpp"
+
+#include <environment/global/physical_constants.hpp>
+
+#include "./relative_orbit_models.hpp"
+#include "./sgp4/sgp4unit.h"  // for getgravconst()
+
+namespace orbit {
+
+RelativeOrbitCarter::RelativeOrbitCarter() {}
+
+RelativeOrbitCarter::~RelativeOrbitCarter() {}
+
+void RelativeOrbitCarter::CalculateInitialInverseMatrix(double gravity_constant_m3_s2, double f_ref_rad, OrbitalElements* reference_oe) {
+  double e = reference_oe->GetEccentricity();
+  double a = reference_oe->GetSemiMajorAxis_m();
+  double h = pow(a * (1.0 - pow(e, 2)) * gravity_constant_m3_s2, 0.5);  // angular momentum
+
+  double E_rad = 2.0 * atan(sqrt((1.0 - e) / (1.0 + e)) * tan(f_ref_rad / 2.0));
+  double k = e * cos(f_ref_rad) + 1.0;
+  // double K1 = pow(1 - e * e, -2.5) * (-1.5 * e * E_rad + (1 + e * e) * sin(E_rad) - e * sin(2. * E_rad) / 4.);
+  double K2 = pow(1.0 - pow(e, 2.0), -2.5) * (0.5 * E_rad - 0.25 * sin(2.0 * E_rad) - e * pow(sin(E_rad), 3.0) / 3.0);
+  double phi1 = sin(f_ref_rad) * k;
+  double phi2 = 2.0 * e * phi1 * (sin(f_ref_rad) / pow(k, 3) - 3.0 * e * K2) - cos(f_ref_rad) / k;
+  double phi3 = 6.0 * e * phi1 * K2 - 2.0 * pow(sin(f_ref_rad), 2.0) / pow(k, 2.0) - pow(cos(f_ref_rad), 2.0) / k - pow(cos(f_ref_rad), 2.0);
+  double phi1_prime = cos(f_ref_rad) * k - e * pow(sin(f_ref_rad), 2.0);
+  double phi2_prime = -6.0 * pow(e, 2.0) * phi1_prime * K2 +
+                      2.0 * e * sin(f_ref_rad) * (2.0 * cos(f_ref_rad) - 3.0 * e * pow(sin(f_ref_rad), 2.0) + 2.0 * e) / pow(k, 3.0) +
+                      sin(f_ref_rad) / pow(k, 2.0);
+  double phi3_prime = 6.0 * e * phi1_prime * K2 - (6.0 * e * pow(sin(f_ref_rad), 3.0) - 4.0 * sin(f_ref_rad) * (e + cos(f_ref_rad))) / pow(k, 3.0) +
+                      0.5 * sin(2 * f_ref_rad) * (2.0 + (2.0 + e * cos(f_ref_rad)) * pow(k, 2.0));
+  double S1 = -cos(f_ref_rad) * (1.0 + 0.5 * e * cos(f_ref_rad));
+  double S2 = 3.0 * e * pow(k, 2.0) * K2 - sin(f_ref_rad) / k;
+  double S3 = -6.0 * pow(k, 2.0) * K2 - (2.0 + k) / 2.0 / k * sin(2.0 * f_ref_rad);
+
+  initial_inverse_matrix_[0][0] = 4.0 * S2 + phi2_prime;
+  initial_inverse_matrix_[0][1] = 0.0;
+  initial_inverse_matrix_[0][2] = 0.0;
+  initial_inverse_matrix_[0][3] = -phi2;
+  initial_inverse_matrix_[0][4] = 2.0 * S2;
+  initial_inverse_matrix_[0][5] = 0.0;
+  initial_inverse_matrix_[1][0] = -2.0;
+  initial_inverse_matrix_[1][1] = 0.0;
+  initial_inverse_matrix_[1][2] = 0.0;
+  initial_inverse_matrix_[1][3] = 0.0;
+  initial_inverse_matrix_[1][4] = -1.0;
+  initial_inverse_matrix_[1][5] = 0.0;
+  initial_inverse_matrix_[2][0] = 0.0;
+  initial_inverse_matrix_[2][1] = 0.0;
+  initial_inverse_matrix_[2][2] = cos(f_ref_rad);
+  initial_inverse_matrix_[2][3] = 0.0;
+  initial_inverse_matrix_[2][4] = 0.0;
+  initial_inverse_matrix_[2][5] = -sin(f_ref_rad);
+  initial_inverse_matrix_[3][0] = -(4.0 * S1 + phi1_prime);
+  initial_inverse_matrix_[3][1] = 0.0;
+  initial_inverse_matrix_[3][2] = 0.0;
+  initial_inverse_matrix_[3][3] = phi1;
+  initial_inverse_matrix_[3][4] = -2.0 * S1;
+  initial_inverse_matrix_[3][5] = 0.0;
+  initial_inverse_matrix_[4][0] = 2.0 * S3 + phi3_prime;
+  initial_inverse_matrix_[4][1] = -1.0;
+  initial_inverse_matrix_[4][2] = 0.0;
+  initial_inverse_matrix_[4][3] = -phi3;
+  initial_inverse_matrix_[4][4] = S3;
+  initial_inverse_matrix_[4][5] = 0.0;
+  initial_inverse_matrix_[5][0] = 0.0;
+  initial_inverse_matrix_[5][1] = 0.0;
+  initial_inverse_matrix_[5][2] = sin(f_ref_rad);
+  initial_inverse_matrix_[5][3] = 0.0;
+  initial_inverse_matrix_[5][4] = 0.0;
+  initial_inverse_matrix_[5][5] = cos(f_ref_rad);
+
+  initial_inverse_matrix_ =
+      initial_inverse_matrix_ * orbit::CalcStateTransformationMatrixLvlhToTschaunerHampel(gravity_constant_m3_s2, e, h, f_ref_rad);
+}
+
+math::Matrix<6, 6> RelativeOrbitCarter::CalculateSTM(double orbit_radius_m, double gravity_constant_m3_s2, double f_ref_rad,
+                                                     OrbitalElements* reference_oe) {
+  math::Matrix<6, 6> stm;
+  double e = reference_oe->GetEccentricity();
+  double a = reference_oe->GetSemiMajorAxis_m();
+  double h = pow(a * (1 - pow(e, 2)) * gravity_constant_m3_s2, 0.5);  // angular momentum
+  double E_rad = 2.0 * atan(sqrt((1.0 - e) / (1.0 + e)) * tan(f_ref_rad / 2.0));
+  double k = e * cos(f_ref_rad) + 1.0;
+  // double K1 = pow(1 - e * e, -2.5) * (-1.5 * e * E_rad + (1 + e * e) * sin(E_rad) - e * sin(2. * E_rad) / 4.);
+  double K2 = pow(1.0 - pow(e, 2.0), -2.5) * (0.5 * E_rad - 0.25 * sin(2.0 * E_rad) - e * pow(sin(E_rad), 3.0) / 3.0);
+  double phi1 = sin(f_ref_rad) * k;
+  double phi2 = 2.0 * e * phi1 * (sin(f_ref_rad) / pow(k, 3) - 3.0 * e * K2) - cos(f_ref_rad) / k;
+  double phi3 = 6.0 * e * phi1 * K2 - 2.0 * pow(sin(f_ref_rad), 2.0) / pow(k, 2.0) - pow(cos(f_ref_rad), 2.0) / k - pow(cos(f_ref_rad), 2.0);
+  double phi1_prime = cos(f_ref_rad) * k - e * pow(sin(f_ref_rad), 2.0);
+  double phi2_prime = -6.0 * pow(e, 2.0) * phi1_prime * K2 +
+                      2.0 * e * sin(f_ref_rad) * (2.0 * cos(f_ref_rad) - 3.0 * e * pow(sin(f_ref_rad), 2.0) + 2.0 * e) / pow(k, 3.0) +
+                      sin(f_ref_rad) / pow(k, 2.0);
+  double phi3_prime = 6.0 * e * phi1_prime * K2 - (6.0 * e * pow(sin(f_ref_rad), 3.0) - 4.0 * sin(f_ref_rad) * (e + cos(f_ref_rad))) / pow(k, 3.0) +
+                      0.5 * sin(2 * f_ref_rad) * (2.0 + (2.0 + e * cos(f_ref_rad)) * pow(k, 2.0));
+  double S1 = -cos(f_ref_rad) * (1.0 + 0.5 * e * cos(f_ref_rad));
+  double S2 = 3.0 * e * pow(k, 2.0) * K2 - sin(f_ref_rad) / k;
+  double S3 = -6.0 * pow(k, 2.0) * K2 - (2.0 + k) / 2.0 / k * sin(2.0 * f_ref_rad);
+
+  stm[0][0] = phi1;
+  stm[0][1] = phi3;
+  stm[0][2] = 0.0;
+  stm[0][3] = phi2;
+  stm[0][4] = 0.0;
+  stm[0][5] = 0.0;
+  stm[1][0] = -2 * S1;
+  stm[1][1] = -S3;
+  stm[1][2] = 0.0;
+  stm[1][3] = -2 * S2;
+  stm[1][4] = -1;
+  stm[1][5] = 0.0;
+  stm[2][0] = 0.0;
+  stm[2][1] = 0.0;
+  stm[2][2] = cos(f_ref_rad);
+  stm[2][3] = 0.0;
+  stm[2][4] = 0.0;
+  stm[2][5] = sin(f_ref_rad);
+  stm[3][0] = phi1_prime;
+  stm[3][1] = phi3_prime;
+  stm[3][2] = 0.0;
+  stm[3][3] = phi2_prime;
+  stm[3][4] = 0.0;
+  stm[3][5] = 0.0;
+  stm[4][0] = -2 * phi1;
+  stm[4][1] = -(2 * phi3 + 1);
+  stm[4][2] = 0.0;
+  stm[4][3] = -2 * phi2;
+  stm[4][4] = 0.0;
+  stm[4][5] = 0.0;
+  stm[5][0] = 0.0;
+  stm[5][1] = 0.0;
+  stm[5][2] = -sin(f_ref_rad);
+  stm[5][3] = 0.0;
+  stm[5][4] = 0.0;
+  stm[5][5] = cos(f_ref_rad);
+  return orbit::CalcStateTransformationMatrixTschaunerHampelToLvlh(gravity_constant_m3_s2, e, h, f_ref_rad) * stm * initial_inverse_matrix_;
+}
+
+}  // namespace orbit

src/math_physics/orbit/relative_orbit_carter.hpp

@@ -0,0 +1,62 @@
+/**
+ * @file relative_orbit_carter.hpp
+ * @brief Functions to calculate Carter's STM for relative orbit
+ */
+
+#ifndef S2E_LIBRARY_ORBIT_RELATIVE_ORBIT_CARTER_HPP_
+#define S2E_LIBRARY_ORBIT_RELATIVE_ORBIT_CARTER_HPP_
+
+#include "../math/matrix.hpp"
+#include "./orbital_elements.hpp"
+
+namespace orbit {
+
+/**
+ * @class RelativeOrbitCarter
+ * @brief Class to calculate Yamanaka-Ankersen relative orbital STM
+ */
+class RelativeOrbitCarter {
+ public:
+  /**
+   * @fn RelativeOrbitCarter
+   * @brief Default Constructor
+   */
+  RelativeOrbitCarter();
+  /**
+   * @fn ~RelativeOrbitCarter
+   * @brief Destructor
+   */
+  ~RelativeOrbitCarter();
+
+  /**
+   * @fn CalculateInitialInverseMatrix
+   * @brief Calculate position and velocity with Kepler orbit propagation
+   * @param [in] gravity_constant_m3_s2: Gravity constant of the center body [m3/s2]
+   * @param [in] f_ref_rad: True anomaly of the reference satellite [rad]
+   * @param [in] reference_oe: Orbital elements of reference satellite
+   */
+  void CalculateInitialInverseMatrix(double gravity_constant_m3_s2, double f_ref_rad, OrbitalElements* reference_oe);
+
+  /**
+   * @fn CalculateSTM
+   * @brief Calculate position and velocity with Kepler orbit propagation
+   * @param [in] orbit_radius_m: Orbit radius [m]
+   * @param [in] gravity_constant_m3_s2: Gravity constant of the center body [m3/s2]
+   * @param [in] f_ref_rad: True anomaly of the reference satellite [rad]
+   * @param [in] reference_oe: Orbital elements of reference satellite
+   */
+  math::Matrix<6, 6> CalculateSTM(double orbit_radius_m, double gravity_constant_m3_s2, double f_ref_rad, OrbitalElements* reference_oe);
+
+  /**
+   * @fn GetInitialInverseMatrix
+   * @brief Return initial inverse matrix
+   */
+  inline const math::Matrix<6, 6> GetInitialInverseMatrix() const { return initial_inverse_matrix_; }
+
+ private:
+  math::Matrix<6, 6> initial_inverse_matrix_{0.0};  //!< Gravity constant of the center body [m3/s2]
+};
+
+}  // namespace orbit
+
+#endif  // S2E_LIBRARY_ORBIT_RELATIVE_ORBIT_MODEL_HPP_

src/math_physics/orbit/relative_orbit_models.cpp

@@ -129,10 +129,90 @@ math::Matrix<6, 6> CalcSabatiniStm(double orbit_radius_m, double gravity_constan
   return stm;
 }
 
-math::Matrix<6, 6> CalcCarterStm(double orbit_radius_m, double gravity_constant_m3_s2, double f_ref_rad, OrbitalElements* reference_oe) {
-  math::Matrix<6, 6> stm;
-  // ここでstmを計算してください
-  return stm;
+math::Matrix<6, 6> CalcStateTransformationMatrixLvlhToTschaunerHampel(const double gravity_constant_m3_s2, const double eccentricity,
+                                                                      const double angular_momentum_kg_m2_s, const double true_anomaly_rad) {
+  math::Matrix<6, 6> transition_matrix;
+  transition_matrix[0][0] = 1.0 + eccentricity * cos(true_anomaly_rad);
+  transition_matrix[0][1] = 0.0;
+  transition_matrix[0][2] = 0.0;
+  transition_matrix[0][3] = 0.0;
+  transition_matrix[0][4] = 0.0;
+  transition_matrix[0][5] = 0.0;
+  transition_matrix[1][0] = 0.0;
+  transition_matrix[1][1] = 1.0 + eccentricity * cos(true_anomaly_rad);
+  transition_matrix[1][2] = 0.0;
+  transition_matrix[1][3] = 0.0;
+  transition_matrix[1][4] = 0.0;
+  transition_matrix[1][5] = 0.0;
+  transition_matrix[2][0] = 0.0;
+  transition_matrix[2][1] = 0.0;
+  transition_matrix[2][2] = 1.0 + eccentricity * cos(true_anomaly_rad);
+  transition_matrix[2][3] = 0.0;
+  transition_matrix[2][4] = 0.0;
+  transition_matrix[2][5] = 0.0;
+  transition_matrix[3][0] = -eccentricity * sin(true_anomaly_rad);
+  transition_matrix[3][1] = 0.0;
+  transition_matrix[3][2] = 0.0;
+  transition_matrix[3][3] = pow(angular_momentum_kg_m2_s, 3.0) / pow(gravity_constant_m3_s2, 2.0) / (1.0 + eccentricity * cos(true_anomaly_rad));
+  transition_matrix[3][4] = 0.0;
+  transition_matrix[3][5] = 0.0;
+  transition_matrix[4][0] = 0.0;
+  transition_matrix[4][1] = -eccentricity * sin(true_anomaly_rad);
+  transition_matrix[4][2] = 0.0;
+  transition_matrix[4][3] = 0.0;
+  transition_matrix[4][4] = pow(angular_momentum_kg_m2_s, 3.0) / pow(gravity_constant_m3_s2, 2.0) / (1.0 + eccentricity * cos(true_anomaly_rad));
+  transition_matrix[4][5] = 0.0;
+  transition_matrix[5][0] = 0.0;
+  transition_matrix[5][1] = 0.0;
+  transition_matrix[5][2] = -eccentricity * sin(true_anomaly_rad);
+  transition_matrix[5][3] = 0.0;
+  transition_matrix[5][4] = 0.0;
+  transition_matrix[5][5] = pow(angular_momentum_kg_m2_s, 3.0) / pow(gravity_constant_m3_s2, 2.0) / (1.0 + eccentricity * cos(true_anomaly_rad));
+
+  return transition_matrix;
+}
+
+math::Matrix<6, 6> CalcStateTransformationMatrixTschaunerHampelToLvlh(const double gravity_constant_m3_s2, const double eccentricity,
+                                                                      const double angular_momentum_kg_m2_s, const double true_anomaly_rad) {
+  math::Matrix<6, 6> transition_matrix;
+  transition_matrix[0][0] = 1.0 / (1.0 + eccentricity * cos(true_anomaly_rad));
+  transition_matrix[0][1] = 0.0;
+  transition_matrix[0][2] = 0.0;
+  transition_matrix[0][3] = 0.0;
+  transition_matrix[0][4] = 0.0;
+  transition_matrix[0][5] = 0.0;
+  transition_matrix[1][0] = 0.0;
+  transition_matrix[1][1] = 1.0 / (1.0 + eccentricity * cos(true_anomaly_rad));
+  transition_matrix[1][2] = 0.0;
+  transition_matrix[1][3] = 0.0;
+  transition_matrix[1][4] = 0.0;
+  transition_matrix[1][5] = 0.0;
+  transition_matrix[2][0] = 0.0;
+  transition_matrix[2][1] = 0.0;
+  transition_matrix[2][2] = 1.0 / (1.0 + eccentricity * cos(true_anomaly_rad));
+  transition_matrix[2][3] = 0.0;
+  transition_matrix[2][4] = 0.0;
+  transition_matrix[2][5] = 0.0;
+  transition_matrix[3][0] = pow(gravity_constant_m3_s2, 2.0) * eccentricity * sin(true_anomaly_rad) / pow(angular_momentum_kg_m2_s, 3.0);
+  transition_matrix[3][1] = 0.0;
+  transition_matrix[3][2] = 0.0;
+  transition_matrix[3][3] = pow(gravity_constant_m3_s2, 2.0) * (1.0 + eccentricity * cos(true_anomaly_rad)) / pow(angular_momentum_kg_m2_s, 3.0);
+  transition_matrix[3][4] = 0.0;
+  transition_matrix[3][5] = 0.0;
+  transition_matrix[4][0] = 0.0;
+  transition_matrix[4][1] = pow(gravity_constant_m3_s2, 2.0) * eccentricity * sin(true_anomaly_rad) / pow(angular_momentum_kg_m2_s, 3.0);
+  transition_matrix[4][2] = 0.0;
+  transition_matrix[4][3] = 0.0;
+  transition_matrix[4][4] = pow(gravity_constant_m3_s2, 2.0) * (1.0 + eccentricity * cos(true_anomaly_rad)) / pow(angular_momentum_kg_m2_s, 3.0);
+  transition_matrix[4][5] = 0.0;
+  transition_matrix[5][0] = 0.0;
+  transition_matrix[5][1] = 0.0;
+  transition_matrix[5][2] = pow(gravity_constant_m3_s2, 2.0) * eccentricity * sin(true_anomaly_rad) / pow(angular_momentum_kg_m2_s, 3.0);
+  transition_matrix[5][3] = 0.0;
+  transition_matrix[5][4] = 0.0;
+  transition_matrix[5][5] = pow(gravity_constant_m3_s2, 2.0) * (1.0 + eccentricity * cos(true_anomaly_rad)) / pow(angular_momentum_kg_m2_s, 3.0);
+
+  return transition_matrix;
 }
 
 }  // namespace orbit

src/math_physics/orbit/relative_orbit_models.hpp

@@ -90,15 +90,28 @@ math::Vector<6> CalcSsCorrectionTerm(double orbit_radius_m, double gravity_const
 math::Matrix<6, 6> CalcSabatiniStm(double orbit_radius_m, double gravity_constant_m3_s2, double elapsed_time_s, OrbitalElements* reference_oe);
 
 /**
- * @fn CalcCarterStm
- * @brief Calculate Carter State Transition Matrix
- * @param [in] orbit_radius_m: Orbit radius [m]
+ * @fn CalcStateTransformationMatrixLvlhToTschaunerHampel
+ * @brief Calculate state tranformation matrix from the state variables in the LVLH frame to the state variables for Tschauner-Hampel equation.
  * @param [in] gravity_constant_m3_s2: Gravity constant of the center body [m3/s2]
- * @param [in] f_ref_rad: True anomaly of the reference satellite [rad]
- * @param [in] reference_oe: Orbital elements of reference satellite
- * @return State Transition Matrix
+ * @param [in] eccentricity: Eccentricity []
+ * @param [in] angular_momentum_kg_m2_s: Angular momentum of the spacecraft [kg*m2/s]
+ * @param [in] true_anomaly_rad: True anomaly of the spacecraft [rad]
+ * @param
+ */
+math::Matrix<6, 6> CalcStateTransformationMatrixLvlhToTschaunerHampel(const double gravity_constant_m3_s2, const double eccentricity,
+                                                                      const double angular_momentum_kg_m2_s, const double true_anomaly_rad);
+
+/**
+ * @fn CalcStateTransformationMatrixTschaunerHampelToLvlh
+ * @brief Calculate state tranformation matrix from the state variables for Tschauner-Hampel equation to the state variables in the LVLH frame.
+ * @param [in] gravity_constant_m3_s2: Gravity constant of the center body [m3/s2]
+ * @param [in] eccentricity: Eccentricity []
+ * @param [in] angular_momentum_kg_m2_s: Angular momentum of the spacecraft [kg*m2/s]
+ * @param [in] true_anomaly_rad: True anomaly of the spacecraft [rad]
+ * @param
  */
-math::Matrix<6, 6> CalcCarterStm(double orbit_radius_m, double gravity_constant_m3_s2, double f_ref_rad, OrbitalElements* reference_oe);
+math::Matrix<6, 6> CalcStateTransformationMatrixTschaunerHampelToLvlh(const double gravity_constant_m3_s2, const double eccentricity,
+                                                                      const double angular_momentum_kg_m2_s, const double true_anomaly_rad);
 
 }  // namespace orbit