implement erfa, iers, math, matrix and time

package.json

@@ -6,11 +6,11 @@
 		"lint": "eslint ."
 	},
 	"devDependencies": {
-        "@types/bun": "latest",
-        "eslint": "9.14.0",
-        "typescript-eslint": "8.14.0"
+		"@types/bun": "latest",
+		"eslint": "9.16.0",
+		"typescript-eslint": "8.18.0"
 	},
 	"peerDependencies": {
-		"typescript": "5.6.3"
+		"typescript": "5.7.2"
 	}
 }

src/erfa.ts

@@ -5,13 +5,13 @@ import { roundToNearestWholeNumber } from './math'
 const DBL_EPSILON = 2.220446049250313e-16
 
 // International Atomic Time, TAI, to Universal Time, UT1.
-export function eraTaiUt1(tai1: number, tai2: number, ut1_minus_tai: number): [number, number] {
-	return [tai1, tai2 + ut1_minus_tai / DAYSEC]
+export function eraTaiUt1(tai1: number, tai2: number, ut1MinusTai: number): [number, number] {
+	return [tai1, tai2 + ut1MinusTai / DAYSEC]
 }
 
 // Universal Time, UT1, to International Atomic Time, TAI.
-export function eraUt1Tai(ut11: number, ut12: number, ut1_minus_tai: number): [number, number] {
-	return [ut11, ut12 - ut1_minus_tai / DAYSEC]
+export function eraUt1Tai(ut11: number, ut12: number, ut1MinusTai: number): [number, number] {
+	return [ut11, ut12 - ut1MinusTai / DAYSEC]
 }
 
 // International Atomic Time, TAI, to Coordinated Universal Time, UTC.

src/iers.ts

@@ -0,0 +1,40 @@
+import type { Angle } from './angle'
+import type { PolarMotion, Time, TimeDelta } from './time'
+
+export type Iers = TimeDelta & PolarMotion
+
+export class IersA implements Iers {
+	delta(time: Time): number {
+		return 0
+	}
+
+	xy(time: Time): [Angle, Angle] {
+		return [0, 0]
+	}
+}
+
+export class IersB implements Iers {
+	delta(time: Time): number {
+		return 0
+	}
+
+	xy(time: Time): [Angle, Angle] {
+		return [0, 0]
+	}
+}
+
+export class IersAB implements Iers {
+	constructor(readonly a: IersA, readonly b: IersB) {}
+
+	delta(time: Time): number {
+		return this.a.delta(time) || this.b.delta(time)
+	}
+
+	xy(time: Time): [Angle, Angle] {
+		return [0, 0]
+	}
+}
+
+export const iersa = new IersA()
+export const iersb = new IersB()
+export const iersab = new IersAB(iersa, iersb)

src/math.ts

@@ -1,9 +1,3 @@
-// Rectangular array of numbers with two rows and two columns.
-export type Matrix2 = [number, number, number, number]
-
-// Rectangular array of numbers with three rows and three columns.
-export type Matrix3 = [number, number, number, number, number, number, number, number, number]
-
 // Adds [a] and [b] exactly, returning the result as two 64-bit floats.
 export function twoSum(a: number, b: number): [number, number] {
 	const x = a + b

src/matrix.test.ts

@@ -0,0 +1,60 @@
+import { expect, test } from 'bun:test'
+import { identity, Mat3, rotX, rotY, rotZ } from './matrix'
+
+test('rotX', () => {
+	const m = new Mat3(2, 3, 2, 3, 2, 3, 3, 4, 5)
+	m.rotX(0.3456789)
+
+	expect(m[0]).toBeCloseTo(2, 12)
+	expect(m[1]).toBeCloseTo(3, 12)
+	expect(m[2]).toBeCloseTo(2, 12)
+	expect(m[3]).toBeCloseTo(3.83904338823561246, 12)
+	expect(m[4]).toBeCloseTo(3.237033249594111899, 12)
+	expect(m[5]).toBeCloseTo(4.516714379005982719, 12)
+	expect(m[6]).toBeCloseTo(1.806030415924501684, 12)
+	expect(m[7]).toBeCloseTo(3.085711545336372503, 12)
+	expect(m[8]).toBeCloseTo(3.687721683977873065, 12)
+
+	expect(identity().rotX(0.3456789)).toEqual(rotX(0.3456789))
+})
+
+test('rotY', () => {
+	const m = new Mat3(2, 3, 2, 3, 2, 3, 3, 4, 5)
+	m.rotY(0.3456789)
+
+	expect(m[0]).toBeCloseTo(0.865184781897815993, 12)
+	expect(m[1]).toBeCloseTo(1.467194920539316554, 12)
+	expect(m[2]).toBeCloseTo(0.1875137911274457342, 12)
+	expect(m[3]).toBeCloseTo(3, 12)
+	expect(m[4]).toBeCloseTo(2, 12)
+	expect(m[5]).toBeCloseTo(3, 12)
+	expect(m[6]).toBeCloseTo(3.50020789285042733, 12)
+	expect(m[7]).toBeCloseTo(4.77988902226229815, 12)
+	expect(m[8]).toBeCloseTo(5.381899160903798712, 12)
+
+	expect(identity().rotY(0.3456789)).toEqual(rotY(0.3456789))
+})
+
+test('rotZ', () => {
+	const m = new Mat3(2, 3, 2, 3, 2, 3, 3, 4, 5)
+	m.rotZ(0.3456789)
+
+	expect(m[0]).toBeCloseTo(2.898197754208926769, 12)
+	expect(m[1]).toBeCloseTo(3.50020789285042733, 12)
+	expect(m[2]).toBeCloseTo(2.898197754208926769, 12)
+	expect(m[3]).toBeCloseTo(2.144865911309686813, 12)
+	expect(m[4]).toBeCloseTo(0.865184781897815993, 12)
+	expect(m[5]).toBeCloseTo(2.144865911309686813, 12)
+	expect(m[6]).toBeCloseTo(3, 12)
+	expect(m[7]).toBeCloseTo(4, 12)
+	expect(m[8]).toBeCloseTo(5, 12)
+
+	expect(identity().rotZ(0.3456789)).toEqual(rotZ(0.3456789))
+})
+
+test('clone', () => {
+	const m = new Mat3(2, 3, 2, 3, 2, 3, 3, 4, 5)
+	const n = m.clone()
+	expect(n).toEqual(m)
+	expect(m === n).toBe(false)
+})

src/matrix.ts

@@ -0,0 +1,117 @@
+import type { Angle } from './angle'
+
+// Rectangular array of numbers with three rows and three columns.
+export class Mat3 extends Float64Array {
+	constructor(a: number = 0, b: number = 0, c: number = 0, d: number = 0, e: number = 0, f: number = 0, g: number = 0, h: number = 0, i: number = 0) {
+		super(9)
+
+		this[0] = a
+		this[1] = b
+		this[2] = c
+		this[3] = d
+		this[4] = e
+		this[5] = f
+		this[6] = g
+		this[7] = h
+		this[8] = i
+	}
+
+	// Clones this matrix.
+	clone() {
+		const matrix = new Mat3()
+		matrix.set(this)
+		return matrix
+	}
+
+	// Rotates this matrix around X axis.
+	rotX(angle: Angle): this {
+		const ca = Math.cos(angle)
+		const sa = Math.sin(angle)
+
+		const m3 = ca * this[3] + sa * this[6]
+		const m4 = ca * this[4] + sa * this[7]
+		const m5 = ca * this[5] + sa * this[8]
+		const m6 = -sa * this[3] + ca * this[6]
+		const m7 = -sa * this[4] + ca * this[7]
+		const m8 = -sa * this[5] + ca * this[8]
+
+		this[3] = m3
+		this[4] = m4
+		this[5] = m5
+		this[6] = m6
+		this[7] = m7
+		this[8] = m8
+
+		return this
+	}
+
+	// Rotates this matrix around Y axis.
+	rotY(angle: Angle): this {
+		const ca = Math.cos(angle)
+		const sa = Math.sin(angle)
+
+		const m0 = ca * this[0] - sa * this[6]
+		const m1 = ca * this[1] - sa * this[7]
+		const m2 = ca * this[2] - sa * this[8]
+		const m6 = sa * this[0] + ca * this[6]
+		const m7 = sa * this[1] + ca * this[7]
+		const m8 = sa * this[2] + ca * this[8]
+
+		this[0] = m0
+		this[1] = m1
+		this[2] = m2
+		this[6] = m6
+		this[7] = m7
+		this[8] = m8
+
+		return this
+	}
+
+	// Rotates this matrix around Z axis.
+	rotZ(angle: Angle): this {
+		const ca = Math.cos(angle)
+		const sa = Math.sin(angle)
+
+		const m0 = ca * this[0] + sa * this[3]
+		const m1 = ca * this[1] + sa * this[4]
+		const m2 = ca * this[2] + sa * this[5]
+		const m3 = -sa * this[0] + ca * this[3]
+		const m4 = -sa * this[1] + ca * this[4]
+		const m5 = -sa * this[2] + ca * this[5]
+
+		this[0] = m0
+		this[1] = m1
+		this[2] = m2
+		this[3] = m3
+		this[4] = m4
+		this[5] = m5
+
+		return this
+	}
+}
+
+// Creates a new Identity Matrix.
+export function identity() {
+	return new Mat3(1, 0, 0, 0, 1, 0, 0, 0, 1)
+}
+
+// Creates a new rotation matrix around X axis.
+export function rotX(angle: Angle): Mat3 {
+	const ca = Math.cos(angle)
+	const sa = Math.sin(angle)
+	return new Mat3(1, 0, 0, 0, ca, sa, 0, -sa, ca)
+}
+
+// Creates a new rotation matrix around Y axis.
+export function rotY(angle: Angle): Mat3 {
+	const ca = Math.cos(angle)
+	const sa = Math.sin(angle)
+	return new Mat3(ca, 0, -sa, 0, 1, 0, sa, 0, ca)
+}
+
+// Creates a new rotation matrix around Z axis.
+export function rotZ(angle: Angle): Mat3 {
+	const ca = Math.cos(angle)
+	const sa = Math.sin(angle)
+	return new Mat3(ca, sa, 0, -sa, ca, 0, 0, 0, 1)
+}

src/time.test.ts

@@ -1,49 +1,72 @@
 import { expect, test } from 'bun:test'
 import { J2000 } from './constants'
-import { time, timeJulian, timeMJD, timeUnix, timeYMDHMS } from './time'
+import { time, timeBesselian, timeGPS, timeJulian, timeMJD, Timescale, timeUnix, timeYMDHMS } from './time'
 
 test('time', () => {
-	const t = time(2449353.623)
+	const t = time(2449353.623, 0, Timescale.TT)
 	expect(t[0]).toBe(2449354)
 	expect(t[1]).toBeCloseTo(-0.377, 3)
+	expect(t[2]).toBe(Timescale.TT)
 })
 
 test('timeUnix', () => {
 	let t = timeUnix(0)
 	expect(t[0]).toBe(2440588)
 	expect(t[1]).toBe(-0.5)
+	expect(t[2]).toBe(Timescale.UTC)
 
-	t = timeUnix(946684800)
+	t = timeUnix(946684800, Timescale.TAI)
 	expect(t[0]).toBe(J2000)
 	expect(t[1]).toBe(-0.5)
+	expect(t[2]).toBe(Timescale.TAI)
 })
 
 test('timeMJD', () => {
-	const t = timeMJD(51544.0)
+	const t = timeMJD(51544.0, Timescale.UT1)
 	expect(t[0]).toBe(J2000)
 	expect(t[1]).toBe(-0.5)
+	expect(t[2]).toBe(Timescale.UT1)
 })
 
 test('timeYMDHMS', () => {
 	let t = timeYMDHMS(2022, 1, 1, 12, 0, 0.0)
 	expect(t[0]).toBe(2459581)
 	expect(t[1]).toBeCloseTo(0, 14)
+	expect(t[2]).toBe(Timescale.UTC)
 
 	t = timeYMDHMS(2023, 6, 1, 23, 59, 59.0)
 	expect(t[0]).toBe(2460097)
 	expect(t[1]).toBeCloseTo(0.49998842592592596, 14)
+	expect(t[2]).toBe(Timescale.UTC)
 
 	t = timeYMDHMS(2024, 1, 1, 0, 0, 0.0)
 	expect(t[0]).toBe(2460311)
 	expect(t[1]).toBeCloseTo(-0.5, 14)
+	expect(t[2]).toBe(Timescale.UTC)
 
-	t = timeYMDHMS(2025, 12, 31, 17, 59, 43.0)
+	t = timeYMDHMS(2025, 12, 31, 17, 59, 43.0, Timescale.TDB)
 	expect(t[0]).toBe(2461041)
 	expect(t[1]).toBeCloseTo(0.24980324074074078, 14)
+	expect(t[2]).toBe(Timescale.TDB)
 })
 
 test('timeJulian', () => {
-	const t = timeJulian(2024.0)
+	const t = timeJulian(2024.0, Timescale.TCG)
 	expect(t[0]).toBe(2460311)
 	expect(t[1]).toBe(0)
+	expect(t[2]).toBe(Timescale.TCG)
+})
+
+test('timeBesselian', () => {
+	const t = timeBesselian(1950.0, Timescale.TCB)
+	expect(t[0]).toBe(2433282)
+	expect(t[1]).toBe(0.42345904977992177)
+	expect(t[2]).toBe(Timescale.TCB)
+})
+
+test('timGPS', () => {
+	const t = timeGPS(630720013)
+	expect(t[0]).toBe(J2000)
+	expect(t[1]).toBe(-0.4996296167373657)
+	expect(t[2]).toBe(Timescale.TAI)
 })