Merge pull request #26 from ironwallaby/master

Optimize phase()

.gitignore

@@ -1 +1,2 @@
+.*swp
 /node_modules

.travis.yml

@@ -1,5 +1,3 @@
 language: node_js
-
 node_js:
-- "4.4"
-- "5.8"
+- "10"

lib/lune.js

@@ -2,10 +2,16 @@
  * This library calculates the current phase of the moon
  * as well as finds the dates of the recent moon phases.
  *
- * Ported from python version found here:
+ * Some functionality is ported from python version found here:
  * https://bazaar.launchpad.net/~keturn/py-moon-phase/trunk/annotate/head:/moon.py
  *
+ * Some functionality is taken from Astronomical Algorithms, 2nd ed.
+ *
+ * Some functionality is taken from the US Naval Observatory, described here:
+ * https://aa.usno.navy.mil/faq/docs/SunApprox.php
+ *
  * Author: Ryan Seys (https://github.com/ryanseys)
+ * Author: Jay LaPorte (https://github.com/ironwallaby)
  */
 
 'use strict'
@@ -20,39 +26,11 @@ const LAST = 3
 const PHASE_MASK = 3
 
 // Astronomical Constants
-// JDN stands for Julian Day Number
-// Angles here are in degrees
-// 1980 January 0.0 in JDN
-// XXX: DateTime(1980).jdn yields 2444239.5 -- which one is right?
-// XXX: even though 2444239.5 is correct for the 1 Jan 1980, 2444238.5 gives
-// better accuracy results... possibly somebody chose all of the below
-// constants based on the wrong epoch?
-const EPOCH = 2444238.5
-
-// Ecliptic longitude of the Sun at epoch 1980.0
-const ECLIPTIC_LONGITUDE_EPOCH = 278.833540
-
-// Ecliptic longitude of the Sun at perigee
-const ECLIPTIC_LONGITUDE_PERIGEE = 282.596403
-
-// Eccentricity of Earth's orbit
-const ECCENTRICITY = 0.016718
-
 // Semi-major axis of Earth's orbit, in kilometers
 const SUN_SMAXIS = 1.49585e8
 
-// Sun's angular size, in degrees, at semi-major axis distance
-const SUN_ANGULAR_SIZE_SMAXIS = 0.533128
-
-// Elements of the Moon's orbit, epoch 1980.0
-// Moon's mean longitude at the epoch
-const MOON_MEAN_LONGITUDE_EPOCH = 64.975464
-
-// Mean longitude of the perigee at the epoch
-const MOON_MEAN_PERIGEE_EPOCH = 349.383063
-
-// Eccentricity of the Moon's orbit
-const MOON_ECCENTRICITY = 0.054900
+// SUN_SMAXIS premultiplied by the angular size of the Sun from the Earth
+const SUN_ANGULAR_SIZE_SMAXIS = SUN_SMAXIS * 0.533128
 
 // Semi-major axis of the Moon's orbit, in kilometers
 const MOON_SMAXIS = 384401.0
@@ -63,10 +41,6 @@ const MOON_ANGULAR_SIZE_SMAXIS = MOON_SMAXIS * 0.5181
 // Synodic month (new Moon to new Moon), in days
 const SYNODIC_MONTH = 29.53058868
 
-function fixangle (a) {
-  return a - 360.0 * Math.floor(a / 360.0)
-}
-
 /**
  * Convert degrees to radians
  * @param  {Number} d Angle in degrees
@@ -81,10 +55,6 @@ function torad (d) {
  * @param  {Number} r Angle in radians
  * @return {Number}   Angle in degrees
  */
-function todeg (r) {
-  return (180.0 / Math.PI) * r
-}
-
 function dsin (d) {
   return Math.sin(torad(d))
 }
@@ -94,98 +64,87 @@ function dcos (d) {
 }
 
 /**
- * Solve the equation of Kepler.
+ * Convert astronomical units to kilometers
+ * @param  {Number} au Distance in astronomical units
+ * @return {Number}    Distance in kilometers
  */
-function kepler (m, ecc) {
-  const epsilon = 1e-6
-
-  m = torad(m)
-  let e = m
-  while (1) {
-    const delta = e - ecc * Math.sin(e) - m
-    e -= delta / (1.0 - ecc * Math.cos(e))
-
-    if (Math.abs(delta) <= epsilon) {
-      break
-    }
-  }
-
-  return e
+function tokm (au) {
+  return 149597870.700 * au
 }
 
 /**
  * Finds the phase information for specific date.
- * @param  {Date} phase_date Date to get phase information of.
- * @return {Object}          Phase data
+ * @param  {Date}   date Date to get phase information of.
+ * @return {Object}      Phase data
  */
-function phase (phase_date) {
-  if (!phase_date) {
-    phase_date = new Date()
+function phase (date) {
+  if (!date) {
+    date = new Date()
   }
-  phase_date = julian.fromDate(phase_date)
-
-  const day = phase_date - EPOCH
-
-  // calculate sun position
-  const sun_mean_anomaly =
-    (360.0 / 365.2422) * day +
-    (ECLIPTIC_LONGITUDE_EPOCH - ECLIPTIC_LONGITUDE_PERIGEE)
-  const sun_true_anomaly =
-    2 * todeg(Math.atan(
-      Math.sqrt((1.0 + ECCENTRICITY) / (1.0 - ECCENTRICITY)) *
-      Math.tan(0.5 * kepler(sun_mean_anomaly, ECCENTRICITY))
-    ))
-  const sun_ecliptic_longitude =
-    ECLIPTIC_LONGITUDE_PERIGEE + sun_true_anomaly
-  const sun_orbital_distance_factor =
-    (1 + ECCENTRICITY * dcos(sun_true_anomaly)) /
-    (1 - ECCENTRICITY * ECCENTRICITY)
-
-  // calculate moon position
-  const moon_mean_longitude =
-    MOON_MEAN_LONGITUDE_EPOCH + 13.1763966 * day
-  const moon_mean_anomaly =
-    moon_mean_longitude - 0.1114041 * day - MOON_MEAN_PERIGEE_EPOCH
-  const moon_evection =
-    1.2739 * dsin(
-      2 * (moon_mean_longitude - sun_ecliptic_longitude) - moon_mean_anomaly
-    )
-  const moon_annual_equation =
-    0.1858 * dsin(sun_mean_anomaly)
-  // XXX: what is the proper name for this value?
-  const moon_mp =
-    moon_mean_anomaly +
-    moon_evection -
-    moon_annual_equation -
-    0.37 * dsin(sun_mean_anomaly)
-  const moon_equation_center_correction =
-    6.2886 * dsin(moon_mp)
-  const moon_corrected_longitude =
-    moon_mean_longitude +
-    moon_evection +
-    moon_equation_center_correction -
-    moon_annual_equation +
-    0.214 * dsin(2.0 * moon_mp)
-  const moon_age =
-    fixangle(
-      moon_corrected_longitude -
-      sun_ecliptic_longitude +
-      0.6583 * dsin(
-        2 * (moon_corrected_longitude - sun_ecliptic_longitude)
-      )
-    )
-  const moon_distance =
-    (MOON_SMAXIS * (1.0 - MOON_ECCENTRICITY * MOON_ECCENTRICITY)) /
-    (1.0 + MOON_ECCENTRICITY * dcos(moon_mp + moon_equation_center_correction))
+
+  // t is the time in "Julian centuries" of 36525 days starting from 2000-01-01
+  // (Note the 0.3 is because the UNIX epoch is on 1970-01-01 and that's 3/10th
+  // of a century before 2000-01-01, which I find cute :3 )
+  const t = date.getTime() * (1 / 3155760000000) - 0.3
+
+  // lunar mean elongation (Astronomical Algorithms, 2nd ed., p. 338)
+  const d = 297.8501921 + t * (445267.1114034 + t * (-0.0018819 + t * ((1 / 545868) + t * (-1 / 113065000))))
+
+  // solar mean anomaly (p. 338)
+  const m = 357.5291092 + t * (35999.0502909 + t * (-0.0001536 + t * (1 / 24490000)))
+
+  // lunar mean anomaly (p. 338)
+  const n = 134.9633964 + t * (477198.8675055 + t * (0.0087414 + t * ((1 / 69699) + t * (-1 / 14712000))))
+
+  // derive sines and cosines necessary for the below calculations
+  const sind = dsin(d)
+  const sinm = dsin(m)
+  const sinn = dsin(n)
+  const cosd = dcos(d)
+  const cosm = dcos(m)
+  const cosn = dcos(n)
+
+  // use trigonometric identities to derive the remainder of the sines and
+  // cosines we need. this reduces us from needing 14 sin/cos calls to only 7,
+  // and makes the lunar distance and solar distance essentially free.
+  // http://mathworld.wolfram.com/Double-AngleFormulas.html
+  // http://mathworld.wolfram.com/TrigonometricAdditionFormulas.html
+  const sin2d = 2 * sind * cosd // sin(2d)
+  const sin2n = 2 * sinn * cosn // sin(2n)
+  const cos2d = 2 * cosd * cosd - 1 // cos(2d)
+  const cos2m = 2 * cosm * cosm - 1 // cos(2m)
+  const cos2n = 2 * cosn * cosn - 1 // cos(2n)
+  const sin2dn = sin2d * cosn - cos2d * sinn // sin(2d - n)
+  const cos2dn = cos2d * cosn + sin2d * sinn // cos(2d - n)
+
+  // lunar phase angle (p. 346)
+  const i = 180 - d - 6.289 * sinn + 2.100 * sinm - 1.274 * sin2dn - 0.658 * sin2d - 0.214 * sin2n - 0.110 * sind
+
+  // fractional illumination (p. 345)
+  const illumination = dcos(i) * 0.5 + 0.5
+
+  // fractional lunar phase
+  let phase = 0.5 - i * (1 / 360)
+  phase -= Math.floor(phase)
+
+  // lunar distance (p. 339-342)
+  // XXX: the book is not clear on how many terms to use for a given level of
+  // accuracy! I used all the easy ones that we already have for free above,
+  // but I imagine this is more than necessary...
+  const moonDistance = 385000.56 - 20905.355 * cosn - 3699.111 * cos2dn - 2955.968 * cos2d - 569.925 * cos2n + 108.743 * cosd
+
+  // solar distance
+  // https://aa.usno.navy.mil/faq/docs/SunApprox.php
+  const sunDistance = tokm(1.00014 - 0.01671 * cosm - 0.00014 * cos2m)
 
   return {
-    phase: (1.0 / 360.0) * moon_age,
-    illuminated: 0.5 * (1.0 - dcos(moon_age)),
-    age: (SYNODIC_MONTH / 360.0) * moon_age,
-    distance: moon_distance,
-    angular_diameter: MOON_ANGULAR_SIZE_SMAXIS / moon_distance,
-    sun_distance: SUN_SMAXIS / sun_orbital_distance_factor,
-    sun_angular_diameter: SUN_ANGULAR_SIZE_SMAXIS * sun_orbital_distance_factor
+    phase: phase,
+    illuminated: illumination,
+    age: phase * SYNODIC_MONTH,
+    distance: moonDistance,
+    angular_diameter: MOON_ANGULAR_SIZE_SMAXIS / moonDistance,
+    sun_distance: sunDistance,
+    sun_angular_diameter: SUN_ANGULAR_SIZE_SMAXIS / sunDistance
   }
 }
 
@@ -201,8 +160,8 @@ function phase (phase_date) {
  */
 function meanphase (sdate, k) {
   // Time in Julian centuries from 1900 January 12 noon UTC
-  const delta_t = (sdate - -2208945600000.0) / 86400000.0
-  const t = delta_t / 36525
+  const delta = (sdate - -2208945600000.0) / 86400000.0
+  const t = delta / 36525
   return 2415020.75933 +
     SYNODIC_MONTH * k +
     (0.0001178 - 0.000000155 * t) * t * t +
@@ -298,7 +257,7 @@ function truephase (k, tphase) {
  * @param  {Date} sdate Date to start hunting from (defaults to current date)
  * @return {Object}     Object containing recent past and future phases
  */
-function phase_hunt (sdate) {
+function phaseHunt (sdate) {
   if (!sdate) {
     sdate = new Date()
   }
@@ -328,7 +287,7 @@ function phase_hunt (sdate) {
   }
 }
 
-function phase_range (start, end, phase) {
+function phaseRange (start, end, phase) {
   start = start.getTime()
   end = end.getTime()
 
@@ -361,5 +320,5 @@ exports.PHASE_FIRST = FIRST
 exports.PHASE_FULL = FULL
 exports.PHASE_LAST = LAST
 exports.phase = phase
-exports.phase_hunt = phase_hunt
-exports.phase_range = phase_range
+exports.phase_hunt = phaseHunt
+exports.phase_range = phaseRange

package.json

@@ -24,9 +24,9 @@
   "main": "lib/lune.js",
   "license": "Apache-2.0",
   "devDependencies": {
-    "chai": "~3.5.0",
-    "mocha": "~2.4.5",
-    "standard": "^6.0.8"
+    "chai": "~4.2.0",
+    "mocha": "~6.1.4",
+    "standard": "^12.0.1"
   },
   "scripts": {
     "test": "standard && mocha --reporter min"

test/index.js

@@ -46,18 +46,32 @@ describe('lune', function () {
   ]
 
   describe('#phase()', function () {
-    it('should return expected values for Feb 17th data', function () {
+    it('should return expected values for 2014-02-17', function () {
       const phase = lune.phase(new Date('2014-02-17T00:00-0500'))
 
       assert.closeTo(phase.phase, 0.568, 0.001)
       assert.closeTo(phase.illuminated, 0.955, 0.001)
       assert.closeTo(phase.age, 16.779, 0.030)
-      assert.closeTo(phase.distance, 396084.5, 384.4)
-      assert.closeTo(phase.angular_diameter, 0.5028, 0.0005)
       assert.closeTo(phase.sun_distance, 147822500, 149600)
       assert.closeTo(phase.sun_angular_diameter, 0.5395, 0.0005)
     })
 
+    // Astronomical Algorithms, 2nd ed., p. 347
+    it('should return expected values for 1992-04-12', function () {
+      const phase = lune.phase(new Date('1992-04-12'))
+
+      assert.closeTo(phase.illuminated, 0.6802, 0.0001)
+      assert.closeTo(phase.distance, 368405, 400)
+      assert.closeTo(phase.angular_diameter, 0.5405, 0.0005)
+    })
+
+    // Astronomical Algorithms, 2nd ed., p. 169
+    it('should return expected values for 1992-10-13', function () {
+      const phase = lune.phase(new Date('1992-10-13'))
+
+      assert.closeTo(phase.sun_distance, 149240000, 10000)
+    })
+
     // http://bazaar.launchpad.net/~keturn/py-moon-phase/trunk/view/head:/moontest.py
     it('should be accurate to astronomical observations', function () {
       let error = 0