Conversion of GCRF to ITRF

Question

I have to convert the $[r,v,a]$ vector from GCRF frame to the ITRF frame (and vice versa), complying with the precession, nutation and polar motion.

How to calculate the rotation matrix? How to transform the velocity and acceleration?

---

• GCRF - Geocentric Celestial Reference Frame (inertial frame)
• ITRF - International Terrestrial Reference frame (Earth fixed)

Found this paper, but still reading to form a single formula.

Answer 1

Having recently gone through the same process, I thought I'd leave the procedure that I found relatively easy to implement here.

Section 5.9 of IERS 2010 technical notes describes 2 methods for performing the conversion from GCRF to ITRF, and state that both agree up to microarcseconds accuracy. I have chosen to implement method 2, which I summarize below (and certainly, as @DavidHammen pointed out in comments, using SOFA routines). The materials used in this implementation are SOFA routines available here and Earth orientation parameters, available for example at CelesTrak. Note also that the SOFA routines specify the time system into which the time should be input (e.g., TT, UT1, etc.). It is important to input it correctly. The time differences for a given date can be inferred from the linked CelesTrak file, which gives the UT1-UTC and TAI-UTC differences in seconds.

1. Calculate a rotation matrix to account for the effects of precession, nutation and frame-bias. Regarding frame-bias, my understanding is that this is to correct for a small misalignment of around 23 milliarcseconds between the mean pole and equinox of the J2000 frame and the GCRS (of which the GCRF is a realization). The procedure described in the IERS note gives two possibilities: to apply the NUT06A SOFA routine to obtain nutation components and then input them into PN06 routine to obtain the so-called nutation-precession-bias (NBP) matrix, which already gives a rotation to account for the 3 effects together; alternatively, to compute each separately with routines BI00 (frame bias), P06E (precession) and NUM06A (nutation), and to then multiply in the order Nutation x Precession x Frame Bias to obtain the NPB matrix. There is in fact a third way, which is to apply the PNM06A routine to obtain directly the NPB matrix. I recommend this, or NUT06A+PN06, to minimize errors. Also, note that you can add corrections to the nutation components before using them to compute the nutation part of the NPB matrix. You can find these, for example, in the CelesTrak file, in the $\delta\Delta\psi$ and $\delta\Delta\epsilon$ columns.
2. Calculate a rotation matrix to account for the rotation of Earth (this makes sense, since being at coordinates X, Y, Z with respect to the surface of Earth at a given time will translate into being at a different point in space depending on the time). You can achieve this by first calculating the Greenwich apparent sidereal time with routine GST06, and then converting it into a rotation matrix corresponding to a rotation along the Z-axis with routine RZ.
3. Calculate a polar motion matrix, the component that @DavidHammen mentioned that you might neglect if you don't require such high precision. If you require it, you will need to get X and Y coordinates for the rotation axis. You can get this from the CelesTrak file as well (columns named X and Y), and then apply routine POM00 to calculate the polar motion matrix. You will notice that this routine requires another variable, called $s'$. This is the so-called locator of the Terrestrial Intermediate Origin on the equator of the Celestial Intermediate Pole, and a decent estimate can be calculated with routine POM00 (this basically just considers the main, secular component. I wonder how much more accuracy getting more components for this parameter would help? But I guess if polar motion is already small, further components will be even more negligible. Still it would be interesting to know for what applications these would be necessary, if any).
4. Multiply all these matrixes together, and then multiply it by your position in GCRF to obtain the position in ITRF. If you want to perform the inverse conversion, transpose it before multiplication.

Note this is only for the position vector. If you want to transform also velocity vectors, you will need to take into account the angular velocity of Earth rotation and the distance to Earth, I guess. I am trying to implement this at the moment myself also.

Answer 2

I don't have the reputation to comment, but I wanted to leave a link to SOFA Tools for Earth Attitude. Section 5 walks through several formulations using the SOFA library. Here's an example implementation from section 5.5. I know this answer is a bit late, but I thought it would still be helpful for future googlers.

#include "sofa.h"

/**
 * Compute the celestial-to-terrestrial rotation matrix (ITRF to GCRF)
 * IAU 2006/2000A, CIO based, using classical angles
 * 
 * This method includes EOP corrections for both the terrestrial pole (Xp, Yp),
 * and celestial pole (dX06, dY06). 
 * 
 * Implements Section 5.5 of “SOFA Tools for Earth Attitude”
 * 
 * References:
 * [1] “SOFA Tools for Earth Attitude”, International Astronomical Union
 *     Standards of Fundamental Astronomy; International Astronomical
 *     Union (13 July 2017)
 * 
 * @param tta First part of two-part Julian date in TT timescale
 * @param ttb Second part of two-part Julian date in TT timescale
 * @param uta First part of two-part Julian date in UT1 timescale 
 * @param utb Second part of two-part Julian date in UT1 timescale
 * @param dx06 Celestial intermediate pole correction (from EOP parameters) [radians]
 * @param dy06 Celestial intermediate pole correction (from EOP parameters) [radians]
 * @param xp Coordinate of the pole (from EOP parameters) [radians]
 * @param yp Coordinate of the pole (from EOP parameters) [radians]
 *
 * @param rc2it (modified by routine) The celestial-terrestrial rotation matrix
 */
void rc2it2006CioClassicalAngles(double tta, double ttb, double uta, double utb, double dx06, double dy06, double xp, double yp, double rc2it[3][3])
{
    double x, y, s;
    double rc2i[3][3];
    double rc2ti[3][3];
    double rpom[3][3];

    // CIP and CIO, IAU 2006/2000A
    iauXys06a(tta, ttb, &x, &y, &s);

    // Add CIP correction
    x += dx06;
    y += dy06;

    // GCRS to CIRS matrix
    iauC2ixys(x, y, s, rc2i);

    // ERA
    auto era = iauEra00(uta, utb);

    // Form celestial-terrestrial matrix (no polar motion yet)
    iauCr(rc2i, rc2ti);
    iauRz(era, rc2ti);

    // Polar motion matrix (TIRS->ITRS, IERS 2003)
    auto sp = iauSp00(tta, ttb);
    iauPom00(xp, yp, sp, rpom);

    // Form celestial-terrestrial matrix (including polar motion)
    iauRxr(rpom, rc2ti, rc2it);
}