Adding moon phases

Signed-off-by: Cédric Foellmi <[email protected]>

src/earth/moon/index.ts

@@ -23,6 +23,7 @@ import {
   getPhaseAngle,
   getPositionAngleOfTheBrightLimb
 } from './details'
+import { getAge, getTimeOfMeanPhase } from './phases'
 
 export const Moon: NaturalMoon = {
   getMeanLongitude,
@@ -44,5 +45,7 @@ export const Moon: NaturalMoon = {
   getPhaseAngle,
   getIlluminatedFraction,
   getEquatorialHorizontalParallax,
-  getPositionAngleOfTheBrightLimb
+  getPositionAngleOfTheBrightLimb,
+  getTimeOfMeanPhase,
+  getAge
 }

src/earth/moon/phases.ts

@@ -0,0 +1,56 @@
+import Decimal from '@/decimal'
+import { Day, JulianDay } from '@/types'
+import { MOON_SYNODIC_PERIOD, MoonPhase } from '@/constants'
+import { fmod } from '@/utils'
+import { getDecimalYear } from '@/dates'
+
+// The value of K must be an integer
+function getK (jd: JulianDay | number): Decimal {
+  const decimalYear = getDecimalYear(jd)
+  const decimalK = new Decimal('12.3685').mul(decimalYear.minus('2000'))
+  return decimalK.isPositive() ? Decimal.floor(decimalK) : Decimal.ceil(decimalK)
+}
+
+function getPhaseK (jd: JulianDay | number, phase: MoonPhase): Decimal {
+  let k = getK(jd)
+  if (phase === MoonPhase.FirstQuarter) {
+    k = k.plus(0.25)
+  } else if (phase == MoonPhase.Full) {
+    k = k.plus(0.5)
+  } else if (phase == MoonPhase.LastQuarter) {
+    k = k.plus(0.75)
+  }
+  return k
+}
+
+/**
+ * The time of a given Moon phase.
+ * Results are already corrected for the Sun's aberration and by the Moon's light-time.
+ * @param {JulianDay} jd The julian day
+ * @param {MoonPhase} phase The requested phase
+ * @return {JulianDay}
+ */
+export function getTimeOfMeanPhase (jd: JulianDay | number, phase: MoonPhase): JulianDay {
+  const k = getPhaseK(jd, phase)
+  const T = k.dividedBy('1236.85')
+  return new Decimal('2451_550.097_66')
+    .plus(new Decimal('29.530_588_861').mul(k))
+    .plus(new Decimal('0.000_154_37').mul(T.pow(2)))
+    .minus(new Decimal('0.000_000_150').mul(T.pow(3)))
+    .plus(new Decimal('0.000_000_000_73').mul(T.pow(4)))
+}
+
+/**
+ * The age of the Moon cycle (0 = New Moon, MOON_SYNODIC_PERIOD/2 = Full Moon).
+ * This is a low-accuracy age of the moon, using the average moon synodic period.
+ * @param {JulianDay} jd The julian day
+ * @return {JulianDay}
+ */
+export function getAge (jd: JulianDay | number): Day {
+  const djd = new Decimal(jd)
+  let jdNewMoon = getTimeOfMeanPhase(djd.minus(MOON_SYNODIC_PERIOD), MoonPhase.New)
+  if (jdNewMoon.greaterThan(djd)) {
+    jdNewMoon = jdNewMoon.minus(MOON_SYNODIC_PERIOD)
+  }
+  return fmod(djd.minus(jdNewMoon), MOON_SYNODIC_PERIOD)
+}

tests/moon.test.js

@@ -1,5 +1,6 @@
 import { constants, Earth, juliandays, Sun, times } from '@'
 import { getDecimalValue } from '@/sexagesimal'
+import { MOON_SYNODIC_PERIOD, MoonPhase } from '@/constants'
 
 describe('moon', () => {
   test('get moon mean longitude', () => {
@@ -53,4 +54,27 @@ describe('moon', () => {
     const k = Earth.Moon.getIlluminatedFraction(jd)
     expect(k.toNumber()).toBeCloseTo(0.6786, 3)
   })
+
+  // See example 49.a, AA p 353.
+  test('get time of new moon', () => {
+    // Mid-february
+    const UTCDate = new Date(Date.UTC(1977, 1, 12))
+    const T = juliandays.getJulianCentury(juliandays.getJulianDay(UTCDate))
+    expect(T.toNumber()).toBeCloseTo(-0.228_81, 4)
+    const newMoonJD = Earth.Moon.getTimeOfMeanPhase(juliandays.getJulianDay(UTCDate), MoonPhase.New)
+    expect(newMoonJD.toNumber()).toBeCloseTo(2443_192.941_02, 5)
+  })
+
+  // See example 49.a, AA p 353.
+  test('get new moon age', () => {
+    const UTCDate = new Date(Date.UTC(1977, 1, 12))
+    const newMoonJD = Earth.Moon.getTimeOfMeanPhase(juliandays.getJulianDay(UTCDate), MoonPhase.New)
+    expect(Earth.Moon.getAge(newMoonJD).toNumber()).toEqual(0)
+    const fqMoonJD = Earth.Moon.getTimeOfMeanPhase(juliandays.getJulianDay(UTCDate), MoonPhase.FirstQuarter)
+    expect(Earth.Moon.getAge(fqMoonJD).toNumber()).toBeCloseTo(MOON_SYNODIC_PERIOD.dividedBy(4).toNumber(), 6)
+    const fullMoonJD = Earth.Moon.getTimeOfMeanPhase(juliandays.getJulianDay(UTCDate), MoonPhase.Full)
+    expect(Earth.Moon.getAge(fullMoonJD).toNumber()).toBeCloseTo(MOON_SYNODIC_PERIOD.dividedBy(2).toNumber(), 5)
+    const lqMoonJD = Earth.Moon.getTimeOfMeanPhase(juliandays.getJulianDay(UTCDate), MoonPhase.LastQuarter)
+    expect(Earth.Moon.getAge(lqMoonJD).toNumber()).toBeCloseTo(MOON_SYNODIC_PERIOD.mul(3).dividedBy(4).toNumber(), 5)
+  })
 })