astronomia/es/perihelion.js

/* eslint standard/no-callback-literal:0 */
/**
 * @copyright 2013 Sonia Keys
 * @copyright 2016 commenthol
 * @license MIT
 * @module perihelion
 */
/**
 * Perihelion: Chapter 38, Planets in Perihelion and Aphelion.
 *
 * Functions Aphelion and Perihelion implement algorithms from the book
 * to return approximate results.
 *
 * For accurate results, Meeus describes the general technique of
 * interpolating from a precise ephemeris but does not give a complete
 * algorithm.  The algorithm implemented here for Aphelion2 and Perihelion2
 * is to start with the approximate result and then crawl along the curve
 * at the specified time resolution until the desired extremum is found.
 * This algorithm slows down as higher accuracy is demanded.  1 day accuracy
 * is generally quick for planets other than Neptune.
 *
 * Meeus doesn't give an algorithm to handle the double extrema of Neptune.
 * The algorithm here is to pick starting points several years either side
 * of the approximate date and follow the slopes inward.  The consequence of
 * starting farther from the extremum is that these functions are particularly
 * slow for Neptune.  They are offered here though as a simple implementation
 * of Meeus's presentation in the book.
 */
import base from './base';
import interp from './interpolation';

/**
 * Planet constants for first argument of Perihelion and Aphelion functions.
 */
var planets = {};
export var mercury = planets.mercury = 0;
export var venus = planets.venus = 1;
export var earth = planets.earth = 2;
export var mars = planets.mars = 3;
export var jupiter = planets.jupiter = 4;
export var saturn = planets.saturn = 5;
export var uranus = planets.uranus = 6;
export var neptune = planets.neptune = 7;
export var embary = planets.embary = 8;

/**
 * Perihelion returns an approximate jde of the perihelion event nearest the given time.
 *
 * @param {perihelion.NAME} p - planet constant from above
 * @param {Number} y - year number indicating a time near the perihelion event.
 * @returns {Number} jde - time of the event
 */
export function perihelion(p, year) {
  return ap(p, year, false, pf);
}

/**
 * Aphelion returns an approximate jde of the aphelion event nearest the given time.
 *
 * @param {perihelion.NAME} p - planet constant from above
 * @param {Number} y - year number indicating a time near the aphelion event.
 * @returns {Number} jde - time of the event
 */
export function aphelion(p, year) {
  return ap(p, year, true, af);
}

var pf = function pf(x) {
  // (x float64)  float64
  return Math.floor(x + 0.5);
};

var af = function af(x) {
  // (x float64)  float64
  return Math.floor(x) + 0.5;
};

var ap = function ap(p, y, a, f) {
  // (p int, y float64, a bool, f func(float64)  float64) float64
  var i = p;
  if (i === embary) {
    i = earth;
  }
  var k = f(ka[i].a * (y - ka[i].b));
  var j = base.horner(k, c[i]);
  if (p === earth) {
    var _c = ep;
    if (a) {
      _c = ea;
    }
    for (var _i = 0; _i < 5; _i++) {
      j += _c[_i] * Math.sin((ec[_i].a + ec[_i].b * k) * Math.PI / 180);
    }
  }
  return j;
};

var ka = [{ a: 4.15201, b: 2000.12 }, // Mercury
{ a: 1.62549, b: 2000.53 }, // ...
{ a: 0.99997, b: 2000.01 }, { a: 0.53166, b: 2001.78 }, { a: 0.0843, b: 2011.2 }, { a: 0.03393, b: 2003.52 }, { a: 0.0119, b: 2051.1 }, // Neptune
{ a: 0.00607, b: 2047.5 // EMBary
}];

var c = [[2451590.257, 87.96934963], [2451738.233, 224.7008188, -0.0000000327], [2451547.507, 365.2596358, 0.0000000156], [2452195.026, 686.9957857, -0.0000001187], [2455636.936, 4332.897065, 0.0001367], [2452830.12, 10764.21676, 0.000827], [2470213.5, 30694.8767, -0.00541], [2468895.1, 60190.33, 0.03429]];

var ec = [{ a: 328.41, b: 132.788585 }, { a: 316.13, b: 584.903153 }, { a: 346.2, b: 450.380738 }, { a: 136.95, b: 659.306737 }, { a: 249.52, b: 329.653368 }];

var ep = [1.278, -0.055, -0.091, -0.056, -0.045];
var ea = [-1.352, 0.061, 0.062, 0.029, 0.031];

/**
 * Perihelion2 returns the perihelion event nearest the given time.
 *
 * @param {planetposition.Planet} planet - VSOP87 planet (EMBary is not allowed)
 * @param {Number} year - (float) decimal year number near the perihelion event
 * @param {Number} precision - desired precision of the time result, in days
 * @param {Function} [cb] - callback function for asynchronous processing `cb([jde, r])`
 * @returns {Array} [jde, r]
 *   {Number} jde - time of the event
 *   {Number} r - the distance of the planet from the Sun in AU.
 */
export function perihelion2(planet, year, precision, cb) {
  return ap2(planets[planet.name], year, precision, planet, false, pf, cb);
}

/**
 * Aphelion2 returns the aphelion event nearest the given time.
 *
 * @param {planetposition.Planet} planet - VSOP87 planet (EMBary is not allowed)
 * @param {Number} year - (float) decimal year number near the perihelion event
 * @param {Number} precision - desired precision of the time result, in days
 * @param {Function} [cb] - callback function for asynchronous processing `cb([jde, r])`
 * @returns {Array} [jde, r]
 *   {Number} jde - time of the event
 *   {Number} r - the distance of the planet from the Sun in AU.
 */
export function aphelion2(planet, year, precision, cb) {
  return ap2(planets[planet.name], year, precision, planet, true, af, cb);
}

if (typeof setImmediate !== 'function') {
  var _setImmediate = setTimeout; // eslint-disable-line no-unused-vars
}

var ap2 = function ap2(p, y, d, v, a, f, cb) {
  var j1 = ap(p, y, a, f);
  if (p !== neptune) {
    return ap2a(j1, d, a, v, cb);
  }
  // handle the double extrema of Neptune
  if (cb) {
    ap2a(j1 - 5000, d, a, v, function (_ref) {
      var j0 = _ref[0],
          r0 = _ref[1];

      ap2a(j1 + 5000, d, a, v, function (_ref2) {
        var j2 = _ref2[0],
            r2 = _ref2[1];

        if (r0 > r2 === a) {
          cb([j0, r0]);
          return;
        }
        cb([j2, r2]);
      });
    });
  } else {
    var _ap2a = ap2a(j1 - 5000, d, a, v),
        j0 = _ap2a[0],
        r0 = _ap2a[1];

    var _ap2a2 = ap2a(j1 + 5000, d, a, v),
        j2 = _ap2a2[0],
        r2 = _ap2a2[1];

    if (r0 > r2 === a) {
      return [j0, r0];
    }
    return [j2, r2];
  }
};

var ap2a = function ap2a(j1, d, a, v, cb) {
  var j0 = j1 - d;
  var j2 = j1 + d;
  var rr = new Array(3);
  rr[1] = v.position2000(j1).range;
  rr[0] = v.position2000(j0).range;
  rr[2] = v.position2000(j2).range;

  function end() {
    var l = new interp.Len3(j0, j2, rr);

    var _l$extremum = l.extremum(),
        jde = _l$extremum[0],
        r = _l$extremum[1];

    return [jde, r];
  }

  function run() {
    if (a) {
      if (rr[1] > rr[0] && rr[1] > rr[2]) {
        cb && cb(end());
        return true;
      }
    } else {
      if (rr[1] < rr[0] && rr[1] < rr[2]) {
        cb && cb(end());
        return true;
      }
    }
    if (rr[0] < rr[2] === a) {
      j0 = j1;
      j1 = j2;
      j2 += d;
      rr[0] = rr[1];
      rr[1] = rr[2];
      rr[2] = v.position2000(j2).range;
    } else {
      j2 = j1;
      j1 = j0;
      j0 -= d;
      rr[2] = rr[1];
      rr[1] = rr[0];
      rr[0] = v.position2000(j0).range;
    }
    if (cb) {
      setImmediate(run, 0);
    }
  }

  if (cb) {
    run();
  } else {
    for (;;) {
      if (run()) {
        return end();
      }
    }
  }
};

export default {
  mercury: mercury,
  venus: venus,
  earth: earth,
  mars: mars,
  jupiter: jupiter,
  saturn: saturn,
  uranus: uranus,
  neptune: neptune,
  embary: embary,
  perihelion: perihelion,
  aphelion: aphelion,
  perihelion2: perihelion2,
  aphelion2: aphelion2
};