implement erfa and time

scripts/astropy/time.py

@@ -17,6 +17,24 @@ def time(scale):
     print('expect(tdb(time)).toMatchTime([{0}, {1:.18f}, Timescale.TDB])'.format(t.tdb.jd1, t.tdb.jd2))
     print('expect(tcb(time)).toMatchTime([{0}, {1:.18f}, Timescale.TCB])'.format(t.tcb.jd1, t.tcb.jd2))
 
+def gast():
+    t = Time('2020-10-07T12:00:00', format='isot', scale='utc')
+    print('GAST: {0:.18f}'.format(t.sidereal_time('apparent', 'tio')))
+
+def gmst():
+    t = Time('2020-10-07T12:00:00', format='isot', scale='utc')
+    print('GMST: {0:.18f}'.format(t.sidereal_time('mean', 'tio')))
+
+def era():
+    t = Time('2020-10-07T12:00:00', format='isot', scale='utc')
+    print('ERA: {0:.18f}'.format(t.earth_rotation_angle('tio')))
+
 match sys.argv[1]:
     case 'ut1' | 'utc' | 'tai' | 'tt' | 'tcg' | 'tdb' | 'tcb':
         time(sys.argv[1])
+    case 'gast':
+        gast()
+    case 'gmst':
+        gmst()
+    case 'era':
+        era()

src/erfa.test.ts

@@ -1,6 +1,6 @@
 import { expect, test } from 'bun:test'
 import { kilometer } from './distance'
-import { eraCalToJd, eraDat, eraDtDb, eraEors, eraEra00, eraFad03, eraFae03, eraFaf03, eraFaju03, eraFal03, eraFalp03, eraFama03, eraFame03, eraFaom03, eraFapa03, eraFasa03, eraFaur03, eraFave03, eraFw2m, eraGst06, eraGst06a, eraJdToCal, eraNut00a, eraNut06a, eraObl06, eraPfw06, eraPnm06a, eraS06, eraSp00, eraTaiTt, eraTaiUt1, eraTaiUtc, eraTcbTdb, eraTcgTt, eraTdbTcb, eraTdbTt, eraTtTai, eraTtTcg, eraTtTdb, eraUt1Tai, eraUt1Utc, eraUtcTai, eraUtcUt1 } from './erfa'
+import { eraCalToJd, eraDat, eraDtDb, eraEors, eraEra00, eraFad03, eraFae03, eraFaf03, eraFaju03, eraFal03, eraFalp03, eraFama03, eraFame03, eraFaom03, eraFapa03, eraFasa03, eraFaur03, eraFave03, eraFw2m, eraGmst06, eraGst06, eraGst06a, eraJdToCal, eraNut00a, eraNut06a, eraObl06, eraPfw06, eraPnm06a, eraS06, eraSp00, eraTaiTt, eraTaiUt1, eraTaiUtc, eraTcbTdb, eraTcgTt, eraTdbTcb, eraTdbTt, eraTtTai, eraTtTcg, eraTtTdb, eraUt1Tai, eraUt1Utc, eraUtcTai, eraUtcUt1 } from './erfa'
 import type { Mat3 } from './matrix'
 
 test('eraTaiUt1', () => {
@@ -128,6 +128,10 @@ test('eraGst06', () => {
 	expect(eraGst06(2453736.0, 0.5, 2453736.0, 0.5, rnpb)).toBeCloseTo(1.754166138018167568, 12)
 })
 
+test('eraGmst06', () => {
+	expect(eraGmst06(2453736.0, 0.5, 2453736.0, 0.5)).toBeCloseTo(1.754174971870091203, 12)
+})
+
 test('eraEra00', () => {
 	expect(eraEra00(2454388.0, 0.5)).toBeCloseTo(0.4022837240028158102, 12)
 })

src/erfa.ts

@@ -443,6 +443,14 @@ export function eraGst06(ut11: number, ut12: number, tt1: number, tt2: number, r
 	return normalize(era - eors)
 }
 
+// Greenwich mean sidereal time (model consistent with IAU 2006 precession).
+export function eraGmst06(ut11: number, ut12: number, tt1: number, tt2: number): Angle {
+	// TT Julian centuries since J2000.0.
+	const t = (tt1 - J2000 + tt2) / DAYSPERJC
+
+	return normalize(eraEra00(ut11, ut12) + arcsec(0.014506 + (4612.156534 + (1.3915817 + (-0.00000044 + (-0.000029956 + -0.0000000368 * t) * t) * t) * t) * t))
+}
+
 // Earth rotation angle (IAU 2000 model).
 export function eraEra00(ut11: number, ut12: number): Angle {
 	const t = ut12 + (ut11 - J2000)
@@ -451,7 +459,7 @@ export function eraEra00(ut11: number, ut12: number): Angle {
 	const f = (ut12 % 1.0) + (ut11 % 1.0)
 
 	// Earth rotation angle at this UT1.
-	return eraAnpm(TAU * (f + 0.779057273264 + 0.00273781191135448 * t))
+	return normalize(TAU * (f + 0.779057273264 + 0.00273781191135448 * t))
 }
 
 // Equation of the origins, given the classical NPB matrix and the quantity [s].

src/time.test.ts

@@ -1,7 +1,8 @@
 import { beforeAll, expect, test, type CustomMatcher } from 'bun:test'
+import { hour } from './angle'
 import { J2000 } from './constants'
 import { iersa, iersb } from './iers'
-import { normalize, tai, tcb, tcg, tdb, tdbMinusTt, tdbMinusTtByFairheadAndBretagnon1990, time, timeBesselian, timeGPS, timeJulian, timeMJD, Timescale, timeUnix, timeYMDHMS, tt, ut1, utc, type Time } from './time'
+import { era, gast, gmst, normalize, tai, tcb, tcg, tdb, tdbMinusTt, tdbMinusTtByFairheadAndBretagnon1990, time, timeBesselian, timeGPS, timeJulian, timeMJD, Timescale, timeUnix, timeYMDHMS, tt, ut1, utc, type Time } from './time'
 
 const toMatchTime: CustomMatcher<Time, never[]> = (actual, expected: Time, precision?: number) => {
 	const b = normalize(expected.day, expected.fraction)
@@ -206,28 +207,28 @@ test('memoize', () => {
 
 	for (let i = 0; i < 10000; i++) {
 		const a = ut1(t)
-		expect(a.tcb).toEqual(t)
-		expect(t.ut1).toEqual(a)
+		expect(a.extra?.tcb).toEqual(t)
+		expect(t.extra?.ut1).toEqual(a)
 
 		const b = utc(t)
-		expect(b.tcb).toEqual(t)
-		expect(t.utc).toEqual(b)
+		expect(b.extra?.tcb).toEqual(t)
+		expect(t.extra?.utc).toEqual(b)
 
 		const c = tai(t)
-		expect(c.tcb).toEqual(t)
-		expect(t.tai).toEqual(c)
+		expect(c.extra?.tcb).toEqual(t)
+		expect(t.extra?.tai).toEqual(c)
 
 		const d = tt(t)
-		expect(d.tcb).toEqual(t)
-		expect(t.tt).toEqual(d)
+		expect(d.extra?.tcb).toEqual(t)
+		expect(t.extra?.tt).toEqual(d)
 
 		const e = tcg(t)
-		expect(e.tcb).toEqual(t)
-		expect(t.tcg).toEqual(e)
+		expect(e.extra?.tcb).toEqual(t)
+		expect(t.extra?.tcg).toEqual(e)
 
 		const f = tdb(t)
-		expect(f.tcb).toEqual(t)
-		expect(t.tdb).toEqual(f)
+		expect(f.extra?.tcb).toEqual(t)
+		expect(t.extra?.tdb).toEqual(f)
 	}
 }, 1000)
 
@@ -244,3 +245,21 @@ test('tdbMinusTtByFairheadAndBretagnon1990', () => {
 
 	expect(tdbMinusTtByFairheadAndBretagnon1990(t0)).toBeCloseTo(tdbMinusTt(t0), 5)
 })
+
+test('gast', () => {
+	const t = timeYMDHMS(2020, 10, 7, 12, 0, 0, Timescale.UTC)
+	expect(gast(t)).toBe(t.extra!.gast!)
+	expect(t.extra?.gast).toBeCloseTo(hour(13.106038262872143463), 15)
+})
+
+test('gmst', () => {
+	const t = timeYMDHMS(2020, 10, 7, 12, 0, 0, Timescale.UTC)
+	expect(gmst(t)).toBe(t.extra!.gmst!)
+	expect(t.extra?.gmst).toBeCloseTo(hour(13.106345240224241522), 15)
+})
+
+test('era', () => {
+	const t = timeYMDHMS(2020, 10, 7, 12, 0, 0, Timescale.UTC)
+	expect(era(t)).toBe(t.extra!.era!)
+	expect(t.extra?.era).toBeCloseTo(hour(13.088607043262001639), 15)
+})