update pck so it directly rotates states in AU, AU/d instead of km, km/s

rahil-makadia · Apr 28, 2024 · 6f6b6f9 · 6f6b6f9
1 parent b4767d6
commit 6f6b6f9
Show file tree

Hide file tree

Showing 4 changed files with 19 additions and 32 deletions.
diff --git a/src/approach.cpp b/src/approach.cpp
@@ -712,20 +712,7 @@ void ImpactParameters::get_impact_parameters(PropSimulation *propSim){
     get_baseBodyFrame(this->centralBodySpiceId, this->t, baseBodyFrame);
     std::vector<std::vector<real>> rotMat(6, std::vector<real>(6));
     get_pck_rotMat("J2000", baseBodyFrame, this->t, propSim->pckEphem, rotMat);
-    std::vector<real> impactRelStateInertial(6);
-    impactRelStateInertial[0] = this->xRel[0]*propSim->consts.du2m/1.0e3L;
-    impactRelStateInertial[1] = this->xRel[1]*propSim->consts.du2m/1.0e3L;
-    impactRelStateInertial[2] = this->xRel[2]*propSim->consts.du2m/1.0e3L;
-    impactRelStateInertial[3] = this->xRel[3]*propSim->consts.duptu2mps/1.0e3L;
-    impactRelStateInertial[4] = this->xRel[4]*propSim->consts.duptu2mps/1.0e3L;
-    impactRelStateInertial[5] = this->xRel[5]*propSim->consts.duptu2mps/1.0e3L;
-    mat_vec_mul(rotMat, impactRelStateInertial, this->xRelBodyFixed);
-    this->xRelBodyFixed[0] *= 1.0e3L/propSim->consts.du2m;
-    this->xRelBodyFixed[1] *= 1.0e3L/propSim->consts.du2m;
-    this->xRelBodyFixed[2] *= 1.0e3L/propSim->consts.du2m;
-    this->xRelBodyFixed[3] *= 1.0e3L/propSim->consts.duptu2mps;
-    this->xRelBodyFixed[4] *= 1.0e3L/propSim->consts.duptu2mps;
-    this->xRelBodyFixed[5] *= 1.0e3L/propSim->consts.duptu2mps;
+    mat_vec_mul(rotMat, this->xRel, this->xRelBodyFixed);
     real x, y, z, lon, lat, dist;
     x = this->xRelBodyFixed[0];
     y = this->xRelBodyFixed[1];

diff --git a/src/pck.cpp b/src/pck.cpp
@@ -679,7 +679,8 @@ void euler313_to_rotMat(const real euler[6], real *rotMat, real *rotMatDot){
     for (int i = 0; i < 3; i++){
         domega[i] = 0.0;
         for (int j = 0; j < 3; j++){
-            domega[i] += solutn[i][j] * euler[j+3];
+            // domega is in per second, convert to per day
+            domega[i] += solutn[i][j] * euler[j+3] * 86400.0;
         }
     }
     real rotDotTimesRotTranspose[3][3] = {

diff --git a/src/simulation.cpp b/src/simulation.cpp
@@ -96,25 +96,18 @@ void get_observer_state(const real &tObsMjd,
     real lon = observerInfo[1];
     real lat = observerInfo[2];
     real rho = observerInfo[3];
-    real bodyFixedX = rho * cos(lat) * cos(lon) / 1.0e3L;
-    real bodyFixedY = rho * cos(lat) * sin(lon) / 1.0e3L;
-    real bodyFixedZ = rho * sin(lat) / 1.0e3L;
+    real bodyFixedX = rho * cos(lat) * cos(lon) / propSim->consts.du2m;
+    real bodyFixedY = rho * cos(lat) * sin(lon) / propSim->consts.du2m;
+    real bodyFixedZ = rho * sin(lat) / propSim->consts.du2m;
     std::vector<real> bodyFixedState = {bodyFixedX, bodyFixedY, bodyFixedZ,
                                             0.0,        0.0,        0.0};
-    std::vector<real> observerStateInertial(6);
-    mat_vec_mul(rotMat, bodyFixedState, observerStateInertial);
-    observerStateInertial[0] *= 1.0e3L / propSim->consts.du2m;
-    observerStateInertial[1] *= 1.0e3L / propSim->consts.du2m;
-    observerStateInertial[2] *= 1.0e3L / propSim->consts.du2m;
-    observerStateInertial[3] *= 1.0e3L / propSim->consts.duptu2mps;
-    observerStateInertial[4] *= 1.0e3L / propSim->consts.duptu2mps;
-    observerStateInertial[5] *= 1.0e3L / propSim->consts.duptu2mps;
-    observerState[0] = baseBodyState[0] + observerStateInertial[0];
-    observerState[1] = baseBodyState[1] + observerStateInertial[1];
-    observerState[2] = baseBodyState[2] + observerStateInertial[2];
-    observerState[3] = baseBodyState[3] + observerStateInertial[3];
-    observerState[4] = baseBodyState[4] + observerStateInertial[4];
-    observerState[5] = baseBodyState[5] + observerStateInertial[5];
+    mat_vec_mul(rotMat, bodyFixedState, observerState);
+    observerState[0] += baseBodyState[0];
+    observerState[1] += baseBodyState[1];
+    observerState[2] += baseBodyState[2];
+    observerState[3] += baseBodyState[3];
+    observerState[4] += baseBodyState[4];
+    observerState[5] += baseBodyState[5];
 }
 
 /** 

diff --git a/tests/cpp/prop/pck_map.cpp b/tests/cpp/prop/pck_map.cpp
@@ -55,6 +55,12 @@ int main(){
     SpiceDouble spiceMat[6][6];
     sxform_c(from.c_str(), to.c_str(), et, spiceMat);
     kclear_c();
+    // convert SPICE frame so it rotates au, au/day instead of km, km/s
+    for (int i = 3; i < 6; i++){
+        for (int j = 0; j < 3; j++){
+            spiceMat[i][j] *= 86400.0;
+        }
+    }
 
     // // print both matrices
     // std::cout << "Spice Matrix:" << std::endl;