FrameTransformations.cpp

Go to the documentation of this file.

//
// Created by vformato on 6/27/23.
//
// Extracted from gbatch, file CC/FrameTrans.C
// Original author: S. Della Torre et al.
// Adapted by V. Formato
//
 
#include <array>
 
#include <fmt/format.h>
#include <fmt/ranges.h>
 
#include <Math/DisplacementVector3D.h>
#include <Math/Rotation3D.h>
#include <Math/RotationX.h>
#include <Math/RotationY.h>
#include <Math/RotationZ.h>
 
#include "Utility/FrameTransformation/FrameTransformations.h"
 
namespace NAIA {
namespace FrameTransformation {
 
double UTCToJulianDate(double utc_time) {
  // convert UTC time (given in unix time seconds) in JD day
 
  // The difference between TAI and UTC is
  // [INTERNATIONAL EARTH ROTATION AND REFERENCE SYSTEMS SERVICE (IERS) - Bulletin C 43]
  // from 2009 January 1, 0h UTC, to 2012 July 1  0h UTC  : UTC-TAI = - 34s
  // from 2012 July 1,    0h UTC, until further notice    : UTC-TAI = - 35s
 
  static constexpr double UTC1Ju2012 = 1341100800; /* unix time of 2012 July 1,    0h UTC*/
  double TAItime;
  // convert UTC Time in TAI time, this Difference is a function of the time
  if (utc_time >= 1483228800)
    TAItime = utc_time - 37;
  else if (utc_time >= 1435708800)
    TAItime = utc_time - 36;
  else if (utc_time >= UTC1Ju2012)
    TAItime = utc_time - 35;
  else
    TAItime = utc_time - 34;
 
  // JD(secconds)=TAI(seconds)+JD_1970*SecondsInADay
  double JDsec = TAItime + (2440587.5 * 86400);
  // convert seconds in days
  return JDsec / 86400;
}
 
double GMSTAngle(double utc_time) {
  // convert from unixtime (UTC) to Julian days
  double time = UTCToJulianDate(utc_time);
 
  // Greenwich mean sidereal time  in radians
 
  // Interval of time, measured in Julian centuries of 36525 days of UT (mean solar day), elapsed since the epoch 2000
  // Jan 1d12hUT
  double T = (time - 2451545.0) / 36525.0;
  double d = (time - 2451545.0);
  // greenwich mean sidereal time in degree (i.e. the Greenwich hour angle of the mean equinox of date)
 
  // GMST from Meeus, J., 2000. Astronomical Algorithms. Willman-Bell, Richmond,VA, 2nd ed. p 84 (eq.11-4) adapted from
  // IDL procedure http://idlastro.gsfc.nasa.gov/ftp/pro/astro/ct2lst.pro
  double GMST_deg = std::fmod((280.46061837 + 360.98564736629 * d + T * T * (0.000387933 - T / 38710000.0)), 360);
 
  double GMSTrad = GMST_deg * pi / 180.; /* greenwich mean sidereal time  in radians*/
  return GMSTrad;
}
 
//
// double degree_to_Rad(double angDeg){
// return angDeg/180.*pi;
// }
// double rad_to_degree(double angRad){
// return angRad*180./pi;
// }
//
// void FT_Cart2Angular(double x, double y, double z, double& r, double& theta, double& phi){
//   /* convert cartesian coordinates into spherical coordinates
//      theta is the latitude (i.e. the angle with the equator plane
//      phi is the aziutal angle*/
//   r=sqrt(x*x+y*y+z*z);
//   phi=atan2(y,x);
//   theta=pi/2.-acos(z/r);
// }
 
CartesianCoo<RefFrame::ECI> GTODToECIPos(CartesianCoo<RefFrame::GTOD> coo, double time) {
  // Get rotation angle for the Greenwich meridian in GTOD ref frame
  float gmt_angle = GMSTAngle(time);
 
  auto rotation = ROOT::Math::RotationZ(gmt_angle);
 
  return CartesianCoo<RefFrame::ECI>{rotation(coo).Coordinates()};
}
 
CartesianCoo<RefFrame::ECI> GTODToECIVel(CartesianCoo<RefFrame::GTOD> vel, CartesianCoo<RefFrame::GTOD> coo,
                                         double time) {
  // Get rotation angle for the Greenwich meridian in GTOD ref frame
  float gmt_angle = GMSTAngle(time);
 
  auto rotation = ROOT::Math::RotationZ(gmt_angle);
 
  CartesianCoo<RefFrame::ECI> result{rotation(vel).Coordinates()};
 
  // Add Earth rotation contribution
  static constexpr double omega = 2.0 * pi / 86400.0;
 
  result.SetCoordinates(result.X() * std::cos(gmt_angle) - result.Y() * std::sin(gmt_angle) +
                            omega * (-coo.Y() * std::cos(gmt_angle) - coo.X() * std::sin(gmt_angle)),
                        result.X() * std::sin(gmt_angle) + result.Y() * std::cos(gmt_angle) +
                            omega * (coo.X() * std::cos(gmt_angle) - coo.Y() * std::sin(gmt_angle)),
                        result.Z());
  return result;
}
 
CartesianCoo<RefFrame::ISSBody> AMSLocalToISSBody(CartesianCoo<RefFrame::AMSLocal> coo) {
  static constexpr double AMS_angle = 12.0001 * pi / 180.0;
 
  auto rotation = ROOT::Math::RotationY(AMS_angle);
 
  CartesianCoo<RefFrame::ISSBody> result{rotation(coo).Coordinates()};
  // TODO: check this transformation: how is ISSBody defined? And then can we express it as a rotation or set of
  //       rotations?
  result.SetCoordinates(-1 * result.Y(), -1 * result.X(), -1 * result.Z());
 
  return result;
}
 
CartesianCoo<RefFrame::ISSLVLH> ISSBodyToLVLH(CartesianCoo<RefFrame::ISSBody> coo, double ISSYaw, double ISSPitch,
                                              double ISSRoll) {
  // rotation along X axis
  auto roll_rotation = ROOT::Math::RotationX(ISSRoll);
  coo = roll_rotation(coo);
 
  // rotation along Y axis
  auto pitch_rotation = ROOT::Math::RotationY(ISSPitch);
  coo = pitch_rotation(coo);
 
  // rotation along Z axis
  auto yaw_rotation = ROOT::Math::RotationZ(ISSYaw);
  coo = yaw_rotation(coo);
 
  return CartesianCoo<RefFrame::ISSLVLH>{coo.Coordinates()};
}
 
CartesianCoo<RefFrame::ECI> ISSLVLHToECI(CartesianCoo<RefFrame::ISSLVLH> coo, CartesianCoo<RefFrame::ECI> ISSECIPos,
                                         CartesianCoo<RefFrame::ECI> ISSECIVel) {
  // 1st step: normalize ISS position and velocity
  ISSECIPos /= std::sqrt(ISSECIPos.mag2());
  ISSECIVel /= std::sqrt(ISSECIVel.mag2());
 
  // 2nd step: create the martix R_m elements
 
  // NOTE: from the original definition:
  //  the matrix of transformation is given by three lines X_v, Y_v, Z_v  of 3 elements each
  // | X_v |   | Y_v x Z_v |
  // | Y_v | = | V_v x R_v | = R_m
  // | Z_v |   | -R_v      |
 
  // NOTE: X_v, Y_v, Z_v vectors must be unitary
  auto Z_v = ISSECIPos * -1.0;
  auto Y_v = ISSECIVel.Cross(
      ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double>, RFrame<RefFrame::ECI>>{ISSECIPos});
  Y_v /= std::sqrt(Y_v.mag2());
  auto X_v = Y_v.Cross(ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double>, RFrame<RefFrame::ECI>>{Z_v});
  X_v /= std::sqrt(X_v.mag2());
 
  // NOTE: this constructor places the three vectors as the columns of the matrix. The original definition is given in
  // terms of the matrix rows so we don't need to transpose it.
  auto R_m = ROOT::Math::Rotation3D{X_v, Y_v, Z_v};
  TMatrixT<double> m{3, 3};
  R_m.GetRotationMatrix(m);
 
  return CartesianCoo<RefFrame::ECI>{R_m(coo).Coordinates()};
}
 
CartesianCoo<RefFrame::GTOD> ECIToGTOD(CartesianCoo<RefFrame::ECI> coo,
                                       double time) { // Get rotation angle for the Greenwich meridian in GTOD ref frame
  float gmt_angle = GMSTAngle(time);
 
  auto rotation = ROOT::Math::RotationZ(-1.0 * gmt_angle);
 
  return CartesianCoo<RefFrame::GTOD>{rotation(coo).Coordinates()};
}
 
} // namespace FrameTransformation
} // namespace NAIA