srp not working

Author: xanthospap

bin/CMakeLists.txt

@@ -3,7 +3,7 @@
 set(BIN_SOURCES
   test_deriv_costg.cpp
   eom.cpp
-  stm.cpp
+  #stm.cpp
 )
 
 # Process each source file and create an executable

bin/config.yaml

@@ -76,6 +76,6 @@ atmospheric-tide:
 
 satellite-attitude:
   satellite: ja3
+  cnes_sat_file: data/ja3mass.txt
   # data_file:
   data_file: data/qua_ja3.csv
-  cnes_sat_file: data/ja3mass.txt

bin/eom.cpp

@@ -265,28 +265,40 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
       ac += dso::iers2010_relativistic_acceleration(y0, rsun);
     }
 
-    /* Solar Radiation Pressure */
-    if (params->matt) {
+    if (params->mCr) {
       const double of =
-          dso::conical_occultation(y0.segment<3>(0), rsun.segment<3>(0));
+          /* dso::conical_occultation(y0.segment<3>(0), rsun.segment<3>(0)); */
+          dso::conical_occultation(y0.segment<3>(0), rsun.segment<3>(0), Eigen::Matrix<double,3,1>::Zero(), ::iers2010::Re) * 
+          dso::conical_occultation(y0.segment<3>(0), rsun.segment<3>(0), rmoon, 1737e3);
+      
       if (of > 0e0) {
-        /* get attitude */
-        if (params->matt->attitude_at(tt, *(params->mattdata))) {
-          fprintf(stderr, "[ERROR] Failed getting attitude!\n");
+        /* Solar Radiation Pressure */
+        if (params->matt) {
+          /* get attitude */
+          if (params->matt->attitude_at(tt, *(params->mattdata))) {
+            fprintf(stderr, "[ERROR] Failed getting attitude!\n");
 #ifdef USE_BOOST
-          throw std::runtime_error(
-              "ERROR. Failed computing derivative [attitude]\n");
+            throw std::runtime_error(
+                "ERROR. Failed computing derivative [attitude]\n");
 #else
-        return 200;
+          return 200;
 #endif
+          }
+          /* we may need (depending on satellite) the satellite-to-sun vector */
+          Eigen::Vector3d sat2sun = -y0.segment<3>(0)+rsun.segment<3>(0);
+          /* compute SRP acceleration */
+          ac += (params->mCr * of) *
+                dso::solar_radiation_pressure(
+                    params->msatmm->rotate_macromodel(
+                        params->mattdata->quaternions(),
+                        params->mattdata->angles(), &sat2sun),
+                    y0.segment<3>(0), rsun.segment<3>(0), params->msatmm->satellite_mass());
+        } else {
+          ac += (params->mCr * of) *
+                dso::solar_radiation_pressure(
+                    params->msatmm->srp_cannoball_area(), y0.segment<3>(0),
+                    rsun.segment<3>(0), params->msatmm->satellite_mass());
         }
-        /* compute SRP acceleration */
-        ac +=
-            (of / params->msatmm->satellite_mass()) *
-            dso::solar_radiation_pressure(params->msatmm->rotate_macromodel(
-                                              params->mattdata->quaternions(),
-                                              params->mattdata->angles()),
-                                          y0.segment<3>(0), rsun.segment<3>(0));
       }
     }
 
@@ -454,7 +466,7 @@ int main(int argc, char *argv[]) {
     }
 
     ++it;
-    if (it > 400)
+    if (it > 600)
       break;
   }
 

include/integration_parameters.hpp

@@ -40,6 +40,8 @@ namespace dso
     dso::attitude_details::MeasuredAttitudeData *mattdata{nullptr};
     /* satellite Macromodel */
     dso::SatelliteMacromodel *msatmm{nullptr};
+    /* dynamic parameter for SRP, i.e. Cr */
+    double mCr;
 
     /* third body gravity
      *                       | Bit Nr. | Default state

src/integration_parameters_from_config.cpp

@@ -230,30 +230,56 @@ dso::IntegrationParameters::from_config(const char *fn,
   }
   params.mdealias = dap;
 
+  /* get the satellite (mandatory) */
+  dso::SATELLITE sat(dso::SATELLITE::UNKNOWN);
+  if (is_nonempty_string(config["satellite-attitude"]["satellite"])) {
+    const auto s = config["satellite-attitude"]["satellite"].as<std::string>();
+    try {
+      sat = dso::translate_satid(s.c_str());
+    } catch (std::exception &e) {
+      fprintf(stderr, e.what());
+      std::string err_msg = "[ERROR] Failed to translate satellite name "
+                            "from "
+                            "config parameters " +
+                            s + " (traceback:" + std::string(__func__) + ")\n";
+      throw std::runtime_error(err_msg);
+    }
+  }
+
   /* Satellite attitude (if model satellite-attitude::data_file is non-empty) */
   dso::SatelliteAttitude *satat = nullptr;
-  dso::SatelliteMacromodel *satmm = nullptr;
   dso::attitude_details::MeasuredAttitudeData *mad = nullptr;
   if (is_nonempty_string(config["satellite-attitude"]["data_file"])) {
     const auto d = config["satellite-attitude"]["data_file"].as<std::string>();
-    const auto s = config["satellite-attitude"]["satellite"].as<std::string>();
     try {
-      const dso::SATELLITE sat = dso::translate_satid(s.c_str());
       satat = new dso::MeasuredAttitude(sat, d.c_str());
       mad = new dso::attitude_details::MeasuredAttitudeData(
           dso::measured_attitude_data_factory(sat));
-      satmm = new dso::SatelliteMacromodel(
-          dso::SatelliteMacromodel::createSatelliteMacromodel(sat));
     } catch (std::exception &e) {
       fprintf(stderr, e.what());
       std::string err_msg = "[ERROR] Failed to construct a SatelliteAttitude "
-                            "and/or SatelliteMacromodel instance from "
+                            "and/or MeasuredAttitudeData instance from "
                             "config parameters " +
-                            s + ", " + d +
+                            d +
                             " (traceback:" + std::string(__func__) + ")\n";
       throw std::runtime_error(err_msg);
     }
+  }
+  params.matt = satat;   // attitude
+  params.mattdata = mad; // attitude data (for retrieving attitude from stream)
 
+  /* create the macromodel (independent of attitude) */
+  dso::SatelliteMacromodel *satmm = nullptr;
+  {
+  try {
+    satmm = new dso::SatelliteMacromodel(
+        dso::SatelliteMacromodel::createSatelliteMacromodel(sat));
+  } catch (std::exception &e) {
+    fprintf(stderr, e.what());
+    std::string err_msg = "[ERROR] Failed to construct a SatelliteMacromodel"
+                          " instance from config parameters (traceback:" + std::string(__func__) + ")\n";
+    throw std::runtime_error(err_msg);
+  }
     if (is_nonempty_string(config["satellite-attitude"]["cnes_sat_file"])) {
       const auto c =
           config["satellite-attitude"]["cnes_sat_file"].as<std::string>();
@@ -265,9 +291,12 @@ dso::IntegrationParameters::from_config(const char *fn,
       }
     }
   }
-  params.matt = satat;
-  params.mattdata = mad;
-  params.msatmm = satmm;
+  params.msatmm = satmm; // macromodel (maybe used or not depending on attitude)
+  
+  /* Dynamic parameters */
+  {
+  params.mCr = config["dynamic-parameters"]["Cr"].as<double>();
+  }
 
   /* return */
   return params;