Heliocentric coordinates from Kepler orbit and time

hcr.js

var moment = require('moment');

module.exports = function(/* orbit, [Date, time String], [format String], [etc] */){

  // Takes an orbit definition of the kind used in GEOF stellar-objects and
  // remaining arguments as moment arguments, returning heliocentric
  // rectangular coordinates in α's unit (probably AU).

  // TODO: Does not yet take into account changes in orbit over time.
  // TODO: Converts between degrees and radians too often, I think.
  // TODO: May not actually need to keep longitudes between 0 and 2π||360.
  // TODO: Might be useful to get spherical coordinates (λ,β,Δ?) too.

  // TODO: Eventually: Rotation, precession, and axial tilt.
  // TODO: Eventually: Can this do satellites?

  // Technique based upon http://www.stargazing.net/kepler/ellipse.html

  var args = Array.prototype.slice.call(arguments)
    , orbit = args.shift();

  var π = Math.PI
    , rad_deg = π / 180 // rad per deg
    , deg_rad = 180 / π // deg per rad
    , ms_d = 24 * 60 * 60 * 1000 // ms per day
    , cos = function(deg){return Math.cos(rad_deg * deg);}
    , sin = function(deg){return Math.sin(rad_deg * deg);}
    , pow = Math.pow
    , t = moment.apply(this, args).diff(moment(orbit._0.epoch), 'milliseconds') / ms_d
    , α = orbit.α
    , e = orbit.e
    , i = orbit.i
    , Ω = orbit.Ω
    , ϖ = orbit.ϖ
    , P = orbit.orbital_period
    , M_0 = orbit._0.L - ϖ
    ;

  // find the mean motion
  var n = 360 / P;

  // find the mean anomaly
  var M = (M_0 + n * t) % 360;
  M = M < 0 ? M + 360 : M;

  // find the true anomaly
  // this is an approximation using equation of the center
  var ν = M + deg_rad * [ (2 * e - pow(e,3)/4) * sin(M)
    + 5/4 * pow(e,2) * sin(2*M)
    + 13/12 * pow(e,3) * sin(3*M) ];

  // find the radius vector
  var r = α * (1 - pow(e,2)) / [1 + e * cos(ν)];

  var λ = (ν + ϖ - Ω) % 360;
  λ = λ < 0 ? λ + 360 : λ;

  var x = r * [cos(Ω) * cos(λ) - sin(Ω) * sin(λ) * cos(i)]
    , y = r * [sin(Ω) * cos(λ) + cos(Ω) * sin(λ) * cos(i)]
    , z = r * [sin(λ) * sin(i)]
    ;

  // var accuracy = (Math.sqrt(Math.pow(x,2) + Math.pow(y,2) + Math.pow(z,2)) - r);

  return {x: x, y: y, z: z};

};