opik total derivatives

rahil-makadia · Jul 17, 2024 · 6c98896 · 6c98896
1 parent 7fbcdbd
commit 6c98896
Show file tree

Hide file tree

Showing 7 changed files with 100 additions and 28 deletions.
diff --git a/grss/prop/prop_utils.py b/grss/prop/prop_utils.py
@@ -360,7 +360,7 @@ def partials_to_ellipse(ca, au2units, n_std, print_ellipse_params, units, analyt
     part_tlin_minus_t = ca.dTLinMinusT + extra_zeros
     part_tlin = part_tlin_minus_t @ stm_tca_tmap
     part_t = ca.dt + extra_zeros
-    # part_tlin += part_t
+    part_tlin += part_t
     part_kizner = np.zeros((3,stm_tca_tmap.shape[0]))
     part_kizner[:2,:] = np.array([ca.kizner.dx+extra_zeros,ca.kizner.dy+extra_zeros])@stm_tca_tmap
     part_kizner[2,:] = part_tlin
@@ -369,6 +369,7 @@ def partials_to_ellipse(ca, au2units, n_std, print_ellipse_params, units, analyt
     part_scaled[2,:] = part_tlin
     part_opik = np.zeros((3,stm_tca_tmap.shape[0]))
     part_opik[:2,:] = np.array([ca.opik.dx+extra_zeros,ca.opik.dy+extra_zeros])@stm_tca_tmap
+    part_opik[:2,:6] += ca.dOpikTotalDerivTerm2
     part_opik[2,:] = part_tlin
     part_mtp = np.zeros((3,stm_tca_tmap.shape[0]))
     part_mtp[:2,:] = np.array([ca.mtp.dx+extra_zeros,ca.mtp.dy+extra_zeros])@stm_tca_tmap

diff --git a/grss/version.txt b/grss/version.txt
@@ -1 +1 @@
-4.1.4
+4.1.5
diff --git a/include/simulation.h b/include/simulation.h
@@ -298,6 +298,7 @@ struct BPlaneParameters {
  * @param mtp Modified Target Plane parameters.
  * @param dTLinMinusT Partial derivatives of the (linearized intersection minus map) time with respect to the state at the close approach.
  * @param dt Partial derivatives of the time of closest approach with respect to the state at the close approach.
+ * @param dOpikTotalDerivTerm2 Partial derivatives of the Öpik total derivative term 2 (non-constant planet velocity).
  */
 class CloseApproachParameters {
    private:
@@ -328,6 +329,7 @@ class CloseApproachParameters {
     BPlaneParameters mtp;
     std::vector<real> dTLinMinusT = std::vector<real>(6, 0.0L);
     std::vector<real> dt = std::vector<real>(6, 0.0L);
+    std::vector< std::vector<real> > dOpikTotalDerivTerm2 = std::vector< std::vector<real> >(2, std::vector<real>(6, 0.0L));
     /**
      * @brief Get the close approach parameters.
      */

diff --git a/src/approach.cpp b/src/approach.cpp
@@ -468,9 +468,9 @@ void CloseApproachParameters::get_ca_parameters(PropSimulation *propSim, const r
     for (size_t k = 0; k < 3; k++) {
         vCentralBodyHelio[k] = xCentralBody[3+k]-xSun[3+k];
     }
-    real xiHat[3], zetaHat[3], vCentralBodyHelioCrossSHat[3];
-    vcross(vCentralBodyHelio, sHat, vCentralBodyHelioCrossSHat);
-    vunit(vCentralBodyHelioCrossSHat, 3, xiHat);
+    real xiHat[3], zetaHat[3], vPlanetCrossSHatVec[3];
+    vcross(vCentralBodyHelio, sHat, vPlanetCrossSHatVec);
+    vunit(vPlanetCrossSHatVec, 3, xiHat);
     vcross(sHat, xiHat, zetaHat);
     for (size_t k = 0; k < 3; k++) {
         zetaHat[k] *= -1;
@@ -655,26 +655,76 @@ void CloseApproachParameters::get_ca_parameters(PropSimulation *propSim, const r
             this->scaled.dy[k] = (this->kizner.dy[k] - this->scaled.y*partial_lambda[k])/this->gravFocusFactor;
         }
 
-        // // Acceleration of the Sun is not included in the calculation
-        // // because it is negligibly small
-        // real accSun[3];
-        // memset(accSun, 0, 3*sizeof(real));
-        // for (size_t k = 0; k < propSim->integParams.nSpice; k++) {
-        //     if (propSim->spiceBodies[k].spiceId == 10) {
-        //         accSun[0] = propSim->spiceBodies[k].acc[0];
-        //         accSun[1] = propSim->spiceBodies[k].acc[1];
-        //         accSun[2] = propSim->spiceBodies[k].acc[2];
-        //         break;
-        //     }
-        // }
-        real **partial_vel_planet = new real*[6];
+        // Acceleration of the Sun is not included in the calculation
+        real accSun[3];
+        for (size_t k = 0; k < propSim->integParams.nSpice; k++) {
+            if (propSim->spiceBodies[k].spiceId == 10) {
+                accSun[0] = propSim->spiceBodies[k].acc[0];
+                accSun[1] = propSim->spiceBodies[k].acc[1];
+                accSun[2] = propSim->spiceBodies[k].acc[2];
+                break;
+            }
+        }
+        real **partial_vel_planet = new real*[3];
+        for (size_t k = 0; k < 3; k++) {
+            partial_vel_planet[k] = new real[6];
+            for (size_t k2 = 0; k2 < 6; k2++) {
+                partial_vel_planet[k][k2] = this->dt[k2]*(accPlanet[k] - accSun[k]);
+            }
+        }
+        // partials of xi and zeta w.r.t planet velocity are needed for total derivative
+        real **partial_vpl_vpl = new real*[3];
+        for (size_t k = 0; k < 3; k++) {
+            partial_vpl_vpl[k] = new real[3];
+            for (size_t k2 = 0; k2 < 3; k2++) {
+                partial_vpl_vpl[k][k2] = 0;
+            }
+            partial_vpl_vpl[k][k] = 1;
+        }
+        real vPlanetCrossSHat;
+        vnorm(vPlanetCrossSHatVec, 3, vPlanetCrossSHat);
+        real **partial_vPlanetCrossSHatVec_vpl = new real*[3];
+        for (size_t k = 0; k < 3; k++) {
+            partial_vPlanetCrossSHatVec_vpl[k] = new real[3];
+            vcross(partial_vpl_vpl[k], sHat, partial_vPlanetCrossSHatVec_vpl[k]);
+        }
+        real *partial_vPlanetCrossSHat_vpl = new real[3];
+        vcross(xiHat, sHat, partial_vPlanetCrossSHat_vpl);
+        for (size_t k = 0; k < 3; k++){
+            partial_vPlanetCrossSHat_vpl[k] *= -1.0;
+        }
+        real **partial_xiHat_vpl = new real*[3];
+        real **partial_zetaHat_vpl = new real*[3];
+        for (size_t k = 0; k < 3; k++) {
+            partial_xiHat_vpl[k] = new real[3];
+            for (size_t k2 = 0; k2 < 3; k2++) {
+                partial_xiHat_vpl[k][k2] =
+                    (vPlanetCrossSHat * partial_vPlanetCrossSHatVec_vpl[k][k2] -
+                     vPlanetCrossSHatVec[k2] * partial_vPlanetCrossSHat_vpl[k]) /
+                    vPlanetCrossSHat / vPlanetCrossSHat;
+            }
+            partial_zetaHat_vpl[k] = new real[3];
+            vcross(sHat, partial_xiHat_vpl[k], partial_zetaHat_vpl[k]);
+            for (size_t k2 = 0; k2 < 3; k2++) {
+                partial_zetaHat_vpl[k][k2] *= -1.0;
+            }
+        }
+        real *partial_xi_vpl = new real[3];
+        real *partial_zeta_vpl = new real[3];
+        for (size_t k = 0; k < 3; k++) {
+            vdot(partial_zetaHat_vpl[k], hVec, 3, partial_xi_vpl[k]);
+            partial_xi_vpl[k] /= vInf;
+            vdot(partial_xiHat_vpl[k], hVec, 3, partial_zeta_vpl[k]);
+            partial_zeta_vpl[k] /= -vInf;
+        }
         for (size_t k = 0; k < 6; k++) {
-            partial_vel_planet[k] = new real[3];
             for (size_t k2 = 0; k2 < 3; k2++) {
-                // partial_vel_planet[k][k2] = this->dt[k]*(accPlanet[k2] - accSun[k2]);
-                partial_vel_planet[k][k2] = this->dt[k]*accPlanet[k2];
+                temp1Vec[k2] = partial_vel_planet[k2][k];
             }
+            vdot(partial_xi_vpl, temp1Vec, 3, this->dOpikTotalDerivTerm2[0][k]);
+            vdot(partial_zeta_vpl, temp1Vec, 3, this->dOpikTotalDerivTerm2[1][k]);
         }
+
         real **partial_xiHat = new real*[6];
         real temp3, temp3Vec[3];
         real temp4, temp4Vec[3];
@@ -723,10 +773,16 @@ void CloseApproachParameters::get_ca_parameters(PropSimulation *propSim, const r
             delete[] partial_rHat[k];
             delete[] partial_sHat[k];
             delete[] partial_tHat[k];
-            delete[] partial_vel_planet[k];
             delete[] partial_xiHat[k];
             delete[] partial_zetaHat[k];
         }
+        for (size_t k = 0; k < 3; k++) {
+            delete[] partial_vel_planet[k];
+            delete[] partial_vpl_vpl[k];
+            delete[] partial_vPlanetCrossSHatVec_vpl[k];
+            delete[] partial_xiHat_vpl[k];
+            delete[] partial_zetaHat_vpl[k];
+        }
         delete[] partial_alpha;
         delete[] partial_vInf;
         delete[] partial_a;
@@ -742,6 +798,13 @@ void CloseApproachParameters::get_ca_parameters(PropSimulation *propSim, const r
         delete[] partial_tHat;
         delete[] partial_lambda;
         delete[] partial_vel_planet;
+        delete[] partial_vpl_vpl;
+        delete[] partial_vPlanetCrossSHatVec_vpl;
+        delete[] partial_vPlanetCrossSHat_vpl;
+        delete[] partial_xiHat_vpl;
+        delete[] partial_zetaHat_vpl;
+        delete[] partial_xi_vpl;
+        delete[] partial_zeta_vpl;
         delete[] partial_xiHat;
         delete[] partial_zetaHat;
     }

diff --git a/src/grss.cpp b/src/grss.cpp
@@ -202,7 +202,8 @@ PYBIND11_MODULE(libgrss, m) {
                        R"mydelimiter(
             Name of the central body.
             )mydelimiter")
-        .def_readwrite("centralBodyIdx", &CloseApproachParameters::centralBodyIdx,
+        .def_readwrite("centralBodyIdx",
+                       &CloseApproachParameters::centralBodyIdx,
                        R"mydelimiter(
             Index of the central body.
             )mydelimiter")
@@ -252,6 +253,11 @@ PYBIND11_MODULE(libgrss, m) {
         .def_readwrite("dt", &CloseApproachParameters::dt, R"mydelimiter(
             Partials of time of periapsis with respect to CA state.
             )mydelimiter")
+        .def_readwrite("dOpikTotalDerivTerm2",
+                       &CloseApproachParameters::dOpikTotalDerivTerm2,
+                       R"mydelimiter(
+            Partial derivatives of the Öpik total derivative term 2 (non-constant planet velocity).
+            )mydelimiter")
         .def("get_ca_parameters", &CloseApproachParameters::get_ca_parameters,
              py::arg("propSim"), py::arg("tMap"), R"mydelimiter(
             Calculate the close approach parameters.

diff --git a/src/pck.cpp b/src/pck.cpp
@@ -24,7 +24,7 @@ void pck_free(PckInfo* bpc) {
         free(bpc->targets);
     }
     munmap(bpc->map, bpc->len);
-    memset(bpc, 0, sizeof(PckInfo));
+    memset((void *)bpc, 0, sizeof(PckInfo));
     free(bpc);
     bpc = nullptr;
 }

diff --git a/src/spk.cpp b/src/spk.cpp
@@ -24,7 +24,7 @@ void spk_free(SpkInfo* bsp) {
         free(bsp->targets);
     }
     munmap(bsp->map, bsp->len);
-    memset(bsp, 0, sizeof(SpkInfo));
+    memset((void *)bsp, 0, sizeof(SpkInfo));
     free(bsp);
     bsp = nullptr;
 }
@@ -227,17 +227,17 @@ void spk_calc(SpkInfo *bsp, double epoch, int spiceId, double *state) {
     U[0] = 0.0;
     U[1] = 0.0;
     U[2] = 4.0;
-    for (size_t p = 2; p < P; p++) {
+    for (int p = 2; p < P; p++) {
         T[p] = 2.0 * z * T[p - 1] - T[p - 2];
         S[p] = 2.0 * z * S[p - 1] + 2.0 * T[p - 1] - S[p - 2];
     }
-    for (size_t p = 3; p < P; p++) {
+    for (int p = 3; p < P; p++) {
         U[p] = 2.0 * z * U[p - 1] + 4.0 * S[p - 1] - U[p - 2];
     }
     const double c = 1.0 / val[1]; // derivative scaling factor from SPICE/spke02.f and chbint.f
     for (size_t i = 0; i < 3; i++) {
         const int b = 2 + i * P;
-        for (size_t p = 0; p < P; p++) {
+        for (int p = 0; p < P; p++) {
             const double v = val[b + p];
             state[i] += v * T[p] / 149597870.7;
             state[i + 3] += v * S[p] * c / 149597870.7 * 86400.0;