I'm writing a program that converts position in the ECI (CIS, EPOCH JD2000.0) coordinate system to WGS 84 (CTS, ECEF) coordinates. To do so, I'm following what's described in this document and implementing C++ code.

My knowledge of astrophysics is very limited, so I find it hard to understand some terms. Especially, I can't seem to understand properly how the Greenwich Mean Sideral Time is computed, and why I always get different values.

For instance (assuming time is 2019-01-01 08:00:00 and Julian Date is 2458484.833333):

  • If I use this on-line GMST calculator, I find that GMST is 14:42:45.3766.
  • If I use the expression used in this MATLAB script (showing the value in time, not degrees) I also get 14:42:45.350.
  • However, if I calculate it myself with the expression given in the above appendix, I get 06:42:45.379163 (with a clear offset of

According to the equations in page A-26, GMST depends upon the number of Julian centuries (Tu). H0 expression is exactly the same than that of CelesTrak for θ(0h). However, I can't really understand what the du truly stands for, or more specifically, why the document shows du → ±0.5, 1.5, 2.5.... Moreover, to find Λ (which is stated to be in seconds of time) the expressions use values in radians and seconds of time indistinctly.

Could somebody shed light into this?

up vote 3 down vote accepted

To do so, I'm following what's described in this document and implementing C++ code.

That is a very old document you are using. It's using the 1982 IAU precession model. There have been multiple updates since then. The latest set of changes are working toward eliminating the concept of Greenwich Sidereal Time. (In software terms, GMST and GAST are deprecated concepts. While still currently supported, support will be dropped in some future release of the Earth rotation model.)


Other than using an outdated version of a deprecated concept, what did you do wrong? The expression you are using to compute GMST is at the top of page A-26: $$\Lambda = H_0 + \Delta H + \omega^{\ast}\,(t-\Delta t)$$ The $\Delta H$ term takes one from GMST to GAST. Dropping that term, the expression for GMST is $$H = H_0 + \omega^{\ast}\,(t-\Delta t)\tag1$$ where $$\begin{align} H_0 = &\,24110.54841 + 8640184.812866 T_u + 0.093104 {T_u}^2 - 6.2\times10^{-6} {T_u}^3\tag 2 \\ &\,\text{is the Greenwich sidereal time at midnight, in seconds of time} \\ \omega^{\ast} = &\,7.2921158553\times10^{-5} + 4.3\times10^{-15} T_u\tag3 \\ &\,\text{is the Earth's sidereal rotation rate, in radians}/\text{second} \\ t\phantom{\Delta}\quad&\text{is UTC time of day, in seconds, and} \\ \Delta t\quad&\text{is the difference between UT1 and UTC.} \end{align}$$

Note that there's a units mismatch in equation (1). Equation (2) yields a result in units of sidereal seconds while equation (3), when multiplied by time of day, yields a result in units of radians. The latter needs to be converted to sidereal seconds. To do this, simply multiple both terms in equation (3) by $\frac{86400}{2\pi}$:

$$\begin{align} \omega^{\ast} = &\,1.00273790935 + 5.9\times10^{-11} T_u\tag4 \\ &\,\text{is the Earth's sidereal rotation rate, in sidereal seconds}/\text{UT second} \end{align}$$

Applying this to 2019 Jan 01 08:00:00 UTC yields the following:

  • The Julian date at midnight UTC 2019 Jan 01 is $2458484.5$, or $2451545+6939.5$. See footnote 1 for an explanation of this convention.
  • Dividing the residual 6939.5 by 36525 yields $T_u$: $T_u = 6939.5 / 36525 = 0.1899931553730322$.
  • Applying equation (2) yields the Greenwich Mean Sidereal Time at midnight: $H_0 = 1665686.527373333$ sidereal seconds.
  • Applying equation (4) and multiplying by $t=8*3600$ seconds yields another 28878.85178960284 sidereal seconds to the above. What about $\Delta t$, aka dUT1? I've ignored it (i.e., treated it as if it was zero.) See footnote 2.
  • 1665686.527373333 + 28878.85178960284 = 52965.37916293584 (modulo 86400). This is the mean sidereal time, in seconds. Converting to hours, minutes and seconds results in 14:42:45.4. The 0.4 seconds is bogus because I ignored dUT1.

Footnotes:

  1. There's a problem with using Julian dates to express time when high precision is of concern. The problem is that adding 2.45 million and change to the time since noon on 1 Jan 2000 results in potential truncation problems. The next representable Julian date after midnight 1 Jan 2000 is 40 microseconds after midnight. I myself try to stay a couple of orders of magnitude away from the quantization level, which means problems for dates at the millisecond level.

    Almost every scientific package that supports Julian dates either uses a pair of double precision numbers or uses some sort of modified Julian date, such as time since 12 Noon 16 Nov 1858 (omits the leading 2.4 million in Julian date), time since 12 noon 23 May 1968 (omits the leading 2.44 million), or time since J2000.

  2. The term $\Delta t$ in equations (2) and (4) is better known as dUT1. This is the difference between UT1 (time using the mean rotating Earth as a clock) and UTC (time that runs at the same rate as an atomic clock, but with leap seconds occasionally inserted to keep UT1 and UTC in sync). Those leap seconds keep $|\text{UT1}-\text{UTC}| < 0.9\,\text{seconds}$. Simply ignore the $\Delta t$ term if that's the case -- but remember that fractions of a second are meaningless.

  • Thanks for your comprehensive answer, it did really help! I finally found the problem that was causing wrong coordinate conversion and, unfortunately for this forum, it has fallen in the area of wrong unit conversion and signed operations. I plan to upload my C++ code in order to contribute a little but I'd love to know what would be the current/appropriate model to use if not this outdated one. Could you kindly provide a reference? – Carles Araguz May 15 '17 at 9:14
  • @CarlesAraguz -- Two references! One is the International Earth Rotation and Reference Systems Service (IERS) Conventions, described in IERS Technical Note 36. Pay particular attention to chapters 4 and 5. The other is the Standards of Fundamental Astronomy (SOFA) open source software maintained by the International Astronomical Union (IAU) at iausofa.org . If you poke, you'll see lots of functions in C and Fortran, plus some "cookbooks".The cookbook of interest is iausofa.org/sofa_pn_c.pdf . – David Hammen May 15 '17 at 9:36

I worked out the GMST using the method described in http://aa.usno.navy.mil/faq/docs/GAST.php. This should be an authoritative source, and claims 0.1 second accuracy.

In pseudo-code:

midnight = floor(2458484.833333) + 0.5            // J0 = 2458484.5
days_since_midnight = 2458484.833333 - 2458484.5  // = 0.333333
hours_since_midnight = 0.333333 * 24              // H = 8.0
days_since_epoch = 2458484.833333 - 2451545.0     // D = 6939.833333
centuries_since_epoch = days_since_epoch / 35625  // T = 0.190002
whole_days_since_epoch = 2458484.5 - 2451545.0    // D0 = 6939.5

GMST = 6.697374558 
     + 0.06570982441908 * whole_days_since_epoch 
     + 1.00273790935 * hours_since_midnight
     + 0.000026 * centuries_since_epoch**2        //=470.712605328

GMST_hours = 470 % 24 // = 14
GMST_minutes =  floor(0.712605328 * 60) //=42(.7563197)
GMST_seconds =  0.7563197*60 // =45.38

This calculation agrees with the online calculator. The method here (or as described at the link) should be easy to implement.

  • James - A minor error your pseudo-code: In your final calculation for GMST, you show the variable hours_since_epoch. This should be hours_since_midnight. – Robert Smith Jun 12 at 21:47
  • @RobertSmith Cured. – James K Jun 13 at 7:47

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.