.. _program_listing_file_environment_body.cc:

Program Listing for File body.cc
================================

|exhale_lsh| :ref:`Return to documentation for file <file_environment_body.cc>` (``environment/body.cc``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "lupnt/environment/body.h"
   
   #include <filesystem>
   #include <fstream>
   #include <iostream>
   #include <sstream>
   #include <string>
   
   #include "lupnt/conversions/frame_converter.h"
   #include "lupnt/core/error.h"
   #include "lupnt/core/file.h"
   
   namespace lupnt {
     namespace {
       void PopulateBodyFromData(Body& body, const BodyData& data) {
         body.id = data.id;
         body.name = data.name;
         body.GM = data.GM;
         body.R = data.R;
         body.omega = data.omega;
         body.fixed_frame = data.fixed_frame;
         body.inertial_frame = data.inertial_frame;
         body.units = data.units;
       }
   
       BodyData ScaleBodyData(BodyData data, const UnitSystem& units) {
         data.GM = units.GravitationalParameter(data.GM.val());
         data.R = units.Length(data.R.val());
         data.omega = units.AngularVelocity(data.omega.val());
         data.units = units;
         return data;
       }
     }  // namespace
   
     Body Body::Moon(int n_max, int m_max, std::string gravity_file) {
       return Body::Moon(SI_UNITS, n_max, m_max, gravity_file);
     }
   
     Body Body::Moon(const UnitSystem& units, int n_max, int m_max, std::string gravity_file) {
       Body moon;
       PopulateBodyFromData(moon, GetBodyData(BodyId::MOON, units));
       moon.use_gravity_field = (n_max > 1 && m_max > 1);
       if (moon.use_gravity_field)
         moon.gravity_field = ReadHarmonicGravityField<Real>(gravity_file, n_max, m_max, true, units);
       return moon;
     }
   
     Body Body::Earth(int n_max, int m_max, std::string gravity_file) {
       return Body::Earth(SI_UNITS, n_max, m_max, gravity_file);
     }
   
     Body Body::Earth(const UnitSystem& units, int n_max, int m_max, std::string gravity_file) {
       Body earth;
       PopulateBodyFromData(earth, GetBodyData(BodyId::EARTH, units));
       earth.use_gravity_field = (n_max > 1 && m_max > 1);
       if (earth.use_gravity_field)
         earth.gravity_field = ReadHarmonicGravityField<Real>(gravity_file, n_max, m_max, true, units);
   
       return earth;
     }
   
     Body Body::Sun() { return Body::Sun(SI_UNITS); }
   
     Body Body::Sun(const UnitSystem& units) {
       Body sun;
       PopulateBodyFromData(sun, GetBodyData(BodyId::SUN, units));
       sun.use_gravity_field = false;
       return sun;
     }
   
     Body Body::Mars(int n_max, int m_max, std::string gravity_file) {
       return Body::Mars(SI_UNITS, n_max, m_max, gravity_file);
     }
   
     Body Body::Mars(const UnitSystem& units, int n_max, int m_max, std::string gravity_file) {
       Body mars;
       PopulateBodyFromData(mars, GetBodyData(BodyId::MARS, units));
       mars.use_gravity_field = (n_max > 1 && m_max > 1);
       if (mars.use_gravity_field)
         mars.gravity_field = ReadHarmonicGravityField<Real>(gravity_file, n_max, m_max, true, units);
       return mars;
     }
   
     Body Body::Venus(int n_max, int m_max, std::string gravity_file) {
       return Body::Venus(SI_UNITS, n_max, m_max, gravity_file);
     }
   
     Body Body::Venus(const UnitSystem& units, int n_max, int m_max, std::string gravity_file) {
       Body venus;
       PopulateBodyFromData(venus, GetBodyData(BodyId::VENUS, units));
       venus.use_gravity_field = (n_max > 1 && m_max > 1);
       if (venus.use_gravity_field)
         venus.gravity_field = ReadHarmonicGravityField<Real>(gravity_file, n_max, m_max, true, units);
       return venus;
     }
   
     Body Body::Jupiter() { return Body::Jupiter(SI_UNITS); }
   
     Body Body::Jupiter(const UnitSystem& units) {
       Body jupiter;
       PopulateBodyFromData(jupiter, GetBodyData(BodyId::JUPITER, units));
       jupiter.use_gravity_field = false;
       return jupiter;
     }
   
     Body Body::Saturn() { return Body::Saturn(SI_UNITS); }
   
     Body Body::Saturn(const UnitSystem& units) {
       Body saturn;
       PopulateBodyFromData(saturn, GetBodyData(BodyId::SATURN, units));
       saturn.use_gravity_field = false;
       return saturn;
     }
   
     Body Body::Uranus() { return Body::Uranus(SI_UNITS); }
   
     Body Body::Uranus(const UnitSystem& units) {
       Body uranus;
       PopulateBodyFromData(uranus, GetBodyData(BodyId::URANUS, units));
       uranus.use_gravity_field = false;
       return uranus;
     }
   
     Body Body::Neptune() { return Body::Neptune(SI_UNITS); }
   
     Body Body::Neptune(const UnitSystem& units) {
       Body neptune;
       PopulateBodyFromData(neptune, GetBodyData(BodyId::NEPTUNE, units));
       neptune.use_gravity_field = false;
       return neptune;
     }
   
     Body Body::CreateBody(BodyId body, int n, int m, std::string gravity_file) {
       return Body::CreateBody(body, SI_UNITS, n, m, gravity_file);
     }
   
     Body Body::CreateBody(BodyId body, const UnitSystem& units, int n, int m,
                           std::string gravity_file) {
       Logger::Debug(fmt::format("Creating body {}", magic_enum::enum_name(body).data()), "Body");
       if (body == BodyId::SUN) return Body::Sun(units);
       if (body == BodyId::EARTH) {
         if (!gravity_file.empty()) return Body::Earth(units, n, m, gravity_file);
         return Body::Earth(units, n, m);
       }
       if (body == BodyId::MOON) {
         if (!gravity_file.empty()) return Body::Moon(units, n, m, gravity_file);
         return Body::Moon(units, n, m);
       }
       if (body == BodyId::MARS) {
         if (!gravity_file.empty()) return Body::Mars(units, n, m, gravity_file);
         return Body::Mars(units, n, m);
       }
       if (body == BodyId::VENUS) {
         if (!gravity_file.empty()) return Body::Venus(units, n, m, gravity_file);
         return Body::Venus(units, n, m);
       }
       if (body == BodyId::JUPITER) return Body::Jupiter(units);
       if (body == BodyId::SATURN) return Body::Saturn(units);
       if (body == BodyId::URANUS) return Body::Uranus(units);
       if (body == BodyId::NEPTUNE) return Body::Neptune(units);
       LUPNT_CHECK(false, "Body not found", "Body");
     }
   
     BodyData GetBodyData(BodyId id) { return GetBodyData(id, SI_UNITS); }
   
     BodyData GetBodyData(BodyId id, const UnitSystem& units) {
       BodyData data;
       switch (id) {
         case BodyId::SUN:
           data = {.id = BodyId::SUN,
                   .name = "SUN",
                   .GM = GM_SUN,
                   .R = R_SUN,
                   .fixed_frame = Frame::ICRF,
                   .inertial_frame = Frame::ICRF,
                   .flattening = SUN_F,
                   .omega = OMEGA_SUN};
           break;
         case BodyId::MERCURY:
           data = {.id = BodyId::MERCURY,
                   .name = "MERCURY",
                   .GM = GM_MERCURY,
                   .R = R_MERCURY,
                   .fixed_frame = Frame::MERCURY_FIXED,
                   .inertial_frame = Frame::MERCURY_CI,
                   .flattening = MERCURY_F,
                   .omega = OMEGA_MERCURY};
           break;
         case BodyId::VENUS:
           data = {.id = BodyId::VENUS,
                   .name = "VENUS",
                   .GM = GM_VENUS,
                   .R = R_VENUS,
                   .fixed_frame = Frame::VENUS_FIXED,
                   .inertial_frame = Frame::VENUS_CI,
                   .flattening = VENUS_F,
                   .omega = OMEGA_VENUS};
           break;
         case BodyId::EARTH:
           data = {.id = BodyId::EARTH,
                   .name = "EARTH",
                   .GM = GM_EARTH,
                   .R = R_EARTH,
                   .fixed_frame = Frame::ITRF,
                   .inertial_frame = Frame::GCRF,
                   .flattening = WGS84_F,
                   .omega = OMEGA_EARTH};
           break;
         case BodyId::MOON:
           data = {.id = BodyId::MOON,
                   .name = "MOON",
                   .GM = GM_MOON,
                   .R = R_MOON,
                   .fixed_frame = Frame::MOON_PA,
                   .inertial_frame = Frame::MOON_CI,
                   .flattening = MOON_F,
                   .omega = OMEGA_MOON};
           break;
         case BodyId::MARS:
           data = {.id = BodyId::MARS,
                   .name = "MARS",
                   .GM = GM_MARS,
                   .R = R_MARS,
                   .fixed_frame = Frame::MARS_FIXED,
                   .inertial_frame = Frame::MARS_CI,
                   .flattening = MARS_F,
                   .omega = OMEGA_MARS};
           break;
         case BodyId::JUPITER:
           data = {.id = BodyId::JUPITER,
                   .name = "JUPITER",
                   .GM = GM_JUPITER,
                   .R = R_JUPITER,
                   .fixed_frame = Frame::JUPITER_FIXED,
                   .inertial_frame = Frame::JUPITER_CI,
                   .flattening = JUPITER_F,
                   .omega = OMEGA_JUPITER};
           break;
         case BodyId::SATURN:
           data = {.id = BodyId::SATURN,
                   .name = "SATURN",
                   .GM = GM_SATURN,
                   .R = R_SATURN,
                   .fixed_frame = Frame::SATURN_FIXED,
                   .inertial_frame = Frame::SATURN_CI,
                   .flattening = SATURN_F,
                   .omega = OMEGA_SATURN};
           break;
         case BodyId::URANUS:
           data = {.id = BodyId::URANUS,
                   .name = "URANUS",
                   .GM = GM_URANUS,
                   .R = R_URANUS,
                   .fixed_frame = Frame::URANUS_FIXED,
                   .inertial_frame = Frame::URANUS_CI,
                   .flattening = URANUS_F,
                   .omega = OMEGA_URANUS};
           break;
         case BodyId::NEPTUNE:
           data = {.id = BodyId::NEPTUNE,
                   .name = "NEPTUNE",
                   .GM = GM_NEPTUNE,
                   .R = R_NEPTUNE,
                   .fixed_frame = Frame::NEPTUNE_FIXED,
                   .inertial_frame = Frame::NEPTUNE_CI,
                   .flattening = NEPUTUNE_F,
                   .omega = OMEGA_NEPTUNE};
           break;
         default: throw std::runtime_error("Body not found");
       }
       return ScaleBodyData(data, units);
     }
   
     double GetBodyRadius(BodyId body) { return GetBodyRadius(body, SI_UNITS); }
   
     double GetBodyRadius(BodyId body, const UnitSystem& units) {
       BodyData data = GetBodyData(body, units);
       return data.R.val();
     }
   
     double GetBodyGM(BodyId body) { return GetBodyGM(body, SI_UNITS); }
   
     double GetBodyGM(BodyId body, const UnitSystem& units) {
       BodyData data = GetBodyData(body, units);
       return data.GM.val();
     }
   
     double GetBodyFlattening(BodyId body) {
       BodyData data = GetBodyData(body);
       return data.flattening.val();
     }
   
     double GetBodyOmega(BodyId body) { return GetBodyOmega(body, SI_UNITS); }
   
     double GetBodyOmega(BodyId body, const UnitSystem& units) {
       BodyData data = GetBodyData(body, units);
       return data.omega.val();
     }
   
     std::string GetBodyName(BodyId body) {
       BodyData data = GetBodyData(body);
       return data.name;
     }
   
     Frame GetInertialFrameName(BodyId body) {
       BodyData data = GetBodyData(body);
       return data.inertial_frame;
     }
   
     Frame GetBodyFixedFrameName(BodyId body) {
       BodyData data = GetBodyData(body);
       return data.fixed_frame;
     }
   
     Body CreateBody(BodyId body, int n, int m) { return CreateBody(body, SI_UNITS, n, m); }
   
     Body CreateBody(BodyId body, const UnitSystem& units, int n, int m) {
       return Body::CreateBody(body, units, n, m);
     }
   
     double kron(int i, int j) { return (i == j) ? 1 : 0; }
   
     double factprod(int n, int m) {
       double f = 1.0;
       for (int i = n - m + 1; i <= n + m; i++) {
         f /= i;
       }
       return f;
     }
   
     template <typename T> GravityField<T> ReadHarmonicGravityField(const std::string& filename, int n,
                                                                    int m, bool normalized) {
       return ReadHarmonicGravityField<T>(filename, n, m, normalized, SI_UNITS);
     }
   
     template <typename T> GravityField<T> ReadHarmonicGravityField(const std::string& filename, int n,
                                                                    int m, bool normalized,
                                                                    const UnitSystem& units) {
       GravityField<T> gravity_field;
       gravity_field.units = units;
       std::filesystem::path filepath = GetFilePath(filename);
       std::ifstream file = OpenFile<std::ifstream>(filepath);
   
       // Read header lines
       std::string line;
       while (std::getline(file, line)) {
         if (line.find("POTFIELD") != std::string::npos) {
           std::string potfield = line.substr(0, 8);
           int n_max_in = std::stoi(line.substr(8, 3));
           int m_max_in = std::stoi(line.substr(11, 3));
   
           std::istringstream iss(line.substr(14));
           double dummy1, GM, r, dummy2;
           iss >> dummy1 >> GM >> r >> dummy2;
           gravity_field.n_max = n_max_in;
           gravity_field.m_max = m_max_in;
           gravity_field.GM = units.GravitationalParameter(GM);
           gravity_field.R = units.Length(r);
           break;
         }
       }
       gravity_field.n = n;
       gravity_field.m = m;
   
       // Initialize Eigen matrices with the specified maxN and maxM
       if (n > gravity_field.n_max || m > gravity_field.m_max)
         throw std::runtime_error("Degree and order exceed maximum values");
   
       gravity_field.CS = Eigen::MatrixXd::Zero(n + 1, m + 1);
       gravity_field.CS(0, 0) = 1.0;  // C00 = 1.0
       // Read coefficient lines
       while (std::getline(file, line)) {
         if (line.find("RECOEF") != std::string::npos) {
           std::string recoef = line.substr(0, 8);
           int n_in = std::stoi(line.substr(8, 3));
           int m_in = std::stoi(line.substr(11, 3));
           std::istringstream iss(line.substr(14));
           double cnm, snm = 0.0;
           iss >> cnm;
           if (n_in >= n + 1) break;
           if (m_in >= m + 1) continue;
   
           if (m_in == 0) {
             double N = (normalized) ? sqrt(2 * n_in + 1) : 1.0;
             gravity_field.CS(n_in, m_in) = N * cnm;
           } else {
             double N = (normalized)
                            ? sqrt((2 - kron(0, m_in)) * (2 * n_in + 1) * factprod(n_in, m_in))
                            : 1.0;
             iss >> snm;
             double C = N * cnm;
             double S = N * snm;
             gravity_field.CS(n_in, m_in) = C;
             gravity_field.CS(m_in - 1, n_in) = S;
           }
         } else if (line.find("END") != std::string::npos) {
           break;
         }
       }
   
       file.close();
       return gravity_field;
     }
   
   }  // namespace lupnt