From b6541273c81d7b11da55a9fa6f6915d42972364f Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 27 Sep 2023 14:14:14 +0200 Subject: [PATCH 01/11] Fix argument name to suit the coding rule --- src/environment/global/celestial_rotation.cpp | 6 +++--- src/environment/global/celestial_rotation.hpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/environment/global/celestial_rotation.cpp b/src/environment/global/celestial_rotation.cpp index bb3e9f1d1..2ab19d920 100644 --- a/src/environment/global/celestial_rotation.cpp +++ b/src/environment/global/celestial_rotation.cpp @@ -131,12 +131,12 @@ void CelestialRotation::InitCelestialRotationAsEarth(const RotationMode rotation } } -void CelestialRotation::Update(const double JulianDate) { - double gmst_rad = gstime(JulianDate); // It is a bit different with 長沢(Nagasawa)'s algorithm. TODO: Check the correctness +void CelestialRotation::Update(const double julian_date) { + double gmst_rad = gstime(julian_date); // It is a bit different with 長沢(Nagasawa)'s algorithm. TODO: Check the correctness if (rotation_mode_ == RotationMode::kFull) { // Compute Julian date for terrestrial time - double jdTT_day = JulianDate + kDtUt1Utc_ * kSec2Day_; // TODO: Check the correctness. Problem is that S2E doesn't have Gregorian calendar. + double jdTT_day = julian_date + kDtUt1Utc_ * kSec2Day_; // TODO: Check the correctness. Problem is that S2E doesn't have Gregorian calendar. // Compute nth power of julian century for terrestrial time the actual unit of tTT_century is [century^(i+1)], i is the index of the array double tTT_century[4]; diff --git a/src/environment/global/celestial_rotation.hpp b/src/environment/global/celestial_rotation.hpp index e2444daba..b483c8a31 100644 --- a/src/environment/global/celestial_rotation.hpp +++ b/src/environment/global/celestial_rotation.hpp @@ -42,9 +42,9 @@ class CelestialRotation { /** * @fn Update * @brief Update rotation - * @param [in] JulianDate: Julian date + * @param [in] julian_date: Julian date */ - void Update(const double JulianDate); + void Update(const double julian_date); /** * @fn GetDcmJ2000ToXcxf From 444ad59e2da335c66507e6c1263e42ee74d70000 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 27 Sep 2023 14:23:35 +0200 Subject: [PATCH 02/11] Add doxygen comment --- src/environment/global/celestial_rotation.cpp | 37 ++++++++++--------- src/environment/global/celestial_rotation.hpp | 34 ++++++++++++++--- 2 files changed, 48 insertions(+), 23 deletions(-) diff --git a/src/environment/global/celestial_rotation.cpp b/src/environment/global/celestial_rotation.cpp index 2ab19d920..65d1c2f92 100644 --- a/src/environment/global/celestial_rotation.cpp +++ b/src/environment/global/celestial_rotation.cpp @@ -138,7 +138,8 @@ void CelestialRotation::Update(const double julian_date) { // Compute Julian date for terrestrial time double jdTT_day = julian_date + kDtUt1Utc_ * kSec2Day_; // TODO: Check the correctness. Problem is that S2E doesn't have Gregorian calendar. - // Compute nth power of julian century for terrestrial time the actual unit of tTT_century is [century^(i+1)], i is the index of the array + // Compute nth power of julian century for terrestrial time. + // The actual unit of tTT_century is [century^(i+1)], i is the index of the array double tTT_century[4]; tTT_century[0] = (jdTT_day - kJulianDateJ2000_) / kDayJulianCentury_; for (int i = 0; i < 3; i++) { @@ -175,40 +176,40 @@ void CelestialRotation::Update(const double julian_date) { } } -libra::Matrix<3, 3> CelestialRotation::AxialRotation(const double GAST_rad) { return libra::MakeRotationMatrixZ(GAST_rad); } +libra::Matrix<3, 3> CelestialRotation::AxialRotation(const double gast_rad) { return libra::MakeRotationMatrixZ(gast_rad); } -libra::Matrix<3, 3> CelestialRotation::Nutation(const double (&tTT_century)[4]) { +libra::Matrix<3, 3> CelestialRotation::Nutation(const double (&t_tt_century)[4]) { // Mean obliquity of the ecliptic epsilon_rad_ = c_epsilon_rad_[0]; for (int i = 0; i < 3; i++) { - epsilon_rad_ += c_epsilon_rad_[i + 1] * tTT_century[i]; + epsilon_rad_ += c_epsilon_rad_[i + 1] * t_tt_century[i]; } // Compute five delauney angles(l=lm,l'=ls,F,D,Ω=O) // Mean anomaly of the moon double lm_rad = c_lm_rad_[0]; for (int i = 0; i < 4; i++) { - lm_rad += c_lm_rad_[i + 1] * tTT_century[i]; + lm_rad += c_lm_rad_[i + 1] * t_tt_century[i]; } // Mean anomaly of the sun double ls_rad = c_ls_rad_[0]; for (int i = 0; i < 4; i++) { - ls_rad += c_ls_rad_[i + 1] * tTT_century[i]; + ls_rad += c_ls_rad_[i + 1] * t_tt_century[i]; } // Mean longitude of the moon - mean longitude of ascending node of the moon double F_rad = c_f_rad_[0]; for (int i = 0; i < 4; i++) { - F_rad += c_f_rad_[i + 1] * tTT_century[i]; + F_rad += c_f_rad_[i + 1] * t_tt_century[i]; } // Mean elogation of the moon from the sun double D_rad = c_d_rad_[0]; for (int i = 0; i < 4; i++) { - D_rad += c_d_rad_[i + 1] * tTT_century[i]; + D_rad += c_d_rad_[i + 1] * t_tt_century[i]; } // Mean longitude of ascending node of the moon double O_rad = c_o_rad_[0]; for (int i = 0; i < 4; i++) { - O_rad += c_o_rad_[i + 1] * tTT_century[i]; + O_rad += c_o_rad_[i + 1] * t_tt_century[i]; } // Additional angles @@ -239,19 +240,19 @@ libra::Matrix<3, 3> CelestialRotation::Nutation(const double (&tTT_century)[4]) return N; } -libra::Matrix<3, 3> CelestialRotation::Precession(const double (&tTT_century)[4]) { +libra::Matrix<3, 3> CelestialRotation::Precession(const double (&t_tt_century)[4]) { // Compute precession angles(zeta, theta, z) double zeta_rad = 0.0; for (int i = 0; i < 3; i++) { - zeta_rad += c_zeta_rad_[i] * tTT_century[i]; + zeta_rad += c_zeta_rad_[i] * t_tt_century[i]; } double theta_rad = 0.0; for (int i = 0; i < 3; i++) { - theta_rad += c_theta_rad_[i] * tTT_century[i]; + theta_rad += c_theta_rad_[i] * t_tt_century[i]; } double z_rad = 0.0; for (int i = 0; i < 3; i++) { - z_rad += c_z_rad_[i] * tTT_century[i]; + z_rad += c_z_rad_[i] * t_tt_century[i]; } // Develop transformation matrix @@ -265,17 +266,17 @@ libra::Matrix<3, 3> CelestialRotation::Precession(const double (&tTT_century)[4] return P; } -libra::Matrix<3, 3> CelestialRotation::PolarMotion(const double Xp, const double Yp) { +libra::Matrix<3, 3> CelestialRotation::PolarMotion(const double x_p, const double y_p) { libra::Matrix<3, 3> W; W[0][0] = 1.0; W[0][1] = 0.0; - W[0][2] = -Xp; + W[0][2] = -x_p; W[1][0] = 0.0; W[1][1] = 1.0; - W[1][2] = -Yp; - W[2][0] = Xp; - W[2][1] = Yp; + W[1][2] = -y_p; + W[2][0] = x_p; + W[2][1] = y_p; W[2][2] = 1.0; return W; diff --git a/src/environment/global/celestial_rotation.hpp b/src/environment/global/celestial_rotation.hpp index b483c8a31..008f58bdf 100644 --- a/src/environment/global/celestial_rotation.hpp +++ b/src/environment/global/celestial_rotation.hpp @@ -98,12 +98,36 @@ class CelestialRotation { */ void InitCelestialRotationAsEarth(const RotationMode rotation_mode, const std::string center_body_name); - // TODO: Add doxygen comments for the private functions and fix argument name + /** + * @fn AxialRotation + * @brief Calculate movement of the coordinate axes due to rotation around the rotation axis + * @param [in] gast_rad: Greenwitch 'Apparent' Sidereal Time [rad] + * @return Rotation matrix + */ + libra::Matrix<3, 3> AxialRotation(const double gast_rad); + + /** + * @fn Nutation + * @brief Calculate movement of the coordinate axes due to Nutation + * @param [in] t_tt_century: nth power of julian century for terrestrial time + * @return Rotation matrix + */ + libra::Matrix<3, 3> Nutation(const double (&t_tt_century)[4]); - libra::Matrix<3, 3> AxialRotation(const double GAST_rad); //!< Movement of the coordinate axes due to rotation around the rotation axis - libra::Matrix<3, 3> Nutation(const double (&tTT_century)[4]); //!< Movement of the coordinate axes due to Nutation - libra::Matrix<3, 3> Precession(const double (&tTT_century)[4]); //!< Movement of the coordinate axes due to Precession - libra::Matrix<3, 3> PolarMotion(const double Xp, const double Yp); //!< Movement of the coordinate axes due to Polar Motion + /** + * @fn Precession + * @brief Calculate movement of the coordinate axes due to Precession + * @param [in] t_tt_century: nth power of julian century for terrestrial time + * @return Rotation matrix + */ + libra::Matrix<3, 3> Precession(const double (&t_tt_century)[4]); + + /** + * @fn PolarMotion + * @brief Calculate movement of the coordinate axes due to Polar Motion + * @note Currently, this function is not used. + */ + libra::Matrix<3, 3> PolarMotion(const double x_p, const double y_p); }; #endif // S2E_ENVIRONMENT_GLOBAL_CELESTIAL_ROTATION_HPP_ From 25d78166a6612a18343ff17b61376cd4aedaaa4a Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 27 Sep 2023 14:32:24 +0200 Subject: [PATCH 03/11] Fix comment --- src/environment/global/celestial_rotation.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/environment/global/celestial_rotation.cpp b/src/environment/global/celestial_rotation.cpp index 65d1c2f92..276290125 100644 --- a/src/environment/global/celestial_rotation.cpp +++ b/src/environment/global/celestial_rotation.cpp @@ -152,8 +152,7 @@ void CelestialRotation::Update(const double julian_date) { libra::Matrix<3, 3> W; // Nutation + Precession P = Precession(tTT_century); - N = Nutation(tTT_century); // epsilon_rad_, d_epsilon_rad_, d_psi_rad_ are - // updated in this procedure + N = Nutation(tTT_century); // epsilon_rad_, d_epsilon_rad_, d_psi_rad_ are updated in this procedure // Axial Rotation double Eq_rad = d_psi_rad_ * cos(epsilon_rad_ + d_epsilon_rad_); // Equation of equinoxes [rad] @@ -185,7 +184,7 @@ libra::Matrix<3, 3> CelestialRotation::Nutation(const double (&t_tt_century)[4]) epsilon_rad_ += c_epsilon_rad_[i + 1] * t_tt_century[i]; } - // Compute five delauney angles(l=lm,l'=ls,F,D,Ω=O) + // Compute five delauney angles(l=lm, l'=ls, F, D, Ω=O) // Mean anomaly of the moon double lm_rad = c_lm_rad_[0]; for (int i = 0; i < 4; i++) { @@ -213,21 +212,21 @@ libra::Matrix<3, 3> CelestialRotation::Nutation(const double (&t_tt_century)[4]) } // Additional angles - double L_rad = F_rad + O_rad; // F + Ω [rad] - double Ld_rad = L_rad - D_rad; // F - D + Ω [rad] + double L_rad = F_rad + O_rad; // F + Ω + double Ld_rad = L_rad - D_rad; // F + Ω - D // Compute luni-solar nutation // Nutation in obliquity d_psi_rad_ = c_d_psi_rad_[0] * sin(O_rad) + c_d_psi_rad_[1] * sin(2 * Ld_rad) + c_d_psi_rad_[2] * sin(2 * O_rad) + - c_d_psi_rad_[3] * sin(2 * L_rad) + c_d_psi_rad_[4] * sin(ls_rad); // [rad] + c_d_psi_rad_[3] * sin(2 * L_rad) + c_d_psi_rad_[4] * sin(ls_rad); d_psi_rad_ = d_psi_rad_ + c_d_psi_rad_[5] * sin(lm_rad) + c_d_psi_rad_[6] * sin(2 * Ld_rad + ls_rad) + c_d_psi_rad_[7] * sin(2 * L_rad + lm_rad) + - c_d_psi_rad_[8] * sin(2 * Ld_rad - ls_rad); // [rad] + c_d_psi_rad_[8] * sin(2 * Ld_rad - ls_rad); // Nutation in longitude d_epsilon_rad_ = c_d_epsilon_rad_[0] * cos(O_rad) + c_d_epsilon_rad_[1] * cos(2 * Ld_rad) + c_d_epsilon_rad_[2] * cos(2 * O_rad) + - c_d_epsilon_rad_[3] * cos(2 * L_rad) + c_d_epsilon_rad_[4] * cos(ls_rad); // [rad] + c_d_epsilon_rad_[3] * cos(2 * L_rad) + c_d_epsilon_rad_[4] * cos(ls_rad); d_epsilon_rad_ = d_epsilon_rad_ + c_d_epsilon_rad_[5] * cos(lm_rad) + c_d_epsilon_rad_[6] * cos(2 * Ld_rad + ls_rad) + - c_d_epsilon_rad_[7] * cos(2 * L_rad + lm_rad) + c_d_epsilon_rad_[8] * cos(2 * Ld_rad - ls_rad); // [rad] + c_d_epsilon_rad_[7] * cos(2 * L_rad + lm_rad) + c_d_epsilon_rad_[8] * cos(2 * Ld_rad - ls_rad); double epsi_mod_rad = epsilon_rad_ + d_epsilon_rad_; libra::Matrix<3, 3> X_epsi_1st = libra::MakeRotationMatrixX(epsilon_rad_); From 521980b9746b7808b60968e8b382a46248307c5a Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 27 Sep 2023 15:10:42 +0200 Subject: [PATCH 04/11] Fix small --- src/environment/global/celestial_rotation.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/environment/global/celestial_rotation.cpp b/src/environment/global/celestial_rotation.cpp index 276290125..2ad43878d 100644 --- a/src/environment/global/celestial_rotation.cpp +++ b/src/environment/global/celestial_rotation.cpp @@ -230,11 +230,11 @@ libra::Matrix<3, 3> CelestialRotation::Nutation(const double (&t_tt_century)[4]) double epsi_mod_rad = epsilon_rad_ + d_epsilon_rad_; libra::Matrix<3, 3> X_epsi_1st = libra::MakeRotationMatrixX(epsilon_rad_); - libra::Matrix<3, 3> Z_dpsi = libra::MakeRotationMatrixZ(-d_psi_rad_); + libra::Matrix<3, 3> Z_d_psi = libra::MakeRotationMatrixZ(-d_psi_rad_); libra::Matrix<3, 3> X_epsi_2nd = libra::MakeRotationMatrixX(-epsi_mod_rad); libra::Matrix<3, 3> N; - N = X_epsi_2nd * Z_dpsi * X_epsi_1st; + N = X_epsi_2nd * Z_d_psi * X_epsi_1st; return N; } From 823771e6de0492efc98ed38ceb7f23bcb2ce8884 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Fri, 29 Sep 2023 14:01:33 +0200 Subject: [PATCH 05/11] Fix typo --- src/environment/global/celestial_rotation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/environment/global/celestial_rotation.hpp b/src/environment/global/celestial_rotation.hpp index 008f58bdf..4bf4c10a0 100644 --- a/src/environment/global/celestial_rotation.hpp +++ b/src/environment/global/celestial_rotation.hpp @@ -101,7 +101,7 @@ class CelestialRotation { /** * @fn AxialRotation * @brief Calculate movement of the coordinate axes due to rotation around the rotation axis - * @param [in] gast_rad: Greenwitch 'Apparent' Sidereal Time [rad] + * @param [in] gast_rad: Greenwich 'Apparent' Sidereal Time [rad] * @return Rotation matrix */ libra::Matrix<3, 3> AxialRotation(const double gast_rad); From 62067d2433ff74d7c1b03cbf3719d85ad11a4579 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Fri, 29 Sep 2023 14:02:43 +0200 Subject: [PATCH 06/11] Fix typo --- src/environment/global/celestial_rotation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/environment/global/celestial_rotation.cpp b/src/environment/global/celestial_rotation.cpp index 2ad43878d..858c4365f 100644 --- a/src/environment/global/celestial_rotation.cpp +++ b/src/environment/global/celestial_rotation.cpp @@ -156,7 +156,7 @@ void CelestialRotation::Update(const double julian_date) { // Axial Rotation double Eq_rad = d_psi_rad_ * cos(epsilon_rad_ + d_epsilon_rad_); // Equation of equinoxes [rad] - double gast_rad = gmst_rad + Eq_rad; // Greenwitch 'Apparent' Sidereal Time [rad] + double gast_rad = gmst_rad + Eq_rad; // Greenwich 'Apparent' Sidereal Time [rad] R = AxialRotation(gast_rad); // Polar motion (is not considered so far, even without polar motion, the result agrees well with the matlab reference) double Xp = 0.0; From 3d0d4f3693180951a7c6d41d3522eb3e8bd80afb Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Fri, 29 Sep 2023 14:11:55 +0200 Subject: [PATCH 07/11] Fix typo --- src/environment/global/celestial_rotation.cpp | 6 +++--- src/environment/global/celestial_rotation.hpp | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/environment/global/celestial_rotation.cpp b/src/environment/global/celestial_rotation.cpp index 858c4365f..c9664f803 100644 --- a/src/environment/global/celestial_rotation.cpp +++ b/src/environment/global/celestial_rotation.cpp @@ -53,7 +53,7 @@ void CelestialRotation::InitCelestialRotationAsEarth(const RotationMode rotation c_epsilon_rad_[2] = -5.9000e-4 * libra::arcsec_to_rad; // [rad/century^2] c_epsilon_rad_[3] = 1.8130e-3 * libra::arcsec_to_rad; // [rad/century^3] - // Coefficients to compute delauney angles + // Coefficients to compute Delaunay angles // The actual unit of the coefficients are [rad/century^i], where i is the index of the array c_lm_rad_[0] = 134.96340251 * libra::deg_to_rad; // [rad] c_lm_rad_[1] = 1717915923.21780000 * libra::arcsec_to_rad; // [rad/century] @@ -184,7 +184,7 @@ libra::Matrix<3, 3> CelestialRotation::Nutation(const double (&t_tt_century)[4]) epsilon_rad_ += c_epsilon_rad_[i + 1] * t_tt_century[i]; } - // Compute five delauney angles(l=lm, l'=ls, F, D, Ω=O) + // Compute five Delaunay angles(l=lm, l'=ls, F, D, Ω=O) // Mean anomaly of the moon double lm_rad = c_lm_rad_[0]; for (int i = 0; i < 4; i++) { @@ -200,7 +200,7 @@ libra::Matrix<3, 3> CelestialRotation::Nutation(const double (&t_tt_century)[4]) for (int i = 0; i < 4; i++) { F_rad += c_f_rad_[i + 1] * t_tt_century[i]; } - // Mean elogation of the moon from the sun + // Mean elongation of the moon from the sun double D_rad = c_d_rad_[0]; for (int i = 0; i < 4; i++) { D_rad += c_d_rad_[i + 1] * t_tt_century[i]; diff --git a/src/environment/global/celestial_rotation.hpp b/src/environment/global/celestial_rotation.hpp index 4bf4c10a0..23e4402ae 100644 --- a/src/environment/global/celestial_rotation.hpp +++ b/src/environment/global/celestial_rotation.hpp @@ -72,11 +72,11 @@ class CelestialRotation { // TODO: Consider to read setting files for these coefficients // TODO: Consider other formats for other planets double c_epsilon_rad_[4]; //!< Coefficients to compute mean obliquity of the ecliptic - double c_lm_rad_[5]; //!< Coefficients to compute delauney angle (l=lm: Mean anomaly of the moon) - double c_ls_rad_[5]; //!< Coefficients to compute delauney angle (l'=ls: Mean anomaly of the sun) - double c_f_rad_[5]; //!< Coefficients to compute delauney angle (F: Mean longitude of the moon - mean longitude of ascending node of the moon) - double c_d_rad_[5]; //!< Coefficients to compute delauney angle (D: Elogation of the moon from the sun) - double c_o_rad_[5]; //!< Coefficients to compute delauney angle (Ω=O: Mean longitude of ascending node of the moon) + double c_lm_rad_[5]; //!< Coefficients to compute Delaunay angle (l=lm: Mean anomaly of the moon) + double c_ls_rad_[5]; //!< Coefficients to compute Delaunay angle (l'=ls: Mean anomaly of the sun) + double c_f_rad_[5]; //!< Coefficients to compute Delaunay angle (F: Mean longitude of the moon - mean longitude of ascending node of the moon) + double c_d_rad_[5]; //!< Coefficients to compute Delaunay angle (D: Elogation of the moon from the sun) + double c_o_rad_[5]; //!< Coefficients to compute Delaunay angle (Ω=O: Mean longitude of ascending node of the moon) double c_d_epsilon_rad_[9]; //!< Coefficients to compute nutation angle (delta-epsilon) double c_d_psi_rad_[9]; //!< Coefficients to compute nutation angle (delta-psi) double c_zeta_rad_[3]; //!< Coefficients to compute precession angle (zeta) From 47a2662ac7723e1c5de4b961108dd093aeb1da3d Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Fri, 29 Sep 2023 14:13:18 +0200 Subject: [PATCH 08/11] Fix typo --- src/environment/global/celestial_rotation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/environment/global/celestial_rotation.hpp b/src/environment/global/celestial_rotation.hpp index 23e4402ae..1f9fb1e75 100644 --- a/src/environment/global/celestial_rotation.hpp +++ b/src/environment/global/celestial_rotation.hpp @@ -75,7 +75,7 @@ class CelestialRotation { double c_lm_rad_[5]; //!< Coefficients to compute Delaunay angle (l=lm: Mean anomaly of the moon) double c_ls_rad_[5]; //!< Coefficients to compute Delaunay angle (l'=ls: Mean anomaly of the sun) double c_f_rad_[5]; //!< Coefficients to compute Delaunay angle (F: Mean longitude of the moon - mean longitude of ascending node of the moon) - double c_d_rad_[5]; //!< Coefficients to compute Delaunay angle (D: Elogation of the moon from the sun) + double c_d_rad_[5]; //!< Coefficients to compute Delaunay angle (D: Elongation of the moon from the sun) double c_o_rad_[5]; //!< Coefficients to compute Delaunay angle (Ω=O: Mean longitude of ascending node of the moon) double c_d_epsilon_rad_[9]; //!< Coefficients to compute nutation angle (delta-epsilon) double c_d_psi_rad_[9]; //!< Coefficients to compute nutation angle (delta-psi) From 0a3612cc8465d542373a4d9c67f7eb790cb2e867 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Sun, 1 Oct 2023 15:18:10 +0200 Subject: [PATCH 09/11] Fix local variables name --- src/environment/global/celestial_rotation.cpp | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/environment/global/celestial_rotation.cpp b/src/environment/global/celestial_rotation.cpp index c9664f803..73233375f 100644 --- a/src/environment/global/celestial_rotation.cpp +++ b/src/environment/global/celestial_rotation.cpp @@ -146,25 +146,25 @@ void CelestialRotation::Update(const double julian_date) { tTT_century[i + 1] = tTT_century[i] * tTT_century[0]; } - libra::Matrix<3, 3> P; - libra::Matrix<3, 3> N; - libra::Matrix<3, 3> R; - libra::Matrix<3, 3> W; + libra::Matrix<3, 3> dcm_precession; + libra::Matrix<3, 3> dcm_nutation; + libra::Matrix<3, 3> dcm_rotation; + libra::Matrix<3, 3> dcm_polar_motion; // Nutation + Precession - P = Precession(tTT_century); - N = Nutation(tTT_century); // epsilon_rad_, d_epsilon_rad_, d_psi_rad_ are updated in this procedure + dcm_precession = Precession(tTT_century); + dcm_nutation = Nutation(tTT_century); // epsilon_rad_, d_epsilon_rad_, d_psi_rad_ are updated in this procedure // Axial Rotation - double Eq_rad = d_psi_rad_ * cos(epsilon_rad_ + d_epsilon_rad_); // Equation of equinoxes [rad] - double gast_rad = gmst_rad + Eq_rad; // Greenwich 'Apparent' Sidereal Time [rad] - R = AxialRotation(gast_rad); + double equinox_rad = d_psi_rad_ * cos(epsilon_rad_ + d_epsilon_rad_); // Equation of equinoxes [rad] + double gast_rad = gmst_rad + equinox_rad; // Greenwich 'Apparent' Sidereal Time [rad] + dcm_rotation = AxialRotation(gast_rad); // Polar motion (is not considered so far, even without polar motion, the result agrees well with the matlab reference) - double Xp = 0.0; - double Yp = 0.0; - W = PolarMotion(Xp, Yp); + double x_p = 0.0; + double y_p = 0.0; + dcm_polar_motion = PolarMotion(x_p, y_p); // Total orientation - dcm_j2000_to_xcxf_ = W * R * N * P; + dcm_j2000_to_xcxf_ = dcm_polar_motion * dcm_rotation * dcm_nutation * dcm_precession; } else if (rotation_mode_ == RotationMode::kSimple) { // In this case, only Axial Rotation is executed, with its argument replaced from G'A'ST to G'M'ST // FIXME: Not suitable when the center body is not the earth From 4d2c3143e933fbd9a22993637f5574136f51b54f Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Sun, 1 Oct 2023 15:25:10 +0200 Subject: [PATCH 10/11] Fix local variables name --- src/environment/global/celestial_rotation.cpp | 92 +++++++++---------- 1 file changed, 46 insertions(+), 46 deletions(-) diff --git a/src/environment/global/celestial_rotation.cpp b/src/environment/global/celestial_rotation.cpp index 73233375f..ca1bb4601 100644 --- a/src/environment/global/celestial_rotation.cpp +++ b/src/environment/global/celestial_rotation.cpp @@ -140,10 +140,10 @@ void CelestialRotation::Update(const double julian_date) { // Compute nth power of julian century for terrestrial time. // The actual unit of tTT_century is [century^(i+1)], i is the index of the array - double tTT_century[4]; - tTT_century[0] = (jdTT_day - kJulianDateJ2000_) / kDayJulianCentury_; + double terrestrial_time_julian_century[4]; + terrestrial_time_julian_century[0] = (jdTT_day - kJulianDateJ2000_) / kDayJulianCentury_; for (int i = 0; i < 3; i++) { - tTT_century[i + 1] = tTT_century[i] * tTT_century[0]; + terrestrial_time_julian_century[i + 1] = terrestrial_time_julian_century[i] * terrestrial_time_julian_century[0]; } libra::Matrix<3, 3> dcm_precession; @@ -151,8 +151,8 @@ void CelestialRotation::Update(const double julian_date) { libra::Matrix<3, 3> dcm_rotation; libra::Matrix<3, 3> dcm_polar_motion; // Nutation + Precession - dcm_precession = Precession(tTT_century); - dcm_nutation = Nutation(tTT_century); // epsilon_rad_, d_epsilon_rad_, d_psi_rad_ are updated in this procedure + dcm_precession = Precession(terrestrial_time_julian_century); + dcm_nutation = Nutation(terrestrial_time_julian_century); // epsilon_rad_, d_epsilon_rad_, d_psi_rad_ are updated in this procedure // Axial Rotation double equinox_rad = d_psi_rad_ * cos(epsilon_rad_ + d_epsilon_rad_); // Equation of equinoxes [rad] @@ -196,47 +196,47 @@ libra::Matrix<3, 3> CelestialRotation::Nutation(const double (&t_tt_century)[4]) ls_rad += c_ls_rad_[i + 1] * t_tt_century[i]; } // Mean longitude of the moon - mean longitude of ascending node of the moon - double F_rad = c_f_rad_[0]; + double f_rad = c_f_rad_[0]; for (int i = 0; i < 4; i++) { - F_rad += c_f_rad_[i + 1] * t_tt_century[i]; + f_rad += c_f_rad_[i + 1] * t_tt_century[i]; } // Mean elongation of the moon from the sun - double D_rad = c_d_rad_[0]; + double d_rad = c_d_rad_[0]; for (int i = 0; i < 4; i++) { - D_rad += c_d_rad_[i + 1] * t_tt_century[i]; + d_rad += c_d_rad_[i + 1] * t_tt_century[i]; } // Mean longitude of ascending node of the moon - double O_rad = c_o_rad_[0]; + double o_rad = c_o_rad_[0]; for (int i = 0; i < 4; i++) { - O_rad += c_o_rad_[i + 1] * t_tt_century[i]; + o_rad += c_o_rad_[i + 1] * t_tt_century[i]; } // Additional angles - double L_rad = F_rad + O_rad; // F + Ω - double Ld_rad = L_rad - D_rad; // F + Ω - D + double l_rad = f_rad + o_rad; // F + Ω + double ld_rad = l_rad - d_rad; // F + Ω - D // Compute luni-solar nutation // Nutation in obliquity - d_psi_rad_ = c_d_psi_rad_[0] * sin(O_rad) + c_d_psi_rad_[1] * sin(2 * Ld_rad) + c_d_psi_rad_[2] * sin(2 * O_rad) + - c_d_psi_rad_[3] * sin(2 * L_rad) + c_d_psi_rad_[4] * sin(ls_rad); - d_psi_rad_ = d_psi_rad_ + c_d_psi_rad_[5] * sin(lm_rad) + c_d_psi_rad_[6] * sin(2 * Ld_rad + ls_rad) + c_d_psi_rad_[7] * sin(2 * L_rad + lm_rad) + - c_d_psi_rad_[8] * sin(2 * Ld_rad - ls_rad); + d_psi_rad_ = c_d_psi_rad_[0] * sin(o_rad) + c_d_psi_rad_[1] * sin(2 * ld_rad) + c_d_psi_rad_[2] * sin(2 * o_rad) + + c_d_psi_rad_[3] * sin(2 * l_rad) + c_d_psi_rad_[4] * sin(ls_rad); + d_psi_rad_ = d_psi_rad_ + c_d_psi_rad_[5] * sin(lm_rad) + c_d_psi_rad_[6] * sin(2 * ld_rad + ls_rad) + c_d_psi_rad_[7] * sin(2 * l_rad + lm_rad) + + c_d_psi_rad_[8] * sin(2 * ld_rad - ls_rad); // Nutation in longitude - d_epsilon_rad_ = c_d_epsilon_rad_[0] * cos(O_rad) + c_d_epsilon_rad_[1] * cos(2 * Ld_rad) + c_d_epsilon_rad_[2] * cos(2 * O_rad) + - c_d_epsilon_rad_[3] * cos(2 * L_rad) + c_d_epsilon_rad_[4] * cos(ls_rad); - d_epsilon_rad_ = d_epsilon_rad_ + c_d_epsilon_rad_[5] * cos(lm_rad) + c_d_epsilon_rad_[6] * cos(2 * Ld_rad + ls_rad) + - c_d_epsilon_rad_[7] * cos(2 * L_rad + lm_rad) + c_d_epsilon_rad_[8] * cos(2 * Ld_rad - ls_rad); + d_epsilon_rad_ = c_d_epsilon_rad_[0] * cos(o_rad) + c_d_epsilon_rad_[1] * cos(2 * ld_rad) + c_d_epsilon_rad_[2] * cos(2 * o_rad) + + c_d_epsilon_rad_[3] * cos(2 * l_rad) + c_d_epsilon_rad_[4] * cos(ls_rad); + d_epsilon_rad_ = d_epsilon_rad_ + c_d_epsilon_rad_[5] * cos(lm_rad) + c_d_epsilon_rad_[6] * cos(2 * ld_rad + ls_rad) + + c_d_epsilon_rad_[7] * cos(2 * l_rad + lm_rad) + c_d_epsilon_rad_[8] * cos(2 * ld_rad - ls_rad); double epsi_mod_rad = epsilon_rad_ + d_epsilon_rad_; - libra::Matrix<3, 3> X_epsi_1st = libra::MakeRotationMatrixX(epsilon_rad_); - libra::Matrix<3, 3> Z_d_psi = libra::MakeRotationMatrixZ(-d_psi_rad_); - libra::Matrix<3, 3> X_epsi_2nd = libra::MakeRotationMatrixX(-epsi_mod_rad); + libra::Matrix<3, 3> x_epsi_1st = libra::MakeRotationMatrixX(epsilon_rad_); + libra::Matrix<3, 3> z_d_psi = libra::MakeRotationMatrixZ(-d_psi_rad_); + libra::Matrix<3, 3> x_epsi_2nd = libra::MakeRotationMatrixX(-epsi_mod_rad); - libra::Matrix<3, 3> N; - N = X_epsi_2nd * Z_d_psi * X_epsi_1st; + libra::Matrix<3, 3> dcm_nutation; + dcm_nutation = x_epsi_2nd * z_d_psi * x_epsi_1st; - return N; + return dcm_nutation; } libra::Matrix<3, 3> CelestialRotation::Precession(const double (&t_tt_century)[4]) { @@ -255,28 +255,28 @@ libra::Matrix<3, 3> CelestialRotation::Precession(const double (&t_tt_century)[4 } // Develop transformation matrix - libra::Matrix<3, 3> Z_zeta = libra::MakeRotationMatrixZ(-zeta_rad); - libra::Matrix<3, 3> Y_theta = libra::MakeRotationMatrixY(theta_rad); - libra::Matrix<3, 3> Z_z = libra::MakeRotationMatrixZ(-z_rad); + libra::Matrix<3, 3> z_zeta = libra::MakeRotationMatrixZ(-zeta_rad); + libra::Matrix<3, 3> y_theta = libra::MakeRotationMatrixY(theta_rad); + libra::Matrix<3, 3> z_z = libra::MakeRotationMatrixZ(-z_rad); - libra::Matrix<3, 3> P; - P = Z_z * Y_theta * Z_zeta; + libra::Matrix<3, 3> dcm_precession; + dcm_precession = z_z * y_theta * z_zeta; - return P; + return dcm_precession; } libra::Matrix<3, 3> CelestialRotation::PolarMotion(const double x_p, const double y_p) { - libra::Matrix<3, 3> W; - - W[0][0] = 1.0; - W[0][1] = 0.0; - W[0][2] = -x_p; - W[1][0] = 0.0; - W[1][1] = 1.0; - W[1][2] = -y_p; - W[2][0] = x_p; - W[2][1] = y_p; - W[2][2] = 1.0; - - return W; + libra::Matrix<3, 3> dcm_polar_motion; + + dcm_polar_motion[0][0] = 1.0; + dcm_polar_motion[0][1] = 0.0; + dcm_polar_motion[0][2] = -x_p; + dcm_polar_motion[1][0] = 0.0; + dcm_polar_motion[1][1] = 1.0; + dcm_polar_motion[1][2] = -y_p; + dcm_polar_motion[2][0] = x_p; + dcm_polar_motion[2][1] = y_p; + dcm_polar_motion[2][2] = 1.0; + + return dcm_polar_motion; } From 8604ccf59909f4e8a81a70463a8cb293dbb0b2cb Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Sun, 1 Oct 2023 15:33:00 +0200 Subject: [PATCH 11/11] Fix local variables name --- src/environment/global/celestial_rotation.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/environment/global/celestial_rotation.cpp b/src/environment/global/celestial_rotation.cpp index ca1bb4601..811f31fe5 100644 --- a/src/environment/global/celestial_rotation.cpp +++ b/src/environment/global/celestial_rotation.cpp @@ -136,12 +136,12 @@ void CelestialRotation::Update(const double julian_date) { if (rotation_mode_ == RotationMode::kFull) { // Compute Julian date for terrestrial time - double jdTT_day = julian_date + kDtUt1Utc_ * kSec2Day_; // TODO: Check the correctness. Problem is that S2E doesn't have Gregorian calendar. + double terrestrial_time_julian_day = julian_date + kDtUt1Utc_ * kSec2Day_; // TODO: Check the correctness. Problem is that S2E doesn't have Gregorian calendar. // Compute nth power of julian century for terrestrial time. // The actual unit of tTT_century is [century^(i+1)], i is the index of the array double terrestrial_time_julian_century[4]; - terrestrial_time_julian_century[0] = (jdTT_day - kJulianDateJ2000_) / kDayJulianCentury_; + terrestrial_time_julian_century[0] = (terrestrial_time_julian_day - kJulianDateJ2000_) / kDayJulianCentury_; for (int i = 0; i < 3; i++) { terrestrial_time_julian_century[i + 1] = terrestrial_time_julian_century[i] * terrestrial_time_julian_century[0]; }