@@ -272,6 +272,11 @@ void pck_calc(PckInfo *bpc, real epoch, int spiceId, real *rotMat,
272272 euler[4 ] = (real)angleDot[1 ];
273273 euler[5 ] = (real)angleDot[0 ];
274274
275+ if (spiceId == 31008 ){
276+ // Moon PA frame is stored relative to J2000 already
277+ euler313_to_rotMat (euler, rotMat, rotMatDot);
278+ return ;
279+ }
275280 real rotMatEclip[3 ][3 ], rotMatEclipDot[3 ][3 ];
276281 euler313_to_rotMat (euler, *rotMatEclip, *rotMatEclipDot);
277282 // We have ecliptic to body rotation matrix and its derivative
@@ -762,7 +767,8 @@ void get_pck_rotMat(const std::string &from, const std::string &to,
762767 std::vector<std::string> fromTo = {from, to};
763768 std::vector<std::string> validBodyFrames = {
764769 " IAU_SUN" " IAU_MERCURY" " IAU_VENUS"
765- " ITRF93" " IAU_EARTH" " IAU_MOON"
770+ " ITRF93" " IAU_EARTH"
771+ " MOON_PA_DE440" " MOON_ME_DE440_ME421" " IAU_MOON"
766772 " IAU_MARS" " IAU_JUPITER" " IAU_SATURN"
767773 " IAU_URANUS" " IAU_NEPTUNE" " IAU_PLUTO"
768774 " IAU_BENNU" " IAU_ITOKAWA" " IAU_GOLEVKA"
@@ -793,18 +799,34 @@ void get_pck_rotMat(const std::string &from, const std::string &to,
793799 iau_to_euler (t0_mjd, fromTo[bodyFrameIdx], euler);
794800 euler313_to_rotMat (euler, *rotMat, *rotMatDot);
795801 } else {
796- if (ephem.histPck == nullptr || ephem.latestPck == nullptr || ephem.predictPck == nullptr ){
797- throw std::invalid_argument (
798- " get_pck_rotMat: PCK kernels are not loaded. Memory map "
799- " the ephemeris using PropSimulation.map_ephemeris() method first."
802+ if (ephem.histPck == nullptr || ephem.latestPck == nullptr ||
803+ ephem.predictPck == nullptr || ephem.moonPck == nullptr ){
804+ throw std::invalid_argument (
805+ " get_pck_rotMat: PCK kernels are not loaded. Memory map "
806+ " the ephemeris using PropSimulation.map_ephemeris() method first."
807+ }
808+ int BinaryFrameSpiceId;
809+ if (fromTo[bodyFrameIdx] == " ITRF93"
810+ BinaryFrameSpiceId = 3000 ;
811+ if (t0_mjd < ephem.histPck ->targets [0 ].beg || t0_mjd > ephem.predictPck ->targets [0 ].end ){
812+ throw std::runtime_error (" get_pck_rotMat: The epoch is outside the range of the binary PCK files."
813+ }
814+ } else if (fromTo[bodyFrameIdx].substr (0 , 5 ) == " MOON_"
815+ BinaryFrameSpiceId = 31008 ;
816+ if (t0_mjd < ephem.moonPck ->targets [0 ].beg || t0_mjd > ephem.moonPck ->targets [0 ].end ){
817+ throw std::runtime_error (" get_pck_rotMat: The epoch is outside the range of the binary PCK files."
818+ }
819+ } else {
820+ throw std::runtime_error (" get_pck_rotMat: The binary body frame is not recognized."
800821 }
801822 // pick the correct PCK file memory map based on the epoch
802823 PckInfo *bpc;
803- if (t0_mjd < ephem.histPck ->targets [0 ].beg || t0_mjd > ephem.predictPck ->targets [0 ].end ){
804- throw std::runtime_error (" get_pck_rotMat: The epoch is outside the range of the binary PCK files."
805- }
806824 const real tSwitchMargin = 0.1 ;
807- if (t0_mjd <= ephem.histPck ->targets [0 ].end -tSwitchMargin){
825+ if (fromTo[bodyFrameIdx].substr (0 , 5 ) == " MOON_"
826+ bpc = ephem.moonPck ;
827+ // std::cout << "Using moonPck" << std::endl;
828+ }
829+ else if (t0_mjd <= ephem.histPck ->targets [0 ].end -tSwitchMargin){
808830 bpc = ephem.histPck ;
809831 // std::cout << "Using histPck" << std::endl;
810832 } else if (t0_mjd <= ephem.latestPck ->targets [0 ].end -tSwitchMargin){
@@ -814,13 +836,25 @@ void get_pck_rotMat(const std::string &from, const std::string &to,
814836 bpc = ephem.predictPck ;
815837 // std::cout << "Using predictPck" << std::endl;
816838 }
817- int BinaryFrameSpiceId;
818- if (fromTo[bodyFrameIdx] == " ITRF93"
819- BinaryFrameSpiceId = 3000 ;
820- } else {
821- throw std::runtime_error (" get_pck_rotMat: The binary body frame is not recognized."
822- }
823839 pck_calc (bpc, t0_mjd, BinaryFrameSpiceId, *rotMat, *rotMatDot);
840+ if (fromTo[bodyFrameIdx] == " MOON_ME_DE440_ME421"
841+ real rotMatCopy[3 ][3 ];
842+ real rotMatDotCopy[3 ][3 ];
843+ for (int i = 0 ; i < 3 ; i++){
844+ for (int j = 0 ; j < 3 ; j++){
845+ rotMatCopy[i][j] = rotMat[i][j];
846+ rotMatDotCopy[i][j] = rotMatDot[i][j];
847+ }
848+ }
849+ // from https://doi.org/10.3847/1538-3881/abd414, eqn. 19
850+ const real pa_to_me[3 ][3 ] = {
851+ {9.9999987311387650e-01 , -3.2895865791419379e-04 , 3.8152120821145725e-04 },
852+ {3.2895919698748533e-04 , 9.9999994589201047e-01 , -1.3502060036227025e-06 },
853+ {-3.8152074340615683e-04 , 1.4757107425872328e-06 , 9.9999992721986974e-01 }
854+ };
855+ mat3_mat3_mul (*pa_to_me, *rotMatCopy, *rotMat);
856+ mat3_mat3_mul (*pa_to_me, *rotMatDotCopy, *rotMatDot);
857+ }
824858 }
825859 // Okay, assemble the state rotation matrix
826860 if (bodyToInertial){
0 commit comments