more

Author: xanthospap

CMakeLists.txt

@@ -11,10 +11,11 @@ project(
 find_package(datetime REQUIRED)
 find_package(iers REQUIRED)
 find_package(geodesy REQUIRED)
+find_package(sysnsats REQUIRED)
 find_package(yaml-cpp REQUIRED)
 
 # Pass the library dependencies to subdirectories
-set(PROJECT_DEPENDENCIES iers geodesy datetime)
+set(PROJECT_DEPENDENCIES sysnsats iers geodesy datetime)
 
 # Define an option for building binaries (defaults to ON)
 option(BUILD_PROGS "Enable building of binaries" ON)

bin/CMakeLists.txt

@@ -1,7 +1,6 @@
 # bin/CMakeLists.txt
 
 set(BIN_SOURCES
-  equations_of_motion.cpp
   test_deriv_costg.cpp
   eom.cpp
   eomboost.cpp

bin/config.yaml

@@ -42,10 +42,9 @@ ocean-pole-tide:
 
 dealiasing:
   # Leave model name empty to not include effect
-  # model: "AOD1B RL06"
-  model:
+  model: "AOD1B RL06"
   # filename of first file is series ...
-  data-file: data/AOD1B_2008-07-03_X_06.asc
+  data-file: data/AOD1B_2023-05-01_X_06.asc
   # ... and folder where subsequent files can be found (same naming convention)
   data-dir: data
   degree: 180
@@ -72,4 +71,8 @@ atmospheric-tide:
     r2: AOD1B_ATM_R2_06.asc
     s1: AOD1B_ATM_S1_06.asc
     s3: AOD1B_ATM_S3_06.asc
-    t3: AOD1B_ATM_T3_06.asc
+    t3: AOD1B_ATM_T3_06.asc
+
+satellite-attitude:
+  satellite: ja3
+  data_file: data/qua_ja3.csv

bin/eom.cpp

@@ -6,6 +6,8 @@
 #include "iers/iau.hpp"
 #include "iers/icgemio.hpp"
 #include "iers/relativity.hpp"
+#include "sysnsats/srp.hpp"
+#include "sysnsats/occultation.hpp"
 #include "integration_parameters.hpp"
 #include "sp3.hpp"
 #include "yaml-cpp/yaml.h"
@@ -14,12 +16,14 @@
 
 constexpr const double GM_Moon = /*0.49028010560e13;*/ 4902.800076e9;
 constexpr const double GM_Sun = /*1.32712442076e20;*/ 132712440040.944e9;
+constexpr const double EVERY_SEC = 180e0;
 
 /* Compute relevant quaternions for the ITRS/GCRS transformation
  */
 int prep_c2i(const dso::MjdEpoch &tai, dso::EopSeries &eops,
              Eigen::Quaterniond &q_c2tirs, Eigen::Quaterniond &q_tirs2i,
-             double *fargs, dso::EopRecord &eopr) noexcept {
+             double *fargs, dso::EopRecord &eopr) noexcept
+{
 
   /* epoch of request in TT */
   const auto tt = tai.tai2tt();
@@ -31,7 +35,8 @@ int prep_c2i(const dso::MjdEpoch &tai, dso::EopSeries &eops,
   dso::xycip06a(tt, Xcip, Ycip, fargs);
 
   /* interpolate EOPs */
-  if (dso::EopSeries::out_of_bounds(eops.interpolate(tt, eopr))) {
+  if (dso::EopSeries::out_of_bounds(eops.interpolate(tt, eopr)))
+  {
     fprintf(stderr, "Failed to interpolate: Epoch is out of bounds!\n");
     return 1;
   }
@@ -96,7 +101,8 @@ int prep_c2i(const dso::MjdEpoch &tai, dso::EopSeries &eops,
 
 int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
           dso::IntegrationParameters *params,
-          Eigen::Ref<Eigen::VectorXd> yp) noexcept {
+          Eigen::Ref<Eigen::VectorXd> yp) noexcept
+{
   /* epoch of request in TT */
   const auto tt =
       (params->t0().add_seconds(dso::FractionalSeconds(tsec))).tai2tt();
@@ -111,14 +117,16 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
 
   /* get Sun position & velocity in ICRF */
   Eigen::Matrix<double, 6, 1> rsun;
-  if (dso::planet_state(dso::Planet::SUN, tt, rsun)) {
+  if (dso::planet_state(dso::Planet::SUN, tt, rsun))
+  {
     fprintf(stderr, "ERROR Failed to compute Sun position!\n");
     assert(1 == 2);
   }
 
   /* get Moon position in ICRF */
   Eigen::Matrix<double, 3, 1> rmoon;
-  if (dso::planet_pos(dso::Planet::MOON, tt, rmoon)) {
+  if (dso::planet_pos(dso::Planet::MOON, tt, rmoon))
+  {
     fprintf(stderr, "ERROR Failed to compute Moon position!\n");
     assert(1 == 2);
   }
@@ -143,7 +151,8 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
   auto acstokes{params->earth_gravity()};
 
   /* add Solid Earth Tides to SH coeffs */
-  if (params->solid_earth_tide()) {
+  if (params->solid_earth_tide())
+  {
     params->solid_earth_tide()->stokes_coeffs(
         tt, tt.tt2ut1(eopr.dut()), q_tirs2i * (q_c2tirs * rmoon),
         q_tirs2i * (q_c2tirs * rsun.segment<3>(0)), fargs);
@@ -152,45 +161,52 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
   }
 
   /* add Ocean Tides to SH coeffs */
-  if (params->ocean_tide()) {
+  if (params->ocean_tide())
+  {
     params->ocean_tide()->stokes_coeffs(tt, tt.tt2ut1(eopr.dut()), fargs);
     acstokes += params->ocean_tide()->stokes_coeffs();
   }
 
   /* add Pole Tide to SH coeffs */
-  if (params->pole_tide()) {
+  if (params->pole_tide())
+  {
     double dC21, dS21;
     params->pole_tide()->stokes_coeffs(tt, eopr.xp(), eopr.yp(), dC21, dS21);
     acstokes.C(2, 1) += dC21;
     acstokes.S(2, 1) += dS21;
   }
 
   /* add Ocean Pole Tide to SH coeffs */
-  if (params->ocean_pole_tide()) {
-    if (params->ocean_pole_tide()->stokes_coeffs(tt, eopr.xp(), eopr.yp())) {
+  if (params->ocean_pole_tide())
+  {
+    if (params->ocean_pole_tide()->stokes_coeffs(tt, eopr.xp(), eopr.yp()))
+    {
       fprintf(stderr, "ERROR Failed computing Stokes Coefficients\n");
       return 1;
     }
     acstokes += params->ocean_pole_tide()->stokes_coeffs();
   }
 
   /* add deAliasing to SH coeffs */
-  if (params->dealias()) {
+  if (params->dealias())
+  {
     /*
     TODO!! WARNING!! The dealias instance should have a function that appends
     the coefficients at a StokesCoeffs instance!
     */
     auto tempstokes{params->earth_gravity()};
     if (params->dealias()->coefficients_at(
-            dso::from_mjdepoch<dso::nanoseconds>(tt), tempstokes)) {
+            dso::from_mjdepoch<dso::nanoseconds>(tt), tempstokes))
+    {
       fprintf(stderr, "Failed interpolating dealiasing coefficients\n");
       return 1;
     }
     acstokes += tempstokes;
   }
 
   /* add atmospheric tides to SH coeffs */
-  if (params->atmospheric_tide()) {
+  if (params->atmospheric_tide())
+  {
     params->atmospheric_tide()->stokes_coeffs(tt, tt.tt2ut1(eopr.dut()), fargs);
     acstokes += params->atmospheric_tide()->stokes_coeffs();
   }
@@ -201,7 +217,8 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
                                   itrs.segment<3>(0), ai, gi,
                                   params->earth_gravity().max_degree(),
                                   params->earth_gravity().max_order(), -1, -1,
-                                  &(params->tw()), &(params->tm()))) {
+                                  &(params->tw()), &(params->tm())))
+  {
     fprintf(stderr, "ERROR Failed computing acceleration/gradient\n");
     return 1;
   }
@@ -215,15 +232,35 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
     ac += dso::iers2010_relativistic_acceleration(y0, rsun);
   }
 
+  /* Solar Radiation Pressure */
+  if (params->matt)
+  {
+    const double of = dso::conical_occultation(y0.segment<3>(0), rsun.segment<3>(0));
+    if (of>0e0) {
+    /* get attitude */
+    if (params->matt->attitude_at(tt, *(params->mattdata)))
+    {
+      fprintf(stderr, "[ERROR] Failed getting attitude!\n");
+      return 9;
+    /* compute SRP acceleration */
+    ac += of * dso::solar_radiation_pressure(
+        params->msatmm->rotate_macromodel(params->mattdata->quaternions(), params->mattdata->angles()),
+        y0.segment<3>(0), rsun.segment<3>(0));
+    }
+    }
+  }
+
   /* form the derivative vector */
   yp.segment<3>(0) = y0.segment<3>(3);
   yp.segment<3>(3) = ac + q_c2tirs.conjugate() * (q_tirs2i.conjugate() * ai);
 
   return 0;
 }
 
-int main(int argc, char *argv[]) {
-  if (argc < 3) {
+int main(int argc, char *argv[])
+{
+  if (argc < 3)
+  {
     fprintf(stderr, "Usage: %s [CONFIG] [SP3] \n", argv[0]);
     return 1;
   }
@@ -233,9 +270,11 @@ int main(int argc, char *argv[]) {
 
   /* choose satellite */
   dso::sp3::SatelliteId sv = sp3.sattellite_vector()[0];
-  if (argc > 3) {
+  if (argc > 3)
+  {
     /* check if the satellite is included in the Sp3 */
-    if (!sp3.has_sv(dso::sp3::SatelliteId(argv[3]))) {
+    if (!sp3.has_sv(dso::sp3::SatelliteId(argv[3])))
+    {
       fprintf(stderr, "Error. Satellite [%s] not included in sp3 file!\n",
               argv[3]);
       return 1;
@@ -245,7 +284,8 @@ int main(int argc, char *argv[]) {
 
   /* get starting epoch in TT */
   auto start_t = sp3.start_epoch();
-  if (!std::strcmp(sp3.time_sys(), "GPS")) {
+  if (!std::strcmp(sp3.time_sys(), "GPS"))
+  {
     start_t = start_t.gps2tai();
   }
 
@@ -267,21 +307,25 @@ int main(int argc, char *argv[]) {
   Eigen::Quaterniond q_c2tirs, q_tirs2i;
   Eigen::Vector3d omega;
   dso::MjdEpoch t0;
+  dso::MjdEpoch tp;
   std::size_t it = 0;
   dso::Sp3DataBlock block;
   int sp3err = 0;
 
-  while (!sp3err) {
+  while (!sp3err)
+  {
     /* get next record from sp3 */
     sp3err = sp3.get_next_data_block(sv, block);
-    if (sp3err > 0) {
+    if (sp3err > 0)
+    {
       printf("Something went wrong ....status = %3d\n", sp3err);
       return 1;
     }
 
     /* time of current block in TAI */
     auto block_tai = block.t;
-    if (!std::strcmp(sp3.time_sys(), "GPS")) {
+    if (!std::strcmp(sp3.time_sys(), "GPS"))
+    {
       block_tai = block_tai.gps2tai();
     }
 
@@ -296,7 +340,8 @@ int main(int argc, char *argv[]) {
     y << block.state[0] * 1e3, block.state[1] * 1e3, block.state[2] * 1e3,
         block.state[4] * 1e-1, block.state[5] * 1e-1, block.state[6] * 1e-1;
 
-    if (!it) {
+    if (!it)
+    {
       /* transform state to GCRS (from ITRS) */
       Eigen::VectorXd yc = Eigen::Matrix<double, 6, 1>::Zero();
       yc.segment<3>(0) =
@@ -306,16 +351,25 @@ int main(int argc, char *argv[]) {
                           omega.cross(q_tirs2i.conjugate() * y.segment<3>(0)));
       y0 = yc;
       t0 = tai;
+      tp = t0;
       params.t0() = tai;
-    } else {
+    }
+    else if (tai.diff<dso::DateTimeDifferenceType::FractionalSeconds>(tp).seconds() >= EVERY_SEC)
+    {
       /* seconds since initial epoch */
       dso::FractionalSeconds sec =
           tai.diff<dso::DateTimeDifferenceType::FractionalSeconds>(t0);
 
       /* integrate from first SP3 record to here - setup initial conditions */
       params.t0() = t0;
       Eigen::VectorXd yc = Eigen::Matrix<double, 6, 1>::Zero();
-      if (dop853.integrate(0e0, sec.seconds(), y0, yc)) {
+      /* note: we are starting from the top here, so reload the attitude stream */
+      if (params.matt->reload()) {
+        fprintf(stderr, "ERROR. Failed reloading attitude stream!\n");
+        return -100;
+      }
+      if (dop853.integrate(0e0, sec.seconds(), y0, yc))
+      {
         fprintf(stderr, "ERROR. Integration failed! sec is %.3f\n",
                 sec.seconds());
         return 1;
@@ -331,12 +385,14 @@ int main(int argc, char *argv[]) {
           "%.9f %.3f %.3f %.3f %.6f %.6f %.6f %.3f %.3f %.3f %.6f %.6f %.6f\n",
           sec.seconds(), y(0), y(1), y(2), y(3), y(4), y(5), yi(0), yi(1),
           yi(2), yi(3), yi(4), yi(5));
+
+      tp = tai;
     }
 
     ++it;
-    if (it > 1000)
+    if (it > 100)
       break;
   }
 
   printf("# Num of epochs parsed/used: %ld\n", it);
-}
+}