Program Listing for File solar_system.cc¶
↰ Return to documentation for file (environment/solar_system.cc)
#include "lupnt/environment/solar_system.h"
#include "lupnt/conversions/time_conversions.h"
#include "lupnt/numerics/math_utils.h"
namespace lupnt {
Real MeanObliquity(Real mjd_tt) {
const Real T = (mjd_tt - MJD_J2000_TT) / 36525.0;
return (23.43929111 - (46.8150 + (0.00059 - 0.001813 * T) * T) * T / 3600.0) * RAD;
}
Mat3 Equatorial2EclipticMatrix(Real mjd_tt) { return RotX(MeanObliquity(mjd_tt)); }
Mat3 PrecessionMatrix(Real mjd_1, Real mjd_2) {
Real T = (mjd_1 - MJD_J2000_TT) / 36525.0;
Real dT = (mjd_2 - mjd_1) / 36525.0;
// Precession angles
Real zeta = ((2306.2181 + (1.39656 - 0.000139 * T) * T)
+ ((0.30188 - 0.000344 * T) + 0.017998 * dT) * dT)
* dT * RAD_ARCSEC;
Real z = zeta + ((0.79280 + 0.000411 * T) + 0.000205 * dT) * dT * dT * RAD_ARCSEC;
Real theta = ((2004.3109 - (0.85330 + 0.000217 * T) * T)
- ((0.42665 + 0.000217 * T) + 0.041833 * dT) * dT)
* dT * RAD_ARCSEC;
// Precession matrix
Mat3 P = RotZ(-z) * RotX(theta) * RotZ(-zeta);
return P;
}
std::pair<Real, Real> NutAngles(Real mjd_tt) {
const Real T = (mjd_tt - MJD_J2000_TT) / 36525.0;
const Real T2 = T * T;
const Real T3 = T2 * T;
const Real rev = 360.0 * 3600.0; // arcsec/revolution
const int N_coeff = 106;
const long C[N_coeff][9] = {
//
// l l' F D Om dpsi *T deps *T #
//
{0, 0, 0, 0, 1, -1719960, -1742, 920250, 89}, // 1
{0, 0, 0, 0, 2, 20620, 2, -8950, 5}, // 2
{-2, 0, 2, 0, 1, 460, 0, -240, 0}, // 3
{2, 0, -2, 0, 0, 110, 0, 0, 0}, // 4
{-2, 0, 2, 0, 2, -30, 0, 10, 0}, // 5
{1, -1, 0, -1, 0, -30, 0, 0, 0}, // 6
{0, -2, 2, -2, 1, -20, 0, 10, 0}, // 7
{2, 0, -2, 0, 1, 10, 0, 0, 0}, // 8
{0, 0, 2, -2, 2, -131870, -16, 57360, -31}, // 9
{0, 1, 0, 0, 0, 14260, -34, 540, -1}, // 10
{0, 1, 2, -2, 2, -5170, 12, 2240, -6}, // 11
{0, -1, 2, -2, 2, 2170, -5, -950, 3}, // 12
{0, 0, 2, -2, 1, 1290, 1, -700, 0}, // 13
{2, 0, 0, -2, 0, 480, 0, 10, 0}, // 14
{0, 0, 2, -2, 0, -220, 0, 0, 0}, // 15
{0, 2, 0, 0, 0, 170, -1, 0, 0}, // 16
{0, 1, 0, 0, 1, -150, 0, 90, 0}, // 17
{0, 2, 2, -2, 2, -160, 1, 70, 0}, // 18
{0, -1, 0, 0, 1, -120, 0, 60, 0}, // 19
{-2, 0, 0, 2, 1, -60, 0, 30, 0}, // 20
{0, -1, 2, -2, 1, -50, 0, 30, 0}, // 21
{2, 0, 0, -2, 1, 40, 0, -20, 0}, // 22
{0, 1, 2, -2, 1, 40, 0, -20, 0}, // 23
{1, 0, 0, -1, 0, -40, 0, 0, 0}, // 24
{2, 1, 0, -2, 0, 10, 0, 0, 0}, // 25
{0, 0, -2, 2, 1, 10, 0, 0, 0}, // 26
{0, 1, -2, 2, 0, -10, 0, 0, 0}, // 27
{0, 1, 0, 0, 2, 10, 0, 0, 0}, // 28
{-1, 0, 0, 1, 1, 10, 0, 0, 0}, // 29
{0, 1, 2, -2, 0, -10, 0, 0, 0}, // 30
{0, 0, 2, 0, 2, -22740, -2, 9770, -5}, // 31
{1, 0, 0, 0, 0, 7120, 1, -70, 0}, // 32
{0, 0, 2, 0, 1, -3860, -4, 2000, 0}, // 33
{1, 0, 2, 0, 2, -3010, 0, 1290, -1}, // 34
{1, 0, 0, -2, 0, -1580, 0, -10, 0}, // 35
{-1, 0, 2, 0, 2, 1230, 0, -530, 0}, // 36
{0, 0, 0, 2, 0, 630, 0, -20, 0}, // 37
{1, 0, 0, 0, 1, 630, 1, -330, 0}, // 38
{-1, 0, 0, 0, 1, -580, -1, 320, 0}, // 39
{-1, 0, 2, 2, 2, -590, 0, 260, 0}, // 40
{1, 0, 2, 0, 1, -510, 0, 270, 0}, // 41
{0, 0, 2, 2, 2, -380, 0, 160, 0}, // 42
{2, 0, 0, 0, 0, 290, 0, -10, 0}, // 43
{1, 0, 2, -2, 2, 290, 0, -120, 0}, // 44
{2, 0, 2, 0, 2, -310, 0, 130, 0}, // 45
{0, 0, 2, 0, 0, 260, 0, -10, 0}, // 46
{-1, 0, 2, 0, 1, 210, 0, -100, 0}, // 47
{-1, 0, 0, 2, 1, 160, 0, -80, 0}, // 48
{1, 0, 0, -2, 1, -130, 0, 70, 0}, // 49
{-1, 0, 2, 2, 1, -100, 0, 50, 0}, // 50
{1, 1, 0, -2, 0, -70, 0, 0, 0}, // 51
{0, 1, 2, 0, 2, 70, 0, -30, 0}, // 52
{0, -1, 2, 0, 2, -70, 0, 30, 0}, // 53
{1, 0, 2, 2, 2, -80, 0, 30, 0}, // 54
{1, 0, 0, 2, 0, 60, 0, 0, 0}, // 55
{2, 0, 2, -2, 2, 60, 0, -30, 0}, // 56
{0, 0, 0, 2, 1, -60, 0, 30, 0}, // 57
{0, 0, 2, 2, 1, -70, 0, 30, 0}, // 58
{1, 0, 2, -2, 1, 60, 0, -30, 0}, // 59
{0, 0, 0, -2, 1, -50, 0, 30, 0}, // 60
{1, -1, 0, 0, 0, 50, 0, 0, 0}, // 61
{2, 0, 2, 0, 1, -50, 0, 30, 0}, // 62
{0, 1, 0, -2, 0, -40, 0, 0, 0}, // 63
{1, 0, -2, 0, 0, 40, 0, 0, 0}, // 64
{0, 0, 0, 1, 0, -40, 0, 0, 0}, // 65
{1, 1, 0, 0, 0, -30, 0, 0, 0}, // 66
{1, 0, 2, 0, 0, 30, 0, 0, 0}, // 67
{1, -1, 2, 0, 2, -30, 0, 10, 0}, // 68
{-1, -1, 2, 2, 2, -30, 0, 10, 0}, // 69
{-2, 0, 0, 0, 1, -20, 0, 10, 0}, // 70
{3, 0, 2, 0, 2, -30, 0, 10, 0}, // 71
{0, -1, 2, 2, 2, -30, 0, 10, 0}, // 72
{1, 1, 2, 0, 2, 20, 0, -10, 0}, // 73
{-1, 0, 2, -2, 1, -20, 0, 10, 0}, // 74
{2, 0, 0, 0, 1, 20, 0, -10, 0}, // 75
{1, 0, 0, 0, 2, -20, 0, 10, 0}, // 76
{3, 0, 0, 0, 0, 20, 0, 0, 0}, // 77
{0, 0, 2, 1, 2, 20, 0, -10, 0}, // 78
{-1, 0, 0, 0, 2, 10, 0, -10, 0}, // 79
{1, 0, 0, -4, 0, -10, 0, 0, 0}, // 80
{-2, 0, 2, 2, 2, 10, 0, -10, 0}, // 81
{-1, 0, 2, 4, 2, -20, 0, 10, 0}, // 82
{2, 0, 0, -4, 0, -10, 0, 0, 0}, // 83
{1, 1, 2, -2, 2, 10, 0, -10, 0}, // 84
{1, 0, 2, 2, 1, -10, 0, 10, 0}, // 85
{-2, 0, 2, 4, 2, -10, 0, 10, 0}, // 86
{-1, 0, 4, 0, 2, 10, 0, 0, 0}, // 87
{1, -1, 0, -2, 0, 10, 0, 0, 0}, // 88
{2, 0, 2, -2, 1, 10, 0, -10, 0}, // 89
{2, 0, 2, 2, 2, -10, 0, 0, 0}, // 90
{1, 0, 0, 2, 1, -10, 0, 0, 0}, // 91
{0, 0, 4, -2, 2, 10, 0, 0, 0}, // 92
{3, 0, 2, -2, 2, 10, 0, 0, 0}, // 93
{1, 0, 2, -2, 0, -10, 0, 0, 0}, // 94
{0, 1, 2, 0, 1, 10, 0, 0, 0}, // 95
{-1, -1, 0, 2, 1, 10, 0, 0, 0}, // 96
{0, 0, -2, 0, 1, -10, 0, 0, 0}, // 97
{0, 0, 2, -1, 2, -10, 0, 0, 0}, // 98
{0, 1, 0, 2, 0, -10, 0, 0, 0}, // 99
{1, 0, -2, -2, 0, -10, 0, 0, 0}, // 100
{0, -1, 2, 0, 1, -10, 0, 0, 0}, // 101
{1, 1, 0, -2, 1, -10, 0, 0, 0}, // 102
{1, 0, -2, 2, 0, -10, 0, 0, 0}, // 103
{2, 0, 0, 2, 0, 10, 0, 0, 0}, // 104
{0, 0, 2, 4, 2, -10, 0, 0, 0}, // 105
{0, 1, 0, 1, 0, 10, 0, 0, 0} // 106
};
// Mean arguments of luni-solar motion
//
// l mean anomaly of the Moon
// l' mean anomaly of the Sun
// F mean argument of latitude
// D mean longitude elongation of the Moon from the Sun
// Om mean longitude of the ascending node
Real l = mod(485866.733 + (1325.0 * rev + 715922.633) * T + 31.310 * T2 + 0.064 * T3, rev);
Real lp = mod(1287099.804 + (99.0 * rev + 1292581.224) * T - 0.577 * T2 - 0.012 * T3, rev);
Real F = mod(335778.877 + (1342.0 * rev + 295263.137) * T - 13.257 * T2 + 0.011 * T3, rev);
Real D = mod(1072261.307 + (1236.0 * rev + 1105601.328) * T - 6.891 * T2 + 0.019 * T3, rev);
Real Om = mod(450160.280 - (5.0 * rev + 482890.539) * T + 7.455 * T2 + 0.008 * T3, rev);
// Nutation in longitude and obliquity [rad]
Real arg;
Real deps = 0, dpsi = 0;
for (int i = 0; i < N_coeff; i++) {
arg = (C[i][0] * l + C[i][1] * lp + C[i][2] * F + C[i][3] * D + C[i][4] * Om) * RAD_ARCSEC;
dpsi += (C[i][5] + C[i][6] * T) * sin(arg);
deps += (C[i][7] + C[i][8] * T) * cos(arg);
};
dpsi = 1.0E-5 * dpsi * RAD_ARCSEC;
deps = 1.0E-5 * deps * RAD_ARCSEC;
return {dpsi, deps};
}
Mat3 NutationMatrix(Real mjd_tt) {
// Mean obliquity of the ecliptic
Real eps = MeanObliquity(mjd_tt);
// Nutation in longitude and obliquity
auto [dpsi, deps] = NutAngles(mjd_tt);
// Transformation from mean to true equator and equinox
Mat3 R = RotX(-eps - deps) * RotZ(-dpsi) * RotX(+eps);
return R;
}
Mat3 NutationMatrixLowPrecision(Real mjd_tt) {
Real T = (mjd_tt - MJD_J2000_TT) / 36525.0;
// Mean arguments of luni-solar motion
// ls mean anomaly of the Sun
// D diff. longitude Moon-Sun
// F mean argument of latitude
// N longit. ascending node
Real ls = TWO_PI * frac(0.993133 + 99.997306 * T);
Real D = TWO_PI * frac(0.827362 + 1236.853087 * T);
Real F = TWO_PI * frac(0.259089 + 1342.227826 * T);
Real N = TWO_PI * frac(0.347346 - 5.372447 * T);
// Nutation angles
Real dpsi = (-17.200 * sin(N) - 1.319 * sin(2 * (F - D + N)) - 0.227 * sin(2 * (F + N))
+ 0.206 * sin(2 * N) + 0.143 * sin(ls))
* RAD_ARCSEC;
Real deps = (+9.203 * cos(N) + 0.574 * cos(2 * (F - D + N)) + 0.098 * cos(2 * (F + N))
- 0.090 * cos(2 * N))
* RAD_ARCSEC;
// Mean obliquity of the ecliptic
Real eps = 0.4090928 - 2.2696E-4 * T;
Mat3 R = RotX(-eps - deps) * RotZ(-dpsi) * RotX(+eps);
return R;
}
Real EquinoxEquation(Real mjd_tt) {
// Nutation in longitude and obliquity
auto [dpsi, reps] = NutAngles(mjd_tt);
// Equation of the equinoxes
return dpsi * cos(MeanObliquity(mjd_tt));
};
Vec3 SunPositionLowPrecision(Real mjd_tt) {
const Real eps = 23.43929111 * RAD; // Obliquity of J2000 ecliptic
const Real T = (mjd_tt - MJD_J2000_TT) / 36525.0; // Julian cent. since J2000
// Mean anomaly, ecliptic longitude and radius
Real M = TWO_PI * frac(0.9931267 + 99.9973583 * T); // [rad]
Real L = TWO_PI
* frac(0.7859444 + M / TWO_PI
+ (6892.0 * sin(M) + 72.0 * sin(2.0 * M)) / 1296.0e3); // [rad]
Real r = 149.619e9 - 2.499e9 * cos(M) - 0.021e9 * cos(2 * M); // [m]
// Equatorial position vector
Vec3 r_sun = RotX(-eps) * Vec3(r * cos(L), r * sin(L), 0.0);
return r_sun;
}
Vec3 MoonPositionLowPrecision(Real mjd_tt) {
const Real eps = 23.43929111 * RAD; // Obliquity of J2000 ecliptic
const Real T = (mjd_tt - MJD_J2000_TT) / 36525.0; // Julian cent. since J2000
// Mean elements of lunar orbit
Real L_0 = frac(0.606433 + 1336.851344 * T); // Mean longitude [rev] w.r.t. J2000 equinox
Real l = TWO_PI * frac(0.374897 + 1325.552410 * T); // Moon's mean anomaly [rad]
Real lp = TWO_PI * frac(0.993133 + 99.997361 * T); // Sun's mean anomaly [rad]
Real D = TWO_PI * frac(0.827361 + 1236.853086 * T); // Diff. long. Moon-Sun [rad]
Real F = TWO_PI * frac(0.259086 + 1342.227825 * T); // Argument of latitude
// Ecliptic longitude (w.r.t. equinox of J2000)
Real dL = +22640 * sin(l) - 4586 * sin(l - 2 * D) + 2370 * sin(2 * D) + 769 * sin(2 * l)
- 668 * sin(lp) - 412 * sin(2 * F) - 212 * sin(2 * l - 2 * D)
- 206 * sin(l + lp - 2 * D) + 192 * sin(l + 2 * D) - 165 * sin(lp - 2 * D)
- 125 * sin(D) - 110 * sin(l + lp) + 148 * sin(l - lp) - 55 * sin(2 * F - 2 * D);
Real L = TWO_PI * frac(L_0 + dL / 1296.0e3); // [rad]
// Ecliptic latitude
Real S = F + (dL + 412 * sin(2 * F) + 541 * sin(lp)) * RAD_ARCSEC;
Real h = F - 2 * D;
Real N = -526 * sin(h) + 44 * sin(l + h) - 31 * sin(-l + h) - 23 * sin(lp + h)
+ 11 * sin(-lp + h) - 25 * sin(-2 * l + F) + 21 * sin(-l + F);
Real B = (18520.0 * sin(S) + N) * RAD_ARCSEC; // [rad]
Real cosB = cos(B);
// Distance [m]
Real R = 385000e3 - 20905e3 * cos(l) - 3699e3 * cos(2 * D - l) - 2956e3 * cos(2 * D)
- 570e3 * cos(2 * l) + 246e3 * cos(2 * l - 2 * D) - 205e3 * cos(lp - 2 * D)
- 171e3 * cos(l + 2 * D) - 152e3 * cos(l + lp - 2 * D);
// Equatorial coordinates
Vec3 r_moon = RotX(-eps) * Vec3(R * cos(L) * cosB, R * sin(L) * cosB, R * sin(B));
return r_moon;
}
//------------------------------------------------------------------------------
//
// GHAMatrix
//
// Purpose:
//
// Transformation from true equator and equinox to Earth equator and
// Greenwich meridian system
//
// Input/Output:
//
// mjd_ut1 Modified Julian Date UT1
// <return> Greenwich Hour Angle matrix
//
//------------------------------------------------------------------------------
Mat3 GreenwichHourAngleMatrix(Real mjd_ut1) {
return RotZ(GreenwichApparentSiderealTime(mjd_ut1));
}
Vec4 PlanetOrientation(BodyId id, Real t_tdb) {
Real d = t_tdb / SECS_DAY; // Interval in days from J2000
Real T = d / 36525.0; // Interval in Julian Centries
Real alpha0, delta0, W, Wdot;
Real M1, M2, M3, M4, M5;
Real Ja, Jb, Jc, Jd, Je;
Real N, Ndot;
switch (id) {
case BodyId::MERCURY:
// Mercury
M1 = 174.7910857 + 4.092335 * d;
M2 = 349.5821714 + 8.184670 * d;
M3 = 164.3732571 + 12.277005 * d;
M4 = 339.1643429 + 16.369340 * d;
M5 = 153.9554286 + 20.461675 * d;
// Calculation
alpha0 = 281.0103 - 0.0328 * T;
delta0 = 61.4155 - 0.0049 * T;
W = 329.5988 + 6.1385108 * d + 0.01067257 * sin(RAD * M1) - 0.00112309 * sin(RAD * M2)
- 0.00011040 * sin(RAD * M3) - 0.00002539 * sin(RAD * M4) - 0.00000571 * sin(RAD * M5);
Wdot = 6.1385025 / SECS_DAY;
break;
case BodyId::VENUS:
alpha0 = 272.76;
delta0 = 67.16;
W = 160.20 - 1.4813688 * d;
Wdot = -1.4813688 / SECS_DAY;
break;
case BodyId::MARS:
alpha0 = 317.269202 - 0.10927547 * T;
alpha0 += 0.000068 * sin(RAD * (198.991226 + 19139.4819985 * T))
+ 0.000238 * sin(RAD * (226.292679 + 38280.8511281 * T))
+ 0.000052 * sin(RAD * (249.663391 + 57420.7251593 * T))
+ 0.000009 * sin(RAD * (266.183510 + 76560.6367950 * T))
+ 0.419057 * sin(RAD * (79.398797 + 0.5042615 * T));
delta0 = 54.432516 - 0.05827105 * T;
delta0 += 0.000051 * cos(RAD * (122.433576 + 19139.9407476 * T))
+ 0.000141 * cos(RAD * (43.058401 + 38280.8753272 * T))
+ 0.000031 * cos(RAD * (57.663379 + 57420.7517205 * T))
+ 0.000005 * cos(RAD * (79.476401 + 76560.6495004 * T))
+ 1.591274 * cos(RAD * (166.325722 + 0.5042615 * T));
W = 176.049863 + 350.891982443297 * d;
W += 0.000145 * sin(RAD * (129.071773 + 19140.0328244 * T))
+ 0.000157 * sin(RAD * (36.352167 + 38281.0473591 * T))
+ 0.000040 * sin(RAD * (56.668646 + 57420.9295360 * T))
+ 0.000001 * sin(RAD * (67.364003 + 76560.2552215 * T))
+ 0.000001 * sin(RAD * (104.792680 + 95700.4387578 * T))
+ 0.584542 * sin(RAD * (95.391654 + 0.5042615 * T));
Wdot = 350.89198226 / SECS_DAY;
break;
case BodyId::JUPITER:
Ja = 99.360714 + 4850.4046 * T;
Jb = 175.895369 + 1191.9605 * T;
Jc = 300.323162 + 262.5475 * T;
Jd = 114.012305 + 6070.2476 * T;
Je = 49.511251 + 64.3000 * T;
alpha0 = 268.056595 - 0.0064997 * T + 0.000117 * sin(RAD * Ja) + 0.000938 * sin(RAD * Jb);
delta0 = 64.495303 + 0.0024137 * T + 0.000050 * cos(RAD * Ja) + 0.000404 * cos(RAD * Jb)
+ 0.000617 * cos(RAD * Jc) - 0.000013 * cos(RAD * Jd) + 0.000926 * cos(RAD * Je);
W = 284.95 + 870.5360000 * d;
Wdot = 870.5366420 / SECS_DAY;
break;
case BodyId::SATURN:
alpha0 = 40.589 - 0.036 * T;
delta0 = 83.537 - 0.004 * T;
W = 38.90 + 810.7939024 * d;
Wdot = 810.793902 / SECS_DAY;
break;
case BodyId::URANUS:
alpha0 = 257.311;
delta0 = -15.175;
W = 203.81 - 501.1600928 * d;
Wdot = -501.1600928 / SECS_DAY;
break;
case BodyId::NEPTUNE:
N = 357.85 + 52.316 * T;
Ndot = 6.0551e-4;
alpha0 = 299.36 + 0.70 * sin(RAD * N);
delta0 = 43.46 - 0.51 * cos(RAD * N);
W = 249.978 + 541.1397757 * d - 0.48 * sin(RAD * N);
Wdot = 536.3128492 - 0.48 * Ndot * cos(RAD * N) / SECS_DAY;
break;
// otherwise throw error
default: throw std::invalid_argument("Invalid planet ID"); break;
}
// Return the three angles
Vec4 angles = {alpha0, delta0, W, Wdot};
// Convert to radians
for (int i = 0; i < 4; i++) {
angles[i] *= RAD;
}
return angles;
}
Mat3 RotPosBodyFixedToInertial(BodyId id, Real t_tdb) {
// Compute the planet orientation
Vec4 angles = PlanetOrientation(id, t_tdb);
Real alpha0 = angles[0];
Real delta0 = angles[1];
Real W = angles[2];
// Compute the rotation matrix
Mat3 Rr_BI = RotZ(alpha0 + PI / 2).transpose() * RotX(PI / 2 - delta0).transpose()
* RotZ(W).transpose();
return Rr_BI;
}
Mat3 RotPosInertialToBodyFixed(BodyId id, Real t_tdb) {
// Compute the planet orientation
Mat3 Rr_BI = RotPosBodyFixedToInertial(id, t_tdb);
Mat3 Rr_IB = Rr_BI.transpose();
return Rr_IB;
}
Mat6 RotPosVelBodyFixedToInertial(BodyId id, Real t_tdb) {
// Compute the planet orientation
Vec4 angles = PlanetOrientation(id, t_tdb);
Real alpha0 = angles[0];
Real delta0 = angles[1];
Real W = angles[2];
Real Wdot = angles[3];
Mat3 dWRotT;
dWRotT << -Wdot * sin(W), -Wdot * cos(W), 0, Wdot * cos(W), -Wdot * sin(W), 0, 0, 0, 0;
// Compute the rotation matrix
Mat3 Rr_BI = RotPosBodyFixedToInertial(id, t_tdb);
Mat3 Rv_BI = RotZ(alpha0 + PI / 2).transpose() * RotX(PI / 2 - delta0).transpose() * dWRotT;
Mat6 Rrv_BI;
// Blocks
Rrv_BI << Rr_BI, Mat3::Zero(), Rv_BI, Rr_BI;
return Rrv_BI;
}
Mat6 RotPosVelInertialToBodyFixed(BodyId id, Real t_tdb) {
// Compute the planet orientation
Vec4 angles = PlanetOrientation(id, t_tdb);
Real alpha0 = angles[0];
Real delta0 = angles[1];
Real W = angles[2];
Real Wdot = angles[3];
Mat3 dWRotT;
dWRotT << -Wdot * sin(W), -Wdot * cos(W), 0, Wdot * cos(W), -Wdot * sin(W), 0, 0, 0, 0;
// Compute the rotation matrix
Mat3 Rr_BI = RotPosBodyFixedToInertial(id, t_tdb);
Mat3 Rv_BI = RotZ(alpha0 + PI / 2).transpose() * RotX(PI / 2 - delta0).transpose() * dWRotT;
Mat3 Rr_IB = Rr_BI.transpose();
Mat3 Rv_IB = Rv_BI.transpose();
Mat6 Rrv_IB;
// Blocks
Rrv_IB << Rr_IB, Mat3::Zero(), Rv_IB, Rr_IB;
return Rrv_IB;
}
} // namespace lupnt