Add lunar gravity field disturbance

ExtLibraries/CMakeLists.txt

@@ -24,3 +24,4 @@ message("ExtLibraries install dir: ${EXT_LIB_DIR}")
 add_subdirectory(nrlmsise00)
 add_subdirectory(cspice)
 add_subdirectory(GeoPotential)
+add_subdirectory(lunar_gravity_field)

ExtLibraries/lunar_gravity_field/CMakeLists.txt

@@ -0,0 +1,30 @@
+project(lunar_gravity_field)
+
+cmake_minimum_required(VERSION 3.18)
+
+include(FetchContent)
+
+set(LUNAR_GRAVITY_FIELD_INSTALL_DIR ${EXT_LIB_DIR}/lunar_gravity_field)
+
+set(LUNAR_GRAVITY_FIELD_URL_BASE https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/shadr)
+set(LUNAR_GRAVITY_FIELD_FILE gggrx_1200a_sha.tab)
+set(LUNAR_GRAVITY_FIELD_TMP_DIR ${CMAKE_CURRENT_BINARY_DIR}/tmp)
+
+function(download_file_without_validate base_url file)
+  message("downloading ${file}")
+  FetchContent_Declare(
+    ${file}
+    SOURCE_DIR ${LUNAR_GRAVITY_FIELD_TMP_DIR}
+    URL ${base_url}/${file}
+    DOWNLOAD_NO_EXTRACT true
+  )
+  FetchContent_MakeAvailable(${file})
+endfunction()
+
+## download table
+download_file_without_validate(${LUNAR_GRAVITY_FIELD_URL_BASE} ${LUNAR_GRAVITY_FIELD_FILE})
+
+## install table
+install(FILES ${LUNAR_GRAVITY_FIELD_TMP_DIR}/${LUNAR_GRAVITY_FIELD_FILE}
+  DESTINATION ${LUNAR_GRAVITY_FIELD_INSTALL_DIR}
+)

data/sample/initialize_files/sample_disturbance.ini

@@ -5,6 +5,13 @@ logging = ENABLE
 degree = 4
 coefficients_file_path = EXT_LIB_DIR_FROM_EXE/GeoPotential/egm96_to360.ascii
 
+[LUNAR_GRAVITY_FIELD]
+// Enable only when the center object is defined as the Moon
+calculation = DISABLE
+logging = ENABLE
+degree = 10
+coefficients_file_path = EXT_LIB_DIR_FROM_EXE/LunarGravityField/gggrx_1200a_sha.tab
+
 
 [MAGNETIC_DISTURBANCE]
 // Enable only when the center object is defined as the Earth

src/disturbances/CMakeLists.txt

@@ -10,6 +10,7 @@ add_library(${PROJECT_NAME} STATIC
   solar_radiation_pressure_disturbance.cpp
   surface_force.cpp
   third_body_gravity.cpp
+  lunar_gravity_field.cpp
   initialize_disturbances.cpp
 )
 include(../../common.cmake)

src/disturbances/disturbances.cpp

@@ -11,6 +11,7 @@
 #include "geopotential.hpp"
 #include "gravity_gradient.hpp"
 #include "initialize_disturbances.hpp"
+#include "lunar_gravity_field.hpp"
 #include "magnetic_disturbance.hpp"
 #include "solar_radiation_pressure_disturbance.hpp"
 #include "third_body_gravity.hpp"
@@ -72,6 +73,11 @@ void Disturbances::InitializeInstances(const SimulationConfiguration* simulation
       new ThirdBodyGravity(InitThirdBodyGravity(initialize_file_name_, simulation_configuration->initialize_base_file_name_));
   disturbances_list_.push_back(third_body_gravity);
 
+  if (global_environment->GetCelestialInformation().GetCenterBodyName() == "MOON") {
+    LunarGravityField* lunar_gravity_field = new LunarGravityField(InitLunarGravityField(initialize_file_name_));
+    disturbances_list_.push_back(lunar_gravity_field);
+  }
+
   if (global_environment->GetCelestialInformation().GetCenterBodyName() != "EARTH") return;
   // Earth only disturbances (TODO: implement disturbances for other center bodies)
   AirDrag* air_dist =

src/disturbances/initialize_disturbances.cpp

@@ -88,6 +88,21 @@ Geopotential InitGeopotential(const std::string initialize_file_path) {
   return geopotential_disturbance;
 }
 
+LunarGravityField InitLunarGravityField(const std::string initialize_file_path) {
+  auto conf = IniAccess(initialize_file_path);
+  const char* section = "LUNAR_GRAVITY_FIELD";
+
+  const int degree = conf.ReadInt(section, "degree");
+  const std::string coefficients_file_path = conf.ReadString(section, "coefficients_file_path");
+
+  const bool is_calc_enable = conf.ReadEnable(section, CALC_LABEL);
+
+  LunarGravityField lunar_gravity_field(degree, coefficients_file_path, is_calc_enable);
+  lunar_gravity_field.is_log_enabled_ = conf.ReadEnable(section, LOG_LABEL);
+
+  return lunar_gravity_field;
+}
+
 ThirdBodyGravity InitThirdBodyGravity(const std::string initialize_file_path, const std::string ini_path_celes) {
   // Generate a list of bodies to be calculated in "CelesInfo"
   auto conf_celes = IniAccess(ini_path_celes);

src/disturbances/initialize_disturbances.hpp

@@ -9,6 +9,7 @@
 #include <disturbances/air_drag.hpp>
 #include <disturbances/geopotential.hpp>
 #include <disturbances/gravity_gradient.hpp>
+#include <disturbances/lunar_gravity_field.hpp>
 #include <disturbances/magnetic_disturbance.hpp>
 #include <disturbances/solar_radiation_pressure_disturbance.hpp>
 #include <disturbances/third_body_gravity.hpp>
@@ -62,6 +63,13 @@ MagneticDisturbance InitMagneticDisturbance(const std::string initialize_file_pa
  */
 Geopotential InitGeopotential(const std::string initialize_file_path);
 
+/**
+ * @fn InitLunarGravityField
+ * @brief Initialize LunarGravityField class with earth gravitational constant
+ * @param [in] initialize_file_path: Initialize file path
+ */
+LunarGravityField InitLunarGravityField(const std::string initialize_file_path);
+
 /**
  * @fn InitThirdBodyGravity
  * @brief Initialize ThirdBodyGravity class with earth gravitational constant

src/disturbances/lunar_gravity_field.cpp

@@ -0,0 +1,133 @@
+/**
+ * @file lunar_gravity_field.cpp
+ * @brief Class to calculate the high-order lunar gravity acceleration
+ */
+
+#include "lunar_gravity_field.hpp"
+
+#include <chrono>
+#include <cmath>
+#include <environment/global/physical_constants.hpp>
+#include <fstream>
+#include <iostream>
+
+#include "../library/logger/log_utility.hpp"
+#include "../library/utilities/macros.hpp"
+
+#define DEBUG_LUNAR_GRAVITY_FIELD
+
+LunarGravityField::LunarGravityField(const int degree, const std::string file_path, const bool is_calculation_enabled)
+    : Disturbance(is_calculation_enabled, false), degree_(degree) {
+  // Initialize
+  acceleration_mcmf_m_s2_ = libra::Vector<3>(0.0);
+  debug_pos_mcmf_m_ = libra::Vector<3>(0.0);
+  debug_pos_mcmf_m_[0] = 2000000;
+  debug_pos_mcmf_m_[1] = 2000000;
+  debug_pos_mcmf_m_[2] = 2000000;
+  // degree
+  if (degree_ > 1200) {
+    degree_ = 1200;
+    std::cout << "Inputted degree of LunarGravityField is too large for GRGM1200A model(limit is 1200)\n";
+    std::cout << "degree of LunarGravityField set as " << degree_ << "\n";
+  } else if (degree_ <= 1) {
+    degree_ = 0;
+  }
+  // coefficients
+  c_.assign(degree_ + 1, std::vector<double>(degree_ + 1, 0.0));
+  s_.assign(degree_ + 1, std::vector<double>(degree_ + 1, 0.0));
+  // For actual GRGM model, c_[0][0] should be 1.0
+  // In S2E, 0 degree term is inside the RK4 orbit calculation
+  c_[0][0] = 0.0;
+  if (degree_ >= 2) {
+    if (!ReadCoefficientsGrgm1200a(file_path)) {
+      degree_ = 0;
+      std::cout << "degree of LunarGravityField set as " << degree_ << "\n";
+    }
+  }
+  // Initialize GravityPotential
+  lunar_potential_ = GravityPotential(degree, c_, s_, gravity_constants_km3_s2_ * 1e9, reference_radius_km_ * 1e3);
+}
+
+bool LunarGravityField::ReadCoefficientsGrgm1200a(std::string file_name) {
+  std::ifstream coeff_file(file_name);
+  if (!coeff_file.is_open()) {
+    std::cerr << "File open error: LunarGravityField\n";
+    return false;
+  }
+
+  // Read header
+  std::string line, cell;
+  getline(coeff_file, cell, ',');
+  reference_radius_km_ = std::stod(cell);
+  getline(coeff_file, cell, ',');
+  gravity_constants_km3_s2_ = std::stod(cell);
+  // next line
+  getline(coeff_file, line);
+
+  size_t num_coeff = ((degree_ + 1) * (degree_ + 2) / 2) - 1;  // -1 for C00
+  for (size_t i = 0; i < num_coeff; i++) {
+    // degree
+    getline(coeff_file, line, ',');
+    int n = std::stoi(line);
+    getline(coeff_file, line, ',');
+    int m = std::stoi(line);
+    // coefficients
+    getline(coeff_file, line, ',');
+    double c_nm_norm = std::stod(line);
+    getline(coeff_file, line, ',');
+    double s_nm_norm = std::stod(line);
+    // next line
+    getline(coeff_file, line);
+
+    c_[n][m] = c_nm_norm;
+    s_[n][m] = s_nm_norm;
+  }
+  return true;
+}
+
+void LunarGravityField::Update(const LocalEnvironment &local_environment, const Dynamics &dynamics) {
+  libra::Vector<3> position_mcmf_m = dynamics.GetOrbit().GetPosition_ecef_m();
+#ifdef DEBUG_LUNAR_GRAVITY_FIELD
+  std::chrono::system_clock::time_point start, end;
+  start = std::chrono::system_clock::now();
+  position_mcmf_m = debug_pos_mcmf_m_;
+#endif
+
+  acceleration_mcmf_m_s2_ = lunar_potential_.CalcAcceleration_xcxf_m_s2(position_mcmf_m);
+#ifdef DEBUG_LUNAR_GRAVITY_FIELD
+  end = std::chrono::system_clock::now();
+  time_ms_ = static_cast<double>(std::chrono::duration_cast<std::chrono::microseconds>(end - start).count() / 1000.0);
+#else
+  UNUSED(time_ms_);
+#endif
+
+  libra::Matrix<3, 3> trans_eci2mcmf_ =
+      local_environment.GetCelestialInformation().GetGlobalInformation().GetEarthRotation().GetDcmJ2000ToXcxf();  // FIXME: use moon rotation
+  libra::Matrix<3, 3> trans_mcmf2eci = trans_eci2mcmf_.Transpose();
+  acceleration_i_m_s2_ = trans_mcmf2eci * acceleration_mcmf_m_s2_;
+}
+
+std::string LunarGravityField::GetLogHeader() const {
+  std::string str_tmp = "";
+
+#ifdef DEBUG_LUNAR_GRAVITY_FIELD
+  str_tmp += WriteVector("lunar_gravity_calculation_position", "mcmf", "m", 3);
+  str_tmp += WriteScalar("lunar_gravity_calculation_time", "ms");
+#endif
+  str_tmp += WriteVector("lunar_gravity_acceleration", "mcmf", "m/s2", 3);
+
+  return str_tmp;
+}
+
+std::string LunarGravityField::GetLogValue() const {
+  std::string str_tmp = "";
+
+#ifdef DEBUG_LUNAR_GRAVITY_FIELD
+  str_tmp += WriteVector(debug_pos_mcmf_m_, 15);
+  str_tmp += WriteScalar(time_ms_);
+#endif
+
+  str_tmp += WriteVector(acceleration_mcmf_m_s2_, 15);
+
+  return str_tmp;
+}

src/disturbances/lunar_gravity_field.hpp

@@ -0,0 +1,85 @@
+/**
+ * @file lunar_gravity_field.hpp
+ * @brief Class to calculate the high-order lunar gravity acceleration
+ */
+
+#ifndef S2E_DISTURBANCES_LUNAR_GRAVITY_FIELD_HPP_
+#define S2E_DISTURBANCES_LUNAR_GRAVITY_FIELD_HPP_
+
+#include <string>
+
+#include "../library/gravity/gravity_potential.hpp"
+#include "../library/math/vector.hpp"
+#include "disturbance.hpp"
+/**
+ * @class LunarGravityField
+ * @brief Class to calculate the high-order earth gravity acceleration
+ */
+class LunarGravityField : public Disturbance {
+ public:
+  /**
+   * @fn Geopotential
+   * @brief Constructor
+   * @param [in] degree: Maximum degree setting to calculate the geo-potential
+   * @param [in] file_path: EGM96 coefficients file path
+   * @param [in] is_calculation_enabled: Calculation flag
+   */
+  LunarGravityField(const int degree, const std::string file_path, const bool is_calculation_enabled = true);
+
+  /**
+   * @fn LunarGravityField
+   * @brief Copy Constructor
+   */
+  LunarGravityField(const LunarGravityField &obj) : Disturbance(obj) {
+    lunar_potential_ = obj.lunar_potential_;
+    reference_radius_km_ = obj.reference_radius_km_;
+    gravity_constants_km3_s2_ = obj.gravity_constants_km3_s2_;
+    degree_ = obj.degree_;
+    c_ = obj.c_;
+    s_ = obj.s_;
+  }
+
+  ~LunarGravityField() {}
+
+  /**
+   * @fn Update
+   * @brief Override Updates function of SimpleDisturbance
+   * @param [in] local_environment: Local environment information
+   * @param [in] dynamics: Dynamics information
+   */
+  virtual void Update(const LocalEnvironment &local_environment, const Dynamics &dynamics);
+
+  // Override ILoggable
+  /**
+   * @fn GetLogHeader
+   * @brief Override GetLogHeader function of ILoggable
+   */
+  virtual std::string GetLogHeader() const;
+  /**
+   * @fn GetLogValue
+   * @brief Override GetLogValue function of ILoggable
+   */
+  virtual std::string GetLogValue() const;
+
+ private:
+  GravityPotential lunar_potential_;
+  double reference_radius_km_;
+  double gravity_constants_km3_s2_;
+  size_t degree_;                       //!< Maximum degree setting to calculate the geo-potential
+  std::vector<std::vector<double>> c_;  //!< Cosine coefficients
+  std::vector<std::vector<double>> s_;  //!< Sine coefficients
+  Vector<3> acceleration_mcmf_m_s2_;    //!< Calculated acceleration in the MCMF(Moon Centered Moon Fixed) frame [m/s2]
+
+  // debug
+  libra::Vector<3> debug_pos_mcmf_m_;  //!< Spacecraft position in MCMF frame [m]
+  double time_ms_ = 0.0;               //!< Calculation time [ms]
+
+  /**
+   * @fn ReadCoefficientsGrgm1200a
+   * @brief Read the lunar gravity field coefficients for the GRGM1200A model
+   * @param [in] file_name: Coefficient file name
+   */
+  bool ReadCoefficientsGrgm1200a(std::string file_name);
+};
+
+#endif  // S2E_DISTURBANCES_LUNAR_GRAVITY_FIELD_HPP_