1 /*
  2  * Astro.js - Static functions to support astronomical calculations
  3  *
  4  * Copyright © 2014-2015, 2018, JEDLSoft
  5  *
  6  * Licensed under the Apache License, Version 2.0 (the "License");
  7  * you may not use this file except in compliance with the License.
  8  * You may obtain a copy of the License at
  9  *
 10  *     http://www.apache.org/licenses/LICENSE-2.0
 11  *
 12  * Unless required by applicable law or agreed to in writing, software
 13  * distributed under the License is distributed on an "AS IS" BASIS,
 14  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 15  *
 16  * See the License for the specific language governing permissions and
 17  * limitations under the License.
 18  */
 19 
 20 // !data astro
 21 
 22 /*
 23  * These routines were derived from a public domain set of JavaScript
 24  * functions for positional astronomy by John Walker of Fourmilab,
 25  * September 1999.
 26  */
 27 
 28 var ilib = require("../index.js");
 29 var Utils = require("./Utils.js");
 30 var SearchUtils = require("./SearchUtils.js");
 31 
 32 var GregorianDate = require("./GregorianDate.js");
 33 var RataDie = require("./RataDie.js");
 34 var GregRataDie = require("./GregRataDie.js");
 35 
 36 var Astro = {};
 37 
 38 /**
 39  * Load in all the data needed for astrological calculations.
 40  *
 41  * @private
 42  * @param {boolean} sync
 43  * @param {*} loadParams
 44  * @param {function(*)|undefined} callback
 45  */
 46 Astro.initAstro = function(sync, loadParams, callback) {
 47     if (!ilib.data.astro) {
 48         Utils.loadData({
 49             object: "Astro",
 50             name: "astro.json", // countries in their own language
 51             locale: "-", // only need to load the root file
 52             nonLocale: true,
 53             sync: sync,
 54             loadParams: loadParams,
 55             callback: ilib.bind(this, function(astroData) {
 56                 /**
 57                  * @type {{
 58                  *      _EquinoxpTerms:Array.<number>,
 59                  *      _JDE0tab1000:Array.<number>,
 60                  *      _JDE0tab2000:Array.<number>,
 61                  *      _deltaTtab:Array.<number>,
 62                  *      _oterms:Array.<number>,
 63                  *      _nutArgMult:Array.<number>,
 64                  *      _nutArgCoeff:Array.<number>,
 65                  *      _nutCoeffA:Array.<number>,
 66                  *      _nutCoeffB:Array.<number>,
 67                  *      _coeff19th:Array.<number>,
 68                  *      _coeff18th:Array.<number>,
 69                  *      _solarLongCoeff:Array.<number>,
 70                  *      _solarLongMultipliers:Array.<number>,
 71                  *      _solarLongAddends:Array.<number>,
 72                  *      _meanMoonCoeff:Array.<number>,
 73                  *      _elongationCoeff:Array.<number>,
 74                  *      _solarAnomalyCoeff:Array.<number>,
 75                  *      _lunarAnomalyCoeff:Array.<number>,
 76                  *      _moonFromNodeCoeff:Array.<number>,
 77                  *      _eCoeff:Array.<number>,
 78                  *      _lunarElongationLongCoeff:Array.<number>,
 79                  *      _solarAnomalyLongCoeff:Array.<number>,
 80                  *      _lunarAnomalyLongCoeff:Array.<number>,
 81                  *      _moonFromNodeLongCoeff:Array.<number>,
 82                  *      _sineCoeff:Array.<number>,
 83                  *      _nmApproxCoeff:Array.<number>,
 84                  *      _nmCapECoeff:Array.<number>,
 85                  *      _nmSolarAnomalyCoeff:Array.<number>,
 86                  *      _nmLunarAnomalyCoeff:Array.<number>,
 87                  *      _nmMoonArgumentCoeff:Array.<number>,
 88                  *      _nmCapOmegaCoeff:Array.<number>,
 89                  *      _nmEFactor:Array.<number>,
 90                  *      _nmSolarCoeff:Array.<number>,
 91                  *      _nmLunarCoeff:Array.<number>,
 92                  *      _nmMoonCoeff:Array.<number>,
 93                  *      _nmSineCoeff:Array.<number>,
 94                  *      _nmAddConst:Array.<number>,
 95                  *      _nmAddCoeff:Array.<number>,
 96                  *      _nmAddFactor:Array.<number>,
 97                  *      _nmExtra:Array.<number>
 98                  *  }}
 99                  */
100                  ilib.data.astro = astroData;
101                 if (callback && typeof(callback) === 'function') {
102                     callback(astroData);
103                 }
104             })
105         });
106     } else {
107         if (callback && typeof(callback) === 'function') {
108             callback(ilib.data.astro);
109         }
110     }
111 };
112 
113 /**
114  * Convert degrees to radians.
115  *
116  * @static
117  * @protected
118  * @param {number} d angle in degrees
119  * @return {number} angle in radians
120  */
121 Astro._dtr = function(d) {
122     return (d * Math.PI) / 180.0;
123 };
124 
125 /**
126  * Convert radians to degrees.
127  *
128  * @static
129  * @protected
130  * @param {number} r angle in radians
131  * @return {number} angle in degrees
132  */
133 Astro._rtd = function(r) {
134     return (r * 180.0) / Math.PI;
135 };
136 
137 /**
138  * Return the cosine of an angle given in degrees.
139  * @static
140  * @protected
141  * @param {number} d angle in degrees
142  * @return {number} cosine of the angle.
143  */
144 Astro._dcos = function(d) {
145     return Math.cos(Astro._dtr(d));
146 };
147 
148 /**
149  * Return the sine of an angle given in degrees.
150  * @static
151  * @protected
152  * @param {number} d angle in degrees
153  * @return {number} sine of the angle.
154  */
155 Astro._dsin = function(d) {
156     return Math.sin(Astro._dtr(d));
157 };
158 
159 /**
160  * Return the tan of an angle given in degrees.
161  * @static
162  * @protected
163  * @param {number} d angle in degrees
164  * @return {number} tan of the angle.
165  */
166 Astro._dtan = function(d) {
167     return Math.tan(Astro._dtr(d));
168 };
169 
170 /**
171  * Range reduce angle in degrees.
172  *
173  * @static
174  * @param {number} a angle to reduce
175  * @return {number} the reduced angle
176  */
177 Astro._fixangle = function(a) {
178     return a - 360.0 * (Math.floor(a / 360.0));
179 };
180 
181 /**
182  * Range reduce angle in radians.
183  *
184  * @static
185  * @protected
186  * @param {number} a angle to reduce
187  * @return {number} the reduced angle
188  */
189 Astro._fixangr = function(a) {
190     return a - (2 * Math.PI) * (Math.floor(a / (2 * Math.PI)));
191 };
192 
193 /**
194  * Determine the Julian Ephemeris Day of an equinox or solstice.  The "which"
195  * argument selects the item to be computed:
196  *
197  * <ul>
198  * <li>0   March equinox
199  * <li>1   June solstice
200  * <li>2   September equinox
201  * <li>3   December solstice
202  * </ul>
203  *
204  * @static
205  * @protected
206  * @param {number} year Gregorian year to calculate for
207  * @param {number} which Which equinox or solstice to calculate
208  */
209 Astro._equinox = function(year, which) {
210     var deltaL, i, j, JDE0, JDE, JDE0tab, S, T, W, Y;
211 
212     /*  Initialize terms for mean equinox and solstices.  We
213         have two sets: one for years prior to 1000 and a second
214         for subsequent years.  */
215 
216     if (year < 1000) {
217         JDE0tab = ilib.data.astro._JDE0tab1000;
218         Y = year / 1000;
219     } else {
220         JDE0tab = ilib.data.astro._JDE0tab2000;
221         Y = (year - 2000) / 1000;
222     }
223 
224     JDE0 = JDE0tab[which][0] + (JDE0tab[which][1] * Y)
225             + (JDE0tab[which][2] * Y * Y) + (JDE0tab[which][3] * Y * Y * Y)
226             + (JDE0tab[which][4] * Y * Y * Y * Y);
227 
228     //document.debug.log.value += "JDE0 = " + JDE0 + "\n";
229 
230     T = (JDE0 - 2451545.0) / 36525;
231     //document.debug.log.value += "T = " + T + "\n";
232     W = (35999.373 * T) - 2.47;
233     //document.debug.log.value += "W = " + W + "\n";
234     deltaL = 1 + (0.0334 * Astro._dcos(W)) + (0.0007 * Astro._dcos(2 * W));
235     //document.debug.log.value += "deltaL = " + deltaL + "\n";
236 
237     //  Sum the periodic terms for time T
238 
239     S = 0;
240     j = 0;
241     for (i = 0; i < 24; i++) {
242         S += ilib.data.astro._EquinoxpTerms[j]
243                 * Astro._dcos(ilib.data.astro._EquinoxpTerms[j + 1] + (ilib.data.astro._EquinoxpTerms[j + 2] * T));
244         j += 3;
245     }
246 
247     //document.debug.log.value += "S = " + S + "\n";
248     //document.debug.log.value += "Corr = " + ((S * 0.00001) / deltaL) + "\n";
249 
250     JDE = JDE0 + ((S * 0.00001) / deltaL);
251 
252     return JDE;
253 };
254 
255 /*
256  * The table of observed Delta T values at the beginning of
257  * years from 1620 through 2014 as found in astro.json is taken from
258  * http://www.staff.science.uu.nl/~gent0113/deltat/deltat.htm
259  * and
260  * ftp://maia.usno.navy.mil/ser7/deltat.data
261  */
262 
263 /**
264  * Determine the difference, in seconds, between dynamical time and universal time.
265  *
266  * @static
267  * @protected
268  * @param {number} year to calculate the difference for
269  * @return {number} difference in seconds between dynamical time and universal time
270  */
271 Astro._deltat = function (year) {
272     var dt, f, i, t;
273 
274     if ((year >= 1620) && (year <= 2014)) {
275         i = Math.floor(year - 1620);
276         f = (year - 1620) - i; /* Fractional part of year */
277         dt = ilib.data.astro._deltaTtab[i] + ((ilib.data.astro._deltaTtab[i + 1] - ilib.data.astro._deltaTtab[i]) * f);
278     } else {
279         t = (year - 2000) / 100;
280         if (year < 948) {
281             dt = 2177 + (497 * t) + (44.1 * t * t);
282         } else {
283             dt = 102 + (102 * t) + (25.3 * t * t);
284             if ((year > 2000) && (year < 2100)) {
285                 dt += 0.37 * (year - 2100);
286             }
287         }
288     }
289     return dt;
290 };
291 
292 /**
293  * Calculate the obliquity of the ecliptic for a given
294  * Julian date.  This uses Laskar's tenth-degree
295  * polynomial fit (J. Laskar, Astronomy and
296  * Astrophysics, Vol. 157, page 68 [1986]) which is
297  * accurate to within 0.01 arc second between AD 1000
298  * and AD 3000, and within a few seconds of arc for
299  * +/-10000 years around AD 2000.  If we're outside the
300  * range in which this fit is valid (deep time) we
301  * simply return the J2000 value of the obliquity, which
302  * happens to be almost precisely the mean.
303  *
304  * @static
305  * @protected
306  * @param {number} jd Julian Day to calculate the obliquity for
307  * @return {number} the obliquity
308  */
309 Astro._obliqeq = function (jd) {
310     var eps, u, v, i;
311 
312      v = u = (jd - 2451545.0) / 3652500.0;
313 
314      eps = 23 + (26 / 60.0) + (21.448 / 3600.0);
315 
316      if (Math.abs(u) < 1.0) {
317          for (i = 0; i < 10; i++) {
318              eps += (ilib.data.astro._oterms[i] / 3600.0) * v;
319              v *= u;
320          }
321      }
322      return eps;
323 };
324 
325 /**
326  * Return the position of the sun.  We return
327  * intermediate values because they are useful in a
328  * variety of other contexts.
329  * @static
330  * @protected
331  * @param {number} jd find the position of sun on this Julian Day
332  * @return {Object} the position of the sun and many intermediate
333  * values
334  */
335 Astro._sunpos = function(jd) {
336     var ret = {},
337         T, T2, T3, Omega, epsilon, epsilon0;
338 
339     T = (jd - 2451545.0) / 36525.0;
340     //document.debug.log.value += "Sunpos.  T = " + T + "\n";
341     T2 = T * T;
342     T3 = T * T2;
343     ret.meanLongitude = Astro._fixangle(280.46646 + 36000.76983 * T + 0.0003032 * T2);
344     //document.debug.log.value += "ret.meanLongitude = " + ret.meanLongitude + "\n";
345     ret.meanAnomaly = Astro._fixangle(357.52911 + (35999.05029 * T) - 0.0001537 * T2 - 0.00000048 * T3);
346     //document.debug.log.value += "ret.meanAnomaly = " + ret.meanAnomaly + "\n";
347     ret.eccentricity = 0.016708634 - 0.000042037 * T - 0.0000001267 * T2;
348     //document.debug.log.value += "e = " + e + "\n";
349     ret.equationOfCenter = ((1.914602 - 0.004817 * T - 0.000014 * T2) * Astro._dsin(ret.meanAnomaly))
350             + ((0.019993 - 0.000101 * T) * Astro._dsin(2 * ret.meanAnomaly))
351             + (0.000289 * Astro._dsin(3 * ret.meanAnomaly));
352     //document.debug.log.value += "ret.equationOfCenter = " + ret.equationOfCenter + "\n";
353     ret.sunLongitude = ret.meanLongitude + ret.equationOfCenter;
354     //document.debug.log.value += "ret.sunLongitude = " + ret.sunLongitude + "\n";
355     //ret.sunAnomaly = ret.meanAnomaly + ret.equationOfCenter;
356     //document.debug.log.value += "ret.sunAnomaly = " + ret.sunAnomaly + "\n";
357     // ret.sunRadius = (1.000001018 * (1 - (ret.eccentricity * ret.eccentricity))) / (1 + (ret.eccentricity * Astro._dcos(ret.sunAnomaly)));
358     //document.debug.log.value += "ret.sunRadius = " + ret.sunRadius + "\n";
359     Omega = 125.04 - (1934.136 * T);
360     //document.debug.log.value += "Omega = " + Omega + "\n";
361     ret.apparentLong = ret.sunLongitude + (-0.00569) + (-0.00478 * Astro._dsin(Omega));
362     //document.debug.log.value += "ret.apparentLong = " + ret.apparentLong + "\n";
363     epsilon0 = Astro._obliqeq(jd);
364     //document.debug.log.value += "epsilon0 = " + epsilon0 + "\n";
365     epsilon = epsilon0 + (0.00256 * Astro._dcos(Omega));
366     //document.debug.log.value += "epsilon = " + epsilon + "\n";
367     //ret.rightAscension = Astro._fixangle(Astro._rtd(Math.atan2(Astro._dcos(epsilon0) * Astro._dsin(ret.sunLongitude), Astro._dcos(ret.sunLongitude))));
368     //document.debug.log.value += "ret.rightAscension = " + ret.rightAscension + "\n";
369     // ret.declination = Astro._rtd(Math.asin(Astro._dsin(epsilon0) * Astro._dsin(ret.sunLongitude)));
370     ////document.debug.log.value += "ret.declination = " + ret.declination + "\n";
371     ret.inclination = Astro._fixangle(23.4392911 - 0.013004167 * T - 0.00000016389 * T2 + 0.0000005036 * T3);
372     ret.apparentRightAscension = Astro._fixangle(Astro._rtd(Math.atan2(Astro._dcos(epsilon) * Astro._dsin(ret.apparentLong), Astro._dcos(ret.apparentLong))));
373     //document.debug.log.value += "ret.apparentRightAscension = " + ret.apparentRightAscension + "\n";
374     //ret.apparentDeclination = Astro._rtd(Math.asin(Astro._dsin(epsilon) * Astro._dsin(ret.apparentLong)));
375     //document.debug.log.value += "ret.apparentDecliation = " + ret.apparentDecliation + "\n";
376 
377     // Angular quantities are expressed in decimal degrees
378     return ret;
379 };
380 
381 /**
382  * Calculate the nutation in longitude, deltaPsi, and obliquity,
383  * deltaEpsilon for a given Julian date jd. Results are returned as an object
384  * giving deltaPsi and deltaEpsilon in degrees.
385  *
386  * @static
387  * @protected
388  * @param {number} jd calculate the nutation of this Julian Day
389  * @return {Object} the deltaPsi and deltaEpsilon of the nutation
390  */
391 Astro._nutation = function(jd) {
392     var i, j,
393         t = (jd - 2451545.0) / 36525.0,
394         t2, t3, to10,
395         ta = [],
396         dp = 0,
397         de = 0,
398         ang,
399         ret = {};
400 
401     t3 = t * (t2 = t * t);
402 
403     /*
404      * Calculate angles. The correspondence between the elements of our array
405      * and the terms cited in Meeus are:
406      *
407      * ta[0] = D ta[0] = M ta[2] = M' ta[3] = F ta[4] = \Omega
408      *
409      */
410 
411     ta[0] = Astro._dtr(297.850363 + 445267.11148 * t - 0.0019142 * t2 + t3 / 189474.0);
412     ta[1] = Astro._dtr(357.52772 + 35999.05034 * t - 0.0001603 * t2 - t3 / 300000.0);
413     ta[2] = Astro._dtr(134.96298 + 477198.867398 * t + 0.0086972 * t2 + t3 / 56250.0);
414     ta[3] = Astro._dtr(93.27191 + 483202.017538 * t - 0.0036825 * t2 + t3 / 327270);
415     ta[4] = Astro._dtr(125.04452 - 1934.136261 * t + 0.0020708 * t2 + t3 / 450000.0);
416 
417     /*
418      * Range reduce the angles in case the sine and cosine functions don't do it
419      * as accurately or quickly.
420      */
421 
422     for (i = 0; i < 5; i++) {
423         ta[i] = Astro._fixangr(ta[i]);
424     }
425 
426     to10 = t / 10.0;
427     for (i = 0; i < 63; i++) {
428         ang = 0;
429         for (j = 0; j < 5; j++) {
430             if (ilib.data.astro._nutArgMult[(i * 5) + j] != 0) {
431                 ang += ilib.data.astro._nutArgMult[(i * 5) + j] * ta[j];
432             }
433         }
434         dp += (ilib.data.astro._nutArgCoeff[(i * 4) + 0] + ilib.data.astro._nutArgCoeff[(i * 4) + 1] * to10) * Math.sin(ang);
435         de += (ilib.data.astro._nutArgCoeff[(i * 4) + 2] + ilib.data.astro._nutArgCoeff[(i * 4) + 3] * to10) * Math.cos(ang);
436     }
437 
438     /*
439      * Return the result, converting from ten thousandths of arc seconds to
440      * radians in the process.
441      */
442 
443     ret.deltaPsi = dp / (3600.0 * 10000.0);
444     ret.deltaEpsilon = de / (3600.0 * 10000.0);
445 
446     return ret;
447 };
448 
449 /**
450  * Returns the equation of time as a fraction of a day.
451  *
452  * @static
453  * @protected
454  * @param {number} jd the Julian Day of the day to calculate for
455  * @return {number} the equation of time for the given day
456  */
457 Astro._equationOfTime = function(jd) {
458     var alpha, deltaPsi, E, epsilon, L0, tau, pos;
459 
460     // 2451545.0 is the Julian day of J2000 epoch
461     // 365250.0 is the number of days in a Julian millenium
462     tau = (jd - 2451545.0) / 365250.0;
463     //console.log("equationOfTime.  tau = " + tau);
464     L0 = 280.4664567 + (360007.6982779 * tau) + (0.03032028 * tau * tau)
465             + ((tau * tau * tau) / 49931)
466             + (-((tau * tau * tau * tau) / 15300))
467             + (-((tau * tau * tau * tau * tau) / 2000000));
468     //console.log("L0 = " + L0);
469     L0 = Astro._fixangle(L0);
470     //console.log("L0 = " + L0);
471     pos = Astro._sunpos(jd);
472     alpha = pos.apparentRightAscension;
473     //console.log("alpha = " + alpha);
474     var nut = Astro._nutation(jd);
475     deltaPsi = nut.deltaPsi;
476     //console.log("deltaPsi = " + deltaPsi);
477     epsilon = Astro._obliqeq(jd) + nut.deltaEpsilon;
478     //console.log("epsilon = " + epsilon);
479     //console.log("L0 - 0.0057183 = " + (L0 - 0.0057183));
480     //console.log("L0 - 0.0057183 - alpha = " + (L0 - 0.0057183 - alpha));
481     //console.log("deltaPsi * cos(epsilon) = " + deltaPsi * Astro._dcos(epsilon));
482 
483     E = L0 - 0.0057183 - alpha + deltaPsi * Astro._dcos(epsilon);
484     // if alpha and L0 are in different quadrants, then renormalize
485     // so that the difference between them is in the right range
486     if (E > 180) {
487         E -= 360;
488     }
489     //console.log("E = " + E);
490     // E = E - 20.0 * (Math.floor(E / 20.0));
491     E = E * 4;
492     //console.log("Efixed = " + E);
493     E = E / (24 * 60);
494     //console.log("Eday = " + E);
495 
496     return E;
497 };
498 
499 /**
500  * @private
501  * @static
502  */
503 Astro._poly = function(x, coefficients) {
504     var result = coefficients[0];
505     var xpow = x;
506     for (var i = 1; i < coefficients.length; i++) {
507         result += coefficients[i] * xpow;
508         xpow *= x;
509     }
510     return result;
511 };
512 
513 /**
514  * Calculate the UTC RD from the local RD given "zone" number of minutes
515  * worth of offset.
516  *
517  * @static
518  * @protected
519  * @param {number} local RD of the locale time, given in any calendar
520  * @param {number} zone number of minutes of offset from UTC for the time zone
521  * @return {number} the UTC equivalent of the local RD
522  */
523 Astro._universalFromLocal = function(local, zone) {
524     return local - zone / 1440;
525 };
526 
527 /**
528  * Calculate the local RD from the UTC RD given "zone" number of minutes
529  * worth of offset.
530  *
531  * @static
532  * @protected
533  * @param {number} local RD of the locale time, given in any calendar
534  * @param {number} zone number of minutes of offset from UTC for the time zone
535  * @return {number} the UTC equivalent of the local RD
536  */
537 Astro._localFromUniversal = function(local, zone) {
538     return local + zone / 1440;
539 };
540 
541 /**
542  * @private
543  * @static
544  * @param {number} c julian centuries of the date to calculate
545  * @return {number} the aberration
546  */
547 Astro._aberration = function(c) {
548     return 9.74e-05 * Astro._dcos(177.63 + 35999.01847999999 * c) - 0.005575;
549 };
550 
551 /**
552  * @private
553  *
554 ilib.data.astro._nutCoeffA = [124.90, -1934.134, 0.002063];
555 ilib.data.astro._nutCoeffB q= [201.11, 72001.5377, 0.00057];
556 */
557 
558 /**
559  * @private
560  * @static
561  * @param {number} c julian centuries of the date to calculate
562  * @return {number} the nutation for the given julian century in radians
563  */
564 Astro._nutation2 = function(c) {
565     var a = Astro._poly(c, ilib.data.astro._nutCoeffA);
566     var b = Astro._poly(c, ilib.data.astro._nutCoeffB);
567     // return -0.0000834 * Astro._dsin(a) - 0.0000064 * Astro._dsin(b);
568     return -0.004778 * Astro._dsin(a) - 0.0003667 * Astro._dsin(b);
569 };
570 
571 /**
572  * @static
573  * @private
574  */
575 Astro._ephemerisCorrection = function(jd) {
576     var year = GregorianDate._calcYear(jd - 1721424.5);
577 
578     if (1988 <= year && year <= 2019) {
579         return (year - 1933) / 86400;
580     }
581 
582     if (1800 <= year && year <= 1987) {
583         var jul1 = new GregRataDie({
584             year: year,
585             month: 7,
586             day: 1,
587             hour: 0,
588             minute: 0,
589             second: 0
590         });
591         // 693596 is the rd of Jan 1, 1900
592         var theta = (jul1.getRataDie() - 693596) / 36525;
593         return Astro._poly(theta, (1900 <= year) ? ilib.data.astro._coeff19th : ilib.data.astro._coeff18th);
594     }
595 
596     if (1620 <= year && year <= 1799) {
597         year -= 1600;
598         return (196.58333 - 4.0675 * year + 0.0219167 * year * year) / 86400;
599     }
600 
601     // 660724 is the rd of Jan 1, 1810
602     var jan1 = new GregRataDie({
603         year: year,
604         month: 1,
605         day: 1,
606         hour: 0,
607         minute: 0,
608         second: 0
609     });
610     // var x = 0.5 + (jan1.getRataDie() - 660724);
611     var x = 0.5 + (jan1.getRataDie() - 660724);
612 
613     return ((x * x / 41048480) - 15) / 86400;
614 };
615 
616 /**
617  * @static
618  * @private
619  */
620 Astro._ephemerisFromUniversal = function(jd) {
621     return jd + Astro._ephemerisCorrection(jd);
622 };
623 
624 /**
625  * @static
626  * @private
627  */
628 Astro._universalFromEphemeris = function(jd) {
629     return jd - Astro._ephemerisCorrection(jd);
630 };
631 
632 /**
633  * @static
634  * @private
635  */
636 Astro._julianCenturies = function(jd) {
637     // 2451545.0 is the Julian day of J2000 epoch
638     // 730119.5 is the Gregorian RD of J2000 epoch
639     // 36525.0 is the number of days in a Julian century
640     return (Astro._ephemerisFromUniversal(jd) - 2451545.0) / 36525.0;
641 };
642 
643 /**
644  * Calculate the solar longitude
645  *
646  * @static
647  * @protected
648  * @param {number} jd julian day of the date to calculate the longitude for
649  * @return {number} the solar longitude in degrees
650  */
651 Astro._solarLongitude = function(jd) {
652     var c = Astro._julianCenturies(jd),
653         longitude = 0,
654         len = ilib.data.astro._solarLongCoeff.length;
655 
656     for (var i = 0; i < len; i++) {
657         longitude += ilib.data.astro._solarLongCoeff[i] *
658             Astro._dsin(ilib.data.astro._solarLongAddends[i] + ilib.data.astro._solarLongMultipliers[i] * c);
659     }
660     longitude *= 5.729577951308232e-06;
661     longitude += 282.77718340000001 + 36000.769537439999 * c;
662     longitude += Astro._aberration(c) + Astro._nutation2(c);
663     return Astro._fixangle(longitude);
664 };
665 
666 /**
667  * @static
668  * @protected
669  * @param {number} jd
670  * @return {number}
671  */
672 Astro._lunarLongitude = function (jd) {
673     var c = Astro._julianCenturies(jd),
674         meanMoon = Astro._fixangle(Astro._poly(c, ilib.data.astro._meanMoonCoeff)),
675         elongation = Astro._fixangle(Astro._poly(c, ilib.data.astro._elongationCoeff)),
676         solarAnomaly = Astro._fixangle(Astro._poly(c, ilib.data.astro._solarAnomalyCoeff)),
677         lunarAnomaly = Astro._fixangle(Astro._poly(c, ilib.data.astro._lunarAnomalyCoeff)),
678         moonNode = Astro._fixangle(Astro._poly(c, ilib.data.astro._moonFromNodeCoeff)),
679         e = Astro._poly(c, ilib.data.astro._eCoeff);
680 
681     var sum = 0;
682     for (var i = 0; i < ilib.data.astro._lunarElongationLongCoeff.length; i++) {
683         var x = ilib.data.astro._solarAnomalyLongCoeff[i];
684 
685         sum += ilib.data.astro._sineCoeff[i] * Math.pow(e, Math.abs(x)) *
686             Astro._dsin(ilib.data.astro._lunarElongationLongCoeff[i] * elongation + x * solarAnomaly +
687                 ilib.data.astro._lunarAnomalyLongCoeff[i] * lunarAnomaly +
688                 ilib.data.astro._moonFromNodeLongCoeff[i] * moonNode);
689     }
690     var longitude = sum / 1000000;
691     var venus = 3958.0 / 1000000 * Astro._dsin(119.75 + c * 131.84899999999999);
692     var jupiter = 318.0 / 1000000 * Astro._dsin(53.090000000000003 + c * 479264.28999999998);
693     var flatEarth = 1962.0 / 1000000 * Astro._dsin(meanMoon - moonNode);
694 
695     return Astro._fixangle(meanMoon + longitude + venus + jupiter + flatEarth + Astro._nutation2(c));
696 };
697 
698 /**
699  * @static
700  * @protected
701  * @param {number} n
702  * @return {number} julian day of the n'th new moon
703  */
704 Astro._newMoonTime = function(n) {
705     var k = n - 24724;
706     var c = k / 1236.8499999999999;
707     var approx = Astro._poly(c, ilib.data.astro._nmApproxCoeff);
708     var capE = Astro._poly(c, ilib.data.astro._nmCapECoeff);
709     var solarAnomaly = Astro._poly(c, ilib.data.astro._nmSolarAnomalyCoeff);
710     var lunarAnomaly = Astro._poly(c, ilib.data.astro._nmLunarAnomalyCoeff);
711     var moonArgument = Astro._poly(c, ilib.data.astro._nmMoonArgumentCoeff);
712     var capOmega = Astro._poly(c, ilib.data.astro._nmCapOmegaCoeff);
713     var correction = -0.00017 * Astro._dsin(capOmega);
714     var i;
715     for (i = 0; i < ilib.data.astro._nmSineCoeff.length; i++) {
716         correction = correction + ilib.data.astro._nmSineCoeff[i] * Math.pow(capE, ilib.data.astro._nmEFactor[i]) *
717         Astro._dsin(ilib.data.astro._nmSolarCoeff[i] * solarAnomaly +
718                 ilib.data.astro._nmLunarCoeff[i] * lunarAnomaly +
719                 ilib.data.astro._nmMoonCoeff[i] * moonArgument);
720     }
721     var additional = 0;
722     for (i = 0; i < ilib.data.astro._nmAddConst.length; i++) {
723         additional = additional + ilib.data.astro._nmAddFactor[i] *
724         Astro._dsin(ilib.data.astro._nmAddConst[i] + ilib.data.astro._nmAddCoeff[i] * k);
725     }
726     var extra = 0.000325 * Astro._dsin(Astro._poly(c, ilib.data.astro._nmExtra));
727     return Astro._universalFromEphemeris(approx + correction + extra + additional + RataDie.gregorianEpoch);
728 };
729 
730 /**
731  * @static
732  * @protected
733  * @param {number} jd
734  * @return {number}
735  */
736 Astro._lunarSolarAngle = function(jd) {
737     var lunar = Astro._lunarLongitude(jd);
738     var solar = Astro._solarLongitude(jd)
739     return Astro._fixangle(lunar - solar);
740 };
741 
742 /**
743  * @static
744  * @protected
745  * @param {number} jd
746  * @return {number}
747  */
748 Astro._newMoonBefore = function (jd) {
749     var phase = Astro._lunarSolarAngle(jd);
750     // 11.450086114414322 is the julian day of the 0th full moon
751     // 29.530588853000001 is the average length of a month
752     var guess = Math.round((jd - 11.450086114414322 - RataDie.gregorianEpoch) / 29.530588853000001 - phase / 360) - 1;
753     var current, last;
754     current = last = Astro._newMoonTime(guess);
755     while (current < jd) {
756         guess++;
757         last = current;
758         current = Astro._newMoonTime(guess);
759     }
760     return last;
761 };
762 
763 /**
764  * @static
765  * @protected
766  * @param {number} jd
767  * @return {number}
768  */
769 Astro._newMoonAtOrAfter = function (jd) {
770     var phase = Astro._lunarSolarAngle(jd);
771     // 11.450086114414322 is the julian day of the 0th full moon
772     // 29.530588853000001 is the average length of a month
773     var guess = Math.round((jd - 11.450086114414322 - RataDie.gregorianEpoch) / 29.530588853000001 - phase / 360);
774     var current;
775     while ((current = Astro._newMoonTime(guess)) < jd) {
776         guess++;
777     }
778     return current;
779 };
780 
781 /**
782  * @static
783  * @protected
784  * @param {number} jd JD to calculate from
785  * @param {number} longitude longitude to seek
786  * @returns {number} the JD of the next time that the solar longitude
787  * is a multiple of the given longitude
788  */
789 Astro._nextSolarLongitude = function(jd, longitude) {
790     var rate = 365.242189 / 360.0;
791     var tau = jd + rate * Astro._fixangle(longitude - Astro._solarLongitude(jd));
792     var start = Math.max(jd, tau - 5.0);
793     var end = tau + 5.0;
794 
795     return SearchUtils.bisectionSearch(0, start, end, 1e-6, function (l) {
796         return 180 - Astro._fixangle(Astro._solarLongitude(l) - longitude);
797     });
798 };
799 
800 /**
801  * Floor the julian day to midnight of the current julian day.
802  *
803  * @static
804  * @protected
805  * @param {number} jd the julian to round
806  * @return {number} the jd floored to the midnight of the julian day
807  */
808 Astro._floorToJD = function(jd) {
809     return Math.floor(jd - 0.5) + 0.5;
810 };
811 
812 /**
813  * Floor the julian day to midnight of the current julian day.
814  *
815  * @static
816  * @protected
817  * @param {number} jd the julian to round
818  * @return {number} the jd floored to the midnight of the julian day
819  */
820 Astro._ceilToJD = function(jd) {
821     return Math.ceil(jd + 0.5) - 0.5;
822 };
823 
824 module.exports = Astro;
825