From f754a6de8289270d678a2a9b98d78eb83373d542 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Wed, 29 Apr 2020 21:53:57 -0400 Subject: [PATCH] Fixed #58 - Solar System Barycenter, Earth/Moon Barycenter. Can now calculate the heliocentric Solar System Barycenter (SSB) and Earth/Moon Barycenter (EMB). Changes made in C, C#, JavaScript and Python: Added new body codes SSB, EMB. Added support for calculating both in HelioVector functions. Verified that all calculations match NOVAS. Verified that all calculations match each other across languages. --- demo/nodejs/culminate.js | 2 +- generate/astro_check.js | 2 +- generate/ctest.c | 5 +- generate/dotnet/csharp_test/csharp_test.cs | 5 +- generate/generate.c | 23 +++- generate/mag_test.js | 13 +-- generate/novas_body.h | 1 + generate/template/astronomy.c | 55 ++++++++- generate/template/astronomy.cs | 70 +++++++++++- generate/template/astronomy.js | 39 ++++++- generate/template/astronomy.py | 37 ++++++- generate/test.py | 5 +- generate/vsop/vsop.h | 15 ++- source/c/README.md | 6 +- source/c/astronomy.c | 55 ++++++++- source/c/astronomy.h | 6 +- source/csharp/README.md | 4 +- source/csharp/astronomy.cs | 70 +++++++++++- source/js/README.md | 2 +- source/js/astronomy.js | 39 ++++++- source/js/astronomy.min.js | 123 +++++++++++---------- source/python/README.md | 4 +- source/python/astronomy.py | 37 ++++++- 23 files changed, 514 insertions(+), 104 deletions(-) diff --git a/demo/nodejs/culminate.js b/demo/nodejs/culminate.js index eec8bda4..23ab16a0 100644 --- a/demo/nodejs/culminate.js +++ b/demo/nodejs/culminate.js @@ -42,7 +42,7 @@ function Demo() { console.log('search : ' + date.toISOString()); for (let body of Astronomy.Bodies) { - if (body !== 'Earth') { + if (body !== 'Earth' && body !== 'EMB' && body !== 'SSB') { let culm = Astronomy.SearchHourAngle(body, observer, 0, date); DisplayEvent(body, culm); } diff --git a/generate/astro_check.js b/generate/astro_check.js index 09ca687f..3fa9b2b6 100644 --- a/generate/astro_check.js +++ b/generate/astro_check.js @@ -17,7 +17,7 @@ while (date.tt < stop.tt) { pos = Astronomy.HelioVector(body, date); console.log(`v ${body} ${pos.t.tt.toFixed(16)} ${pos.x.toFixed(16)} ${pos.y.toFixed(16)} ${pos.z.toFixed(16)}`); - if (body !== 'Earth') { + if (body !== 'Earth' && body !== 'EMB' && body !== 'SSB') { j2000 = Astronomy.Equator(body, date, observer, false, false); ofdate = Astronomy.Equator(body, date, observer, true, true); hor = Astronomy.Horizon(date, observer, ofdate.ra, ofdate.dec); diff --git a/generate/ctest.c b/generate/ctest.c index 89c497f6..f8c74b70 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -259,7 +259,8 @@ static int AstroCheck(void) static const astro_body_t bodylist[] = /* match the order in the JavaScript unit test */ { BODY_SUN, BODY_MERCURY, BODY_VENUS, BODY_EARTH, BODY_MARS, - BODY_JUPITER, BODY_SATURN, BODY_URANUS, BODY_NEPTUNE, BODY_PLUTO + BODY_JUPITER, BODY_SATURN, BODY_URANUS, BODY_NEPTUNE, BODY_PLUTO, + BODY_SSB, BODY_EMB }; static int nbodies = sizeof(bodylist) / sizeof(bodylist[0]); @@ -283,7 +284,7 @@ static int AstroCheck(void) CHECK_VECTOR(pos, Astronomy_HelioVector(body, time)); fprintf(outfile, "v %s %0.16lf %0.16lf %0.16lf %0.16lf\n", Astronomy_BodyName(body), pos.t.tt, pos.x, pos.y, pos.z); - if (body != BODY_EARTH) + if (body != BODY_EARTH && body != BODY_EMB && body != BODY_SSB) { CHECK_EQU(j2000, Astronomy_Equator(body, &time, observer, EQUATOR_J2000, NO_ABERRATION)); CHECK_EQU(ofdate, Astronomy_Equator(body, &time, observer, EQUATOR_OF_DATE, ABERRATION)); diff --git a/generate/dotnet/csharp_test/csharp_test.cs b/generate/dotnet/csharp_test/csharp_test.cs index be5ccbd2..fec8c158 100644 --- a/generate/dotnet/csharp_test/csharp_test.cs +++ b/generate/dotnet/csharp_test/csharp_test.cs @@ -103,7 +103,8 @@ namespace csharp_test var bodylist = new Body[] { Body.Sun, Body.Mercury, Body.Venus, Body.Earth, Body.Mars, - Body.Jupiter, Body.Saturn, Body.Uranus, Body.Neptune, Body.Pluto + Body.Jupiter, Body.Saturn, Body.Uranus, Body.Neptune, Body.Pluto, + Body.SSB, Body.EMB }; var observer = new Observer(29.0, -81.0, 10.0); @@ -121,7 +122,7 @@ namespace csharp_test { pos = Astronomy.HelioVector(body, time); outfile.WriteLine("v {0} {1} {2} {3} {4}", body, pos.t.tt.ToString("G17"), pos.x.ToString("G17"), pos.y.ToString("G17"), pos.z.ToString("G17")); - if (body != Body.Earth) + if (body != Body.Earth && body != Body.SSB && body != Body.EMB) { j2000 = Astronomy.Equator(body, time, observer, EquatorEpoch.J2000, Aberration.None); ofdate = Astronomy.Equator(body, time, observer, EquatorEpoch.OfDate, Aberration.Corrected); diff --git a/generate/generate.c b/generate/generate.c index 536c456e..7b02613c 100644 --- a/generate/generate.c +++ b/generate/generate.c @@ -204,7 +204,17 @@ static int NovasBodyPos(double jd, int body, double pos[3]) jed[0] = jd; jed[1] = 0.0; - if (body == BODY_EARTH || body == BODY_MOON) + if (body == BODY_SSB) + { + /* + The NOVAS state() function does everything with respect to + the Solar System Barycenter (SSB) already. + Below we will get the Sun's position relative to the SSB. + Negating that will result in the SSB position with respect to the Sun. + */ + pos[0] = pos[1] = pos[2] = 0.0; + } + else if (body == BODY_EARTH || body == BODY_MOON) { double factor; @@ -913,9 +923,7 @@ static int PositionArcminError(int body, double jd, const double a[3], const dou return 0; } - /* Exclude Sun and Geocentric Moon (GM) from error estimates. */ - /* Can also use Earth position rather than EMB position. */ - if (body < BODY_MERCURY || body > BODY_PLUTO) + if ((body != BODY_SSB) && (body < BODY_MERCURY || body > BODY_PLUTO)) { fprintf(stderr, "PositionArcminError(FATAL): Invalid body = %d\n", body); return 1; @@ -1520,6 +1528,7 @@ static vsop_body_t LookupBody(const char *name) if (!strcmp(name, "Sun")) return BODY_SUN; if (!strcmp(name, "EMB")) return BODY_EMB; if (!strcmp(name, "GM")) return BODY_GM; + if (!strcmp(name, "SSB")) return BODY_SSB; return BODY_INVALID; } @@ -1564,9 +1573,11 @@ static int PrintErrorStats(vsop_body_t body, const char *tag, const error_stat_t case VSOP_SATURN: name = "Saturn"; break; case VSOP_URANUS: name = "Uranus"; break; case VSOP_NEPTUNE: name = "Neptune"; break; - case 8: name = "Pluto"; break; - case 9: name = "Moon"; break; + case VSOP_PLUTO: name = "Pluto"; break; + case VSOP_GM: name = "GM"; break; case VSOP_SUN: name = "Sun"; break; + case VSOP_MOON: name = "Moon"; break; + case VSOP_SSB: name = "SSB"; break; default: name = ""; break; } printf("STATS(%-8s %s): count= %6d max= %10.8lf rms= %10.8lf\n", name, tag, stats->count, stats->max, sqrt(stats->sum / stats->count)); diff --git a/generate/mag_test.js b/generate/mag_test.js index da337435..77bee39b 100644 --- a/generate/mag_test.js +++ b/generate/mag_test.js @@ -140,12 +140,11 @@ function TestMaxMag(filename, body) { function Test() { let all_passed = true; - for (let body of Astronomy.Bodies) { - if (body !== 'Earth' && body !== 'Saturn') { - const data = LoadMagnitudeData(`magnitude/${body}.txt`); - if (!CheckMagnitudeData(body, data)) - all_passed = false; - } + const bodies = ['Sun', 'Moon', 'Mercury', 'Venus', 'Mars', 'Jupiter', 'Uranus', 'Neptune', 'Pluto']; + for (let body of bodies) { + const data = LoadMagnitudeData(`magnitude/${body}.txt`); + if (!CheckMagnitudeData(body, data)) + all_passed = false; } if (!TestSaturn()) @@ -159,4 +158,4 @@ function Test() { Test(); console.log('mag_test: success'); -process.exit(0); \ No newline at end of file +process.exit(0); diff --git a/generate/novas_body.h b/generate/novas_body.h index 058edbaa..571b11d4 100644 --- a/generate/novas_body.h +++ b/generate/novas_body.h @@ -48,5 +48,6 @@ /* The following are extensions supported by my own code... */ #define BODY_EARTH 11 #define BODY_MOON 12 +#define BODY_SSB 13 /* Solar System Barycenter */ #endif /* __DDC_NOVAS_H */ diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 03137286..de524070 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -62,6 +62,12 @@ static const double REFRACTION_NEAR_HORIZON = 34.0 / 60.0; /* degrees of refra static const double SUN_RADIUS_AU = 4.6505e-3; static const double MOON_RADIUS_AU = 1.15717e-5; static const double ASEC180 = 180.0 * 60.0 * 60.0; /* arcseconds per 180 degrees (or pi radians) */ +static const double EARTH_MOON_MASS_RATIO = 81.30056; +static const double SUN_MASS = 333054.25318; /* Sun's mass relative to Earth. */ +static const double JUPITER_MASS = 317.84997; /* Jupiter's mass relative to Earth. */ +static const double SATURN_MASS = 95.16745; /* Saturn's mass relative to Earth. */ +static const double URANUS_MASS = 14.53617; /* Uranus's mass relative to Earth. */ +static const double NEPTUNE_MASS = 17.14886; /* Neptune's mass relative to Earth. */ /** @cond DOXYGEN_SKIP */ #define ARRAYSIZE(x) (sizeof(x) / sizeof(x[0])) @@ -133,13 +139,15 @@ const char *Astronomy_BodyName(astro_body_t body) case BODY_PLUTO: return "Pluto"; case BODY_SUN: return "Sun"; case BODY_MOON: return "Moon"; + case BODY_EMB: return "EMB"; + case BODY_SSB: return "SSB"; default: return ""; } } /** * @brief Returns the #astro_body_t value corresponding to the given English name. - * @param name One of the following strings: Sun, Moon, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto. + * @param name One of the following strings: Sun, Moon, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto, EMB, SSB. * @return If `name` is one of the strings (case-sensitive) listed above, the returned value is the corresponding #astro_body_t value, otherwise it is `BODY_INVALID`. */ astro_body_t Astronomy_BodyCode(const char *name) @@ -157,6 +165,8 @@ astro_body_t Astronomy_BodyCode(const char *name) if (!strcmp(name, "Pluto")) return BODY_PLUTO; if (!strcmp(name, "Sun")) return BODY_SUN; if (!strcmp(name, "Moon")) return BODY_MOON; + if (!strcmp(name, "EMB")) return BODY_EMB; + if (!strcmp(name, "SSB")) return BODY_SSB; } return BODY_INVALID; } @@ -1524,6 +1534,34 @@ static astro_vector_t CalcChebyshev(const astro_cheb_record_t model[], int nrecs /*------------------ end of generated code ------------------*/ +static void AdjustBarycenter(astro_vector_t *ssb, astro_time_t time, astro_body_t body, double pmass) +{ + astro_vector_t planet; + double shift; + + shift = pmass / (pmass + SUN_MASS); + planet = CalcVsop(&vsop[body], time); + ssb->x += shift * planet.x; + ssb->y += shift * planet.y; + ssb->z += shift * planet.z; +} + +static astro_vector_t CalcSolarSystemBarycenter(astro_time_t time) +{ + astro_vector_t ssb; + + ssb.status = ASTRO_SUCCESS; + ssb.t = time; + ssb.x = ssb.y = ssb.z = 0.0; + + AdjustBarycenter(&ssb, time, BODY_JUPITER, JUPITER_MASS); + AdjustBarycenter(&ssb, time, BODY_SATURN, SATURN_MASS); + AdjustBarycenter(&ssb, time, BODY_URANUS, URANUS_MASS); + AdjustBarycenter(&ssb, time, BODY_NEPTUNE, NEPTUNE_MASS); + + return ssb; +} + /** * @brief Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. * @@ -1540,7 +1578,9 @@ static astro_vector_t CalcChebyshev(const astro_cheb_record_t model[], int nrecs * the `status` field inside the returned #astro_vector_t for `ASTRO_SUCCESS` (success) * or any other value (failure) before trusting the resulting vector. * - * @param body A body for which to calculate a heliocentric position: the Sun, Moon, or any of the planets. + * @param body + * A body for which to calculate a heliocentric position: the Sun, Moon, any of the planets, + * the Solar System Barycenter (SSB), or the Earth Moon Barycenter (EMB). * @param time The date and time for which to calculate the position. * @return A heliocentric position vector of the center of the given body. */ @@ -1579,6 +1619,17 @@ astro_vector_t Astronomy_HelioVector(astro_body_t body, astro_time_t time) vector.z += earth.z; return vector; + case BODY_EMB: + vector = Astronomy_GeoMoon(time); + earth = CalcEarth(time); + vector.x = earth.x + (vector.x / (1.0 + EARTH_MOON_MASS_RATIO)); + vector.y = earth.y + (vector.y / (1.0 + EARTH_MOON_MASS_RATIO)); + vector.z = earth.z + (vector.z / (1.0 + EARTH_MOON_MASS_RATIO)); + return vector; + + case BODY_SSB: + return CalcSolarSystemBarycenter(time); + default: return VecError(ASTRO_INVALID_BODY, time); } diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index 10eb25ec..0052fa78 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -122,6 +122,16 @@ namespace CosineKitty /// The Earth's natural satellite, the Moon. /// Moon, + + /// + /// The Earth/Moon Barycenter. + /// + EMB, + + /// + /// The Solar System Barycenter. + /// + SSB, } /// @@ -1285,6 +1295,12 @@ $ASTRO_ADDSOL() internal const double MOON_RADIUS_AU = 1.15717e-5; private const double ASEC180 = 180.0 * 60.0 * 60.0; /* arcseconds per 180 degrees (or pi radians) */ private const double AU_PER_PARSEC = (ASEC180 / Math.PI); /* exact definition of how many AU = one parsec */ + private const double EARTH_MOON_MASS_RATIO = 81.30056; + private const double SUN_MASS = 333054.25318; /* Sun's mass relative to Earth. */ + private const double JUPITER_MASS = 317.84997; /* Jupiter's mass relative to Earth. */ + private const double SATURN_MASS = 95.16745; /* Saturn's mass relative to Earth. */ + private const double URANUS_MASS = 14.53617; /* Uranus's mass relative to Earth. */ + private const double NEPTUNE_MASS = 17.14886; /* Neptune's mass relative to Earth. */ internal static double LongitudeOffset(double diff) { @@ -1971,6 +1987,32 @@ $ASTRO_CSHARP_CHEBYSHEV(8); return new AstroVector(mpos2.x, mpos2.y, mpos2.z, time); } + private static AstroVector BarycenterContrib(AstroTime time, Body body, double pmass) + { + double shift = pmass / (pmass + SUN_MASS); + AstroVector p = CalcVsop(vsop[(int)body], time); + return new AstroVector( + shift * p.x, + shift * p.y, + shift * p.z, + time + ); + } + + private static AstroVector CalcSolarSystemBarycenter(AstroTime time) + { + AstroVector j = BarycenterContrib(time, Body.Jupiter, JUPITER_MASS); + AstroVector s = BarycenterContrib(time, Body.Saturn, SATURN_MASS); + AstroVector u = BarycenterContrib(time, Body.Uranus, URANUS_MASS); + AstroVector n = BarycenterContrib(time, Body.Neptune, NEPTUNE_MASS); + return new AstroVector( + j.x + s.x + u.x + n.x, + j.y + s.y + u.y + n.y, + j.z + s.z + u.z + n.z, + time + ); + } + /// /// Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. /// @@ -1986,11 +2028,13 @@ $ASTRO_CSHARP_CHEBYSHEV(8); /// If given an invalid value for `body`, or the body is `Body.Pluto` and the `time` is outside /// the year range 1700..2200, this function will throw an `ArgumentException`. /// - /// A body for which to calculate a heliocentric position: the Sun, Moon, or any of the planets. + /// A body for which to calculate a heliocentric position: the Sun, Moon, EMB, SSB, or any of the planets. /// The date and time for which to calculate the position. /// A heliocentric position vector of the center of the given body. public static AstroVector HelioVector(Body body, AstroTime time) { + AstroVector earth, geomoon; + switch (body) { case Body.Sun: @@ -2009,6 +2053,30 @@ $ASTRO_CSHARP_CHEBYSHEV(8); case Body.Pluto: return CalcChebyshev(cheb_8, time); + case Body.Moon: + geomoon = GeoMoon(time); + earth = CalcEarth(time); + return new AstroVector( + earth.x + geomoon.x, + earth.y + geomoon.y, + earth.z + geomoon.z, + time + ); + + case Body.EMB: + geomoon = GeoMoon(time); + earth = CalcEarth(time); + double denom = 1.0 + EARTH_MOON_MASS_RATIO; + return new AstroVector( + earth.x + (geomoon.x / denom), + earth.y + (geomoon.y / denom), + earth.z + (geomoon.z / denom), + time + ); + + case Body.SSB: + return CalcSolarSystemBarycenter(time); + default: throw new ArgumentException(string.Format("Invalid body: {0}", body)); } diff --git a/generate/template/astronomy.js b/generate/template/astronomy.js index b7bf1bbc..e56d6c36 100644 --- a/generate/template/astronomy.js +++ b/generate/template/astronomy.js @@ -66,6 +66,12 @@ const SOLAR_DAYS_PER_SIDEREAL_DAY = 0.9972695717592592; const SUN_RADIUS_AU = 4.6505e-3; const MOON_RADIUS_AU = 1.15717e-5; const REFRACTION_NEAR_HORIZON = 34 / 60; // degrees of refractive "lift" seen for objects near horizon +const EARTH_MOON_MASS_RATIO = 81.30056; +const SUN_MASS = 333054.25318; /* Sun's mass relative to Earth. */ +const JUPITER_MASS = 317.84997; /* Jupiter's mass relative to Earth. */ +const SATURN_MASS = 95.16745; /* Saturn's mass relative to Earth. */ +const URANUS_MASS = 14.53617; /* Uranus's mass relative to Earth. */ +const NEPTUNE_MASS = 17.14886; /* Neptune's mass relative to Earth. */ let ob2000; // lazy-evaluated mean obliquity of the ecliptic at J2000, in radians let cos_ob2000; let sin_ob2000; @@ -214,7 +220,9 @@ Astronomy.Bodies = [ 'Saturn', 'Uranus', 'Neptune', - 'Pluto' + 'Pluto', + 'SSB', // Solar System Barycenter + 'EMB' // Earth/Moon Barycenter ]; const Planet = { @@ -1539,6 +1547,23 @@ function CalcChebyshev(model, time) { throw `Cannot extrapolate Chebyshev model for given Terrestrial Time: ${time.tt}`; } +function AdjustBarycenter(ssb, time, body, pmass) { + const shift = pmass / (pmass + SUN_MASS); + const planet = CalcVsop(vsop[body], time); + ssb.x += shift * planet.x; + ssb.y += shift * planet.y; + ssb.z += shift * planet.z; +} + +function CalcSolarSystemBarycenter(time) { + const ssb = new Vector(0.0, 0.0, 0.0, time); + AdjustBarycenter(ssb, time, 'Jupiter', JUPITER_MASS); + AdjustBarycenter(ssb, time, 'Saturn', SATURN_MASS); + AdjustBarycenter(ssb, time, 'Uranus', URANUS_MASS); + AdjustBarycenter(ssb, time, 'Neptune', NEPTUNE_MASS); + return ssb; +} + /** * Calculates heliocentric (i.e., with respect to the center of the Sun) * Cartesian coordinates in the J2000 equatorial system of a celestial @@ -1548,7 +1573,8 @@ function CalcChebyshev(model, time) { * One of the strings * `"Sun"`, `"Moon"`, `"Mercury"`, `"Venus"`, * `"Earth"`, `"Mars"`, `"Jupiter"`, `"Saturn"`, - * `"Uranus"`, `"Neptune"`, or `"Pluto"`. + * `"Uranus"`, `"Neptune"`, `"Pluto"`, + * `"SSB"`, or `"EMB"`. * * @param {(Date | number | Astronomy.AstroTime)} date * The date and time for which the body's position is to be calculated. @@ -1571,6 +1597,15 @@ Astronomy.HelioVector = function(body, date) { var m = Astronomy.GeoMoon(time); return new Vector(e.x+m.x, e.y+m.y, e.z+m.z, time); } + if (body === 'EMB') { + const e = CalcVsop(vsop.Earth, time); + const m = Astronomy.GeoMoon(time); + const denom = 1.0 + EARTH_MOON_MASS_RATIO; + return new Vector(e.x+(m.x/denom), e.y+(m.y/denom), e.z+(m.z/denom), time); + } + if (body === 'SSB') { + return CalcSolarSystemBarycenter(time); + } throw `Astronomy.HelioVector: Unknown body "${body}"`; }; diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 42184cd1..f98f7967 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -59,6 +59,13 @@ _SUN_RADIUS_AU = 4.6505e-3 _MOON_RADIUS_AU = 1.15717e-5 _ASEC180 = 180.0 * 60.0 * 60.0 _AU_PER_PARSEC = _ASEC180 / math.pi +_EARTH_MOON_MASS_RATIO = 81.30056 +_SUN_MASS = 333054.25318 # Sun's mass relative to Earth. +_JUPITER_MASS = 317.84997 # Jupiter's mass relative to Earth. +_SATURN_MASS = 95.16745 # Saturn's mass relative to Earth. +_URANUS_MASS = 14.53617 # Uranus's mass relative to Earth. +_NEPTUNE_MASS = 17.14886 # Neptune's mass relative to Earth. + def _LongitudeOffset(diff): offset = diff @@ -121,6 +128,8 @@ class Body(enum.IntEnum): Pluto: The planet Pluto. Sun: The Sun. Moon: The Earth's moon. + EMB: The Earth/Moon Barycenter. + SSB: The Solar System Barycenter. """ Invalid = -1 Mercury = 0 @@ -134,6 +143,8 @@ class Body(enum.IntEnum): Pluto = 8 Sun = 9 Moon = 10 + EMB = 11 + SSB = 12 def BodyCode(name): """Finds the Body enumeration value, given the name of a body. @@ -1226,6 +1237,21 @@ def Search(func, context, t1, t2, dt_tolerance_seconds): # END Search #---------------------------------------------------------------------------- +def _AdjustBarycenter(ssb, time, body, pmass): + shift = pmass / (pmass + _SUN_MASS) + planet = _CalcVsop(_vsop[body], time) + ssb.x += shift * planet.x + ssb.y += shift * planet.y + ssb.z += shift * planet.z + +def _CalcSolarSystemBarycenter(time): + ssb = Vector(0.0, 0.0, 0.0, time) + _AdjustBarycenter(ssb, time, Body.Jupiter, _JUPITER_MASS) + _AdjustBarycenter(ssb, time, Body.Saturn, _SATURN_MASS) + _AdjustBarycenter(ssb, time, Body.Uranus, _URANUS_MASS) + _AdjustBarycenter(ssb, time, Body.Neptune, _NEPTUNE_MASS) + return ssb + def HelioVector(body, time): """Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. @@ -1244,7 +1270,7 @@ def HelioVector(body, time): ---------- body : Body The celestial body whose heliocentric position is to be calculated: - The Sun, Moon, or any of the planets. + The Sun, Moon, EMB, SSB, or any of the planets. time : Time The time at which to calculate the heliocentric position. @@ -1268,6 +1294,15 @@ def HelioVector(body, time): m = GeoMoon(time) return Vector(e.x+m.x, e.y+m.y, e.z+m.z, time) + if body == Body.EMB: + e = _CalcEarth(time) + m = GeoMoon(time) + d = 1.0 + _EARTH_MOON_MASS_RATIO + return Vector(e.x+(m.x/d), e.y+(m.y/d), e.z+(m.z/d), time) + + if body == Body.SSB: + return _CalcSolarSystemBarycenter(time) + raise InvalidBodyError() diff --git a/generate/test.py b/generate/test.py index a1c0ea0f..859f0897 100755 --- a/generate/test.py +++ b/generate/test.py @@ -111,7 +111,8 @@ def Test_AstroCheck(printflag): bodylist = [ astronomy.Body.Sun, astronomy.Body.Moon, astronomy.Body.Mercury, astronomy.Body.Venus, astronomy.Body.Earth, astronomy.Body.Mars, astronomy.Body.Jupiter, astronomy.Body.Saturn, - astronomy.Body.Uranus, astronomy.Body.Neptune, astronomy.Body.Pluto + astronomy.Body.Uranus, astronomy.Body.Neptune, astronomy.Body.Pluto, + astronomy.Body.SSB, astronomy.Body.EMB ] while time.tt < stop.tt: @@ -121,7 +122,7 @@ def Test_AstroCheck(printflag): pos = astronomy.HelioVector(body, time) if printflag: print('v {} {:0.16f} {:0.16f} {:0.16f} {:0.16f}'.format(name, pos.t.tt, pos.x, pos.y, pos.z)) - if body != astronomy.Body.Earth: + if body != astronomy.Body.Earth and body != astronomy.Body.EMB and body != astronomy.Body.SSB: j2000 = astronomy.Equator(body, time, observer, False, False) ofdate = astronomy.Equator(body, time, observer, True, True) hor = astronomy.Horizon(time, observer, ofdate.ra, ofdate.dec, astronomy.Refraction.Airless) diff --git a/generate/vsop/vsop.h b/generate/vsop/vsop.h index 6d9af10d..a2c9c575 100644 --- a/generate/vsop/vsop.h +++ b/generate/vsop/vsop.h @@ -51,6 +51,7 @@ vsop_formula_t; typedef enum /* values created for compatibility with NOVAS; these are *NOT* the body codes used by VSOP! */ { VSOP_INVALID_BODY = -1, + VSOP_MERCURY = 0, VSOP_VENUS = 1, VSOP_EMB = 2, /* Earth/Moon barycenter, to match NOVAS */ @@ -59,13 +60,21 @@ typedef enum /* values created for compatibility with NOVAS; these are *NOT* VSOP_SATURN = 5, VSOP_URANUS = 6, VSOP_NEPTUNE = 7, - /* NOTE: VSOP does not provide Pluto or Moon. */ + + /* Remaining values are not supported by VSOP. */ + /* Placeholders for NOVAS. */ + VSOP_PLUTO = 8, + VSOP_GM = 9, VSOP_SUN = 10, - VSOP_EARTH = 11, /* weird value so as not to conflict with NOVAS body values */ + + /* Extensions beyond what NOVAS supports. */ + VSOP_EARTH = 11, + VSOP_MOON = 12, + VSOP_SSB = 13 } vsop_body_t; -#define VSOP_BODY_LIMIT 12 +#define VSOP_BODY_LIMIT 14 typedef struct { diff --git a/source/c/README.md b/source/c/README.md index ce9e063f..fcfa11fc 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -224,7 +224,7 @@ This function calculates the angular separation between the given body and the S | Type | Parameter | Description | | --- | --- | --- | -| `const char *` | `name` | One of the following strings: Sun, Moon, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto. | +| `const char *` | `name` | One of the following strings: Sun, Moon, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto, EMB, SSB. | @@ -526,7 +526,7 @@ If given an invalid value for `body`, or the body is `BODY_PLUTO` and the `time` | Type | Parameter | Description | | --- | --- | --- | -| [`astro_body_t`](#astro_body_t) | `body` | A body for which to calculate a heliocentric position: the Sun, Moon, or any of the planets. | +| [`astro_body_t`](#astro_body_t) | `body` | A body for which to calculate a heliocentric position: the Sun, Moon, any of the planets, the Solar System Barycenter (SSB), or the Earth Moon Barycenter (EMB). | | [`astro_time_t`](#astro_time_t) | `time` | The date and time for which to calculate the position. | @@ -1861,6 +1861,8 @@ Aberration correction is useful to improve accuracy of coordinates of apparent l | `BODY_PLUTO` | Pluto | | `BODY_SUN` | Sun | | `BODY_MOON` | Moon | +| `BODY_EMB` | Earth/Moon Barycenter | +| `BODY_SSB` | Solar System Barycenter | diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 7fa52146..66c8b3cc 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -62,6 +62,12 @@ static const double REFRACTION_NEAR_HORIZON = 34.0 / 60.0; /* degrees of refra static const double SUN_RADIUS_AU = 4.6505e-3; static const double MOON_RADIUS_AU = 1.15717e-5; static const double ASEC180 = 180.0 * 60.0 * 60.0; /* arcseconds per 180 degrees (or pi radians) */ +static const double EARTH_MOON_MASS_RATIO = 81.30056; +static const double SUN_MASS = 333054.25318; /* Sun's mass relative to Earth. */ +static const double JUPITER_MASS = 317.84997; /* Jupiter's mass relative to Earth. */ +static const double SATURN_MASS = 95.16745; /* Saturn's mass relative to Earth. */ +static const double URANUS_MASS = 14.53617; /* Uranus's mass relative to Earth. */ +static const double NEPTUNE_MASS = 17.14886; /* Neptune's mass relative to Earth. */ /** @cond DOXYGEN_SKIP */ #define ARRAYSIZE(x) (sizeof(x) / sizeof(x[0])) @@ -133,13 +139,15 @@ const char *Astronomy_BodyName(astro_body_t body) case BODY_PLUTO: return "Pluto"; case BODY_SUN: return "Sun"; case BODY_MOON: return "Moon"; + case BODY_EMB: return "EMB"; + case BODY_SSB: return "SSB"; default: return ""; } } /** * @brief Returns the #astro_body_t value corresponding to the given English name. - * @param name One of the following strings: Sun, Moon, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto. + * @param name One of the following strings: Sun, Moon, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto, EMB, SSB. * @return If `name` is one of the strings (case-sensitive) listed above, the returned value is the corresponding #astro_body_t value, otherwise it is `BODY_INVALID`. */ astro_body_t Astronomy_BodyCode(const char *name) @@ -157,6 +165,8 @@ astro_body_t Astronomy_BodyCode(const char *name) if (!strcmp(name, "Pluto")) return BODY_PLUTO; if (!strcmp(name, "Sun")) return BODY_SUN; if (!strcmp(name, "Moon")) return BODY_MOON; + if (!strcmp(name, "EMB")) return BODY_EMB; + if (!strcmp(name, "SSB")) return BODY_SSB; } return BODY_INVALID; } @@ -2815,6 +2825,34 @@ static astro_vector_t CalcChebyshev(const astro_cheb_record_t model[], int nrecs /*------------------ end of generated code ------------------*/ +static void AdjustBarycenter(astro_vector_t *ssb, astro_time_t time, astro_body_t body, double pmass) +{ + astro_vector_t planet; + double shift; + + shift = pmass / (pmass + SUN_MASS); + planet = CalcVsop(&vsop[body], time); + ssb->x += shift * planet.x; + ssb->y += shift * planet.y; + ssb->z += shift * planet.z; +} + +static astro_vector_t CalcSolarSystemBarycenter(astro_time_t time) +{ + astro_vector_t ssb; + + ssb.status = ASTRO_SUCCESS; + ssb.t = time; + ssb.x = ssb.y = ssb.z = 0.0; + + AdjustBarycenter(&ssb, time, BODY_JUPITER, JUPITER_MASS); + AdjustBarycenter(&ssb, time, BODY_SATURN, SATURN_MASS); + AdjustBarycenter(&ssb, time, BODY_URANUS, URANUS_MASS); + AdjustBarycenter(&ssb, time, BODY_NEPTUNE, NEPTUNE_MASS); + + return ssb; +} + /** * @brief Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. * @@ -2831,7 +2869,9 @@ static astro_vector_t CalcChebyshev(const astro_cheb_record_t model[], int nrecs * the `status` field inside the returned #astro_vector_t for `ASTRO_SUCCESS` (success) * or any other value (failure) before trusting the resulting vector. * - * @param body A body for which to calculate a heliocentric position: the Sun, Moon, or any of the planets. + * @param body + * A body for which to calculate a heliocentric position: the Sun, Moon, any of the planets, + * the Solar System Barycenter (SSB), or the Earth Moon Barycenter (EMB). * @param time The date and time for which to calculate the position. * @return A heliocentric position vector of the center of the given body. */ @@ -2870,6 +2910,17 @@ astro_vector_t Astronomy_HelioVector(astro_body_t body, astro_time_t time) vector.z += earth.z; return vector; + case BODY_EMB: + vector = Astronomy_GeoMoon(time); + earth = CalcEarth(time); + vector.x = earth.x + (vector.x / (1.0 + EARTH_MOON_MASS_RATIO)); + vector.y = earth.y + (vector.y / (1.0 + EARTH_MOON_MASS_RATIO)); + vector.z = earth.z + (vector.z / (1.0 + EARTH_MOON_MASS_RATIO)); + return vector; + + case BODY_SSB: + return CalcSolarSystemBarycenter(time); + default: return VecError(ASTRO_INVALID_BODY, time); } diff --git a/source/c/astronomy.h b/source/c/astronomy.h index 703be79c..f49fbf79 100644 --- a/source/c/astronomy.h +++ b/source/c/astronomy.h @@ -209,12 +209,14 @@ typedef enum BODY_NEPTUNE, /**< Neptune */ BODY_PLUTO, /**< Pluto */ BODY_SUN, /**< Sun */ - BODY_MOON /**< Moon */ + BODY_MOON, /**< Moon */ + BODY_EMB, /**< Earth/Moon Barycenter */ + BODY_SSB /**< Solar System Barycenter */ } astro_body_t; #define MIN_BODY BODY_MERCURY /**< Minimum valid astro_body_t value; useful for iteration. */ -#define MAX_BODY BODY_MOON /**< Maximum valid astro_body_t value; useful for iteration. */ +#define MAX_BODY BODY_SSB /**< Maximum valid astro_body_t value; useful for iteration. */ #define MIN_YEAR 1700 /**< Minimum year value supported by Astronomy Engine. */ #define MAX_YEAR 2200 /**< Maximum year value supported by Astronomy Engine. */ diff --git a/source/csharp/README.md b/source/csharp/README.md index 3e8b82a9..8567da8c 100644 --- a/source/csharp/README.md +++ b/source/csharp/README.md @@ -327,7 +327,7 @@ the year range 1700..2200, this function will throw an `ArgumentException`. | Type | Parameter | Description | | --- | --- | --- | -| [`Body`](#Body) | `body` | A body for which to calculate a heliocentric position: the Sun, Moon, or any of the planets. | +| [`Body`](#Body) | `body` | A body for which to calculate a heliocentric position: the Sun, Moon, EMB, SSB, or any of the planets. | | [`AstroTime`](#AstroTime) | `time` | The date and time for which to calculate the position. | **Returns:** A heliocentric position vector of the center of the given body. @@ -1391,6 +1391,8 @@ according to the historical and predictive Delta-T model provided by the | `Pluto` | The planet Pluto. | | `Sun` | The Sun. | | `Moon` | The Earth's natural satellite, the Moon. | +| `EMB` | The Earth/Moon Barycenter. | +| `SSB` | The Solar System Barycenter. | --- diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index ff21025d..f29719dd 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -122,6 +122,16 @@ namespace CosineKitty /// The Earth's natural satellite, the Moon. /// Moon, + + /// + /// The Earth/Moon Barycenter. + /// + EMB, + + /// + /// The Solar System Barycenter. + /// + SSB, } /// @@ -1390,6 +1400,12 @@ namespace CosineKitty internal const double MOON_RADIUS_AU = 1.15717e-5; private const double ASEC180 = 180.0 * 60.0 * 60.0; /* arcseconds per 180 degrees (or pi radians) */ private const double AU_PER_PARSEC = (ASEC180 / Math.PI); /* exact definition of how many AU = one parsec */ + private const double EARTH_MOON_MASS_RATIO = 81.30056; + private const double SUN_MASS = 333054.25318; /* Sun's mass relative to Earth. */ + private const double JUPITER_MASS = 317.84997; /* Jupiter's mass relative to Earth. */ + private const double SATURN_MASS = 95.16745; /* Saturn's mass relative to Earth. */ + private const double URANUS_MASS = 14.53617; /* Uranus's mass relative to Earth. */ + private const double NEPTUNE_MASS = 17.14886; /* Neptune's mass relative to Earth. */ internal static double LongitudeOffset(double diff) { @@ -3267,6 +3283,32 @@ namespace CosineKitty return new AstroVector(mpos2.x, mpos2.y, mpos2.z, time); } + private static AstroVector BarycenterContrib(AstroTime time, Body body, double pmass) + { + double shift = pmass / (pmass + SUN_MASS); + AstroVector p = CalcVsop(vsop[(int)body], time); + return new AstroVector( + shift * p.x, + shift * p.y, + shift * p.z, + time + ); + } + + private static AstroVector CalcSolarSystemBarycenter(AstroTime time) + { + AstroVector j = BarycenterContrib(time, Body.Jupiter, JUPITER_MASS); + AstroVector s = BarycenterContrib(time, Body.Saturn, SATURN_MASS); + AstroVector u = BarycenterContrib(time, Body.Uranus, URANUS_MASS); + AstroVector n = BarycenterContrib(time, Body.Neptune, NEPTUNE_MASS); + return new AstroVector( + j.x + s.x + u.x + n.x, + j.y + s.y + u.y + n.y, + j.z + s.z + u.z + n.z, + time + ); + } + /// /// Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. /// @@ -3282,11 +3324,13 @@ namespace CosineKitty /// If given an invalid value for `body`, or the body is `Body.Pluto` and the `time` is outside /// the year range 1700..2200, this function will throw an `ArgumentException`. /// - /// A body for which to calculate a heliocentric position: the Sun, Moon, or any of the planets. + /// A body for which to calculate a heliocentric position: the Sun, Moon, EMB, SSB, or any of the planets. /// The date and time for which to calculate the position. /// A heliocentric position vector of the center of the given body. public static AstroVector HelioVector(Body body, AstroTime time) { + AstroVector earth, geomoon; + switch (body) { case Body.Sun: @@ -3305,6 +3349,30 @@ namespace CosineKitty case Body.Pluto: return CalcChebyshev(cheb_8, time); + case Body.Moon: + geomoon = GeoMoon(time); + earth = CalcEarth(time); + return new AstroVector( + earth.x + geomoon.x, + earth.y + geomoon.y, + earth.z + geomoon.z, + time + ); + + case Body.EMB: + geomoon = GeoMoon(time); + earth = CalcEarth(time); + double denom = 1.0 + EARTH_MOON_MASS_RATIO; + return new AstroVector( + earth.x + (geomoon.x / denom), + earth.y + (geomoon.y / denom), + earth.z + (geomoon.z / denom), + time + ); + + case Body.SSB: + return CalcSolarSystemBarycenter(time); + default: throw new ArgumentException(string.Format("Invalid body: {0}", body)); } diff --git a/source/js/README.md b/source/js/README.md index eea5b2ed..145e327a 100644 --- a/source/js/README.md +++ b/source/js/README.md @@ -708,7 +708,7 @@ body at a specified time. The position is not corrected for light travel time or | Param | Type | Description | | --- | --- | --- | -| body | string | One of the strings `"Sun"`, `"Moon"`, `"Mercury"`, `"Venus"`, `"Earth"`, `"Mars"`, `"Jupiter"`, `"Saturn"`, `"Uranus"`, `"Neptune"`, or `"Pluto"`. | +| body | string | One of the strings `"Sun"`, `"Moon"`, `"Mercury"`, `"Venus"`, `"Earth"`, `"Mars"`, `"Jupiter"`, `"Saturn"`, `"Uranus"`, `"Neptune"`, `"Pluto"`, `"SSB"`, or `"EMB"`. | | date | Date \| number \| [AstroTime](#Astronomy.AstroTime) | The date and time for which the body's position is to be calculated. | diff --git a/source/js/astronomy.js b/source/js/astronomy.js index a1bf9e15..e23d3634 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -66,6 +66,12 @@ const SOLAR_DAYS_PER_SIDEREAL_DAY = 0.9972695717592592; const SUN_RADIUS_AU = 4.6505e-3; const MOON_RADIUS_AU = 1.15717e-5; const REFRACTION_NEAR_HORIZON = 34 / 60; // degrees of refractive "lift" seen for objects near horizon +const EARTH_MOON_MASS_RATIO = 81.30056; +const SUN_MASS = 333054.25318; /* Sun's mass relative to Earth. */ +const JUPITER_MASS = 317.84997; /* Jupiter's mass relative to Earth. */ +const SATURN_MASS = 95.16745; /* Saturn's mass relative to Earth. */ +const URANUS_MASS = 14.53617; /* Uranus's mass relative to Earth. */ +const NEPTUNE_MASS = 17.14886; /* Neptune's mass relative to Earth. */ let ob2000; // lazy-evaluated mean obliquity of the ecliptic at J2000, in radians let cos_ob2000; let sin_ob2000; @@ -214,7 +220,9 @@ Astronomy.Bodies = [ 'Saturn', 'Uranus', 'Neptune', - 'Pluto' + 'Pluto', + 'SSB', // Solar System Barycenter + 'EMB' // Earth/Moon Barycenter ]; const Planet = { @@ -2609,6 +2617,23 @@ function CalcChebyshev(model, time) { throw `Cannot extrapolate Chebyshev model for given Terrestrial Time: ${time.tt}`; } +function AdjustBarycenter(ssb, time, body, pmass) { + const shift = pmass / (pmass + SUN_MASS); + const planet = CalcVsop(vsop[body], time); + ssb.x += shift * planet.x; + ssb.y += shift * planet.y; + ssb.z += shift * planet.z; +} + +function CalcSolarSystemBarycenter(time) { + const ssb = new Vector(0.0, 0.0, 0.0, time); + AdjustBarycenter(ssb, time, 'Jupiter', JUPITER_MASS); + AdjustBarycenter(ssb, time, 'Saturn', SATURN_MASS); + AdjustBarycenter(ssb, time, 'Uranus', URANUS_MASS); + AdjustBarycenter(ssb, time, 'Neptune', NEPTUNE_MASS); + return ssb; +} + /** * Calculates heliocentric (i.e., with respect to the center of the Sun) * Cartesian coordinates in the J2000 equatorial system of a celestial @@ -2618,7 +2643,8 @@ function CalcChebyshev(model, time) { * One of the strings * `"Sun"`, `"Moon"`, `"Mercury"`, `"Venus"`, * `"Earth"`, `"Mars"`, `"Jupiter"`, `"Saturn"`, - * `"Uranus"`, `"Neptune"`, or `"Pluto"`. + * `"Uranus"`, `"Neptune"`, `"Pluto"`, + * `"SSB"`, or `"EMB"`. * * @param {(Date | number | Astronomy.AstroTime)} date * The date and time for which the body's position is to be calculated. @@ -2641,6 +2667,15 @@ Astronomy.HelioVector = function(body, date) { var m = Astronomy.GeoMoon(time); return new Vector(e.x+m.x, e.y+m.y, e.z+m.z, time); } + if (body === 'EMB') { + const e = CalcVsop(vsop.Earth, time); + const m = Astronomy.GeoMoon(time); + const denom = 1.0 + EARTH_MOON_MASS_RATIO; + return new Vector(e.x+(m.x/denom), e.y+(m.y/denom), e.z+(m.z/denom), time); + } + if (body === 'SSB') { + return CalcSolarSystemBarycenter(time); + } throw `Astronomy.HelioVector: Unknown body "${body}"`; }; diff --git a/source/js/astronomy.min.js b/source/js/astronomy.min.js index 5c6b79d9..7eba19d0 100644 --- a/source/js/astronomy.min.js +++ b/source/js/astronomy.min.js @@ -28,35 +28,35 @@ SOFTWARE. */ var $jscomp=$jscomp||{};$jscomp.scope={};$jscomp.arrayIteratorImpl=function(d){var u=0;return function(){return uMath.abs(c))throw"AngleBetween: first vector is too short.";var e=b.x*b.x+b.y*b.y+b.z*b.z;if(1E-8>Math.abs(e))throw"AngleBetween: second vector is too short.";a=(a.x*b.x+a.y*b.y+a.z*b.z)/Math.sqrt(c*e);return-1>=a?180:1<=a?0:57.29577951308232*Math.acos(a)}function w(a){a:{var b=a+51544.5;if(b<=r[0].mjd)b=r[0].dt;else if(b>=r[r.length-1].mjd)b=r[r.length-1].dt;else{for(var c=0,e=r.length-2;c<=e;){var d= -c+e>>1;if(br[d+1].mjd)c=d+1;else{b=r[d].dt+(b-r[d].mjd)/(r[d+1].mjd-r[d].mjd)*(r[d+1].dt-r[d].dt);break a}}throw"Could not find Delta-T value for MJD="+b;}}return a+b/86400}function E(a){a=a.tt/36525;return(((((-4.34E-8*a-5.76E-7)*a+.0020034)*a-1.831E-4)*a-46.836769)*a+84381.406)/3600}function P(a){var b;if(!Q||1E-6=b;++b)0!==a[b]&&h(f.x,f.y,e(N,a[b],b),e(r,a[b],b),function(a,b){return f.x=a,f.y=b});return f}function f(a,b,c,e,d,f,g,h){d=m(d,f,g,h);q+=a*d.y;v+=b*d.y;A+=c*d.x;B+=e*d.x}++J.calcmoon;a=a.tt/36525;var n,p,q,v,N=c(-6,6,1,4),r=c(-6,6,1,4);var t=a*a;var A=v=q=0;var B=3422.7;var w=l(.19833+.05611*a);var z=l(.27869+.04508*a);var x=l(.16827-.36903*a);var E=l(.34734-5.37261* -a),G=l(.10498-5.37899*a),C=l(.42681-.41855*a),M=l(.14943-5.37511*a);var F=.84*w+.31*z+14.27*x+7.26*E+.28*G+.24*C;var K=2.94*w+.31*z+14.27*x+9.34*E+1.12*G+.83*C;var H=-6.4*w-1.89*C;x=.21*w+.31*z+14.27*x-88.7*E-15.3*G+.24*C-1.86*M;z=F-H;w=-3.332E-6*l(.59734-5.37261*a)-5.39E-7*l(.35498-5.37899*a)-6.4E-8*l(.39943-5.37511*a);F=I*u(.60643382+1336.85522467*a-3.13E-6*t)+F/L;K=I*u(.37489701+1325.55240982*a+2.565E-5*t)+K/L;H=I*u(.99312619+99.99735956*a-4.4E-7*t)+H/L;x=I*u(.25909118+1342.2278298*a-8.92E-6*t)+ -x/L;t=I*u(.82736186+1236.85308708*a-3.97E-6*t)+z/L;for(n=1;4>=n;++n){switch(n){case 1:var D=K;var y=4;var O=1.000002208;break;case 2:D=H;y=3;O=.997504612-.002495388*a;break;case 3:D=x;y=4;O=1.000002708+139.978*w;break;case 4:D=t,y=6,O=1}d(0,n,1);d(1,n,Math.cos(D)*O);k(0,n,0);k(1,n,Math.sin(D)*O);for(p=2;p<=y;++p)h(e(N,p-1,n),e(r,p-1,n),e(N,1,n),e(r,1,n),function(a,b){return d(p,n,a),k(p,n,b)});for(p=1;p<=y;++p)d(-p,n,e(N,p,n)),k(-p,n,-e(r,p,n))}f(13.902,14.06,-.001,.2607,0,0,0,4);f(.403,-4.01,.394, +$jscomp.defineProperty=$jscomp.ASSUME_ES5||"function"==typeof Object.defineProperties?Object.defineProperty:function(d,u,B){d!=Array.prototype&&d!=Object.prototype&&(d[u]=B.value)};$jscomp.getGlobal=function(d){return"undefined"!=typeof window&&window===d?d:"undefined"!=typeof global&&null!=global?global:d};$jscomp.global=$jscomp.getGlobal(this); +$jscomp.polyfill=function(d,u,B,w){if(u){B=$jscomp.global;d=d.split(".");for(w=0;wMath.abs(c))throw"AngleBetween: first vector is too short.";var e=b.x*b.x+b.y*b.y+b.z*b.z;if(1E-8>Math.abs(e))throw"AngleBetween: second vector is too short.";a=(a.x*b.x+a.y*b.y+a.z*b.z)/Math.sqrt(c*e);return-1>=a?180:1<=a?0:57.29577951308232*Math.acos(a)}function w(a){a:{var b=a+51544.5;if(b<=r[0].mjd)b=r[0].dt;else if(b>=r[r.length-1].mjd)b=r[r.length-1].dt;else{for(var c=0,e=r.length-2;c<=e;){var d= +c+e>>1;if(br[d+1].mjd)c=d+1;else{b=r[d].dt+(b-r[d].mjd)/(r[d+1].mjd-r[d].mjd)*(r[d+1].dt-r[d].dt);break a}}throw"Could not find Delta-T value for MJD="+b;}}return a+b/86400}function F(a){a=a.tt/36525;return(((((-4.34E-8*a-5.76E-7)*a+.0020034)*a-1.831E-4)*a-46.836769)*a+84381.406)/3600}function P(a){var b;if(!Q||1E-6=b;++b)0!==a[b]&&h(f.x,f.y,e(N,a[b],b),e(t,a[b],b),function(a,b){return f.x=a,f.y=b});return f}function f(a,b,c,e,d,f,g,h){d=m(d,f,g,h);q+=a*d.y;v+=b*d.y;y+=c*d.x;B+=e*d.x}++J.calcmoon;a=a.tt/36525;var n,p,q,v,N=c(-6,6,1,4),t=c(-6,6,1,4);var r=a*a;var y=v=q=0;var B=3422.7;var w=l(.19833+.05611*a);var A=l(.27869+.04508*a);var x=l(.16827-.36903*a);var C=l(.34734-5.37261* +a),F=l(.10498-5.37899*a),D=l(.42681-.41855*a),M=l(.14943-5.37511*a);var G=.84*w+.31*A+14.27*x+7.26*C+.28*F+.24*D;var K=2.94*w+.31*A+14.27*x+9.34*C+1.12*F+.83*D;var H=-6.4*w-1.89*D;x=.21*w+.31*A+14.27*x-88.7*C-15.3*F+.24*D-1.86*M;A=G-H;w=-3.332E-6*l(.59734-5.37261*a)-5.39E-7*l(.35498-5.37899*a)-6.4E-8*l(.39943-5.37511*a);G=I*u(.60643382+1336.85522467*a-3.13E-6*r)+G/L;K=I*u(.37489701+1325.55240982*a+2.565E-5*r)+K/L;H=I*u(.99312619+99.99735956*a-4.4E-7*r)+H/L;x=I*u(.25909118+1342.2278298*a-8.92E-6*r)+ +x/L;r=I*u(.82736186+1236.85308708*a-3.97E-6*r)+A/L;for(n=1;4>=n;++n){switch(n){case 1:var E=K;var z=4;var O=1.000002208;break;case 2:E=H;z=3;O=.997504612-.002495388*a;break;case 3:E=x;z=4;O=1.000002708+139.978*w;break;case 4:E=r,z=6,O=1}d(0,n,1);d(1,n,Math.cos(E)*O);k(0,n,0);k(1,n,Math.sin(E)*O);for(p=2;p<=z;++p)h(e(N,p-1,n),e(t,p-1,n),e(N,1,n),e(t,1,n),function(a,b){return d(p,n,a),k(p,n,b)});for(p=1;p<=z;++p)d(-p,n,e(N,p,n)),k(-p,n,-e(t,p,n))}f(13.902,14.06,-.001,.2607,0,0,0,4);f(.403,-4.01,.394, .0023,0,0,0,3);f(2369.912,2373.36,.601,28.2333,0,0,0,2);f(-125.154,-112.79,-.725,-.9781,0,0,0,1);f(1.979,6.98,-.445,.0433,1,0,0,4);f(191.953,192.72,.029,3.0861,1,0,0,2);f(-8.466,-13.51,.455,-.1093,1,0,0,1);f(22639.5,22609.07,.079,186.5398,1,0,0,0);f(18.609,3.59,-.094,.0118,1,0,0,-1);f(-4586.465,-4578.13,-.077,34.3117,1,0,0,-2);f(3.215,5.44,.192,-.0386,1,0,0,-3);f(-38.428,-38.64,.001,.6008,1,0,0,-4);f(-.393,-1.43,-.092,.0086,1,0,0,-6);f(-.289,-1.59,.123,-.0053,0,1,0,4);f(-24.42,-25.1,.04,-.3,0,1,0, 2);f(18.023,17.93,.007,.1494,0,1,0,1);f(-668.146,-126.98,-1.302,-.3997,0,1,0,0);f(.56,.32,-.001,-.0037,0,1,0,-1);f(-165.145,-165.06,.054,1.9178,0,1,0,-2);f(-1.877,-6.46,-.416,.0339,0,1,0,-4);f(.213,1.02,-.074,.0054,2,0,0,4);f(14.387,14.78,-.017,.2833,2,0,0,2);f(-.586,-1.2,.054,-.01,2,0,0,1);f(769.016,767.96,.107,10.1657,2,0,0,0);f(1.75,2.01,-.018,.0155,2,0,0,-1);f(-211.656,-152.53,5.679,-.3039,2,0,0,-2);f(1.225,.91,-.03,-.0088,2,0,0,-3);f(-30.773,-34.07,-.308,.3722,2,0,0,-4);f(-.57,-1.4,-.074,.0109, 2,0,0,-6);f(-2.921,-11.75,.787,-.0484,1,1,0,2);f(1.267,1.52,-.022,.0164,1,1,0,1);f(-109.673,-115.18,.461,-.949,1,1,0,0);f(-205.962,-182.36,2.056,1.4437,1,1,0,-2);f(.233,.36,.012,-.0025,1,1,0,-3);f(-4.391,-9.66,-.471,.0673,1,1,0,-4);f(.283,1.53,-.111,.006,1,-1,0,4);f(14.577,31.7,-1.54,.2302,1,-1,0,2);f(147.687,138.76,.679,1.1528,1,-1,0,0);f(-1.089,.55,.021,0,1,-1,0,-1);f(28.475,23.59,-.443,-.2257,1,-1,0,-2);f(-.276,-.38,-.006,-.0036,1,-1,0,-3);f(.636,2.27,.146,-.0102,1,-1,0,-4);f(-.189,-1.68,.131, -.0028,0,2,0,2);f(-7.486,-.66,-.037,-.0086,0,2,0,0);f(-8.096,-16.35,-.74,.0918,0,2,0,-2);f(-5.741,-.04,0,-9E-4,0,0,2,2);f(.255,0,0,0,0,0,2,1);f(-411.608,-.2,0,-.0124,0,0,2,0);f(.584,.84,0,.0071,0,0,2,-1);f(-55.173,-52.14,0,-.1052,0,0,2,-2);f(.254,.25,0,-.0017,0,0,2,-3);f(.025,-1.67,0,.0031,0,0,2,-4);f(1.06,2.96,-.166,.0243,3,0,0,2);f(36.124,50.64,-1.3,.6215,3,0,0,0);f(-13.193,-16.4,.258,-.1187,3,0,0,-2);f(-1.187,-.74,.042,.0074,3,0,0,-4);f(-.293,-.31,-.002,.0046,3,0,0,-6);f(-.29,-1.45,.116,-.0051, 2,1,0,2);f(-7.649,-10.56,.259,-.1038,2,1,0,0);f(-8.627,-7.59,.078,-.0192,2,1,0,-2);f(-2.74,-2.54,.022,.0324,2,1,0,-4);f(1.181,3.32,-.212,.0213,2,-1,0,2);f(9.703,11.67,-.151,.1268,2,-1,0,0);f(-.352,-.37,.001,-.0028,2,-1,0,-1);f(-2.494,-1.17,-.003,-.0017,2,-1,0,-2);f(.36,.2,-.012,-.0043,2,-1,0,-4);f(-1.167,-1.25,.008,-.0106,1,2,0,0);f(-7.412,-6.12,.117,.0484,1,2,0,-2);f(-.311,-.65,-.032,.0044,1,2,0,-4);f(.757,1.82,-.105,.0112,1,-2,0,2);f(2.58,2.32,.027,.0196,1,-2,0,0);f(2.533,2.4,-.014,-.0212,1,-2, 0,-2);f(-.344,-.57,-.025,.0036,0,3,0,-2);f(-.992,-.02,0,0,1,0,2,2);f(-45.099,-.02,0,-.001,1,0,2,0);f(-.179,-9.52,0,-.0833,1,0,2,-2);f(-.301,-.33,0,.0014,1,0,2,-4);f(-6.382,-3.37,0,-.0481,1,0,-2,2);f(39.528,85.13,0,-.7136,1,0,-2,0);f(9.366,.71,0,-.0112,1,0,-2,-2);f(.202,.02,0,0,1,0,-2,-4);f(.415,.1,0,.0013,0,1,2,0);f(-2.152,-2.26,0,-.0066,0,1,2,-2);f(-1.44,-1.3,0,.0014,0,1,-2,2);f(.384,-.04,0,0,0,1,-2,-2);f(1.938,3.6,-.145,.0401,4,0,0,0);f(-.952,-1.58,.052,-.013,4,0,0,-2);f(-.551,-.94,.032,-.0097, -3,1,0,0);f(-.482,-.57,.005,-.0045,3,1,0,-2);f(.681,.96,-.026,.0115,3,-1,0,0);f(-.297,-.27,.002,-9E-4,2,2,0,-2);f(.254,.21,-.003,0,2,-2,0,-2);f(-.25,-.22,.004,.0014,1,3,0,-2);f(-3.996,0,0,4E-4,2,0,2,0);f(.557,-.75,0,-.009,2,0,2,-2);f(-.459,-.38,0,-.0053,2,0,-2,2);f(-1.298,.74,0,4E-4,2,0,-2,0);f(.538,1.14,0,-.0141,2,0,-2,-2);f(.263,.02,0,0,1,1,2,0);f(.426,.07,0,-6E-4,1,1,-2,-2);f(-.304,.03,0,3E-4,1,-1,2,0);f(-.372,-.19,0,-.0027,1,-1,-2,2);f(.418,0,0,0,0,0,4,0);f(-.33,-.04,0,0,3,0,2,0);y=-526.069*m(0, -0,1,-2).y;y+=-3.352*m(0,0,1,-4).y;y+=44.297*m(1,0,1,-2).y;y+=-6*m(1,0,1,-4).y;y+=20.599*m(-1,0,1,0).y;y+=-30.598*m(-1,0,1,-2).y;y+=-24.649*m(-2,0,1,0).y;y+=-2*m(-2,0,1,-2).y;y+=-22.571*m(0,1,1,-2).y;y+=10.985*m(0,-1,1,-2).y;q+=.82*l(.7736-62.5512*a)+.31*l(.0466-125.1025*a)+.35*l(.5785-25.1042*a)+.66*l(.4591+1335.8075*a)+.64*l(.313-91.568*a)+1.14*l(.148+1331.2898*a)+.21*l(.5918+1056.5859*a)+.44*l(.5784+1322.8595*a)+.24*l(.2275-5.7374*a)+.28*l(.2965+2.6929*a)+.33*l(.3132+6.3368*a);a=x+v/L;a=(1.000002708+ -139.978*w)*(18518.511+1.189+A)*Math.sin(a)-6.24*Math.sin(3*a)+y;return{geo_eclip_lon:I*u((F+q/L)/I),geo_eclip_lat:Math.PI/648E3*a,distance_au:4.263520978299708E-5*L/(.999953253*B)}}function R(a,b,c){a=X(a,c);return[a.rot[0][0]*b[0]+a.rot[1][0]*b[1]+a.rot[2][0]*b[2],a.rot[0][1]*b[0]+a.rot[1][1]*b[1]+a.rot[2][1]*b[2],a.rot[0][2]*b[0]+a.rot[1][2]*b[1]+a.rot[2][2]*b[2]]}function X(a,b){var c=84381.406;if(0!==a&&0!==b)throw"One of (tt1, tt2) must be 0.";a=(b-a)/36525;0===b&&(a=-a);var e=((((3.337E-7*a- +3,1,0,0);f(-.482,-.57,.005,-.0045,3,1,0,-2);f(.681,.96,-.026,.0115,3,-1,0,0);f(-.297,-.27,.002,-9E-4,2,2,0,-2);f(.254,.21,-.003,0,2,-2,0,-2);f(-.25,-.22,.004,.0014,1,3,0,-2);f(-3.996,0,0,4E-4,2,0,2,0);f(.557,-.75,0,-.009,2,0,2,-2);f(-.459,-.38,0,-.0053,2,0,-2,2);f(-1.298,.74,0,4E-4,2,0,-2,0);f(.538,1.14,0,-.0141,2,0,-2,-2);f(.263,.02,0,0,1,1,2,0);f(.426,.07,0,-6E-4,1,1,-2,-2);f(-.304,.03,0,3E-4,1,-1,2,0);f(-.372,-.19,0,-.0027,1,-1,-2,2);f(.418,0,0,0,0,0,4,0);f(-.33,-.04,0,0,3,0,2,0);z=-526.069*m(0, +0,1,-2).y;z+=-3.352*m(0,0,1,-4).y;z+=44.297*m(1,0,1,-2).y;z+=-6*m(1,0,1,-4).y;z+=20.599*m(-1,0,1,0).y;z+=-30.598*m(-1,0,1,-2).y;z+=-24.649*m(-2,0,1,0).y;z+=-2*m(-2,0,1,-2).y;z+=-22.571*m(0,1,1,-2).y;z+=10.985*m(0,-1,1,-2).y;q+=.82*l(.7736-62.5512*a)+.31*l(.0466-125.1025*a)+.35*l(.5785-25.1042*a)+.66*l(.4591+1335.8075*a)+.64*l(.313-91.568*a)+1.14*l(.148+1331.2898*a)+.21*l(.5918+1056.5859*a)+.44*l(.5784+1322.8595*a)+.24*l(.2275-5.7374*a)+.28*l(.2965+2.6929*a)+.33*l(.3132+6.3368*a);a=x+v/L;a=(1.000002708+ +139.978*w)*(18518.511+1.189+y)*Math.sin(a)-6.24*Math.sin(3*a)+z;return{geo_eclip_lon:I*u((G+q/L)/I),geo_eclip_lat:Math.PI/648E3*a,distance_au:4.263520978299708E-5*L/(.999953253*B)}}function R(a,b,c){a=Y(a,c);return[a.rot[0][0]*b[0]+a.rot[1][0]*b[1]+a.rot[2][0]*b[2],a.rot[0][1]*b[0]+a.rot[1][1]*b[1]+a.rot[2][1]*b[2],a.rot[0][2]*b[0]+a.rot[1][2]*b[1]+a.rot[2][2]*b[2]]}function Y(a,b){var c=84381.406;if(0!==a&&0!==b)throw"One of (tt1, tt2) must be 0.";a=(b-a)/36525;0===b&&(a=-a);var e=((((3.337E-7*a- 4.67E-7)*a-.00772503)*a+.0512623)*a-.025754)*a+c;c*=4.84813681109536E-6;var d=((((-9.51E-8*a+1.32851E-4)*a-.00114045)*a-1.0790069)*a+5038.481507)*a*4.84813681109536E-6;e*=4.84813681109536E-6;var k=((((-5.6E-8*a+1.70663E-4)*a-.00121197)*a-2.3814292)*a+10.556403)*a*4.84813681109536E-6;a=Math.sin(c);c=Math.cos(c);var h=Math.sin(-d);d=Math.cos(-d);var l=Math.sin(-e);var m=Math.cos(-e);var f=Math.sin(k);var n=Math.cos(k);k=n*d-h*f*m;e=n*h*c+f*m*d*c-a*f*l;var p=n*h*a+f*m*d*a+c*f*l;var q=-f*d-h*n*m;var v= --f*h*c+n*m*d*c-a*n*l;f=-f*h*a+n*m*d*a+c*n*l;h*=l;n=-l*d*c-a*m;a=-l*d*a+m*c;return 0===b?new z([[k,e,p],[q,v,f],[h,n,a]]):new z([[k,q,h],[e,v,n],[p,f,a]])}function S(a){var b=a.tt/36525,c=15*P(a).ee;a=(.779057273264+.00273781191135448*a.ut+a.ut%1)%1*360;0>a&&(a+=360);b=((c+.014506+((((-3.68E-8*b-2.9956E-5)*b-4.4E-7)*b+1.3915817)*b+4612.156534)*b)/3600+a)%360/15;0>b&&(b+=24);return b}function Y(a,b,c){a=Z(a,b);return[a.rot[0][0]*c[0]+a.rot[1][0]*c[1]+a.rot[2][0]*c[2],a.rot[0][1]*c[0]+a.rot[1][1]*c[1]+ -a.rot[2][1]*c[2],a.rot[0][2]*c[0]+a.rot[1][2]*c[1]+a.rot[2][2]*c[2]]}function Z(a,b){a=P(a);var c=.017453292519943295*a.mobl,e=.017453292519943295*a.tobl,d=4.84813681109536E-6*a.dpsi;a=Math.cos(c);c=Math.sin(c);var k=Math.cos(e),h=Math.sin(e);e=Math.cos(d);var l=Math.sin(d);d=-l*a;var m=-l*c,f=l*k,n=e*a*k+c*h,p=e*c*k-a*h;l*=h;var q=e*a*h-c*k;a=e*c*h+a*k;return 0===b?new z([[e,f,l],[d,n,q],[m,p,a]]):new z([[e,d,m],[f,n,p],[l,q,a]])}function la(a){if(!(a instanceof Array)||3!==a.length)return!1;for(var b= -0;3>b;++b){if(!(a[b]instanceof Array)||3!==a[b].length)return!1;for(var c=0;3>c;++c)if("number"!==typeof a[b][c])return!1}return!0}function ca(a){var b=a[0]*a[0]+a[1]*a[1],c=Math.sqrt(b+a[2]*a[2]);if(0===b){if(0===a[2])throw"Indeterminate sky coordinates";return 0>a[2]?{ra:0,dec:-90,dist:c}:{ra:0,dec:90,dist:c}}var e=Math.atan2(a[1],a[0])/.2617993877991494;0>e&&(e+=24);return new da(e,Math.atan2(a[2],Math.sqrt(b))/.017453292519943295,c)}function G(a,b){var c=.017453292519943295*a;a=Math.cos(c);c= -Math.sin(c);return[a*b[0]+c*b[1]+0*b[2],-c*b[0]+a*b[1]+0*b[2],0*b[0]+0*b[1]+1*b[2]]}function ea(a,b,c,e,d){var g=b*e+c*d;b=-b*d+c*e;c=Math.sqrt(a*a+g*g);e=0;0e&&(e+=360));return new ma(a,g,b,57.29577951308232*Math.atan2(b,c),e)}function K(a,b){var c=1,e=0;a=$jscomp.makeIterator(a);for(var d=a.next();!d.done;d=a.next()){var k=0;d=$jscomp.makeIterator(d.value);for(var h=d.next();!h.done;h=d.next())h=h.value,k+=h[0]*Math.cos(h[1]+b*h[2]);e+=c*k;c*=b}return e} -function C(a,b){var c=b.tt/365250,e=[];a=$jscomp.makeIterator(a);for(var d=a.next();!d.done;d=a.next())e.push(K(d.value,c));c=e[2]*Math.cos(e[1]);e=[c*Math.cos(e[0]),c*Math.sin(e[0]),e[2]*Math.sin(e[1])];return new t(e[0]+4.4036E-7*e[1]-1.90919E-7*e[2],-4.79966E-7*e[0]+.917482137087*e[1]-.397776982902*e[2],.397776982902*e[1]+.917482137087*e[2],b)}function na(a,b,c,e,d){var g=(d+c)/2-e;c=(d-c)/2;if(0==g){if(0==c)return null;e=-e/c;if(-1>e||1=e)return null;d=Math.sqrt(e); -e=(-c+d)/(2*g);d=(-c-d)/(2*g);if(-1<=e&&1>=e){if(-1<=d&&1>=d)return null}else if(-1<=d&&1>=d)e=d;else return null}return{x:e,t:a+e*b,df_dt:(2*g*e+c)/b}}function H(a){for(;-180>=a;)a+=360;for(;180h;++h){var l=b.AddDays(h*c);l=d*K(B.Neptune[2],l.tt/365250);if(0==h||l>k)g=h,k=l}b=b.AddDays((g-1)*c);c*=2}}function oa(a){var b=a.AddDays(-30/360*x.Neptune.OrbitalPeriod),c=a.AddDays(.75*x.Neptune.OrbitalPeriod),d=b,g=b,k=-1,h=-1;c=(c.ut-b.ut)/99;for(var l=0;100>l;++l){var m=b.AddDays(l*c),f=K(B.Neptune[2],m.tt/365250);0===l?h=k=f:(f>h&&(h=f,g=m),f=a.tt)return g.time.tt>=a.tt&&g.time.tt=a.tt)return g;throw"Internal error: failed to find Neptune apsis.";}function ha(a){a=360-a;360<=a?a-=360:0>a&&(a+=360);return a}var aa=new Date("2000-01-01T12:00:00Z"),I=2*Math.PI,L=180/Math.PI*3600,pa=-.17-5*Math.log10(648E3/Math.PI),qa=34/60,U,ia,ja,V=function(){this.calcmoon=this.lunar_apsis_iter=this.lunar_apsis_calls=this.longitude_iter=this.longitude_search=this.search=this.search_func= -0};V.prototype.Clone=function(){var a=new V;a.search_func=this.search_func;a.search=this.search;a.longitude_search=this.longitude_search;a.longitude_iter=this.longitude_iter;a.lunar_apsis_calls=this.lunar_apsis_calls;a.lunar_apsis_iter=this.lunar_apsis_iter;a.calcmoon=this.calcmoon;return a};var J=new V;d.GetPerformanceMetrics=function(){return J.Clone()};d.ResetPerformanceMetrics=function(){J=new V};d.Bodies="Sun Moon Mercury Venus Earth Mars Jupiter Saturn Uranus Neptune Pluto".split(" ");var x= -{Mercury:{OrbitalPeriod:87.969},Venus:{OrbitalPeriod:224.701},Earth:{OrbitalPeriod:365.256},Mars:{OrbitalPeriod:686.98},Jupiter:{OrbitalPeriod:4332.589},Saturn:{OrbitalPeriod:10759.22},Uranus:{OrbitalPeriod:30685.4},Neptune:{OrbitalPeriod:60189},Pluto:{OrbitalPeriod:90560}},B={Mercury:[[[[4.40250710144,0,0],[.40989414977,1.48302034195,26087.9031415742],[.050462942,4.47785489551,52175.8062831484],[.00855346844,1.16520322459,78263.70942472259],[.00165590362,4.11969163423,104351.61256629678],[3.4561897E-4, +-f*h*c+n*m*d*c-a*n*l;f=-f*h*a+n*m*d*a+c*n*l;h*=l;n=-l*d*c-a*m;a=-l*d*a+m*c;return 0===b?new A([[k,e,p],[q,v,f],[h,n,a]]):new A([[k,q,h],[e,v,n],[p,f,a]])}function S(a){var b=a.tt/36525,c=15*P(a).ee;a=(.779057273264+.00273781191135448*a.ut+a.ut%1)%1*360;0>a&&(a+=360);b=((c+.014506+((((-3.68E-8*b-2.9956E-5)*b-4.4E-7)*b+1.3915817)*b+4612.156534)*b)/3600+a)%360/15;0>b&&(b+=24);return b}function Z(a,b,c){a=aa(a,b);return[a.rot[0][0]*c[0]+a.rot[1][0]*c[1]+a.rot[2][0]*c[2],a.rot[0][1]*c[0]+a.rot[1][1]*c[1]+ +a.rot[2][1]*c[2],a.rot[0][2]*c[0]+a.rot[1][2]*c[1]+a.rot[2][2]*c[2]]}function aa(a,b){a=P(a);var c=.017453292519943295*a.mobl,e=.017453292519943295*a.tobl,d=4.84813681109536E-6*a.dpsi;a=Math.cos(c);c=Math.sin(c);var k=Math.cos(e),h=Math.sin(e);e=Math.cos(d);var l=Math.sin(d);d=-l*a;var m=-l*c,f=l*k,n=e*a*k+c*h,p=e*c*k-a*h;l*=h;var q=e*a*h-c*k;a=e*c*h+a*k;return 0===b?new A([[e,f,l],[d,n,q],[m,p,a]]):new A([[e,d,m],[f,n,p],[l,q,a]])}function ma(a){if(!(a instanceof Array)||3!==a.length)return!1;for(var b= +0;3>b;++b){if(!(a[b]instanceof Array)||3!==a[b].length)return!1;for(var c=0;3>c;++c)if("number"!==typeof a[b][c])return!1}return!0}function da(a){var b=a[0]*a[0]+a[1]*a[1],c=Math.sqrt(b+a[2]*a[2]);if(0===b){if(0===a[2])throw"Indeterminate sky coordinates";return 0>a[2]?{ra:0,dec:-90,dist:c}:{ra:0,dec:90,dist:c}}var e=Math.atan2(a[1],a[0])/.2617993877991494;0>e&&(e+=24);return new ea(e,Math.atan2(a[2],Math.sqrt(b))/.017453292519943295,c)}function D(a,b){var c=.017453292519943295*a;a=Math.cos(c);c= +Math.sin(c);return[a*b[0]+c*b[1]+0*b[2],-c*b[0]+a*b[1]+0*b[2],0*b[0]+0*b[1]+1*b[2]]}function fa(a,b,c,e,d){var g=b*e+c*d;b=-b*d+c*e;c=Math.sqrt(a*a+g*g);e=0;0e&&(e+=360));return new na(a,g,b,57.29577951308232*Math.atan2(b,c),e)}function K(a,b){var c=1,e=0;a=$jscomp.makeIterator(a);for(var d=a.next();!d.done;d=a.next()){var k=0;d=$jscomp.makeIterator(d.value);for(var h=d.next();!h.done;h=d.next())h=h.value,k+=h[0]*Math.cos(h[1]+b*h[2]);e+=c*k;c*=b}return e} +function x(a,b){var c=b.tt/365250,e=[];a=$jscomp.makeIterator(a);for(var d=a.next();!d.done;d=a.next())e.push(K(d.value,c));c=e[2]*Math.cos(e[1]);e=[c*Math.cos(e[0]),c*Math.sin(e[0]),e[2]*Math.sin(e[1])];return new t(e[0]+4.4036E-7*e[1]-1.90919E-7*e[2],-4.79966E-7*e[0]+.917482137087*e[1]-.397776982902*e[2],.397776982902*e[1]+.917482137087*e[2],b)}function T(a,b,c,e){e/=e+333054.25318;b=x(y[c],b);a.x+=e*b.x;a.y+=e*b.y;a.z+=e*b.z}function oa(a,b,c,e,d){var g=(d+c)/2-e;c=(d-c)/2;if(0==g){if(0==c)return null; +e=-e/c;if(-1>e||1=e)return null;d=Math.sqrt(e);e=(-c+d)/(2*g);d=(-c-d)/(2*g);if(-1<=e&&1>=e){if(-1<=d&&1>=d)return null}else if(-1<=d&&1>=d)e=d;else return null}return{x:e,t:a+e*b,df_dt:(2*g*e+c)/b}}function H(a){for(;-180>=a;)a+=360;for(;180h;++h){var l=b.AddDays(h*c);l=e*K(y.Neptune[2],l.tt/365250);if(0==h||l>k)d=h,k=l}b=b.AddDays((d-1)*c);c*=2}}function pa(a){var b=a.AddDays(-30/360*C.Neptune.OrbitalPeriod),c=a.AddDays(.75*C.Neptune.OrbitalPeriod),d=b,g=b,k=-1,h=-1;c=(c.ut-b.ut)/99;for(var l=0;100>l;++l){var m=b.AddDays(l*c),f=K(y.Neptune[2], +m.tt/365250);0===l?h=k=f:(f>h&&(h=f,g=m),f=a.tt)return g.time.tt>=a.tt&&g.time.tt=a.tt)return g;throw"Internal error: failed to find Neptune apsis.";}function ia(a){a=360-a;360<=a?a-=360:0>a&&(a+=360);return a}var ba=new Date("2000-01-01T12:00:00Z"),I=2*Math.PI,L=180/Math.PI*3600,qa=-.17-5*Math.log10(648E3/Math.PI),ra=34/60,V,ja,ka,W=function(){this.calcmoon=this.lunar_apsis_iter=this.lunar_apsis_calls= +this.longitude_iter=this.longitude_search=this.search=this.search_func=0};W.prototype.Clone=function(){var a=new W;a.search_func=this.search_func;a.search=this.search;a.longitude_search=this.longitude_search;a.longitude_iter=this.longitude_iter;a.lunar_apsis_calls=this.lunar_apsis_calls;a.lunar_apsis_iter=this.lunar_apsis_iter;a.calcmoon=this.calcmoon;return a};var J=new W;d.GetPerformanceMetrics=function(){return J.Clone()};d.ResetPerformanceMetrics=function(){J=new W};d.Bodies="Sun Moon Mercury Venus Earth Mars Jupiter Saturn Uranus Neptune Pluto SSB EMB".split(" "); +var C={Mercury:{OrbitalPeriod:87.969},Venus:{OrbitalPeriod:224.701},Earth:{OrbitalPeriod:365.256},Mars:{OrbitalPeriod:686.98},Jupiter:{OrbitalPeriod:4332.589},Saturn:{OrbitalPeriod:10759.22},Uranus:{OrbitalPeriod:30685.4},Neptune:{OrbitalPeriod:60189},Pluto:{OrbitalPeriod:90560}},y={Mercury:[[[[4.40250710144,0,0],[.40989414977,1.48302034195,26087.9031415742],[.050462942,4.47785489551,52175.8062831484],[.00855346844,1.16520322459,78263.70942472259],[.00165590362,4.11969163423,104351.61256629678],[3.4561897E-4, .77930768443,130439.51570787099],[7.583476E-5,3.71348404924,156527.41884944518]],[[26087.90313685529,0,0],[.01131199811,6.21874197797,26087.9031415742],[.00292242298,3.04449355541,52175.8062831484],[7.5775081E-4,6.08568821653,78263.70942472259],[1.9676525E-4,2.80965111777,104351.61256629678]]],[[[.11737528961,1.98357498767,26087.9031415742],[.02388076996,5.03738959686,52175.8062831484],[.01222839532,3.14159265359,0],[.0054325181,1.79644363964,78263.70942472259],[.0012977877,4.83232503958,104351.61256629678], [3.1866927E-4,1.58088495658,130439.51570787099],[7.963301E-5,4.60972126127,156527.41884944518]],[[.00274646065,3.95008450011,26087.9031415742],[9.9737713E-4,3.14159265359,0]]],[[[.39528271651,0,0],[.07834131818,6.19233722598,26087.9031415742],[.00795525558,2.95989690104,52175.8062831484],[.00121281764,6.01064153797,78263.70942472259],[2.1921969E-4,2.77820093972,104351.61256629678],[4.354065E-5,5.82894543774,130439.51570787099]],[[.0021734774,4.65617158665,26087.9031415742],[4.4141826E-4,1.42385544001, 52175.8062831484]]]],Venus:[[[[3.17614666774,0,0],[.01353968419,5.59313319619,10213.285546211],[8.9891645E-4,5.30650047764,20426.571092422],[5.477194E-5,4.41630661466,7860.4193924392],[3.455741E-5,2.6996444782,11790.6290886588],[2.372061E-5,2.99377542079,3930.2096962196],[1.317168E-5,5.18668228402,26.2983197998],[1.664146E-5,4.25018630147,1577.3435424478],[1.438387E-5,4.15745084182,9683.5945811164],[1.200521E-5,6.15357116043,30639.856638633]],[[10213.28554621638,0,0],[9.5617813E-4,2.4640651111,10213.285546211], @@ -94,7 +94,7 @@ fa(1,g.AddDays(-2*c),4*c);if(b.time.tt>=a.tt)return g.time.tt>=a.tt&&g.time.ttf&&(f+=360),360<=f&& -(f-=360));n=57.29577951308232*Math.atan2(p,n);p=e;if(g&&(e=n,g=d.Refraction(g,90-n),n-=g,0l;++l)g.push((b[l]-e*a[l])/q*c+a[l]*p);p=Math.sqrt(g[0]*g[0]+g[1]*g[1]);0c&&(c+=24),24<=c&&(c-=24)):c=0;p=57.29577951308232*Math.atan2(g[2],p)}return new ra(f,90-n,c,p)};var sa=function(a,b,c){this.latitude= -a;this.longitude=b;this.height=c};d.MakeObserver=function(a,b,c){return new sa(a,b,c||0)};d.SunPosition=function(a){a=d.MakeTime(a).AddDays(-.005775518331089121);var b=C(B.Earth,a);b=R(0,[-b.x,-b.y,-b.z],a.tt);b=Y(a,0,b);a=.017453292519943295*P(a).tobl;return ea(b[0],b[1],b[2],Math.cos(a),Math.sin(a))};d.Equator=function(a,b,c,e,g){b=d.MakeTime(b);var k=S(b),h=.017453292519943295*c.latitude,l=Math.sin(h);h=Math.cos(h);var m=1/Math.sqrt(h*h+.9933056020041345*l*l),f=c.height/1E3,n=6378.1366*m+f;c=.017453292519943295* -(15*k+c.longitude);c=Y(b,-1,[n*h*Math.cos(c)/1.4959787069098932E8,n*h*Math.sin(c)/1.4959787069098932E8,(6335.438815127603*m+f)*l/1.4959787069098932E8]);c=R(b.tt,c,0);a=d.GeoVector(a,b,g);a=[a.x-c[0],a.y-c[1],a.z-c[2]];if(!e)return ca(a);e=R(0,a,b.tt);e=Y(b,0,e);return ca(e)};d.Ecliptic=function(a,b,c){void 0===U&&(U=.017453292519943295*P(d.MakeTime(aa)).mobl,ia=Math.cos(U),ja=Math.sin(U));return ea(a,b,c,ia,ja)};d.GeoMoon=function(a){a=d.MakeTime(a);var b=F(a),c=b.distance_au*Math.cos(b.geo_eclip_lat); -b=[c*Math.cos(b.geo_eclip_lon),c*Math.sin(b.geo_eclip_lon),b.distance_au*Math.sin(b.geo_eclip_lat)];var e=.017453292519943295*E(a);c=Math.cos(e);e=Math.sin(e);b=R(a.tt,[b[0],b[1]*c-b[2]*e,b[1]*e+b[2]*c],0);return new t(b[0],b[1],b[2],a)};d.HelioVector=function(a,b){b=d.MakeTime(b);if(a in B)return C(B[a],b);if(a in ka){a:{var c=$jscomp.makeIterator(ka[a]);for(a=c.next();!a.done;a=c.next()){a=a.value;var e=a.tt;var g=a.tt+a.ndays;e=(2*b.tt-(g+e))/(g-e);if(-1<=e&&1>=e){var k=[];for(g=0;3>g;++g){var h= -1;var l=a.coeff[0][g];var m=e;l+=a.coeff[1][g]*m;for(c=2;cl;++l){g=d.HelioVector(a,h);c&&(e=C(B.Earth,h));k=new t(g.x-e.x,g.y-e.y,g.z-e.z,b);if("Sun"===a)return k;var m=b.AddDays(-k.Length()/173.1446326846693);g=Math.abs(m.tt-h.tt);if(1E-9>g)return k;h=m}throw"Light-travel time solver did not converge: dt="+g;};d.Search=function(a,b,c,e){function g(b){++J.search_func; -return a(b)}var k=Math.abs((e&&e.dt_tolerance_seconds||1)/86400);++J.search;var h=e&&e.init_f1||g(b),l=e&&e.init_f2||g(c),m,f=0;e=e&&e.iter_limit||20;for(var n=!0;;){if(++f>e)throw"Excessive iteration in Search()";var p=new M(b.ut+.5*(c.ut-b.ut)),q=p.ut-b.ut;if(Math.abs(q)(q.ut-b.ut)*(q.ut- -c.ut)&&0>(r.ut-b.ut)*(r.ut-c.ut))){v=g(q);var u=g(r);if(0>v&&0<=u){h=v;l=u;b=q;c=r;m=t;n=!1;continue}}}}if(0>h&&0<=m)c=p,l=m;else if(0>m&&0<=l)b=p,h=m;else return null}};d.SearchSunLongitude=function(a,b,c){b=d.MakeTime(b);c=b.AddDays(c);return d.Search(function(b){b=d.SunPosition(b);return H(b.elon-a)},b,c)};d.LongitudeFromSun=function(a,b){if("Earth"===a)throw"The Earth does not have a longitude as seen from itself.";b=d.MakeTime(b);a=d.GeoVector(a,b,!0);a=d.Ecliptic(a.x,a.y,a.z);b=d.GeoVector("Sun", -b,!0);b=d.Ecliptic(b.x,b.y,b.z);for(b=a.elon-b.elon;0>b;)b+=360;for(;360<=b;)b-=360;return b};d.AngleFromSun=function(a,b){if("Earth"==a)throw"The Earth does not have an angle as seen from itself.";var c=d.GeoVector("Sun",b,!0);a=d.GeoVector(a,b,!0);return A(c,a)};d.EclipticLongitude=function(a,b){if("Sun"===a)throw"Cannot calculate heliocentric longitude of the Sun.";a=d.HelioVector(a,b);return d.Ecliptic(a.x,a.y,a.z).elon};var ta=function(a,b,c,d,g,k,h,l){this.time=a;this.mag=b;this.phase_angle= -c;this.phase_fraction=(1+Math.cos(.017453292519943295*c))/2;this.helio_dist=d;this.geo_dist=g;this.gc=k;this.hc=h;this.ring_tilt=l};d.Illumination=function(a,b){if("Earth"===a)throw"The illumination of the Earth is not defined.";var c=d.MakeTime(b),e=C(B.Earth,c);if("Sun"===a){var g=new t(-e.x,-e.y,-e.z,c);b=new t(0,0,0,c);e=0}else"Moon"===a?(g=d.GeoMoon(c),b=new t(e.x+g.x,e.y+g.y,e.z+g.z,c)):(b=d.HelioVector(a,b),g=new t(b.x-e.x,b.y-e.y,b.z-e.z,c)),e=A(g,b);var k=g.Length(),h=b.Length(),l=null;if("Sun"=== -a)a=pa+5*Math.log10(k);else if("Moon"===a){a=.017453292519943295*e;var m=a*a;a=-12.717+1.49*Math.abs(a)+.0431*m*m;a+=5*Math.log10(k/.002573570052980638*h)}else if("Saturn"===a)a=e,l=d.Ecliptic(g.x,g.y,g.z),m=.017453292519943295*l.elat,l=Math.asin(Math.sin(m)*Math.cos(.4897393881096089)-Math.cos(m)*Math.sin(.4897393881096089)*Math.sin(.017453292519943295*l.elon-.017453292519943295*(169.51+3.82E-5*c.tt))),m=Math.sin(Math.abs(l)),a=-9+.044*a+m*(-2.6+1.2*m)+5*Math.log10(h*k),l*=57.29577951308232;else{var f= -m=0,n=0;switch(a){case "Mercury":a=-.6;m=4.98;f=-4.88;n=3.02;break;case "Venus":163.6>e?(a=-4.47,m=1.03,f=.57,n=.13):(a=.98,m=-1.02);break;case "Mars":a=-1.52;m=1.6;break;case "Jupiter":a=-9.4;m=.5;break;case "Uranus":a=-7.19;m=.25;break;case "Neptune":a=-6.87;break;case "Pluto":a=-1;m=4;break;default:throw"VisualMagnitude: unsupported body "+a;}var p=e/100;a=a+p*(m+p*(f+p*n))+5*Math.log10(h*k)}return new ta(c,a,e,h,k,g,b,l)};d.SearchRelativeLongitude=function(a,b,c){function e(c){var e=d.EclipticLongitude(a, -c);c=d.EclipticLongitude("Earth",c);return H(k*(c-e)-b)}var g=x[a];if(!g)throw"Cannot search relative longitude because body is not a planet: "+a;if("Earth"===a)throw"Cannot search relative longitude for the Earth (it is always 0)";var k=g.OrbitalPeriod>x.Earth.OrbitalPeriod?1:-1;g=D(a);c=d.MakeTime(c);var h=e(c);0l;++l){++J.longitude_iter;var m=-h/360*g;c=c.AddDays(m);if(1>86400*Math.abs(m))return c;m=h;h=e(c);30>Math.abs(m)&&m!==h&&(m/=m-h,.5m&&(g*=m))}throw"Relative longitude search failed to converge for "+a+" near "+c.toString()+" (error_angle = "+h+").";};d.MoonPhase=function(a){return d.LongitudeFromSun("Moon",a)};d.SearchMoonPhase=function(a,b,c){function e(b){b=d.MoonPhase(b);return H(b-a)}b=d.MakeTime(b);var g=e(b);0c)return null;c=Math.min(c,k+.9);g=b.AddDays(g);b=b.AddDays(c);return d.Search(e,g,b)};var ua=function(a,b){this.quarter=a;this.time=b};d.SearchMoonQuarter=function(a){var b= -d.MoonPhase(a);b=(Math.floor(b/90)+1)%4;return(a=d.SearchMoonPhase(90*b,a,10))&&new ua(b,a)};d.NextMoonQuarter=function(a){a=new Date(a.time.date.getTime()+5184E5);return d.SearchMoonQuarter(a)};d.SearchRiseSet=function(a,b,c,e,g){function k(e){var f=d.Equator(a,e,b,!0,!0);e=d.Horizon(e,b,f.ra,f.dec).altitude+h/f.dist*57.29577951308232+qa;return c*e}var h={Sun:.0046505,Moon:1.15717E-5}[a]||0;if("Earth"===a)throw"Cannot find rise or set time of the Earth.";if(1===c){var l=12;var m=0}else if(-1===c)l= -0,m=12;else throw"Astronomy.SearchRiseSet: Invalid direction parameter "+c+" -- must be +1 or -1";e=d.MakeTime(e);var f=k(e);var n;if(0=f&&0=e.ut+g)return null;p=f.time;f=k(f.time);n=k(q.time)}};var va=function(a,b){this.time=a;this.hor=b};d.SearchHourAngle= -function(a,b,c,e){e=d.MakeTime(e);var g=0;if("Earth"===a)throw"Cannot search for hour angle of the Earth.";if(0>c||24<=c)throw"Invalid hour angle "+c;for(;;){++g;var k=S(e),h=d.Equator(a,e,b,!0,!0);k=(c+h.ra-b.longitude/15-k)%24;1===g?0>k&&(k+=24):-12>k?k+=24:123600*Math.abs(k))return a=d.Horizon(e,b,h.ra,h.dec,"normal"),new va(e,a);e=e.AddDays(k/24*.9972695717592592)}};var wa=function(a,b,c,d){this.mar_equinox=a;this.jun_solstice=b;this.sep_equinox=c;this.dec_solstice=d};d.Seasons= -function(a){function b(b,c,e){c=new Date(Date.UTC(a,c-1,e));b=d.SearchSunLongitude(b,c,4);if(!b)throw"Cannot find season change near "+c.toISOString();return b}a instanceof Date&&(a=a.getUTCFullYear());var c=b(0,3,19),e=b(90,6,19),g=b(180,9,21),k=b(270,12,20);return new wa(c,e,g,k)};var xa=function(a,b,c,d){this.time=a;this.visibility=b;this.elongation=c;this.ecliptic_separation=d};d.Elongation=function(a,b){b=d.MakeTime(b);var c=d.LongitudeFromSun(a,b);if(180=++g;){var k=d.EclipticLongitude(a,b),h=d.EclipticLongitude("Earth",b);h=H(k-h);var l=k=void 0,m=l=k=void 0;h>=-e.s1&&h<+e.s1?(m=0,k=+e.s1,l=+e.s2):h>=+e.s2|| -h<-e.s2?(m=0,k=-e.s2,l=-e.s1):0<=h?(m=-D(a)/4,k=+e.s1,l=+e.s2):(m=-D(a)/4,k=-e.s2,l=-e.s1);h=b.AddDays(m);k=d.SearchRelativeLongitude(a,k,h);l=d.SearchRelativeLongitude(a,l,k);h=c(k);if(0<=h)throw"SearchMaxElongation: internal error: m1 = "+h;m=c(l);if(0>=m)throw"SearchMaxElongation: internal error: m2 = "+m;h=d.Search(c,k,l,{init_f1:h,init_f2:m,dt_tolerance_seconds:10});if(!h)throw"SearchMaxElongation: failed search iter "+g+" (t1="+k.toString()+", t2="+l.toString()+")";if(h.tt>=b.tt)return d.Elongation(a, -h);b=l.AddDays(1)}throw"SearchMaxElongation: failed to find event after 2 tries.";};d.SearchPeakMagnitude=function(a,b){function c(b){var c=b.AddDays(-.005);b=b.AddDays(.005);c=d.Illumination(a,c).mag;return(d.Illumination(a,b).mag-c)/.01}if("Venus"!==a)throw"SearchPeakMagnitude currently works for Venus only.";b=d.MakeTime(b);for(var e=0;2>=++e;){var g=d.EclipticLongitude(a,b),k=d.EclipticLongitude("Earth",b);k=H(g-k);var h=g=void 0,l=h=g=void 0;-10<=k&&10>k?(l=0,g=10,h=30):30<=k||-30>k?(l=0,g=-30, -h=-10):0<=k?(l=-D(a)/4,g=10,h=30):(l=-D(a)/4,g=-30,h=-10);k=b.AddDays(l);g=d.SearchRelativeLongitude(a,g,k);h=d.SearchRelativeLongitude(a,h,g);k=c(g);if(0<=k)throw"SearchPeakMagnitude: internal error: m1 = "+k;l=c(h);if(0>=l)throw"SearchPeakMagnitude: internal error: m2 = "+l;k=d.Search(c,g,h,{init_f1:k,init_f2:l,dt_tolerance_seconds:10});if(!k)throw"SearchPeakMagnitude: failed search iter "+e+" (t1="+g.toString()+", t2="+h.toString()+")";if(k.tt>=b.tt)return d.Illumination(a,k);b=h.AddDays(1)}throw"SearchPeakMagnitude: failed to find event after 2 tries."; -};var T=function(a,b,c){this.time=a;this.kind=b;this.dist_au=c;this.dist_km=1.4959787069098932E8*c};d.SearchLunarApsis=function(a){function b(a){var b=a.AddDays(-5E-4);a=a.AddDays(5E-4);b=F(b).distance_au;return(F(a).distance_au-b)/.001}function c(a){return-b(a)}a=d.MakeTime(a);var e=b(a);++J.lunar_apsis_calls;for(var g=0;59.061176>5*g;++g){++J.lunar_apsis_iter;var k=a.AddDays(5),h=b(k);if(0>=e*h){if(0>e||0h){a=d.Search(c,a,k,{init_f1:-e,init_f2:-h});if(null==a)throw"SearchLunarApsis INTERNAL ERROR: apogee search failed!";e=F(a).distance_au;return new T(a,1,e)}throw"SearchLunarApsis INTERNAL ERROR: cannot classify apsis event!";}a=k;e=h}throw"SearchLunarApsis INTERNAL ERROR: could not find apsis within 2 synodic months of start date.";};d.NextLunarApsis=function(a){var b=d.SearchLunarApsis(a.time.AddDays(11));if(1!==b.kind+a.kind)throw"NextLunarApsis INTERNAL ERROR: did not find alternating apogee/perigee: prev="+ -a.kind+" @ "+a.time.toString()+", next="+b.kind+" @ "+b.time.toString();return b};d.SearchPlanetApsis=function(a,b){function c(b){var c=b.AddDays(-5E-4);b=b.AddDays(5E-4);c=d.HelioDistance(a,c);return(d.HelioDistance(a,b)-c)/.001}function e(a){return-c(a)}if("Neptune"===a)return oa(b);for(var g=x[a].OrbitalPeriod,k=g/6,h=c(b),l=0;l*k<2*g;++l){var m=b.AddDays(k),f=c(m);if(0>=h*f){g=k=void 0;if(0>h||0f)k=e,g=1;else throw"Internal error with slopes in SearchPlanetApsis";b=d.Search(k, -b,m,1);if(null==b)throw"Failed to find slope transition in planetary apsis search.";h=d.HelioDistance(a,b);return new T(b,g,h)}b=m;h=f}throw"Internal error: should have found planetary apsis within 2 orbital periods.";};d.NextPlanetApsis=function(a,b){if(0!==b.kind&&1!==b.kind)throw"Invalid apsis kind: "+b.kind;var c=b.time.AddDays(.25*x[a].OrbitalPeriod);a=d.SearchPlanetApsis(a,c);if(1!==a.kind+b.kind)throw"Internal error: previous apsis was "+b.kind+", but found "+a.kind+" for next apsis.";return a}; -d.InverseRotation=function(a){return new z([[a.rot[0][0],a.rot[1][0],a.rot[2][0]],[a.rot[0][1],a.rot[1][1],a.rot[2][1]],[a.rot[0][2],a.rot[1][2],a.rot[2][2]]])};d.CombineRotation=function(a,b){return new z([[b.rot[0][0]*a.rot[0][0]+b.rot[1][0]*a.rot[0][1]+b.rot[2][0]*a.rot[0][2],b.rot[0][1]*a.rot[0][0]+b.rot[1][1]*a.rot[0][1]+b.rot[2][1]*a.rot[0][2],b.rot[0][2]*a.rot[0][0]+b.rot[1][2]*a.rot[0][1]+b.rot[2][2]*a.rot[0][2]],[b.rot[0][0]*a.rot[1][0]+b.rot[1][0]*a.rot[1][1]+b.rot[2][0]*a.rot[1][2],b.rot[0][1]* -a.rot[1][0]+b.rot[1][1]*a.rot[1][1]+b.rot[2][1]*a.rot[1][2],b.rot[0][2]*a.rot[1][0]+b.rot[1][2]*a.rot[1][1]+b.rot[2][2]*a.rot[1][2]],[b.rot[0][0]*a.rot[2][0]+b.rot[1][0]*a.rot[2][1]+b.rot[2][0]*a.rot[2][2],b.rot[0][1]*a.rot[2][0]+b.rot[1][1]*a.rot[2][1]+b.rot[2][1]*a.rot[2][2],b.rot[0][2]*a.rot[2][0]+b.rot[1][2]*a.rot[2][1]+b.rot[2][2]*a.rot[2][2]]])};d.VectorFromSphere=function(a,b){var c=.017453292519943295*a.lat,d=.017453292519943295*a.lon,g=a.dist*Math.cos(c);return new t(g*Math.cos(d),g*Math.sin(d), -a.dist*Math.sin(c),b)};d.VectorFromEquator=function(a,b){return d.VectorFromSphere(new W(a.dec,15*a.ra,a.dist),b)};d.EquatorFromVector=function(a){a=d.SphereFromVector(a);return new da(a.lon/15,a.lat,a.dist)};d.SphereFromVector=function(a){var b=a.x*a.x+a.y*a.y,c=Math.sqrt(b+a.z*a.z);if(0===b){if(0===a.z)throw"Zero-length vector not allowed.";var d=0;a=0>a.z?-90:90}else d=57.29577951308232*Math.atan2(a.y,a.x),0>d&&(d+=360),a=57.29577951308232*Math.atan2(a.z,Math.sqrt(b));return new W(a,d,c)};d.HorizonFromVector= -function(a,b){a=d.SphereFromVector(a);a.lon=ha(a.lon);a.lat+=d.Refraction(b,a.lat);return a};d.VectorFromHorizon=function(a,b,c){var e=ha(a.lon);c=a.lat+d.InverseRefraction(c,a.lat);a=new W(c,e,a.dist);return d.VectorFromSphere(a,b)};d.Refraction=function(a,b){if(-90>b||90c&&(c=-1);c=1.02/Math.tan(.017453292519943295*(c+10.3/(c+5.11)))/60;"normal"===a&&-1>b&&(c*=(b+90)/89)}else c=0;return c};d.InverseRefraction=function(a,b){if(-90>b||90Math.abs(e))return c-b;c-=e}};d.RotateVector=function(a,b){return new t(a.rot[0][0]*b.x+a.rot[1][0]*b.y+a.rot[2][0]*b.z,a.rot[0][1]*b.x+a.rot[1][1]*b.y+a.rot[2][1]*b.z,a.rot[0][2]*b.x+a.rot[1][2]*b.y+a.rot[2][2]*b.z,b.t)};d.Rotation_EQJ_ECL=function(){return new z([[1,0,0],[0,.9174821430670688,-.3977769691083922],[0,.3977769691083922,.9174821430670688]])};d.Rotation_ECL_EQJ=function(){return new z([[1,0,0],[0,.9174821430670688, -.3977769691083922],[0,-.3977769691083922,.9174821430670688]])};d.Rotation_EQJ_EQD=function(a){var b=X(0,a.tt);a=Z(a,0);return d.CombineRotation(b,a)};d.Rotation_EQD_EQJ=function(a){var b=Z(a,1);a=X(a.tt,0);return d.CombineRotation(b,a)};d.Rotation_EQD_HOR=function(a,b){var c=Math.sin(.017453292519943295*b.latitude),d=Math.cos(.017453292519943295*b.latitude),g=Math.sin(.017453292519943295*b.longitude),k=Math.cos(.017453292519943295*b.longitude);b=[d*k,d*g,c];c=[-c*k,-c*g,d];g=[g,-k,0];a=-15*S(a);b= -G(a,b);c=G(a,c);a=G(a,g);return new z([[c[0],a[0],b[0]],[c[1],a[1],b[1]],[c[2],a[2],b[2]]])};d.Rotation_HOR_EQD=function(a,b){a=d.Rotation_EQD_HOR(a,b);return d.InverseRotation(a)};d.Rotation_HOR_EQJ=function(a,b){b=d.Rotation_HOR_EQD(a,b);a=d.Rotation_EQD_EQJ(a);return d.CombineRotation(b,a)};d.Rotation_EQJ_HOR=function(a,b){a=d.Rotation_HOR_EQJ(a,b);return d.InverseRotation(a)};d.Rotation_EQD_ECL=function(a){a=d.Rotation_EQD_EQJ(a);var b=d.Rotation_EQJ_ECL();return d.CombineRotation(a,b)};d.Rotation_ECL_EQD= -function(a){a=d.Rotation_EQD_ECL(a);return d.InverseRotation(a)};d.Rotation_ECL_HOR=function(a,b){var c=d.Rotation_ECL_EQD(a);a=d.Rotation_EQD_HOR(a,b);return d.CombineRotation(c,a)};d.Rotation_HOR_ECL=function(a,b){a=d.Rotation_ECL_HOR(a,b);return d.InverseRotation(a)}})("undefined"===typeof exports?this.Astronomy={}:exports); +647,0,4]},{nals:[-1,1,0,1,1],cls:[1314,0,0,-700,0,0]},{nals:[0,-2,2,-2,1],cls:[-1283,0,0,672,0,0]},{nals:[1,0,2,2,1],cls:[-1331,0,8,663,0,4]},{nals:[-2,0,2,2,2],cls:[1383,0,-2,-594,0,-2]},{nals:[-1,0,0,0,2],cls:[1405,0,4,-610,0,2]},{nals:[1,1,2,-2,2],cls:[1290,0,0,-556,0,0]}],Q,t=function(a,b,c,d){this.x=a;this.y=b;this.z=c;this.t=d};t.prototype.Length=function(){return Math.sqrt(this.x*this.x+this.y*this.y+this.z*this.z)};var X=function(a,b,c){this.lat=a;this.lon=b;this.dist=c};d.MakeSpherical=function(a, +b,c){return new X(a,b,c)};var ea=function(a,b,c){this.ra=a;this.dec=b;this.dist=c},A=function(a){this.rot=a};d.MakeRotation=function(a){if(!ma(a))throw"Argument must be a [3][3] array of numbers";return new A(a)};var sa=function(a,b,c,d){this.azimuth=a;this.altitude=b;this.ra=c;this.dec=d},na=function(a,b,c,d,g){this.ex=a;this.ey=b;this.ez=c;this.elat=d;this.elon=g};d.Horizon=function(a,b,c,e,g){a=d.MakeTime(a);var k=Math.sin(.017453292519943295*b.latitude),h=Math.cos(.017453292519943295*b.latitude), +l=Math.sin(.017453292519943295*b.longitude),m=Math.cos(.017453292519943295*b.longitude);b=Math.sin(.017453292519943295*e);var f=Math.cos(.017453292519943295*e),n=Math.sin(.2617993877991494*c),p=Math.cos(.2617993877991494*c),q=[h*m,h*l,k];k=[-k*m,-k*l,h];l=[l,-m,0];h=-15*S(a);a=D(h,q);q=D(h,k);l=D(h,l);b=[f*p,f*n,b];n=b[0]*a[0]+b[1]*a[1]+b[2]*a[2];q=b[0]*q[0]+b[1]*q[1]+b[2]*q[2];l=b[0]*l[0]+b[1]*l[1]+b[2]*l[2];p=Math.sqrt(q*q+l*l);f=0;0f&&(f+=360),360<=f&& +(f-=360));n=57.29577951308232*Math.atan2(p,n);p=e;if(g&&(e=n,g=d.Refraction(g,90-n),n-=g,0l;++l)g.push((b[l]-e*a[l])/q*c+a[l]*p);p=Math.sqrt(g[0]*g[0]+g[1]*g[1]);0c&&(c+=24),24<=c&&(c-=24)):c=0;p=57.29577951308232*Math.atan2(g[2],p)}return new sa(f,90-n,c,p)};var ta=function(a,b,c){this.latitude= +a;this.longitude=b;this.height=c};d.MakeObserver=function(a,b,c){return new ta(a,b,c||0)};d.SunPosition=function(a){a=d.MakeTime(a).AddDays(-.005775518331089121);var b=x(y.Earth,a);b=R(0,[-b.x,-b.y,-b.z],a.tt);b=Z(a,0,b);a=.017453292519943295*P(a).tobl;return fa(b[0],b[1],b[2],Math.cos(a),Math.sin(a))};d.Equator=function(a,b,c,e,g){b=d.MakeTime(b);var k=S(b),h=.017453292519943295*c.latitude,l=Math.sin(h);h=Math.cos(h);var m=1/Math.sqrt(h*h+.9933056020041345*l*l),f=c.height/1E3,n=6378.1366*m+f;c=.017453292519943295* +(15*k+c.longitude);c=Z(b,-1,[n*h*Math.cos(c)/1.4959787069098932E8,n*h*Math.sin(c)/1.4959787069098932E8,(6335.438815127603*m+f)*l/1.4959787069098932E8]);c=R(b.tt,c,0);a=d.GeoVector(a,b,g);a=[a.x-c[0],a.y-c[1],a.z-c[2]];if(!e)return da(a);e=R(0,a,b.tt);e=Z(b,0,e);return da(e)};d.Ecliptic=function(a,b,c){void 0===V&&(V=.017453292519943295*P(d.MakeTime(ba)).mobl,ja=Math.cos(V),ka=Math.sin(V));return fa(a,b,c,ja,ka)};d.GeoMoon=function(a){a=d.MakeTime(a);var b=G(a),c=b.distance_au*Math.cos(b.geo_eclip_lat); +b=[c*Math.cos(b.geo_eclip_lon),c*Math.sin(b.geo_eclip_lon),b.distance_au*Math.sin(b.geo_eclip_lat)];var e=.017453292519943295*F(a);c=Math.cos(e);e=Math.sin(e);b=R(a.tt,[b[0],b[1]*c-b[2]*e,b[1]*e+b[2]*c],0);return new t(b[0],b[1],b[2],a)};d.HelioVector=function(a,b){b=d.MakeTime(b);if(a in y)return x(y[a],b);if(a in la){a:{var c=$jscomp.makeIterator(la[a]);for(a=c.next();!a.done;a=c.next()){a=a.value;var e=a.tt;var g=a.tt+a.ndays;e=(2*b.tt-(g+e))/(g-e);if(-1<=e&&1>=e){var k=[];for(g=0;3>g;++g){var h= +1;var l=a.coeff[0][g];var m=e;l+=a.coeff[1][g]*m;for(c=2;cl;++l){g=d.HelioVector(a,h);c&&(e=x(y.Earth, +h));k=new t(g.x-e.x,g.y-e.y,g.z-e.z,b);if("Sun"===a)return k;var m=b.AddDays(-k.Length()/173.1446326846693);g=Math.abs(m.tt-h.tt);if(1E-9>g)return k;h=m}throw"Light-travel time solver did not converge: dt="+g;};d.Search=function(a,b,c,e){function g(b){++J.search_func;return a(b)}var k=Math.abs((e&&e.dt_tolerance_seconds||1)/86400);++J.search;var h=e&&e.init_f1||g(b),l=e&&e.init_f2||g(c),m,f=0;e=e&&e.iter_limit||20;for(var n=!0;;){if(++f>e)throw"Excessive iteration in Search()";var p=new M(b.ut+.5* +(c.ut-b.ut)),q=p.ut-b.ut;if(Math.abs(q)(q.ut-b.ut)*(q.ut-c.ut)&&0>(r.ut-b.ut)*(r.ut-c.ut))){v=g(q);var u=g(r);if(0>v&&0<=u){h=v;l=u;b=q;c=r;m=t;n=!1;continue}}}}if(0>h&&0<=m)c=p,l=m;else if(0>m&&0<=l)b=p,h=m;else return null}};d.SearchSunLongitude=function(a,b,c){b=d.MakeTime(b);c=b.AddDays(c); +return d.Search(function(b){b=d.SunPosition(b);return H(b.elon-a)},b,c)};d.LongitudeFromSun=function(a,b){if("Earth"===a)throw"The Earth does not have a longitude as seen from itself.";b=d.MakeTime(b);a=d.GeoVector(a,b,!0);a=d.Ecliptic(a.x,a.y,a.z);b=d.GeoVector("Sun",b,!0);b=d.Ecliptic(b.x,b.y,b.z);for(b=a.elon-b.elon;0>b;)b+=360;for(;360<=b;)b-=360;return b};d.AngleFromSun=function(a,b){if("Earth"==a)throw"The Earth does not have an angle as seen from itself.";var c=d.GeoVector("Sun",b,!0);a=d.GeoVector(a, +b,!0);return B(c,a)};d.EclipticLongitude=function(a,b){if("Sun"===a)throw"Cannot calculate heliocentric longitude of the Sun.";a=d.HelioVector(a,b);return d.Ecliptic(a.x,a.y,a.z).elon};var ua=function(a,b,c,d,g,k,h,l){this.time=a;this.mag=b;this.phase_angle=c;this.phase_fraction=(1+Math.cos(.017453292519943295*c))/2;this.helio_dist=d;this.geo_dist=g;this.gc=k;this.hc=h;this.ring_tilt=l};d.Illumination=function(a,b){if("Earth"===a)throw"The illumination of the Earth is not defined.";var c=d.MakeTime(b), +e=x(y.Earth,c);if("Sun"===a){var g=new t(-e.x,-e.y,-e.z,c);b=new t(0,0,0,c);e=0}else"Moon"===a?(g=d.GeoMoon(c),b=new t(e.x+g.x,e.y+g.y,e.z+g.z,c)):(b=d.HelioVector(a,b),g=new t(b.x-e.x,b.y-e.y,b.z-e.z,c)),e=B(g,b);var k=g.Length(),h=b.Length(),l=null;if("Sun"===a)a=qa+5*Math.log10(k);else if("Moon"===a){a=.017453292519943295*e;var m=a*a;a=-12.717+1.49*Math.abs(a)+.0431*m*m;a+=5*Math.log10(k/.002573570052980638*h)}else if("Saturn"===a)a=e,l=d.Ecliptic(g.x,g.y,g.z),m=.017453292519943295*l.elat,l=Math.asin(Math.sin(m)* +Math.cos(.4897393881096089)-Math.cos(m)*Math.sin(.4897393881096089)*Math.sin(.017453292519943295*l.elon-.017453292519943295*(169.51+3.82E-5*c.tt))),m=Math.sin(Math.abs(l)),a=-9+.044*a+m*(-2.6+1.2*m)+5*Math.log10(h*k),l*=57.29577951308232;else{var f=m=0,n=0;switch(a){case "Mercury":a=-.6;m=4.98;f=-4.88;n=3.02;break;case "Venus":163.6>e?(a=-4.47,m=1.03,f=.57,n=.13):(a=.98,m=-1.02);break;case "Mars":a=-1.52;m=1.6;break;case "Jupiter":a=-9.4;m=.5;break;case "Uranus":a=-7.19;m=.25;break;case "Neptune":a= +-6.87;break;case "Pluto":a=-1;m=4;break;default:throw"VisualMagnitude: unsupported body "+a;}var p=e/100;a=a+p*(m+p*(f+p*n))+5*Math.log10(h*k)}return new ua(c,a,e,h,k,g,b,l)};d.SearchRelativeLongitude=function(a,b,c){function e(c){var e=d.EclipticLongitude(a,c);c=d.EclipticLongitude("Earth",c);return H(k*(c-e)-b)}var g=C[a];if(!g)throw"Cannot search relative longitude because body is not a planet: "+a;if("Earth"===a)throw"Cannot search relative longitude for the Earth (it is always 0)";var k=g.OrbitalPeriod> +C.Earth.OrbitalPeriod?1:-1;g=E(a);c=d.MakeTime(c);var h=e(c);0l;++l){++J.longitude_iter;var m=-h/360*g;c=c.AddDays(m);if(1>86400*Math.abs(m))return c;m=h;h=e(c);30>Math.abs(m)&&m!==h&&(m/=m-h,.5m&&(g*=m))}throw"Relative longitude search failed to converge for "+a+" near "+c.toString()+" (error_angle = "+h+").";};d.MoonPhase=function(a){return d.LongitudeFromSun("Moon",a)};d.SearchMoonPhase=function(a,b,c){function e(b){b=d.MoonPhase(b);return H(b- +a)}b=d.MakeTime(b);var g=e(b);0c)return null;c=Math.min(c,k+.9);g=b.AddDays(g);b=b.AddDays(c);return d.Search(e,g,b)};var va=function(a,b){this.quarter=a;this.time=b};d.SearchMoonQuarter=function(a){var b=d.MoonPhase(a);b=(Math.floor(b/90)+1)%4;return(a=d.SearchMoonPhase(90*b,a,10))&&new va(b,a)};d.NextMoonQuarter=function(a){a=new Date(a.time.date.getTime()+5184E5);return d.SearchMoonQuarter(a)};d.SearchRiseSet=function(a,b,c,e,g){function k(e){var f= +d.Equator(a,e,b,!0,!0);e=d.Horizon(e,b,f.ra,f.dec).altitude+h/f.dist*57.29577951308232+ra;return c*e}var h={Sun:.0046505,Moon:1.15717E-5}[a]||0;if("Earth"===a)throw"Cannot find rise or set time of the Earth.";if(1===c){var l=12;var m=0}else if(-1===c)l=0,m=12;else throw"Astronomy.SearchRiseSet: Invalid direction parameter "+c+" -- must be +1 or -1";e=d.MakeTime(e);var f=k(e);var n;if(0= +f&&0=e.ut+g)return null;p=f.time;f=k(f.time);n=k(q.time)}};var wa=function(a,b){this.time=a;this.hor=b};d.SearchHourAngle=function(a,b,c,e){e=d.MakeTime(e);var g=0;if("Earth"===a)throw"Cannot search for hour angle of the Earth.";if(0>c||24<=c)throw"Invalid hour angle "+c;for(;;){++g;var k=S(e),h=d.Equator(a,e,b,!0,!0);k=(c+h.ra-b.longitude/15-k)%24;1===g?0>k&& +(k+=24):-12>k?k+=24:123600*Math.abs(k))return a=d.Horizon(e,b,h.ra,h.dec,"normal"),new wa(e,a);e=e.AddDays(k/24*.9972695717592592)}};var xa=function(a,b,c,d){this.mar_equinox=a;this.jun_solstice=b;this.sep_equinox=c;this.dec_solstice=d};d.Seasons=function(a){function b(b,c,e){c=new Date(Date.UTC(a,c-1,e));b=d.SearchSunLongitude(b,c,4);if(!b)throw"Cannot find season change near "+c.toISOString();return b}a instanceof Date&&(a=a.getUTCFullYear());var c=b(0,3,19),e=b(90,6,19),g=b(180, +9,21),k=b(270,12,20);return new xa(c,e,g,k)};var ya=function(a,b,c,d){this.time=a;this.visibility=b;this.elongation=c;this.ecliptic_separation=d};d.Elongation=function(a,b){b=d.MakeTime(b);var c=d.LongitudeFromSun(a,b);if(180=++g;){var k=d.EclipticLongitude(a,b),h=d.EclipticLongitude("Earth",b);h=H(k-h);var l=k=void 0,m=l=k=void 0;h>=-e.s1&&h<+e.s1?(m=0,k=+e.s1,l=+e.s2):h>=+e.s2||h<-e.s2?(m=0,k=-e.s2,l=-e.s1):0<=h?(m=-E(a)/4,k=+e.s1,l=+e.s2):(m=-E(a)/4,k=-e.s2,l=-e.s1);h=b.AddDays(m);k=d.SearchRelativeLongitude(a,k,h);l=d.SearchRelativeLongitude(a,l,k);h=c(k);if(0<=h)throw"SearchMaxElongation: internal error: m1 = "+ +h;m=c(l);if(0>=m)throw"SearchMaxElongation: internal error: m2 = "+m;h=d.Search(c,k,l,{init_f1:h,init_f2:m,dt_tolerance_seconds:10});if(!h)throw"SearchMaxElongation: failed search iter "+g+" (t1="+k.toString()+", t2="+l.toString()+")";if(h.tt>=b.tt)return d.Elongation(a,h);b=l.AddDays(1)}throw"SearchMaxElongation: failed to find event after 2 tries.";};d.SearchPeakMagnitude=function(a,b){function c(b){var c=b.AddDays(-.005);b=b.AddDays(.005);c=d.Illumination(a,c).mag;return(d.Illumination(a,b).mag- +c)/.01}if("Venus"!==a)throw"SearchPeakMagnitude currently works for Venus only.";b=d.MakeTime(b);for(var e=0;2>=++e;){var g=d.EclipticLongitude(a,b),k=d.EclipticLongitude("Earth",b);k=H(g-k);var h=g=void 0,l=h=g=void 0;-10<=k&&10>k?(l=0,g=10,h=30):30<=k||-30>k?(l=0,g=-30,h=-10):0<=k?(l=-E(a)/4,g=10,h=30):(l=-E(a)/4,g=-30,h=-10);k=b.AddDays(l);g=d.SearchRelativeLongitude(a,g,k);h=d.SearchRelativeLongitude(a,h,g);k=c(g);if(0<=k)throw"SearchPeakMagnitude: internal error: m1 = "+k;l=c(h);if(0>=l)throw"SearchPeakMagnitude: internal error: m2 = "+ +l;k=d.Search(c,g,h,{init_f1:k,init_f2:l,dt_tolerance_seconds:10});if(!k)throw"SearchPeakMagnitude: failed search iter "+e+" (t1="+g.toString()+", t2="+h.toString()+")";if(k.tt>=b.tt)return d.Illumination(a,k);b=h.AddDays(1)}throw"SearchPeakMagnitude: failed to find event after 2 tries.";};var U=function(a,b,c){this.time=a;this.kind=b;this.dist_au=c;this.dist_km=1.4959787069098932E8*c};d.SearchLunarApsis=function(a){function b(a){var b=a.AddDays(-5E-4);a=a.AddDays(5E-4);b=G(b).distance_au;return(G(a).distance_au- +b)/.001}function c(a){return-b(a)}a=d.MakeTime(a);var e=b(a);++J.lunar_apsis_calls;for(var g=0;59.061176>5*g;++g){++J.lunar_apsis_iter;var k=a.AddDays(5),h=b(k);if(0>=e*h){if(0>e||0h){a=d.Search(c,a,k,{init_f1:-e,init_f2:-h});if(null==a)throw"SearchLunarApsis INTERNAL ERROR: apogee search failed!";e=G(a).distance_au;return new U(a,1, +e)}throw"SearchLunarApsis INTERNAL ERROR: cannot classify apsis event!";}a=k;e=h}throw"SearchLunarApsis INTERNAL ERROR: could not find apsis within 2 synodic months of start date.";};d.NextLunarApsis=function(a){var b=d.SearchLunarApsis(a.time.AddDays(11));if(1!==b.kind+a.kind)throw"NextLunarApsis INTERNAL ERROR: did not find alternating apogee/perigee: prev="+a.kind+" @ "+a.time.toString()+", next="+b.kind+" @ "+b.time.toString();return b};d.SearchPlanetApsis=function(a,b){function c(b){var c=b.AddDays(-5E-4); +b=b.AddDays(5E-4);c=d.HelioDistance(a,c);return(d.HelioDistance(a,b)-c)/.001}function e(a){return-c(a)}if("Neptune"===a)return pa(b);for(var g=C[a].OrbitalPeriod,k=g/6,h=c(b),l=0;l*k<2*g;++l){var m=b.AddDays(k),f=c(m);if(0>=h*f){g=k=void 0;if(0>h||0f)k=e,g=1;else throw"Internal error with slopes in SearchPlanetApsis";b=d.Search(k,b,m,1);if(null==b)throw"Failed to find slope transition in planetary apsis search.";h=d.HelioDistance(a,b);return new U(b,g,h)}b=m;h=f}throw"Internal error: should have found planetary apsis within 2 orbital periods."; +};d.NextPlanetApsis=function(a,b){if(0!==b.kind&&1!==b.kind)throw"Invalid apsis kind: "+b.kind;var c=b.time.AddDays(.25*C[a].OrbitalPeriod);a=d.SearchPlanetApsis(a,c);if(1!==a.kind+b.kind)throw"Internal error: previous apsis was "+b.kind+", but found "+a.kind+" for next apsis.";return a};d.InverseRotation=function(a){return new A([[a.rot[0][0],a.rot[1][0],a.rot[2][0]],[a.rot[0][1],a.rot[1][1],a.rot[2][1]],[a.rot[0][2],a.rot[1][2],a.rot[2][2]]])};d.CombineRotation=function(a,b){return new A([[b.rot[0][0]* +a.rot[0][0]+b.rot[1][0]*a.rot[0][1]+b.rot[2][0]*a.rot[0][2],b.rot[0][1]*a.rot[0][0]+b.rot[1][1]*a.rot[0][1]+b.rot[2][1]*a.rot[0][2],b.rot[0][2]*a.rot[0][0]+b.rot[1][2]*a.rot[0][1]+b.rot[2][2]*a.rot[0][2]],[b.rot[0][0]*a.rot[1][0]+b.rot[1][0]*a.rot[1][1]+b.rot[2][0]*a.rot[1][2],b.rot[0][1]*a.rot[1][0]+b.rot[1][1]*a.rot[1][1]+b.rot[2][1]*a.rot[1][2],b.rot[0][2]*a.rot[1][0]+b.rot[1][2]*a.rot[1][1]+b.rot[2][2]*a.rot[1][2]],[b.rot[0][0]*a.rot[2][0]+b.rot[1][0]*a.rot[2][1]+b.rot[2][0]*a.rot[2][2],b.rot[0][1]* +a.rot[2][0]+b.rot[1][1]*a.rot[2][1]+b.rot[2][1]*a.rot[2][2],b.rot[0][2]*a.rot[2][0]+b.rot[1][2]*a.rot[2][1]+b.rot[2][2]*a.rot[2][2]]])};d.VectorFromSphere=function(a,b){var c=.017453292519943295*a.lat,d=.017453292519943295*a.lon,g=a.dist*Math.cos(c);return new t(g*Math.cos(d),g*Math.sin(d),a.dist*Math.sin(c),b)};d.VectorFromEquator=function(a,b){return d.VectorFromSphere(new X(a.dec,15*a.ra,a.dist),b)};d.EquatorFromVector=function(a){a=d.SphereFromVector(a);return new ea(a.lon/15,a.lat,a.dist)};d.SphereFromVector= +function(a){var b=a.x*a.x+a.y*a.y,c=Math.sqrt(b+a.z*a.z);if(0===b){if(0===a.z)throw"Zero-length vector not allowed.";var d=0;a=0>a.z?-90:90}else d=57.29577951308232*Math.atan2(a.y,a.x),0>d&&(d+=360),a=57.29577951308232*Math.atan2(a.z,Math.sqrt(b));return new X(a,d,c)};d.HorizonFromVector=function(a,b){a=d.SphereFromVector(a);a.lon=ia(a.lon);a.lat+=d.Refraction(b,a.lat);return a};d.VectorFromHorizon=function(a,b,c){var e=ia(a.lon);c=a.lat+d.InverseRefraction(c,a.lat);a=new X(c,e,a.dist);return d.VectorFromSphere(a, +b)};d.Refraction=function(a,b){if(-90>b||90c&&(c=-1);c=1.02/Math.tan(.017453292519943295*(c+10.3/(c+5.11)))/60;"normal"===a&&-1>b&&(c*=(b+90)/89)}else c=0;return c};d.InverseRefraction=function(a,b){if(-90>b||90Math.abs(e))return c-b;c-=e}};d.RotateVector=function(a,b){return new t(a.rot[0][0]*b.x+a.rot[1][0]*b.y+a.rot[2][0]*b.z,a.rot[0][1]*b.x+a.rot[1][1]*b.y+a.rot[2][1]* +b.z,a.rot[0][2]*b.x+a.rot[1][2]*b.y+a.rot[2][2]*b.z,b.t)};d.Rotation_EQJ_ECL=function(){return new A([[1,0,0],[0,.9174821430670688,-.3977769691083922],[0,.3977769691083922,.9174821430670688]])};d.Rotation_ECL_EQJ=function(){return new A([[1,0,0],[0,.9174821430670688,.3977769691083922],[0,-.3977769691083922,.9174821430670688]])};d.Rotation_EQJ_EQD=function(a){var b=Y(0,a.tt);a=aa(a,0);return d.CombineRotation(b,a)};d.Rotation_EQD_EQJ=function(a){var b=aa(a,1);a=Y(a.tt,0);return d.CombineRotation(b, +a)};d.Rotation_EQD_HOR=function(a,b){var c=Math.sin(.017453292519943295*b.latitude),d=Math.cos(.017453292519943295*b.latitude),g=Math.sin(.017453292519943295*b.longitude),k=Math.cos(.017453292519943295*b.longitude);b=[d*k,d*g,c];c=[-c*k,-c*g,d];g=[g,-k,0];a=-15*S(a);b=D(a,b);c=D(a,c);a=D(a,g);return new A([[c[0],a[0],b[0]],[c[1],a[1],b[1]],[c[2],a[2],b[2]]])};d.Rotation_HOR_EQD=function(a,b){a=d.Rotation_EQD_HOR(a,b);return d.InverseRotation(a)};d.Rotation_HOR_EQJ=function(a,b){b=d.Rotation_HOR_EQD(a, +b);a=d.Rotation_EQD_EQJ(a);return d.CombineRotation(b,a)};d.Rotation_EQJ_HOR=function(a,b){a=d.Rotation_HOR_EQJ(a,b);return d.InverseRotation(a)};d.Rotation_EQD_ECL=function(a){a=d.Rotation_EQD_EQJ(a);var b=d.Rotation_EQJ_ECL();return d.CombineRotation(a,b)};d.Rotation_ECL_EQD=function(a){a=d.Rotation_EQD_ECL(a);return d.InverseRotation(a)};d.Rotation_ECL_HOR=function(a,b){var c=d.Rotation_ECL_EQD(a);a=d.Rotation_EQD_HOR(a,b);return d.CombineRotation(c,a)};d.Rotation_HOR_ECL=function(a,b){a=d.Rotation_ECL_HOR(a, +b);return d.InverseRotation(a)}})("undefined"===typeof exports?this.Astronomy={}:exports); diff --git a/source/python/README.md b/source/python/README.md index fe4f19ea..44d58275 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -504,6 +504,8 @@ two cases applies to a particular apsis event. | `Pluto` | The planet Pluto. | | `Sun` | The Sun. | | `Moon` | The Earth's moon. | +| `EMB` | The Earth/Moon Barycenter. | +| `SSB` | The Solar System Barycenter. | --- @@ -893,7 +895,7 @@ the year range 1700..2200, this function raise an exception. | Type | Parameter | Description | | --- | --- | --- | -| [`Body`](#Body) | `body` | The celestial body whose heliocentric position is to be calculated: The Sun, Moon, or any of the planets. | +| [`Body`](#Body) | `body` | The celestial body whose heliocentric position is to be calculated: The Sun, Moon, EMB, SSB, or any of the planets. | | [`Time`](#Time) | `time` | The time at which to calculate the heliocentric position. | ### Returns: [`Vector`](#Vector) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 232e58ee..70d96794 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -59,6 +59,13 @@ _SUN_RADIUS_AU = 4.6505e-3 _MOON_RADIUS_AU = 1.15717e-5 _ASEC180 = 180.0 * 60.0 * 60.0 _AU_PER_PARSEC = _ASEC180 / math.pi +_EARTH_MOON_MASS_RATIO = 81.30056 +_SUN_MASS = 333054.25318 # Sun's mass relative to Earth. +_JUPITER_MASS = 317.84997 # Jupiter's mass relative to Earth. +_SATURN_MASS = 95.16745 # Saturn's mass relative to Earth. +_URANUS_MASS = 14.53617 # Uranus's mass relative to Earth. +_NEPTUNE_MASS = 17.14886 # Neptune's mass relative to Earth. + def _LongitudeOffset(diff): offset = diff @@ -121,6 +128,8 @@ class Body(enum.IntEnum): Pluto: The planet Pluto. Sun: The Sun. Moon: The Earth's moon. + EMB: The Earth/Moon Barycenter. + SSB: The Solar System Barycenter. """ Invalid = -1 Mercury = 0 @@ -134,6 +143,8 @@ class Body(enum.IntEnum): Pluto = 8 Sun = 9 Moon = 10 + EMB = 11 + SSB = 12 def BodyCode(name): """Finds the Body enumeration value, given the name of a body. @@ -3336,6 +3347,21 @@ def Search(func, context, t1, t2, dt_tolerance_seconds): # END Search #---------------------------------------------------------------------------- +def _AdjustBarycenter(ssb, time, body, pmass): + shift = pmass / (pmass + _SUN_MASS) + planet = _CalcVsop(_vsop[body], time) + ssb.x += shift * planet.x + ssb.y += shift * planet.y + ssb.z += shift * planet.z + +def _CalcSolarSystemBarycenter(time): + ssb = Vector(0.0, 0.0, 0.0, time) + _AdjustBarycenter(ssb, time, Body.Jupiter, _JUPITER_MASS) + _AdjustBarycenter(ssb, time, Body.Saturn, _SATURN_MASS) + _AdjustBarycenter(ssb, time, Body.Uranus, _URANUS_MASS) + _AdjustBarycenter(ssb, time, Body.Neptune, _NEPTUNE_MASS) + return ssb + def HelioVector(body, time): """Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. @@ -3354,7 +3380,7 @@ def HelioVector(body, time): ---------- body : Body The celestial body whose heliocentric position is to be calculated: - The Sun, Moon, or any of the planets. + The Sun, Moon, EMB, SSB, or any of the planets. time : Time The time at which to calculate the heliocentric position. @@ -3378,6 +3404,15 @@ def HelioVector(body, time): m = GeoMoon(time) return Vector(e.x+m.x, e.y+m.y, e.z+m.z, time) + if body == Body.EMB: + e = _CalcEarth(time) + m = GeoMoon(time) + d = 1.0 + _EARTH_MOON_MASS_RATIO + return Vector(e.x+(m.x/d), e.y+(m.y/d), e.z+(m.z/d), time) + + if body == Body.SSB: + return _CalcSolarSystemBarycenter(time) + raise InvalidBodyError()