From f9ef46c5cc1dba9d2b0b1f7042e52f90a2c6e6f2 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sun, 14 Jun 2020 14:55:52 -0400 Subject: [PATCH] Implemented Python version of transit search functions. --- generate/template/astronomy.js | 10 -- generate/template/astronomy.py | 162 +++++++++++++++++++++++++++++++++ generate/test.py | 59 ++++++++++++ source/js/README.md | 7 -- source/js/astronomy.js | 10 -- source/js/astronomy.min.js | 62 ++++++------- source/python/README.md | 63 +++++++++++++ source/python/astronomy.py | 162 +++++++++++++++++++++++++++++++++ 8 files changed, 477 insertions(+), 58 deletions(-) diff --git a/generate/template/astronomy.js b/generate/template/astronomy.js index 25bf0c5a..e2bd383d 100644 --- a/generate/template/astronomy.js +++ b/generate/template/astronomy.js @@ -4995,13 +4995,6 @@ Astronomy.NextLocalSolarEclipse = function(prevEclipseTime, observer) { * A transit is when Mercury or Venus passes between the Sun and Earth so that * the other planet is seen in silhouette against the Sun. * - * The `start` field reports the moment in time when the planet first becomes - * visible against the Sun in its background. - * The `peak` field reports when the planet is most aligned with the Sun, - * as seen from the Earth. - * The `finish` field reports the last moment when the planet is visible - * against the Sun in its background. - * * The calculations are performed from the point of view of a geocentric observer. * * @class @@ -5040,7 +5033,6 @@ function PlanetShadowBoundary(time, body, planet_radius_km, direction) { function PlanetTransitBoundary(body, planet_radius_km, t1, t2, direction) { // Search for the time the planet's penumbra begins/ends making contact with the center of the Earth. - // context = new SearchContext_PlanetShadowBoundary(body, planet_radius_km, direction); const tx = Astronomy.Search(time => PlanetShadowBoundary(time, body, planet_radius_km, direction), t1, t2, 1.0); if (tx == null) throw 'Planet transit boundary search failed'; @@ -5104,8 +5096,6 @@ Astronomy.SearchTransit = function(body, startTime) { const shadow = PeakPlanetShadow(body, planet_radius_km, conj); if (shadow.r < shadow.p) { // does the planet's penumbra touch the Earth's center? - var transit = new TransitInfo(); - // Find the beginning and end of the penumbral contact. const time_before = shadow.time.AddDays(-dt_days); const start = PlanetTransitBoundary(body, planet_radius_km, time_before, shadow.time, -1.0); diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index e5eec68f..e00d91cc 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -4183,6 +4183,20 @@ def _LocalMoonShadow(time, observer): return _CalcShadow(_MOON_MEAN_RADIUS_KM, time, o, m) +def _PlanetShadow(body, planet_radius_km, time): + # Calculate light-travel-corrected vector from Earth to planet. + g = GeoVector(body, time, False) + # Calculate light-travel-corrected vector from Earth to Sun. + e = GeoVector(Body.Sun, time, False) + # Deduce light-travel-corrected vector from Sun to planet. + p = Vector(g.x - e.x, g.y - e.y, g.z - e.z, time) + # Calcluate Earth's position from the planet's point of view. + e.x = -g.x + e.y = -g.y + e.z = -g.z + return _CalcShadow(planet_radius_km, time, e, p) + + def _ShadowDistanceSlope(shadowfunc, time): dt = 1.0 / 86400.0 t1 = time.AddDays(-dt) @@ -4192,6 +4206,14 @@ def _ShadowDistanceSlope(shadowfunc, time): return (shadow2.r - shadow1.r) / dt +def _PlanetShadowSlope(context, time): + (body, planet_radius_km) = context + dt = 1.0 / 86400.0 + shadow1 = _PlanetShadow(body, planet_radius_km, time.AddDays(-dt)) + shadow2 = _PlanetShadow(body, planet_radius_km, time.AddDays(+dt)) + return (shadow2.r - shadow1.r) / dt + + def _PeakEarthShadow(search_center_time): window = 0.03 # initial search window, in days, before/after given time t1 = search_center_time.AddDays(-window) @@ -4216,6 +4238,14 @@ def _PeakLocalMoonShadow(search_center_time, observer): tx = Search(_ShadowDistanceSlope, lambda time: _LocalMoonShadow(time, observer), t1, t2, 1.0) return _LocalMoonShadow(tx, observer) +def _PeakPlanetShadow(body, planet_radius_km, search_center_time): + # Search for when the body's shadow is closest to the center of the Earth. + window = 1.0 # days before/after inferior conjunction to search for minimum shadow distance. + t1 = search_center_time.AddDays(-window) + t2 = search_center_time.AddDays(+window) + tx = Search(_PlanetShadowSlope, (body, planet_radius_km), t1, t2, 1.0) + return _PlanetShadow(body, planet_radius_km, tx) + class _ShadowDiffContext: def __init__(self, radius_limit, direction): self.radius_limit = radius_limit @@ -4796,3 +4826,135 @@ def NextLocalSolarEclipse(prevEclipseTime, observer): """ startTime = prevEclipseTime.AddDays(10.0) return SearchLocalSolarEclipse(startTime, observer) + + +class TransitInfo: + """Information about a transit of Mercury or Venus, as seen from the Earth. + + Returned by #SearchTransit or #NextTransit to report + information about a transit of Mercury or Venus. + A transit is when Mercury or Venus passes between the Sun and Earth so that + the other planet is seen in silhouette against the Sun. + + The calculations are performed from the point of view of a geocentric observer. + + Attributes + ---------- + start : Time + The date and time at the beginning of the transit. + This is the moment the planet first becomes visible against the Sun in its background. + + peak : Time + When the planet is most aligned with the Sun, as seen from the Earth. + + finish : Time + The date and time at the end of the transit. + This is the moment the planet is last seen against the Sun in its background. + + separation : float + The minimum angular separation, in arcminutes, between the centers of the Sun and the planet. + This angle pertains to the time stored in `peak`. + """ + def __init__(self, start, peak, finish, separation): + self.start = start + self.peak = peak + self.finish = finish + self.separation = separation + + +def _PlanetShadowBoundary(context, time): + (body, planet_radius_km, direction) = context + shadow = _PlanetShadow(body, planet_radius_km, time) + return direction * (shadow.r - shadow.p) + + +def _PlanetTransitBoundary(body, planet_radius_km, t1, t2, direction): + # Search for the time the planet's penumbra begins/ends making contact with the center of the Earth. + # context = new SearchContext_PlanetShadowBoundary(body, planet_radius_km, direction); + tx = Search(_PlanetShadowBoundary, (body, planet_radius_km, direction), t1, t2, 1.0) + if tx is None: + raise Error('Planet transit boundary search failed') + return tx + + +def SearchTransit(body, startTime): + """Searches for the first transit of Mercury or Venus after a given date. + + Finds the first transit of Mercury or Venus after a specified date. + A transit is when an inferior planet passes between the Sun and the Earth + so that the silhouette of the planet is visible against the Sun in the background. + To continue the search, pass the `finish` time in the returned structure to + #NextTransit. + + Parameters + ---------- + body : Body + The planet whose transit is to be found. Must be `Body.Mercury` or `Body.Venus`. + startTime : Time + The date and time for starting the search for a transit. + + Returns + ------- + TransitInfo + """ + threshold_angle = 0.4 # maximum angular separation to attempt transit calculation + dt_days = 1.0 + + # Validate the planet and find its mean radius. + if body == Body.Mercury: + planet_radius_km = 2439.7 + elif body == Body.Venus: + planet_radius_km = 6051.8 + else: + raise InvalidBodyError() + + search_time = startTime + while True: + # Search for the next inferior conjunction of the given planet. + # This is the next time the Earth and the other planet have the same + # ecliptic longitude as seen from the Sun. + conj = SearchRelativeLongitude(body, 0.0, search_time) + + # Calculate the angular separation between the body and the Sun at this time. + conj_separation = AngleFromSun(body, conj) + + if conj_separation < threshold_angle: + # The planet's angular separation from the Sun is small enough + # to consider it a transit candidate. + # Search for the moment when the line passing through the Sun + # and planet are closest to the Earth's center. + shadow = _PeakPlanetShadow(body, planet_radius_km, conj) + if shadow.r < shadow.p: # does the planet's penumbra touch the Earth's center? + # Find the beginning and end of the penumbral contact. + time_before = shadow.time.AddDays(-dt_days) + start = _PlanetTransitBoundary(body, planet_radius_km, time_before, shadow.time, -1.0) + time_after = shadow.time.AddDays(+dt_days) + finish = _PlanetTransitBoundary(body, planet_radius_km, shadow.time, time_after, +1.0) + min_separation = 60.0 * AngleFromSun(body, shadow.time) + return TransitInfo(start, shadow.time, finish, min_separation) + + # This inferior conjunction was not a transit. Try the next inferior conjunction. + search_time = conj.AddDays(10.0) + + + +def NextTransit(body, prevTransitTime): + """Searches for another transit of Mercury or Venus. + + After calling #SearchTransit to find a transit of Mercury or Venus, + this function finds the next transit after that. + Keep calling this function as many times as you want to keep finding more transits. + + Parameters + ---------- + body : Body + The planet whose transit is to be found. Must be `Body.Mercury` or `Body.Venus`. + prevTransitTime : Time + A date and time near the previous transit. + + Returns + ------- + TransitInfo + """ + startTime = prevTransitTime.AddDays(100.0) + return SearchTransit(body, startTime) \ No newline at end of file diff --git a/generate/test.py b/generate/test.py index 909a7c4a..d25bd43e 100755 --- a/generate/test.py +++ b/generate/test.py @@ -1479,6 +1479,64 @@ def LocalSolarEclipse(): #----------------------------------------------------------------------------------------------------------- +def TransitFile(body, filename, limit_minutes, limit_sep): + lnum = 0 + max_minutes = 0 + max_sep = 0 + with open(filename, 'rt') as infile: + transit = astronomy.SearchTransit(body, astronomy.Time.Make(1600, 1, 1, 0, 0, 0)) + for line in infile: + lnum += 1 + token = line.strip().split() + # 22:17 1881-11-08T00:57Z 03:38 3.8633 + if len(token) != 4: + print('PY TransitFile({} line {}): bad data format.'.format(filename, lnum)) + return 1 + + textp = token[1] + text1 = textp[0:11] + token[0] + 'Z' + text2 = textp[0:11] + token[2] + 'Z' + timep = astronomy.Time.Parse(textp) + time1 = astronomy.Time.Parse(text1) + time2 = astronomy.Time.Parse(text2) + separation = float(token[3]) + + # If the start time is after the peak time, it really starts on the previous day. + if time1.ut > timep.ut: + time1 = time1.AddDays(-1.0) + + # If the finish time is before the peak time, it really starts on the following day. + if time2.ut < timep.ut: + time2 = time2.AddDays(+1.0) + + diff_start = (24.0 * 60.0) * vabs(time1.ut - transit.start.ut ) + diff_peak = (24.0 * 60.0) * vabs(timep.ut - transit.peak.ut ) + diff_finish = (24.0 * 60.0) * vabs(time2.ut - transit.finish.ut) + diff_sep = vabs(separation - transit.separation) + max_minutes = vmax(max_minutes, diff_start) + max_minutes = vmax(max_minutes, diff_peak) + max_minutes = vmax(max_minutes, diff_finish) + if max_minutes > limit_minutes: + print('PY TransitFile({} line {}): EXCESSIVE TIME ERROR = {} minutes.'.format(filename, lnum, max_minutes)) + return 1 + max_sep = vmax(max_sep, diff_sep) + if max_sep > limit_sep: + print('PY TransitFile({} line {}): EXCESSIVE SEPARATION ERROR = {} arcminutes.'.format(filename, lnum, max_sep)) + return 1 + transit = astronomy.NextTransit(body, transit.finish) + print('PY TransitFile({}): PASS - verified {}, max minutes = {}, max sep arcmin = {}'.format(filename, lnum, max_minutes, max_sep)) + return 0 + + +def Transit(): + if 0 != TransitFile(astronomy.Body.Mercury, 'eclipse/mercury.txt', 10.710, 0.2121): + return 1 + if 0 != TransitFile(astronomy.Body.Venus, 'eclipse/venus.txt', 9.109, 0.6772): + return 1 + return 0 + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { 'constellation': Constellation, 'elongation': Elongation, @@ -1495,6 +1553,7 @@ UnitTests = { 'rotation': Rotation, 'seasons': Seasons, 'time': AstroTime, + 'transit': Transit, } #----------------------------------------------------------------------------------------------------------- diff --git a/source/js/README.md b/source/js/README.md index b4ce9ea4..094283f0 100644 --- a/source/js/README.md +++ b/source/js/README.md @@ -640,13 +640,6 @@ information about a transit of Mercury or Venus. A transit is when Mercury or Venus passes between the Sun and Earth so that the other planet is seen in silhouette against the Sun. -The `start` field reports the moment in time when the planet first becomes -visible against the Sun in its background. -The `peak` field reports when the planet is most aligned with the Sun, -as seen from the Earth. -The `finish` field reports the last moment when the planet is visible -against the Sun in its background. - The calculations are performed from the point of view of a geocentric observer. **Properties** diff --git a/source/js/astronomy.js b/source/js/astronomy.js index 316d6ab2..b8ff99ba 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -6425,13 +6425,6 @@ Astronomy.NextLocalSolarEclipse = function(prevEclipseTime, observer) { * A transit is when Mercury or Venus passes between the Sun and Earth so that * the other planet is seen in silhouette against the Sun. * - * The `start` field reports the moment in time when the planet first becomes - * visible against the Sun in its background. - * The `peak` field reports when the planet is most aligned with the Sun, - * as seen from the Earth. - * The `finish` field reports the last moment when the planet is visible - * against the Sun in its background. - * * The calculations are performed from the point of view of a geocentric observer. * * @class @@ -6470,7 +6463,6 @@ function PlanetShadowBoundary(time, body, planet_radius_km, direction) { function PlanetTransitBoundary(body, planet_radius_km, t1, t2, direction) { // Search for the time the planet's penumbra begins/ends making contact with the center of the Earth. - // context = new SearchContext_PlanetShadowBoundary(body, planet_radius_km, direction); const tx = Astronomy.Search(time => PlanetShadowBoundary(time, body, planet_radius_km, direction), t1, t2, 1.0); if (tx == null) throw 'Planet transit boundary search failed'; @@ -6534,8 +6526,6 @@ Astronomy.SearchTransit = function(body, startTime) { const shadow = PeakPlanetShadow(body, planet_radius_km, conj); if (shadow.r < shadow.p) { // does the planet's penumbra touch the Earth's center? - var transit = new TransitInfo(); - // Find the beginning and end of the penumbral contact. const time_before = shadow.time.AddDays(-dt_days); const start = PlanetTransitBoundary(body, planet_radius_km, time_before, shadow.time, -1.0); diff --git a/source/js/astronomy.min.js b/source/js/astronomy.min.js index e0153058..55e4c991 100644 --- a/source/js/astronomy.min.js +++ b/source/js/astronomy.min.js @@ -49,17 +49,17 @@ 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]]}functi c=Math.cos(c);var l=Math.sin(-e);e=Math.cos(-e);var g=Math.sin(-d);var m=Math.cos(-d);var f=Math.sin(k);var n=Math.cos(k);k=n*e-l*f*m;d=n*l*c+f*m*e*c-a*f*g;var p=n*l*a+f*m*e*a+c*f*g;var q=-f*e-l*n*m;var r=-f*l*c+n*m*e*c-a*n*g;f=-f*l*a+n*m*e*a+c*n*g;l*=g;n=-g*e*c-a*m;a=-g*e*a+m*c;return 0===b?new w([[k,d,p],[q,r,f],[l,n,a]]):new w([[k,q,l],[d,r,n],[p,f,a]])}function L(a){var b=a.tt/36525,c=15*N(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=N(a);var c=.017453292519943295*a.mobl,d=.017453292519943295*a.tobl,e=4.84813681109536E-6*a.dpsi;a=Math.cos(c);c=Math.sin(c);var k=Math.cos(d),l=Math.sin(d);d=Math.cos(e);var g=Math.sin(e);e=-g*a;var m=-g* c,f=g*k,n=d*a*k+c*l,p=d*c*k-a*l;g*=l;var q=d*a*l-c*k;a=d*c*l+a*k;return 0===b?new w([[d,f,g],[e,n,q],[m,p,a]]):new w([[d,e,m],[f,n,p],[g,q,a]])}function ka(a,b){var c=L(a),d=.017453292519943295*b.latitude,e=Math.sin(d);d=Math.cos(d);var k=1/Math.sqrt(d*d+.9933056020041345*e*e),l=b.height/1E3,g=6378.1366*k+l;b=.017453292519943295*(15*c+b.longitude);e=Z(a,-1,[g*d*Math.cos(b)/1.4959787069098932E8,g*d*Math.sin(b)/1.4959787069098932E8,(6335.438815127603*k+l)*e/1.4959787069098932E8]);return P(a.tt,e,0)} -function za(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 la(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 d=Math.atan2(a[1],a[0])/.2617993877991494;0>d&&(d+=24);return new ba(d,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 ma(a,b,c,d,e){var h=b*d+c*e;b=-b*e+c*d;c=Math.sqrt(a*a+h*h);d=0;0d&&(d+=360));return new Aa(a,h,b,57.29577951308232*Math.atan2(b,c),d)}function G(a,b){var c=1,d=0;a=$jscomp.makeIterator(a);for(var e=a.next();!e.done;e=a.next()){var k=0;e=$jscomp.makeIterator(e.value);for(var l=e.next();!l.done;l=e.next())l= -l.value,k+=l[0]*Math.cos(l[1]+b*l[2]);d+=c*k;c*=b}return d}function A(a,b){var c=b.tt/365250,d=[];a=$jscomp.makeIterator(a);for(var e=a.next();!e.done;e=a.next())d.push(G(e.value,c));c=d[2]*Math.cos(d[1]);d=[c*Math.cos(d[0]),c*Math.sin(d[0]),d[2]*Math.sin(d[1])];return new r(d[0]+4.4036E-7*d[1]-1.90919E-7*d[2],-4.79966E-7*d[0]+.917482137087*d[1]-.397776982902*d[2],.397776982902*d[1]+.917482137087*d[2],b)}function Q(a,b,c,d){d/=d+333054.25318;b=A(t[c],b);a.x+=d*b.x;a.y+=d*b.y;a.z+=d*b.z}function Ba(a, +function ya(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 la(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 d=Math.atan2(a[1],a[0])/.2617993877991494;0>d&&(d+=24);return new ba(d,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 ma(a,b,c,d,e){var h=b*d+c*e;b=-b*e+c*d;c=Math.sqrt(a*a+h*h);d=0;0d&&(d+=360));return new za(a,h,b,57.29577951308232*Math.atan2(b,c),d)}function G(a,b){var c=1,d=0;a=$jscomp.makeIterator(a);for(var e=a.next();!e.done;e=a.next()){var k=0;e=$jscomp.makeIterator(e.value);for(var l=e.next();!l.done;l=e.next())l= +l.value,k+=l[0]*Math.cos(l[1]+b*l[2]);d+=c*k;c*=b}return d}function A(a,b){var c=b.tt/365250,d=[];a=$jscomp.makeIterator(a);for(var e=a.next();!e.done;e=a.next())d.push(G(e.value,c));c=d[2]*Math.cos(d[1]);d=[c*Math.cos(d[0]),c*Math.sin(d[0]),d[2]*Math.sin(d[1])];return new r(d[0]+4.4036E-7*d[1]-1.90919E-7*d[2],-4.79966E-7*d[0]+.917482137087*d[1]-.397776982902*d[2],.397776982902*d[1]+.917482137087*d[2],b)}function Q(a,b,c,d){d/=d+333054.25318;b=A(t[c],b);a.x+=d*b.x;a.y+=d*b.y;a.z+=d*b.z}function Aa(a, b,c,d,e){var h=(e+c)/2-d;c=(e-c)/2;if(0==h){if(0==c)return null;d=-d/c;if(-1>d||1=d)return null;e=Math.sqrt(d);d=(-c+e)/(2*h);e=(-c-e)/(2*h);if(-1<=d&&1>=d){if(-1<=e&&1>=e)return null}else if(-1<=e&&1>=e)d=e;else return null}return{x:d,t:a+d*b,df_dt:(2*h*d+c)/b}}function J(a){for(;-180>=a;)a+=360;for(;180l;++l){var g=b.AddDays(l*c);g=d*G(t.Neptune[2],g.tt/365250);if(0==l||g>k)e=l,k=g}b=b.AddDays((e-1)*c);c*=2}}function Ca(a){var b=a.AddDays(-30/360*x.Neptune.OrbitalPeriod),c=a.AddDays(.75*x.Neptune.OrbitalPeriod),d=b,e=b,k=-1,l=-1;c=(c.ut-b.ut)/ -99;for(var g=0;100>g;++g){var m=b.AddDays(g*c),f=G(t.Neptune[2],m.tt/365250);0===g?l=k=f:(f>l&&(l=f,e=m),f=a.tt)return e.time.tt>=a.tt&&e.time.tt=a.tt)return e;throw"Internal error: failed to find Neptune apsis.";}function oa(a){a=360-a;360<=a?a-=360:0>a&&(a+=360);return a}function E(a,b,c,d){var e=(d.x*c.x+d.y*c.y+d.z*c.z)/(d.x*d.x+d.y*d.y+d.z*d.z),k=e*d.x-c.x,l=e*d.y-c.y,g=e*d.z-c.z;return new Da(b, +if(!b)throw"Not a valid planet name: "+a;a=x.Earth.OrbitalPeriod;return Math.abs(a/(a/b.OrbitalPeriod-1))}function na(a,b,c){for(var d=1===a?1:-1;;){c/=9;if(c<1/1440)return b=b.AddDays(c/2),d=G(t.Neptune[2],b.tt/365250),new R(b,a,d);for(var e=-1,k=0,l=0;10>l;++l){var g=b.AddDays(l*c);g=d*G(t.Neptune[2],g.tt/365250);if(0==l||g>k)e=l,k=g}b=b.AddDays((e-1)*c);c*=2}}function Ba(a){var b=a.AddDays(-30/360*x.Neptune.OrbitalPeriod),c=a.AddDays(.75*x.Neptune.OrbitalPeriod),d=b,e=b,k=-1,l=-1;c=(c.ut-b.ut)/ +99;for(var g=0;100>g;++g){var m=b.AddDays(g*c),f=G(t.Neptune[2],m.tt/365250);0===g?l=k=f:(f>l&&(l=f,e=m),f=a.tt)return e.time.tt>=a.tt&&e.time.tt=a.tt)return e;throw"Internal error: failed to find Neptune apsis.";}function oa(a){a=360-a;360<=a?a-=360:0>a&&(a+=360);return a}function E(a,b,c,d){var e=(d.x*c.x+d.y*c.y+d.z*c.z)/(d.x*d.x+d.y*d.y+d.z*d.z),k=e*d.x-c.x,l=e*d.y-c.y,g=e*d.z-c.z;return new Ca(b, e,1.4959787069098932E8*Math.sqrt(k*k+l*l+g*g),695700-(1+e)*(695700-a),-695700+(1+e)*(695700+a),c,d)}function S(a){var b=A(t.Earth,a),c=e.GeoMoon(a);return E(6459,a,c,b)}function pa(a){var b=A(t.Earth,a),c=e.GeoMoon(a),d=new r(-c.x,-c.y,-c.z,c.t);c.x+=b.x;c.y+=b.y;c.z+=b.z;return E(1737.4,a,d,c)}function ca(a,b){var c=ka(a,b);b=A(t.Earth,a);var d=e.GeoMoon(a);c=new r(c[0]-d.x,c[1]-d.y,c[2]-d.z,a);d.x+=b.x;d.y+=b.y;d.z+=b.z;return E(1737.4,a,c,d)}function T(a,b,c){a=e.GeoVector(a,c,!1);var d=e.GeoVector("Sun", -c,!1),h=new r(a.x-d.x,a.y-d.y,a.z-d.z,c);d.x=-a.x;d.y=-a.y;d.z=-a.z;return E(b,c,d,h)}function da(a,b){var c=1/86400,d=b.AddDays(-c);b=b.AddDays(+c);d=a(d);return(a(b).r-d.r)/c}function Ea(a){var b=a.AddDays(-.03);a=a.AddDays(.03);b=e.Search(function(a){return da(S,a)},b,a,1);return S(b)}function Fa(a){var b=a.AddDays(-.03);a=a.AddDays(.03);b=e.Search(function(a){return da(pa,a)},b,a,1);return pa(b)}function Ga(a,b,c){var d=c.AddDays(-1);c=c.AddDays(1);d=e.Search(function(c){var d=1/86400,e=T(a,b, -c.AddDays(-d));return(T(a,b,c.AddDays(+d)).r-e.r)/d},d,c,1);return T(a,b,d)}function Ha(a,b){function c(a){return ca(a,b)}var d=a.AddDays(-.2),h=a.AddDays(.2);d=e.Search(function(a){return da(c,a)},d,h,1);if(!d)throw"PeakLocalMoonShadow: search failure for search_center_time = "+a;return ca(d,b)}function ea(a,b,c){var d=c/1440;c=a.AddDays(-d);d=a.AddDays(+d);c=e.Search(function(a){return-(S(a).r-b)},c,a,1);a=e.Search(function(a){return+(S(a).r-b)},a,d,1);if(null===c||null===a)throw"Failed to find shadow semiduration"; -return 720*(a.ut-c.ut)}function qa(a){return a.p-a.r}function ra(a){return Math.abs(a.k)-a.r}function U(a,b,c,d,h){d=e.Search(function(d){d=ca(d,a);return b*c(d)},d,h,1);if(null==d)throw"Local eclipse transition search failed.";return sa(a,d)}function sa(a,b){var c=e.Equator("Sun",b,a,!0,!0);a=e.Horizon(b,a,c.ra,c.dec,"normal").altitude;return new Ia(b,a)}function ta(a,b,c,d,h){c=e.Search(function(c){c=T(a,b,c);return h*(c.r-c.p)},c,d,1);if(null==c)throw"Planet transit boundary search failed";return c} -var fa=new Date("2000-01-01T12:00:00Z"),I=2*Math.PI,K=180/Math.PI*3600,Ja=-.17-5*Math.log10(648E3/Math.PI),Ka=34/60,V,ua,va;e.Bodies="Sun Moon Mercury Venus Earth Mars Jupiter Saturn Uranus Neptune Pluto SSB EMB".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}}, +c,!1),h=new r(a.x-d.x,a.y-d.y,a.z-d.z,c);d.x=-a.x;d.y=-a.y;d.z=-a.z;return E(b,c,d,h)}function da(a,b){var c=1/86400,d=b.AddDays(-c);b=b.AddDays(+c);d=a(d);return(a(b).r-d.r)/c}function Da(a){var b=a.AddDays(-.03);a=a.AddDays(.03);b=e.Search(function(a){return da(S,a)},b,a,1);return S(b)}function Ea(a){var b=a.AddDays(-.03);a=a.AddDays(.03);b=e.Search(function(a){return da(pa,a)},b,a,1);return pa(b)}function Fa(a,b,c){var d=c.AddDays(-1);c=c.AddDays(1);d=e.Search(function(c){var d=1/86400,e=T(a,b, +c.AddDays(-d));return(T(a,b,c.AddDays(+d)).r-e.r)/d},d,c,1);return T(a,b,d)}function Ga(a,b){function c(a){return ca(a,b)}var d=a.AddDays(-.2),h=a.AddDays(.2);d=e.Search(function(a){return da(c,a)},d,h,1);if(!d)throw"PeakLocalMoonShadow: search failure for search_center_time = "+a;return ca(d,b)}function ea(a,b,c){var d=c/1440;c=a.AddDays(-d);d=a.AddDays(+d);c=e.Search(function(a){return-(S(a).r-b)},c,a,1);a=e.Search(function(a){return+(S(a).r-b)},a,d,1);if(null===c||null===a)throw"Failed to find shadow semiduration"; +return 720*(a.ut-c.ut)}function qa(a){return a.p-a.r}function ra(a){return Math.abs(a.k)-a.r}function U(a,b,c,d,h){d=e.Search(function(d){d=ca(d,a);return b*c(d)},d,h,1);if(null==d)throw"Local eclipse transition search failed.";return sa(a,d)}function sa(a,b){var c=e.Equator("Sun",b,a,!0,!0);a=e.Horizon(b,a,c.ra,c.dec,"normal").altitude;return new Ha(b,a)}function ta(a,b,c,d,h){c=e.Search(function(c){c=T(a,b,c);return h*(c.r-c.p)},c,d,1);if(null==c)throw"Planet transit boundary search failed";return c} +var fa=new Date("2000-01-01T12:00:00Z"),I=2*Math.PI,K=180/Math.PI*3600,Ia=-.17-5*Math.log10(648E3/Math.PI),Ja=34/60,V,ua,va;e.Bodies="Sun Moon Mercury Venus Earth Mars Jupiter Saturn Uranus Neptune Pluto SSB EMB".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}}, t={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], @@ -124,33 +124,33 @@ cls:[-59641,-11,149,25543,-11,66]},{nals:[1,0,2,0,1],cls:[-51613,-42,129,26366,0 0,0,-2,1],cls:[-4940,-11,-21,2720,0,-9]},{nals:[-1,-1,0,2,0],cls:[7350,0,-8,-51,0,4]},{nals:[2,0,0,-2,1],cls:[4065,0,6,-2206,0,1]},{nals:[1,0,0,2,0],cls:[6579,0,-24,-199,0,2]},{nals:[0,1,2,-2,1],cls:[3579,0,5,-1900,0,1]},{nals:[1,-1,0,0,0],cls:[4725,0,-6,-41,0,3]},{nals:[-2,0,2,0,2],cls:[-3075,0,-2,1313,0,-1]},{nals:[3,0,2,0,2],cls:[-2904,0,15,1233,0,7]},{nals:[0,-1,0,2,0],cls:[4348,0,-10,-81,0,2]},{nals:[1,-1,2,0,2],cls:[-2878,0,8,1232,0,4]},{nals:[0,0,0,1,0],cls:[-4230,0,5,-20,0,-2]},{nals:[-1, -1,2,2,2],cls:[-2819,0,7,1207,0,3]},{nals:[-1,0,2,0,0],cls:[-4056,0,5,40,0,-2]},{nals:[0,-1,2,2,2],cls:[-2647,0,11,1129,0,5]},{nals:[-2,0,0,0,1],cls:[-2294,0,-10,1266,0,-4]},{nals:[1,1,2,0,2],cls:[2481,0,-7,-1062,0,-3]},{nals:[2,0,0,0,1],cls:[2179,0,-2,-1129,0,-2]},{nals:[-1,1,0,1,0],cls:[3276,0,1,-9,0,0]},{nals:[1,1,0,0,0],cls:[-3389,0,5,35,0,-2]},{nals:[1,0,2,0,0],cls:[3339,0,-13,-107,0,1]},{nals:[-1,0,2,-2,1],cls:[-1987,0,-6,1073,0,-2]},{nals:[1,0,0,0,2],cls:[-1981,0,0,854,0,0]},{nals:[-1,0,0, 1,0],cls:[4026,0,-353,-553,0,-139]},{nals:[0,0,2,1,2],cls:[1660,0,-5,-710,0,-2]},{nals:[-1,0,2,4,2],cls:[-1521,0,9,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]}],O;e.CalcMoonCount=0;var r=function(a,b,c,d){this.x=a;this.y=b;this.z=c;this.t=d};r.prototype.Length=function(){return Math.sqrt(this.x* -this.x+this.y*this.y+this.z*this.z)};var W=function(a,b,c){this.lat=a;this.lon=b;this.dist=c};e.MakeSpherical=function(a,b,c){return new W(a,b,c)};var ba=function(a,b,c){this.ra=a;this.dec=b;this.dist=c},w=function(a){this.rot=a};e.MakeRotation=function(a){if(!za(a))throw"Argument must be a [3][3] array of numbers";return new w(a)};var La=function(a,b,c,d){this.azimuth=a;this.altitude=b;this.ra=c;this.dec=d},Aa=function(a,b,c,d,e){this.ex=a;this.ey=b;this.ez=c;this.elat=d;this.elon=e};e.Horizon=function(a, +this.x+this.y*this.y+this.z*this.z)};var W=function(a,b,c){this.lat=a;this.lon=b;this.dist=c};e.MakeSpherical=function(a,b,c){return new W(a,b,c)};var ba=function(a,b,c){this.ra=a;this.dec=b;this.dist=c},w=function(a){this.rot=a};e.MakeRotation=function(a){if(!ya(a))throw"Argument must be a [3][3] array of numbers";return new w(a)};var Ka=function(a,b,c,d){this.azimuth=a;this.altitude=b;this.ra=c;this.dec=d},za=function(a,b,c,d,e){this.ex=a;this.ey=b;this.ez=c;this.elat=d;this.elon=e};e.Horizon=function(a, b,c,d,h){a=e.MakeTime(a);var k=Math.sin(.017453292519943295*b.latitude),l=Math.cos(.017453292519943295*b.latitude),g=Math.sin(.017453292519943295*b.longitude),m=Math.cos(.017453292519943295*b.longitude);b=Math.sin(.017453292519943295*d);var f=Math.cos(.017453292519943295*d),n=Math.sin(.2617993877991494*c),p=Math.cos(.2617993877991494*c),q=[l*m,l*g,k];k=[-k*m,-k*g,l];g=[g,-m,0];l=-15*L(a);a=D(l,q);q=D(l,k);g=D(l,g);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];g=b[0]* g[0]+b[1]*g[1]+b[2]*g[2];p=Math.sqrt(q*q+g*g);f=0;0f&&(f+=360),360<=f&&(f-=360));n=57.29577951308232*Math.atan2(p,n);p=d;if(h&&(d=n,h=e.Refraction(h,90-n),n-=h,0g;++g)h.push((b[g]-d*a[g])/q*c+a[g]*p);p=Math.sqrt(h[0]*h[0]+h[1]*h[1]);0c&&(c+=24), -24<=c&&(c-=24)):c=0;p=57.29577951308232*Math.atan2(h[2],p)}return new La(f,90-n,c,p)};var Ma=function(a,b,c){this.latitude=a;this.longitude=b;this.height=c};e.MakeObserver=function(a,b,c){return new Ma(a,b,c||0)};e.SunPosition=function(a){a=e.MakeTime(a).AddDays(-.005775518331089121);var b=A(t.Earth,a);b=P(0,[-b.x,-b.y,-b.z],a.tt);b=Z(a,0,b);a=.017453292519943295*N(a).tobl;return ma(b[0],b[1],b[2],Math.cos(a),Math.sin(a))};e.Equator=function(a,b,c,d,h){b=e.MakeTime(b);c=ka(b,c);a=e.GeoVector(a,b, +24<=c&&(c-=24)):c=0;p=57.29577951308232*Math.atan2(h[2],p)}return new Ka(f,90-n,c,p)};var La=function(a,b,c){this.latitude=a;this.longitude=b;this.height=c};e.MakeObserver=function(a,b,c){return new La(a,b,c||0)};e.SunPosition=function(a){a=e.MakeTime(a).AddDays(-.005775518331089121);var b=A(t.Earth,a);b=P(0,[-b.x,-b.y,-b.z],a.tt);b=Z(a,0,b);a=.017453292519943295*N(a).tobl;return ma(b[0],b[1],b[2],Math.cos(a),Math.sin(a))};e.Equator=function(a,b,c,d,h){b=e.MakeTime(b);c=ka(b,c);a=e.GeoVector(a,b, h);a=[a.x-c[0],a.y-c[1],a.z-c[2]];if(!d)return la(a);d=P(0,a,b.tt);d=Z(b,0,d);return la(d)};e.Ecliptic=function(a,b,c){void 0===V&&(V=.017453292519943295*N(e.MakeTime(fa)).mobl,ua=Math.cos(V),va=Math.sin(V));return ma(a,b,c,ua,va)};e.GeoMoon=function(a){a=e.MakeTime(a);var b=z(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 d=.017453292519943295*F(a);c=Math.cos(d);d=Math.sin(d);b=P(a.tt,[b[0],b[1]* c-b[2]*d,b[1]*d+b[2]*c],0);return new r(b[0],b[1],b[2],a)};e.HelioVector=function(a,b){b=e.MakeTime(b);if(a in t)return A(t[a],b);if(a in wa){a:{var c=$jscomp.makeIterator(wa[a]);for(a=c.next();!a.done;a=c.next()){a=a.value;var d=a.tt;var h=a.tt+a.ndays;d=(2*b.tt-(h+d))/(h-d);if(-1<=d&&1>=d){var k=[];for(h=0;3>h;++h){var l=1;var g=a.coeff[0][h];var m=d;g+=a.coeff[1][h]*m;for(c=2;cg;++g){h=e.HelioVector(a,l);c&&(d=A(t.Earth,l));k=new r(h.x-d.x,h.y-d.y,h.z-d.z,b);var m=b.AddDays(-k.Length()/173.1446326846693);h=Math.abs(m.tt-l.tt);if(1E-9>h)return k;l=m}throw"Light-travel time solver did not converge: dt="+h;};e.Search=function(a,b,c,d){var h=Math.abs((d&&d.dt_tolerance_seconds|| -1)/86400),k=d&&d.init_f1||a(b),l=d&&d.init_f2||a(c),g,m=0;d=d&&d.iter_limit||20;for(var f=!0;;){if(++m>d)throw"Excessive iteration in Search()";var n=new C(b.ut+.5*(c.ut-b.ut)),p=n.ut-b.ut;if(Math.abs(p)(p.ut-b.ut)*(p.ut-c.ut)&&0>(r.ut-b.ut)*(r.ut-c.ut))){q=a(p);var u=a(r);if(0>q&&0<=u){k= +1)/86400),k=d&&d.init_f1||a(b),l=d&&d.init_f2||a(c),g,m=0;d=d&&d.iter_limit||20;for(var f=!0;;){if(++m>d)throw"Excessive iteration in Search()";var n=new C(b.ut+.5*(c.ut-b.ut)),p=n.ut-b.ut;if(Math.abs(p)(p.ut-b.ut)*(p.ut-c.ut)&&0>(r.ut-b.ut)*(r.ut-c.ut))){q=a(p);var u=a(r);if(0>q&&0<=u){k= q;l=u;b=p;c=r;g=t;f=!1;continue}}}}if(0>k&&0<=g)c=n,l=g;else if(0>g&&0<=l)b=n,k=g;else return null}};e.SearchSunLongitude=function(a,b,c){b=e.MakeTime(b);c=b.AddDays(c);return e.Search(function(b){b=e.SunPosition(b);return J(b.elon-a)},b,c)};e.LongitudeFromSun=function(a,b){if("Earth"===a)throw"The Earth does not have a longitude as seen from itself.";b=e.MakeTime(b);a=e.GeoVector(a,b,!1);a=e.Ecliptic(a.x,a.y,a.z);b=e.GeoVector("Sun",b,!1);b=e.Ecliptic(b.x,b.y,b.z);for(b=a.elon-b.elon;0>b;)b+=360; -for(;360<=b;)b-=360;return b};e.AngleFromSun=function(a,b){if("Earth"==a)throw"The Earth does not have an angle as seen from itself.";var c=e.GeoVector("Sun",b,!0);a=e.GeoVector(a,b,!0);return v(c,a)};e.EclipticLongitude=function(a,b){if("Sun"===a)throw"Cannot calculate heliocentric longitude of the Sun.";a=e.HelioVector(a,b);return e.Ecliptic(a.x,a.y,a.z).elon};var Na=function(a,b,c,d,e,k,l,g){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=e;this.gc=k;this.hc=l;this.ring_tilt=g};e.Illumination=function(a,b){if("Earth"===a)throw"The illumination of the Earth is not defined.";var c=e.MakeTime(b),d=A(t.Earth,c);if("Sun"===a){var h=new r(-d.x,-d.y,-d.z,c);b=new r(0,0,0,c);d=0}else"Moon"===a?(h=e.GeoMoon(c),b=new r(d.x+h.x,d.y+h.y,d.z+h.z,c)):(b=e.HelioVector(a,b),h=new r(b.x-d.x,b.y-d.y,b.z-d.z,c)),d=v(h,b);var k=h.Length(),l=b.Length(),g=null;if("Sun"===a)a=Ja+5*Math.log10(k);else if("Moon"===a){a=.017453292519943295*d; +for(;360<=b;)b-=360;return b};e.AngleFromSun=function(a,b){if("Earth"==a)throw"The Earth does not have an angle as seen from itself.";var c=e.GeoVector("Sun",b,!0);a=e.GeoVector(a,b,!0);return v(c,a)};e.EclipticLongitude=function(a,b){if("Sun"===a)throw"Cannot calculate heliocentric longitude of the Sun.";a=e.HelioVector(a,b);return e.Ecliptic(a.x,a.y,a.z).elon};var Ma=function(a,b,c,d,e,k,l,g){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=e;this.gc=k;this.hc=l;this.ring_tilt=g};e.Illumination=function(a,b){if("Earth"===a)throw"The illumination of the Earth is not defined.";var c=e.MakeTime(b),d=A(t.Earth,c);if("Sun"===a){var h=new r(-d.x,-d.y,-d.z,c);b=new r(0,0,0,c);d=0}else"Moon"===a?(h=e.GeoMoon(c),b=new r(d.x+h.x,d.y+h.y,d.z+h.z,c)):(b=e.HelioVector(a,b),h=new r(b.x-d.x,b.y-d.y,b.z-d.z,c)),d=v(h,b);var k=h.Length(),l=b.Length(),g=null;if("Sun"===a)a=Ia+5*Math.log10(k);else if("Moon"===a){a=.017453292519943295*d; var m=a*a;a=-12.717+1.49*Math.abs(a)+.0431*m*m;a+=5*Math.log10(k/.002573570052980638*l)}else if("Saturn"===a)a=d,g=e.Ecliptic(h.x,h.y,h.z),m=.017453292519943295*g.elat,g=Math.asin(Math.sin(m)*Math.cos(.4897393881096089)-Math.cos(m)*Math.sin(.4897393881096089)*Math.sin(.017453292519943295*g.elon-.017453292519943295*(169.51+3.82E-5*c.tt))),m=Math.sin(Math.abs(g)),a=-9+.044*a+m*(-2.6+1.2*m)+5*Math.log10(l*k),g*=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>d?(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=d/100;a=a+p*(m+p*(f+p*n))+5*Math.log10(l*k)}return new Na(c,a,d,l,k,h,b,g)};e.SearchRelativeLongitude=function(a,b,c){function d(c){var d=e.EclipticLongitude(a,c);c=e.EclipticLongitude("Earth",c);return J(k*(c- +break;case "Venus":163.6>d?(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=d/100;a=a+p*(m+p*(f+p*n))+5*Math.log10(l*k)}return new Ma(c,a,d,l,k,h,b,g)};e.SearchRelativeLongitude=function(a,b,c){function d(c){var d=e.EclipticLongitude(a,c);c=e.EclipticLongitude("Earth",c);return J(k*(c- d)-b)}var h=x[a];if(!h)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=h.OrbitalPeriod>x.Earth.OrbitalPeriod?1:-1;h=H(a);c=e.MakeTime(c);var l=d(c);0g;++g){var m=-l/360*h;c=c.AddDays(m);if(1>86400*Math.abs(m))return c;m=l;l=d(c);30>Math.abs(m)&&m!==l&&(m/=m-l,.5m&&(h*=m))}throw"Relative longitude search failed to converge for "+a+" near "+c.toString()+ -" (error_angle = "+l+").";};e.MoonPhase=function(a){return e.LongitudeFromSun("Moon",a)};e.SearchMoonPhase=function(a,b,c){function d(b){b=e.MoonPhase(b);return J(b-a)}b=e.MakeTime(b);var h=d(b);0c)return null;c=Math.min(c,k+.9);h=b.AddDays(h);b=b.AddDays(c);return e.Search(d,h,b)};var Oa=function(a,b){this.quarter=a;this.time=b};e.SearchMoonQuarter=function(a){var b=e.MoonPhase(a);b=(Math.floor(b/90)+1)%4;return(a=e.SearchMoonPhase(90*b,a,10))&&new Oa(b, -a)};e.NextMoonQuarter=function(a){a=new Date(a.time.date.getTime()+5184E5);return e.SearchMoonQuarter(a)};e.SearchRiseSet=function(a,b,c,d,h){function k(d){var f=e.Equator(a,d,b,!0,!0);d=e.Horizon(d,b,f.ra,f.dec).altitude+l/f.dist*57.29577951308232+Ka;return c*d}var l={Sun:.0046504672612422675,Moon:1.1618480877914597E-5}[a]||0;if("Earth"===a)throw"Cannot find rise or set time of the Earth.";if(1===c){var g=12;var m=0}else if(-1===c)g=0,m=12;else throw"Astronomy.SearchRiseSet: Invalid direction parameter "+ -c+" -- must be +1 or -1";d=e.MakeTime(d);var f=k(d);var n;if(0=f&&0=d.ut+h)return null;p=f.time;f=k(f.time);n=k(q.time)}};var Pa=function(a,b){this.time=a;this.hor=b};e.SearchHourAngle=function(a,b,c,d){d=e.MakeTime(d);var h=0;if("Earth"===a)throw"Cannot search for hour angle of the Earth."; -if(0>c||24<=c)throw"Invalid hour angle "+c;for(;;){++h;var k=L(d),l=e.Equator(a,d,b,!0,!0);k=(c+l.ra-b.longitude/15-k)%24;1===h?0>k&&(k+=24):-12>k?k+=24:123600*Math.abs(k))return a=e.Horizon(d,b,l.ra,l.dec,"normal"),new Pa(d,a);d=d.AddDays(k/24*.9972695717592592)}};var Qa=function(a,b,c,d){this.mar_equinox=a;this.jun_solstice=b;this.sep_equinox=c;this.dec_solstice=d};e.Seasons=function(a){function b(b,c,d){c=new Date(Date.UTC(a,c-1,d));b=e.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),d=b(90,6,19),h=b(180,9,21),k=b(270,12,20);return new Qa(c,d,h,k)};var Ra=function(a,b,c,d){this.time=a;this.visibility=b;this.elongation=c;this.ecliptic_separation=d};e.Elongation=function(a,b){b=e.MakeTime(b);var c=e.LongitudeFromSun(a,b);if(180c)return null;c=Math.min(c,k+.9);h=b.AddDays(h);b=b.AddDays(c);return e.Search(d,h,b)};var Na=function(a,b){this.quarter=a;this.time=b};e.SearchMoonQuarter=function(a){var b=e.MoonPhase(a);b=(Math.floor(b/90)+1)%4;return(a=e.SearchMoonPhase(90*b,a,10))&&new Na(b, +a)};e.NextMoonQuarter=function(a){a=new Date(a.time.date.getTime()+5184E5);return e.SearchMoonQuarter(a)};e.SearchRiseSet=function(a,b,c,d,h){function k(d){var f=e.Equator(a,d,b,!0,!0);d=e.Horizon(d,b,f.ra,f.dec).altitude+l/f.dist*57.29577951308232+Ja;return c*d}var l={Sun:.0046504672612422675,Moon:1.1618480877914597E-5}[a]||0;if("Earth"===a)throw"Cannot find rise or set time of the Earth.";if(1===c){var g=12;var m=0}else if(-1===c)g=0,m=12;else throw"Astronomy.SearchRiseSet: Invalid direction parameter "+ +c+" -- must be +1 or -1";d=e.MakeTime(d);var f=k(d);var n;if(0=f&&0=d.ut+h)return null;p=f.time;f=k(f.time);n=k(q.time)}};var Oa=function(a,b){this.time=a;this.hor=b};e.SearchHourAngle=function(a,b,c,d){d=e.MakeTime(d);var h=0;if("Earth"===a)throw"Cannot search for hour angle of the Earth."; +if(0>c||24<=c)throw"Invalid hour angle "+c;for(;;){++h;var k=L(d),l=e.Equator(a,d,b,!0,!0);k=(c+l.ra-b.longitude/15-k)%24;1===h?0>k&&(k+=24):-12>k?k+=24:123600*Math.abs(k))return a=e.Horizon(d,b,l.ra,l.dec,"normal"),new Oa(d,a);d=d.AddDays(k/24*.9972695717592592)}};var Pa=function(a,b,c,d){this.mar_equinox=a;this.jun_solstice=b;this.sep_equinox=c;this.dec_solstice=d};e.Seasons=function(a){function b(b,c,d){c=new Date(Date.UTC(a,c-1,d));b=e.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),d=b(90,6,19),h=b(180,9,21),k=b(270,12,20);return new Pa(c,d,h,k)};var Qa=function(a,b,c,d){this.time=a;this.visibility=b;this.elongation=c;this.ecliptic_separation=d};e.Elongation=function(a,b){b=e.MakeTime(b);var c=e.LongitudeFromSun(a,b);if(180=++h;){var k=e.EclipticLongitude(a,b),l=e.EclipticLongitude("Earth",b);l=J(k-l);var g=k=void 0,m=g=k=void 0;l>=-d.s1&&l<+d.s1?(m=0,k=+d.s1,g=+d.s2):l>=+d.s2||l<-d.s2?(m=0,k=-d.s2,g=-d.s1):0<=l?(m=-H(a)/4,k=+d.s1,g=+d.s2):(m=-H(a)/4,k=-d.s2,g=-d.s1);l=b.AddDays(m);k=e.SearchRelativeLongitude(a, k,l);g=e.SearchRelativeLongitude(a,g,k);l=c(k);if(0<=l)throw"SearchMaxElongation: internal error: m1 = "+l;m=c(g);if(0>=m)throw"SearchMaxElongation: internal error: m2 = "+m;l=e.Search(c,k,g,{init_f1:l,init_f2:m,dt_tolerance_seconds:10});if(!l)throw"SearchMaxElongation: failed search iter "+h+" (t1="+k.toString()+", t2="+g.toString()+")";if(l.tt>=b.tt)return e.Elongation(a,l);b=g.AddDays(1)}throw"SearchMaxElongation: failed to find event after 2 tries.";};e.SearchPeakMagnitude=function(a,b){function c(b){var c= b.AddDays(-.005);b=b.AddDays(.005);c=e.Illumination(a,c).mag;return(e.Illumination(a,b).mag-c)/.01}if("Venus"!==a)throw"SearchPeakMagnitude currently works for Venus only.";b=e.MakeTime(b);for(var d=0;2>=++d;){var h=e.EclipticLongitude(a,b),k=e.EclipticLongitude("Earth",b);k=J(h-k);var l=h=void 0,g=l=h=void 0;-10<=k&&10>k?(g=0,h=10,l=30):30<=k||-30>k?(g=0,h=-30,l=-10):0<=k?(g=-H(a)/4,h=10,l=30):(g=-H(a)/4,h=-30,l=-10);k=b.AddDays(g);h=e.SearchRelativeLongitude(a,h,k);l=e.SearchRelativeLongitude(a, l,h);k=c(h);if(0<=k)throw"SearchPeakMagnitude: internal error: m1 = "+k;g=c(l);if(0>=g)throw"SearchPeakMagnitude: internal error: m2 = "+g;k=e.Search(c,h,l,{init_f1:k,init_f2:g,dt_tolerance_seconds:10});if(!k)throw"SearchPeakMagnitude: failed search iter "+d+" (t1="+h.toString()+", t2="+l.toString()+")";if(k.tt>=b.tt)return e.Illumination(a,k);b=l.AddDays(1)}throw"SearchPeakMagnitude: failed to find event after 2 tries.";};var R=function(a,b,c){this.time=a;this.kind=b;this.dist_au=c;this.dist_km= 1.4959787069098932E8*c};e.SearchLunarApsis=function(a){function b(a){var b=a.AddDays(-5E-4);a=a.AddDays(5E-4);b=z(b).distance_au;return(z(a).distance_au-b)/.001}function c(a){return-b(a)}a=e.MakeTime(a);for(var d=b(a),h=0;59.061176>5*h;++h){var k=a.AddDays(5),l=b(k);if(0>=d*l){if(0>d||0l){a=e.Search(c,a,k,{init_f1:-d,init_f2:-l});if(null== a)throw"SearchLunarApsis INTERNAL ERROR: apogee search failed!";d=z(a).distance_au;return new R(a,1,d)}throw"SearchLunarApsis INTERNAL ERROR: cannot classify apsis event!";}a=k;d=l}throw"SearchLunarApsis INTERNAL ERROR: could not find apsis within 2 synodic months of start date.";};e.NextLunarApsis=function(a){var b=e.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};e.SearchPlanetApsis=function(a,b){function c(b){var c=b.AddDays(-5E-4);b=b.AddDays(5E-4);c=e.HelioDistance(a,c);return(e.HelioDistance(a,b)-c)/.001}function d(a){return-c(a)}if("Neptune"===a)return Ca(b);for(var h=x[a].OrbitalPeriod,k=h/6,l=c(b),g=0;g*k<2*h;++g){var m=b.AddDays(k),f=c(m);if(0>=l*f){h=k=void 0;if(0>l||0f)k=d,h=1;else throw"Internal error with slopes in SearchPlanetApsis";b=e.Search(k,b,m,1);if(null==b)throw"Failed to find slope transition in planetary apsis search."; +b.kind+" @ "+b.time.toString();return b};e.SearchPlanetApsis=function(a,b){function c(b){var c=b.AddDays(-5E-4);b=b.AddDays(5E-4);c=e.HelioDistance(a,c);return(e.HelioDistance(a,b)-c)/.001}function d(a){return-c(a)}if("Neptune"===a)return Ba(b);for(var h=x[a].OrbitalPeriod,k=h/6,l=c(b),g=0;g*k<2*h;++g){var m=b.AddDays(k),f=c(m);if(0>=l*f){h=k=void 0;if(0>l||0f)k=d,h=1;else throw"Internal error with slopes in SearchPlanetApsis";b=e.Search(k,b,m,1);if(null==b)throw"Failed to find slope transition in planetary apsis search."; l=e.HelioDistance(a,b);return new R(b,h,l)}b=m;l=f}throw"Internal error: should have found planetary apsis within 2 orbital periods.";};e.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=e.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};e.InverseRotation=function(a){return new w([[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]]])};e.CombineRotation=function(a,b){return new w([[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]]])};e.VectorFromSphere=function(a,b){var c=.017453292519943295*a.lat,d=.017453292519943295*a.lon,e=a.dist*Math.cos(c);return new r(e*Math.cos(d),e*Math.sin(d),a.dist*Math.sin(c),b)};e.VectorFromEquator=function(a,b){return e.VectorFromSphere(new W(a.dec, @@ -159,10 +159,10 @@ return a};e.VectorFromHorizon=function(a,b,c){var d=oa(a.lon);c=a.lat+e.InverseR b;c-=d}};e.RotateVector=function(a,b){return new r(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)};e.Rotation_EQJ_ECL=function(){return new w([[1,0,0],[0,.9174821430670688,-.3977769691083922],[0,.3977769691083922,.9174821430670688]])};e.Rotation_ECL_EQJ=function(){return new w([[1,0,0],[0,.9174821430670688,.3977769691083922],[0,-.3977769691083922,.9174821430670688]])};e.Rotation_EQJ_EQD=function(a){var b= Y(0,a.tt);a=aa(a,0);return e.CombineRotation(b,a)};e.Rotation_EQD_EQJ=function(a){var b=aa(a,1);a=Y(a.tt,0);return e.CombineRotation(b,a)};e.Rotation_EQD_HOR=function(a,b){var c=Math.sin(.017453292519943295*b.latitude),d=Math.cos(.017453292519943295*b.latitude),e=Math.sin(.017453292519943295*b.longitude),k=Math.cos(.017453292519943295*b.longitude);b=[d*k,d*e,c];c=[-c*k,-c*e,d];e=[e,-k,0];a=-15*L(a);b=D(a,b);c=D(a,c);a=D(a,e);return new w([[c[0],a[0],b[0]],[c[1],a[1],b[1]],[c[2],a[2],b[2]]])};e.Rotation_HOR_EQD= function(a,b){a=e.Rotation_EQD_HOR(a,b);return e.InverseRotation(a)};e.Rotation_HOR_EQJ=function(a,b){b=e.Rotation_HOR_EQD(a,b);a=e.Rotation_EQD_EQJ(a);return e.CombineRotation(b,a)};e.Rotation_EQJ_HOR=function(a,b){a=e.Rotation_HOR_EQJ(a,b);return e.InverseRotation(a)};e.Rotation_EQD_ECL=function(a){a=e.Rotation_EQD_EQJ(a);var b=e.Rotation_EQJ_ECL();return e.CombineRotation(a,b)};e.Rotation_ECL_EQD=function(a){a=e.Rotation_EQD_ECL(a);return e.InverseRotation(a)};e.Rotation_ECL_HOR=function(a,b){var c= -e.Rotation_ECL_EQD(a);a=e.Rotation_EQD_HOR(a,b);return e.CombineRotation(c,a)};e.Rotation_HOR_ECL=function(a,b){a=e.Rotation_ECL_HOR(a,b);return e.InverseRotation(a)};var Sa=[["And","Andromeda"],["Ant","Antila"],["Aps","Apus"],["Aql","Aquila"],["Aqr","Aquarius"],["Ara","Ara"],["Ari","Aries"],["Aur","Auriga"],["Boo","Bootes"],["Cae","Caelum"],["Cam","Camelopardis"],["Cap","Capricornus"],["Car","Carina"],["Cas","Cassiopeia"],["Cen","Centaurus"],["Cep","Cepheus"],["Cet","Cetus"],["Cha","Chamaeleon"], +e.Rotation_ECL_EQD(a);a=e.Rotation_EQD_HOR(a,b);return e.CombineRotation(c,a)};e.Rotation_HOR_ECL=function(a,b){a=e.Rotation_ECL_HOR(a,b);return e.InverseRotation(a)};var Ra=[["And","Andromeda"],["Ant","Antila"],["Aps","Apus"],["Aql","Aquila"],["Aqr","Aquarius"],["Ara","Ara"],["Ari","Aries"],["Aur","Auriga"],["Boo","Bootes"],["Cae","Caelum"],["Cam","Camelopardis"],["Cap","Capricornus"],["Car","Carina"],["Cas","Cassiopeia"],["Cen","Centaurus"],["Cep","Cepheus"],["Cet","Cetus"],["Cha","Chamaeleon"], ["Cir","Circinus"],["CMa","Canis Major"],["CMi","Canis Minor"],["Cnc","Cancer"],["Col","Columba"],["Com","Coma Berenices"],["CrA","Corona Australis"],["CrB","Corona Borealis"],["Crt","Crater"],["Cru","Crux"],["Crv","Corvus"],["CVn","Canes Venatici"],["Cyg","Cygnus"],["Del","Delphinus"],["Dor","Dorado"],["Dra","Draco"],["Equ","Equuleus"],["Eri","Eridanus"],["For","Fornax"],["Gem","Gemini"],["Gru","Grus"],["Her","Hercules"],["Hor","Horologium"],["Hya","Hydra"],["Hyi","Hydrus"],["Ind","Indus"],["Lac", "Lacerta"],["Leo","Leo"],["Lep","Lepus"],["Lib","Libra"],["LMi","Leo Minor"],["Lup","Lupus"],["Lyn","Lynx"],["Lyr","Lyra"],["Men","Mensa"],["Mic","Microscopium"],["Mon","Monoceros"],["Mus","Musca"],["Nor","Norma"],["Oct","Octans"],["Oph","Ophiuchus"],["Ori","Orion"],["Pav","Pavo"],["Peg","Pegasus"],["Per","Perseus"],["Phe","Phoenix"],["Pic","Pictor"],["PsA","Pisces Austrinus"],["Psc","Pisces"],["Pup","Puppis"],["Pyx","Pyxis"],["Ret","Reticulum"],["Scl","Sculptor"],["Sco","Scorpius"],["Sct","Scutum"], -["Ser","Serpens"],["Sex","Sextans"],["Sge","Sagitta"],["Sgr","Sagittarius"],["Tau","Taurus"],["Tel","Telescopium"],["TrA","Triangulum Australe"],["Tri","Triangulum"],["Tuc","Tucana"],["UMa","Ursa Major"],["UMi","Ursa Minor"],["Vel","Vela"],["Vir","Virgo"],["Vol","Volans"],["Vul","Vulpecula"]],Ta=[[83,0,8640,2112],[83,2880,5220,2076],[83,7560,8280,2068],[83,6480,7560,2064],[15,0,2880,2040],[10,3300,3840,1968],[15,0,1800,1920],[10,3840,5220,1920],[83,6300,6480,1920],[33,7260,7560,1920],[15,0,1263,1848], +["Ser","Serpens"],["Sex","Sextans"],["Sge","Sagitta"],["Sgr","Sagittarius"],["Tau","Taurus"],["Tel","Telescopium"],["TrA","Triangulum Australe"],["Tri","Triangulum"],["Tuc","Tucana"],["UMa","Ursa Major"],["UMi","Ursa Minor"],["Vel","Vela"],["Vir","Virgo"],["Vol","Volans"],["Vul","Vulpecula"]],Sa=[[83,0,8640,2112],[83,2880,5220,2076],[83,7560,8280,2068],[83,6480,7560,2064],[15,0,2880,2040],[10,3300,3840,1968],[15,0,1800,1920],[10,3840,5220,1920],[83,6300,6480,1920],[33,7260,7560,1920],[15,0,1263,1848], [10,4140,4890,1848],[83,5952,6300,1800],[15,7260,7440,1800],[10,2868,3300,1764],[33,3300,4080,1764],[83,4680,5952,1680],[13,1116,1230,1632],[33,7350,7440,1608],[33,4080,4320,1596],[15,0,120,1584],[83,5040,5640,1584],[15,8490,8640,1584],[33,4320,4860,1536],[33,4860,5190,1512],[15,8340,8490,1512],[10,2196,2520,1488],[33,7200,7350,1476],[15,7393.2,7416,1462],[10,2520,2868,1440],[82,2868,3030,1440],[33,7116,7200,1428],[15,7200,7393.2,1428],[15,8232,8340,1418],[13,0,876,1404],[33,6990,7116,1392],[13,612, 687,1380],[13,876,1116,1368],[10,1116,1140,1368],[15,8034,8232,1350],[10,1800,2196,1344],[82,5052,5190,1332],[33,5190,6990,1332],[10,1140,1200,1320],[15,7968,8034,1320],[15,7416,7908,1316],[13,0,612,1296],[50,2196,2340,1296],[82,4350,4860,1272],[33,5490,5670,1272],[15,7908,7968,1266],[10,1200,1800,1260],[13,8232,8400,1260],[33,5670,6120,1236],[62,735,906,1212],[33,6120,6564,1212],[13,0,492,1200],[62,492,600,1200],[50,2340,2448,1200],[13,8400,8640,1200],[82,4860,5052,1164],[13,0,402,1152],[13,8490, 8640,1152],[39,6543,6564,1140],[33,6564,6870,1140],[30,6870,6900,1140],[62,600,735,1128],[82,3030,3300,1128],[13,60,312,1104],[82,4320,4350,1080],[50,2448,2652,1068],[30,7887,7908,1056],[30,7875,7887,1050],[30,6900,6984,1044],[82,3300,3660,1008],[82,3660,3882,960],[8,5556,5670,960],[39,5670,5880,960],[50,3330,3450,954],[0,0,906,882],[62,906,924,882],[51,6969,6984,876],[62,1620,1689,864],[30,7824,7875,864],[44,7875,7920,864],[7,2352,2652,852],[50,2652,2790,852],[0,0,720,840],[44,7920,8214,840],[44, @@ -176,12 +176,12 @@ e.Rotation_ECL_EQD(a);a=e.Rotation_EQD_HOR(a,b);return e.CombineRotation(c,a)};e -1032],[35,1230,1392,-1056],[71,5911.5,6420,-1092],[24,6420,6900,-1092],[76,6900,7320,-1092],[53,7320,7680,-1092],[35,1080,1230,-1104],[9,1620,1740,-1116],[49,5520,5640,-1152],[63,0,840,-1156],[35,960,1080,-1176],[40,1470,1536,-1176],[9,1536,1620,-1176],[38,7680,7920,-1200],[67,2160,2880,-1218],[84,2880,2940,-1218],[35,870,960,-1224],[40,1380,1470,-1224],[63,0,660,-1236],[12,2160,2220,-1260],[84,2940,3042,-1272],[40,1260,1380,-1276],[32,1380,1440,-1276],[63,0,570,-1284],[35,780,870,-1296],[64,1620, 1800,-1296],[49,5418,5520,-1296],[84,3042,3180,-1308],[12,2220,2340,-1320],[14,4260,4620,-1320],[49,5100,5418,-1320],[56,5418,5520,-1320],[32,1440,1560,-1356],[84,3180,3960,-1356],[14,3960,4050,-1356],[5,6300,6480,-1368],[78,6480,7320,-1368],[38,7920,8400,-1368],[40,1152,1260,-1380],[64,1800,1980,-1380],[12,2340,2460,-1392],[63,0,480,-1404],[35,480,780,-1404],[63,8400,8640,-1404],[32,1560,1650,-1416],[56,5520,5911.5,-1440],[43,7320,7680,-1440],[64,1980,2160,-1464],[18,5460,5520,-1464],[5,5911.5,5970, -1464],[18,5370,5460,-1526],[5,5970,6030,-1526],[64,2160,2460,-1536],[12,2460,3252,-1536],[14,4050,4260,-1536],[27,4260,4620,-1536],[14,4620,5232,-1536],[18,4860,4920,-1560],[5,6030,6060,-1560],[40,780,1152,-1620],[69,1152,1650,-1620],[18,5310,5370,-1620],[5,6060,6300,-1620],[60,6300,6480,-1620],[81,7920,8400,-1620],[32,1650,2370,-1680],[18,4920,5310,-1680],[79,5310,6120,-1680],[81,0,480,-1800],[42,1260,1650,-1800],[86,2370,3252,-1800],[12,3252,4050,-1800],[55,4050,4920,-1800],[60,6480,7680,-1800], -[43,7680,8400,-1800],[81,8400,8640,-1800],[81,270,480,-1824],[42,0,1260,-1980],[17,2760,4920,-1980],[2,4920,6480,-1980],[52,1260,2760,-2040],[57,0,8640,-2160]],ha,xa,Ua=function(a,b,c,d){this.symbol=a;this.name=b;this.ra1875=c;this.dec1875=d};e.Constellation=function(a,b){if(-90>b||90a&&(a+=24);ha||(ha=e.Rotation_EQJ_EQD(new C(-45655.74141261017)),xa=new C(0));a=new ba(a,b,1);a=e.VectorFromEquator(a,xa);a=e.RotateVector(ha,a);a=e.EquatorFromVector(a); -b=10/240;for(var c=b/15,d=$jscomp.makeIterator(Ta),h=d.next();!h.done;h=d.next()){h=h.value;var k=h[1]*c,l=h[2]*c;if(h[3]*b<=a.dec&&k<=a.ra&&a.rab;++b){var c=e.SearchMoonPhase(180,a,40);if(null===c)throw"Cannot find full moon.";a=57.29577951308232*z(c).geo_eclip_lat;if(1.8>Math.abs(a)&&(a=Ea(c),a.rb;++b){var c=e.SearchMoonPhase(0,a,40);if(null===c)throw"Cannot find new moon";a=57.29577951308232*z(c).geo_eclip_lat;if(1.8>Math.abs(a)&&(a=Fa(c),a.rb||90a&&(a+=24);ha||(ha=e.Rotation_EQJ_EQD(new C(-45655.74141261017)),xa=new C(0));a=new ba(a,b,1);a=e.VectorFromEquator(a,xa);a=e.RotateVector(ha,a);a=e.EquatorFromVector(a); +b=10/240;for(var c=b/15,d=$jscomp.makeIterator(Sa),h=d.next();!h.done;h=d.next()){h=h.value;var k=h[1]*c,l=h[2]*c;if(h[3]*b<=a.dec&&k<=a.ra&&a.rab;++b){var c=e.SearchMoonPhase(180,a,40);if(null===c)throw"Cannot find full moon.";a=57.29577951308232*z(c).geo_eclip_lat;if(1.8>Math.abs(a)&&(a=Da(c),a.rb;++b){var c=e.SearchMoonPhase(0,a,40);if(null===c)throw"Cannot find new moon";a=57.29577951308232*z(c).geo_eclip_lat;if(1.8>Math.abs(a)&&(a=Ea(c),a.r=d?d+=360:180a.r)throw"Unexpected shadow distance from geoid intersection = "+a.r;k=.014Math.abs(c)){var d=Ha(a,b);if(d.re.AngleFromSun(a,d)&&(b=Ga(a,c,d),b.ra.r)throw"Unexpected shadow distance from geoid intersection = "+a.r;k=.014Math.abs(c)){var d=Ga(a,b);if(d.re.AngleFromSun(a,d)&&(b=Fa(a,c,d),b.r +### class TransitInfo + +**Information about a transit of Mercury or Venus, as seen from the Earth.** + +Returned by [`SearchTransit`](#SearchTransit) or [`NextTransit`](#NextTransit) to report +information about a transit of Mercury or Venus. +A transit is when Mercury or Venus passes between the Sun and Earth so that +the other planet is seen in silhouette against the Sun. +The calculations are performed from the point of view of a geocentric observer. +peak : Time + When the planet is most aligned with the Sun, as seen from the Earth. +finish : Time + The date and time at the end of the transit. + This is the moment the planet is last seen against the Sun in its background. +separation : float + The minimum angular separation, in arcminutes, between the centers of the Sun and the planet. + This angle pertains to the time stored in `peak`. + +| Type | Attribute | Description | +| --- | --- | --- | +| [`Time`](#Time) | `start` | The date and time at the beginning of the transit. This is the moment the planet first becomes visible against the Sun in its background. | + +--- + ### class Vector @@ -1422,6 +1447,24 @@ See [`SearchPlanetApsis`](#SearchPlanetApsis) for more details. --- + +### NextTransit(body, prevTransitTime) + +**Searches for another transit of Mercury or Venus.** + +After calling [`SearchTransit`](#SearchTransit) to find a transit of Mercury or Venus, +this function finds the next transit after that. +Keep calling this function as many times as you want to keep finding more transits. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Body`](#Body) | `body` | The planet whose transit is to be found. Must be `Body.Mercury` or `Body.Venus`. | +| [`Time`](#Time) | `prevTransitTime` | A date and time near the previous transit. | + +### Returns: [`TransitInfo`](#TransitInfo) + +--- + ### RefractionAngle(refraction, altitude) @@ -2119,6 +2162,26 @@ limitDays : float --- + +### SearchTransit(body, startTime) + +**Searches for the first transit of Mercury or Venus after a given date.** + +Finds the first transit of Mercury or Venus after a specified date. +A transit is when an inferior planet passes between the Sun and the Earth +so that the silhouette of the planet is visible against the Sun in the background. +To continue the search, pass the `finish` time in the returned structure to +[`NextTransit`](#NextTransit). + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Body`](#Body) | `body` | The planet whose transit is to be found. Must be `Body.Mercury` or `Body.Venus`. | +| [`Time`](#Time) | `startTime` | The date and time for starting the search for a transit. | + +### Returns: [`TransitInfo`](#TransitInfo) + +--- + ### Seasons(year) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index b561eade..5dbd737b 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -6653,6 +6653,20 @@ def _LocalMoonShadow(time, observer): return _CalcShadow(_MOON_MEAN_RADIUS_KM, time, o, m) +def _PlanetShadow(body, planet_radius_km, time): + # Calculate light-travel-corrected vector from Earth to planet. + g = GeoVector(body, time, False) + # Calculate light-travel-corrected vector from Earth to Sun. + e = GeoVector(Body.Sun, time, False) + # Deduce light-travel-corrected vector from Sun to planet. + p = Vector(g.x - e.x, g.y - e.y, g.z - e.z, time) + # Calcluate Earth's position from the planet's point of view. + e.x = -g.x + e.y = -g.y + e.z = -g.z + return _CalcShadow(planet_radius_km, time, e, p) + + def _ShadowDistanceSlope(shadowfunc, time): dt = 1.0 / 86400.0 t1 = time.AddDays(-dt) @@ -6662,6 +6676,14 @@ def _ShadowDistanceSlope(shadowfunc, time): return (shadow2.r - shadow1.r) / dt +def _PlanetShadowSlope(context, time): + (body, planet_radius_km) = context + dt = 1.0 / 86400.0 + shadow1 = _PlanetShadow(body, planet_radius_km, time.AddDays(-dt)) + shadow2 = _PlanetShadow(body, planet_radius_km, time.AddDays(+dt)) + return (shadow2.r - shadow1.r) / dt + + def _PeakEarthShadow(search_center_time): window = 0.03 # initial search window, in days, before/after given time t1 = search_center_time.AddDays(-window) @@ -6686,6 +6708,14 @@ def _PeakLocalMoonShadow(search_center_time, observer): tx = Search(_ShadowDistanceSlope, lambda time: _LocalMoonShadow(time, observer), t1, t2, 1.0) return _LocalMoonShadow(tx, observer) +def _PeakPlanetShadow(body, planet_radius_km, search_center_time): + # Search for when the body's shadow is closest to the center of the Earth. + window = 1.0 # days before/after inferior conjunction to search for minimum shadow distance. + t1 = search_center_time.AddDays(-window) + t2 = search_center_time.AddDays(+window) + tx = Search(_PlanetShadowSlope, (body, planet_radius_km), t1, t2, 1.0) + return _PlanetShadow(body, planet_radius_km, tx) + class _ShadowDiffContext: def __init__(self, radius_limit, direction): self.radius_limit = radius_limit @@ -7266,3 +7296,135 @@ def NextLocalSolarEclipse(prevEclipseTime, observer): """ startTime = prevEclipseTime.AddDays(10.0) return SearchLocalSolarEclipse(startTime, observer) + + +class TransitInfo: + """Information about a transit of Mercury or Venus, as seen from the Earth. + + Returned by #SearchTransit or #NextTransit to report + information about a transit of Mercury or Venus. + A transit is when Mercury or Venus passes between the Sun and Earth so that + the other planet is seen in silhouette against the Sun. + + The calculations are performed from the point of view of a geocentric observer. + + Attributes + ---------- + start : Time + The date and time at the beginning of the transit. + This is the moment the planet first becomes visible against the Sun in its background. + + peak : Time + When the planet is most aligned with the Sun, as seen from the Earth. + + finish : Time + The date and time at the end of the transit. + This is the moment the planet is last seen against the Sun in its background. + + separation : float + The minimum angular separation, in arcminutes, between the centers of the Sun and the planet. + This angle pertains to the time stored in `peak`. + """ + def __init__(self, start, peak, finish, separation): + self.start = start + self.peak = peak + self.finish = finish + self.separation = separation + + +def _PlanetShadowBoundary(context, time): + (body, planet_radius_km, direction) = context + shadow = _PlanetShadow(body, planet_radius_km, time) + return direction * (shadow.r - shadow.p) + + +def _PlanetTransitBoundary(body, planet_radius_km, t1, t2, direction): + # Search for the time the planet's penumbra begins/ends making contact with the center of the Earth. + # context = new SearchContext_PlanetShadowBoundary(body, planet_radius_km, direction); + tx = Search(_PlanetShadowBoundary, (body, planet_radius_km, direction), t1, t2, 1.0) + if tx is None: + raise Error('Planet transit boundary search failed') + return tx + + +def SearchTransit(body, startTime): + """Searches for the first transit of Mercury or Venus after a given date. + + Finds the first transit of Mercury or Venus after a specified date. + A transit is when an inferior planet passes between the Sun and the Earth + so that the silhouette of the planet is visible against the Sun in the background. + To continue the search, pass the `finish` time in the returned structure to + #NextTransit. + + Parameters + ---------- + body : Body + The planet whose transit is to be found. Must be `Body.Mercury` or `Body.Venus`. + startTime : Time + The date and time for starting the search for a transit. + + Returns + ------- + TransitInfo + """ + threshold_angle = 0.4 # maximum angular separation to attempt transit calculation + dt_days = 1.0 + + # Validate the planet and find its mean radius. + if body == Body.Mercury: + planet_radius_km = 2439.7 + elif body == Body.Venus: + planet_radius_km = 6051.8 + else: + raise InvalidBodyError() + + search_time = startTime + while True: + # Search for the next inferior conjunction of the given planet. + # This is the next time the Earth and the other planet have the same + # ecliptic longitude as seen from the Sun. + conj = SearchRelativeLongitude(body, 0.0, search_time) + + # Calculate the angular separation between the body and the Sun at this time. + conj_separation = AngleFromSun(body, conj) + + if conj_separation < threshold_angle: + # The planet's angular separation from the Sun is small enough + # to consider it a transit candidate. + # Search for the moment when the line passing through the Sun + # and planet are closest to the Earth's center. + shadow = _PeakPlanetShadow(body, planet_radius_km, conj) + if shadow.r < shadow.p: # does the planet's penumbra touch the Earth's center? + # Find the beginning and end of the penumbral contact. + time_before = shadow.time.AddDays(-dt_days) + start = _PlanetTransitBoundary(body, planet_radius_km, time_before, shadow.time, -1.0) + time_after = shadow.time.AddDays(+dt_days) + finish = _PlanetTransitBoundary(body, planet_radius_km, shadow.time, time_after, +1.0) + min_separation = 60.0 * AngleFromSun(body, shadow.time) + return TransitInfo(start, shadow.time, finish, min_separation) + + # This inferior conjunction was not a transit. Try the next inferior conjunction. + search_time = conj.AddDays(10.0) + + + +def NextTransit(body, prevTransitTime): + """Searches for another transit of Mercury or Venus. + + After calling #SearchTransit to find a transit of Mercury or Venus, + this function finds the next transit after that. + Keep calling this function as many times as you want to keep finding more transits. + + Parameters + ---------- + body : Body + The planet whose transit is to be found. Must be `Body.Mercury` or `Body.Venus`. + prevTransitTime : Time + A date and time near the previous transit. + + Returns + ------- + TransitInfo + """ + startTime = prevTransitTime.AddDays(100.0) + return SearchTransit(body, startTime) \ No newline at end of file