Upgrading mercury with constants and Decimal.js

Signed-off-by: Cédric Foellmi <cedric@onekiloparsec.dev>

src/planets/mercury/base.ts

@@ -1,28 +1,31 @@
-import { JulianDay, PlanetConstants } from '@/types'
+import { JulianDay } from '@/types'
 import { getFractionalYear } from '@/dates'
+import Decimal from 'decimal.js'
 
 // The value of K must be an integer
-function getK (jd: JulianDay): number {
+function getK (jd: JulianDay | number): Decimal {
  const decimalYear = getFractionalYear(jd)
- return Math.floor(4.15201 * (decimalYear - 2000.12))
+ return Decimal.floor(new Decimal(4.15201).mul(decimalYear.minus(2000.12)))
 }
 
 /**
  * Aphelion (time of passage at the closest point to the Sun)
  * @param {JulianDay} jd The julian day
  * @returns {JulianDay}
  */
-export function getAphelion (jd: JulianDay): JulianDay {
- const kdash = getK(jd) + 0.5
- return 2451590.257 + 87.96934963 * kdash
+export function getAphelion (jd: JulianDay| number): JulianDay {
+ const kdash = getK(jd).plus(0.5)
+ return new Decimal(2451590.257)
+ .plus(new Decimal(87.96934963).mul(kdash))
 }
 
 /**
  * Perihelion (time of passage at the farthest point to the Sun)
  * @param {JulianDay} jd The julian day
  * @returns {JulianDay}
  */
-export function getPerihelion (jd: JulianDay): JulianDay {
+export function getPerihelion (jd: JulianDay| number): JulianDay {
  const k = getK(jd)
- return 2451590.257 + 87.96934963 * k
+ return new Decimal(2451590.257)
+ .plus(new Decimal(87.96934963).mul(k))
 }

src/planets/mercury/coefficients.ts

@@ -1,3 +1,5 @@
+import Decimal from 'decimal.js'
+
 export const g_L0MercuryCoefficients =
  [
  [440250710, 0, 0],
@@ -39,7 +41,7 @@ export const g_L0MercuryCoefficients =
  [118, 2.781, 77204.327],
  [106, 4.206, 19804.827]
  ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
  })
 
 export const g_L1MercuryCoefficients =
@@ -61,7 +63,7 @@ export const g_L1MercuryCoefficients =
  [28, 3.04, 51066.43],
  [27, 5.09, 234791.13]
  ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
  })
 
 export const g_L2MercuryCoefficients =
@@ -77,7 +79,7 @@ export const g_L2MercuryCoefficients =
  [15, 4.63, 1109.38],
  [12, 0.79, 208703.23]
  ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
  })
 
 export const g_L3MercuryCoefficients = [
@@ -90,7 +92,7 @@ export const g_L3MercuryCoefficients = [
  [7, 5.82, 156527.42],
  [3, 2.57, 182615.32]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_L4MercuryCoefficients = [
@@ -101,13 +103,13 @@ export const g_L4MercuryCoefficients = [
  [1, 4.50, 104351.61],
  [1, 1.27, 130439.52]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_L5MercuryCoefficients = [
  [1, 3.14, 0]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_B0MercuryCoefficients = [
@@ -126,7 +128,7 @@ export const g_B0MercuryCoefficients = [
  [121, 1.813, 53285.185],
  [100, 5.657, 20426.571]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_B1MercuryCoefficients = [
@@ -142,7 +144,7 @@ export const g_B1MercuryCoefficients = [
  [28, 0.29, 27197.28],
  [26, 5.98, 234791.13]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_B2MercuryCoefficients = [
@@ -156,7 +158,7 @@ export const g_B2MercuryCoefficients = [
  [18, 4.67, 182615.32],
  [7, 1.43, 208703.23]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_B3MercuryCoefficients = [
@@ -168,14 +170,14 @@ export const g_B3MercuryCoefficients = [
  [3, 3.12, 130439.52],
  [2, 6.27, 156527.42]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_B4MercuryCoefficients = [
  [4, 1.75, 26087.90],
  [1, 3.14, 0]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_R0MercuryCoefficients = [
@@ -193,7 +195,7 @@ export const g_R0MercuryCoefficients = [
  [142, 6.253, 24978.525],
  [100, 3.734, 21535.950]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_R1MercuryCoefficients = [
@@ -206,7 +208,7 @@ export const g_R1MercuryCoefficients = [
  [153, 1.061, 156527.419],
  [39, 4.11, 182615.32]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_R2MercuryCoefficients = [
@@ -218,7 +220,7 @@ export const g_R2MercuryCoefficients = [
  [22, 3.14, 0],
  [13, 5.80, 156527.42]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
 })
 
 export const g_R3MercuryCoefficients = [
@@ -228,5 +230,5 @@ export const g_R3MercuryCoefficients = [
  [5, 4.44, 104351.61],
  [2, 1.21, 130439.52]
 ].map((a) => {
- return { A: a[0], B: a[1], C: a[2] }
-})
+ return { A: new Decimal(a[0]), B: new Decimal(a[1]), C: new Decimal(a[2]) }
+})

src/planets/mercury/constants.ts

@@ -1,46 +1,66 @@
-import { PlanetConstants, PlanetOrbitalElements } from '@/types'
+import { AstronomicalUnit, Degree, LengthArray, PlanetConstants, PlanetOrbitalElements } from '@/types'
+import Decimal from 'decimal.js'
 
 /**
  * Planet constants, copied from the JPL
  * Reference: https://ssd.jpl.nasa.gov/?planet_phys_par
  */
 export const constants: PlanetConstants = {
- equatorialRadius: 2440.53,
- meanRadius: 2439.4,
- mass: 0.330114,
- bulkDensity: 5.4291,
- siderealRotationPeriod: 58.6462,
- siderealOrbitPeriod: 0.2408467,
- visualMagnitude: -0.60,
- geometricAlbedo: 0.106,
- equatorialGravity: 3.70,
- escapeVelocity: 4.25
+ equatorialRadius: new Decimal(2440.53),
+ meanRadius: new Decimal(2439.4),
+ mass: new Decimal(0.330114),
+ bulkDensity: new Decimal(5.4291),
+ siderealRotationPeriod: new Decimal(58.6462),
+ siderealOrbitPeriod: new Decimal(0.2408467),
+ visualMagnitude: new Decimal(-0.60),
+ geometricAlbedo: new Decimal(0.106),
+ equatorialGravity: new Decimal(3.70),
+ escapeVelocity: new Decimal(4.25)
 }
 
 
-
 /**
  * Orbital Elements for the mean equinox of date
  * Reference: Astronomical Algorithms, J. Meus, pp. 212-213 (Table 31.A).
  */
 export const orbitalElements: PlanetOrbitalElements = {
- meanLongitude: [181.979801, 58519.2130302, 0.00031014, 0.000000015],
- semiMajorAxis: 0.723329820,
- eccentricity: [0.00677192, -0.000047765, 0.0000000981, 0.00000000046],
- inclination: [3.394662, 0.0010037, -0.00000088, -0.000000007],
- longitudeOfAscendingNode: [76.679920, 0.9011206, 0.00040618, -0.000000093],
- longitudeOfPerihelion: [131.563703, 1.4022288, -0.00107618, -0.000005678]
+ meanLongitude: [252.250_906, 149474.072_2491, 0.000_303_50, 0.000_000_018]
+ .map(v => new Decimal(v)) as LengthArray<Degree, 4>,
+
+ semiMajorAxis: [0.387_098_310, 0, 0, 0]
+ .map(v => new Decimal(v)) as LengthArray<AstronomicalUnit, 4>,
+
+ eccentricity: [0.205_631_75, 0.000_020_407, -0.000_000_0283, -0.000_000_000_18]
+ .map(v => new Decimal(v)) as LengthArray<Degree, 4>,
+
+ inclination: [7.004_986, 0.001_8215, -0.000_018_10, 0.000_000_056]
+ .map(v => new Decimal(v)) as LengthArray<Degree, 4>,
+
+ longitudeOfAscendingNode: [48.330_893, 1.186_1883, 0.000_175_42, 0.000_000_215]
+ .map(v => new Decimal(v)) as LengthArray<Degree, 4>,
+
+ longitudeOfPerihelion: [77.456_119, 1.556_4776, 0.000_295_44, 0.000_000_009]
+ .map(v => new Decimal(v)) as LengthArray<Degree, 4>
 }
 
 /**
  * Orbital Elements for the standard equinox J2000
  * Reference: Astronomical Algorithms, J. Meus, pp. 212-213 (Table 31.B).
  */
 export const orbitalElementsJ2000: PlanetOrbitalElements = {
- meanLongitude: [181.979801, 58517.8156760, 0.00000165, -0.000000002],
+ meanLongitude: [252.250_906, 149472.674_6358, -0.000_005_36, 0.000_000_002]
+ .map(v => new Decimal(v)) as LengthArray<Degree, 4>,
+
  semiMajorAxis: orbitalElements.semiMajorAxis,
+
  eccentricity: orbitalElements.eccentricity,
- inclination: [3.394662, -0.0008568, -0.00003244, 0.00000009],
- longitudeOfAscendingNode: [76.679920, -0.2780134, -0.00014257, -0.000000164],
- longitudeOfPerihelion: [131.563703, 0.0048746, -0.00138467, -0.000005695]
+
+ inclination: [7.004_986, -0.005_9516, 0.000_000_80, 0.000_000_043]
+ .map(v => new Decimal(v)) as LengthArray<Degree, 4>,
+
+ longitudeOfAscendingNode: [48.330_893, -0.125_4227, -0.000_088_33, -0.000_000_200]
+ .map(v => new Decimal(v)) as LengthArray<Degree, 4>,
+
+ longitudeOfPerihelion: [77.456_119, 0.158_8643, -0.000_013_42, -0.000_000_007]
+ .map(v => new Decimal(v)) as LengthArray<Degree, 4>
 }