From 4f842627dac5a6f0379003c813d2d6ea2eac0894 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 13 Jun 2020 13:45:59 -0400 Subject: [PATCH] Fixed mistake in GeoVector(SUN): we do need to correct for light-travel time. To be consistent, when calculating the geocentric position of the Sun, we do need to correct for light travel time just like we would for any other object. This reduces the maximum time error for predicting transits from 25 minutes to 11 minutes. Also had to disable aberration when calculating moon phases (longitude from Sun) in order to keep a good fit with test data. --- demo/c/test/culminate_correct.txt | 2 +- demo/c/test/positions_correct.txt | 2 +- demo/c/test/riseset_correct.txt | 4 ++-- demo/csharp/test/culminate_correct.txt | 2 +- demo/csharp/test/positions_correct.txt | 2 +- demo/csharp/test/riseset_correct.txt | 4 ++-- demo/nodejs/test/culminate_correct.txt | 2 +- demo/nodejs/test/positions_correct.txt | 2 +- demo/nodejs/test/riseset_correct.txt | 4 ++-- demo/python/test/culminate_correct.txt | 2 +- demo/python/test/positions_correct.txt | 2 +- demo/python/test/riseset_correct.txt | 4 ++-- generate/template/astronomy.c | 29 ++++++++++---------------- generate/template/astronomy.cs | 9 ++------ generate/template/astronomy.js | 7 ++----- generate/template/astronomy.py | 8 ++----- source/c/astronomy.c | 29 ++++++++++---------------- source/csharp/astronomy.cs | 9 ++------ source/js/astronomy.js | 7 ++----- source/js/astronomy.min.js | 6 +++--- source/python/astronomy.py | 8 ++----- 21 files changed, 53 insertions(+), 91 deletions(-) diff --git a/demo/c/test/culminate_correct.txt b/demo/c/test/culminate_correct.txt index 06419b9b..0f21f547 100644 --- a/demo/c/test/culminate_correct.txt +++ b/demo/c/test/culminate_correct.txt @@ -1,5 +1,5 @@ search : 2015-02-28 00:00:00 UTC -Sun : 2015-02-28 18:12:32 UTC altitude= 52.13 azimuth=180.00 +Sun : 2015-02-28 18:12:31 UTC altitude= 52.13 azimuth=180.00 Moon : 2015-02-28 01:58:15 UTC altitude= 77.88 azimuth=180.00 Mercury : 2015-02-28 16:31:18 UTC altitude= 42.57 azimuth=180.00 Venus : 2015-02-28 20:03:56 UTC altitude= 63.24 azimuth=180.00 diff --git a/demo/c/test/positions_correct.txt b/demo/c/test/positions_correct.txt index b8250b6c..1197cc3e 100644 --- a/demo/c/test/positions_correct.txt +++ b/demo/c/test/positions_correct.txt @@ -1,7 +1,7 @@ UTC date = 2018-11-30T17:55:07.234Z BODY RA DEC AZ ALT -Sun 16.43 -21.68 180.90 22.72 +Sun 16.43 -21.68 180.91 22.72 Moon 11.30 7.81 266.73 14.07 Mercury 15.92 -18.41 189.05 25.50 Venus 13.81 -9.79 224.05 23.85 diff --git a/demo/c/test/riseset_correct.txt b/demo/c/test/riseset_correct.txt index d0fa1e6c..c5887676 100644 --- a/demo/c/test/riseset_correct.txt +++ b/demo/c/test/riseset_correct.txt @@ -1,5 +1,5 @@ search : 2019-06-14 21:45:24 UTC -sunrise : 2019-06-15 10:12:44 UTC -sunset : 2019-06-15 01:48:01 UTC +sunrise : 2019-06-15 10:12:42 UTC +sunset : 2019-06-15 01:48:00 UTC moonrise : 2019-06-14 23:02:43 UTC moonset : 2019-06-15 09:12:25 UTC diff --git a/demo/csharp/test/culminate_correct.txt b/demo/csharp/test/culminate_correct.txt index 04795519..d19fc5c3 100644 --- a/demo/csharp/test/culminate_correct.txt +++ b/demo/csharp/test/culminate_correct.txt @@ -1,5 +1,5 @@ search : 2015-02-28T00:00:00.000Z -Sun : 2015-02-28T18:12:32.334Z altitude=52.13 azimuth=180.00 +Sun : 2015-02-28T18:12:31.045Z altitude=52.13 azimuth=180.00 Moon : 2015-02-28T01:58:15.631Z altitude=77.88 azimuth=180.00 Mercury : 2015-02-28T16:31:18.940Z altitude=42.57 azimuth=180.00 Venus : 2015-02-28T20:03:56.711Z altitude=63.24 azimuth=180.00 diff --git a/demo/csharp/test/positions_correct.txt b/demo/csharp/test/positions_correct.txt index b8250b6c..1197cc3e 100644 --- a/demo/csharp/test/positions_correct.txt +++ b/demo/csharp/test/positions_correct.txt @@ -1,7 +1,7 @@ UTC date = 2018-11-30T17:55:07.234Z BODY RA DEC AZ ALT -Sun 16.43 -21.68 180.90 22.72 +Sun 16.43 -21.68 180.91 22.72 Moon 11.30 7.81 266.73 14.07 Mercury 15.92 -18.41 189.05 25.50 Venus 13.81 -9.79 224.05 23.85 diff --git a/demo/csharp/test/riseset_correct.txt b/demo/csharp/test/riseset_correct.txt index 9977e966..9c80491a 100644 --- a/demo/csharp/test/riseset_correct.txt +++ b/demo/csharp/test/riseset_correct.txt @@ -1,5 +1,5 @@ search : 2019-06-14T21:45:25.000Z -sunrise : 2019-06-15T10:12:44.103Z -sunset : 2019-06-15T01:48:01.922Z +sunrise : 2019-06-15T10:12:42.721Z +sunset : 2019-06-15T01:48:00.374Z moonrise : 2019-06-14T23:02:43.343Z moonset : 2019-06-15T09:12:25.416Z diff --git a/demo/nodejs/test/culminate_correct.txt b/demo/nodejs/test/culminate_correct.txt index 93a8886b..bac22fae 100644 --- a/demo/nodejs/test/culminate_correct.txt +++ b/demo/nodejs/test/culminate_correct.txt @@ -1,5 +1,5 @@ search : 2015-02-28T00:00:00.000Z -Sun : 2015-02-28T18:12:32.333Z altitude= 52.13 azimuth= 180.00 +Sun : 2015-02-28T18:12:31.044Z altitude= 52.13 azimuth= 180.00 Moon : 2015-02-28T01:58:15.630Z altitude= 77.88 azimuth= 180.00 Mercury : 2015-02-28T16:31:18.940Z altitude= 42.57 azimuth= 180.00 Venus : 2015-02-28T20:03:56.711Z altitude= 63.24 azimuth= 180.00 diff --git a/demo/nodejs/test/positions_correct.txt b/demo/nodejs/test/positions_correct.txt index b8250b6c..1197cc3e 100644 --- a/demo/nodejs/test/positions_correct.txt +++ b/demo/nodejs/test/positions_correct.txt @@ -1,7 +1,7 @@ UTC date = 2018-11-30T17:55:07.234Z BODY RA DEC AZ ALT -Sun 16.43 -21.68 180.90 22.72 +Sun 16.43 -21.68 180.91 22.72 Moon 11.30 7.81 266.73 14.07 Mercury 15.92 -18.41 189.05 25.50 Venus 13.81 -9.79 224.05 23.85 diff --git a/demo/nodejs/test/riseset_correct.txt b/demo/nodejs/test/riseset_correct.txt index 5337b85c..04982d2d 100644 --- a/demo/nodejs/test/riseset_correct.txt +++ b/demo/nodejs/test/riseset_correct.txt @@ -1,5 +1,5 @@ search : 2018-11-30T17:55:07.234Z -sunrise : 2018-12-01T13:22:52.850Z -sunset : 2018-11-30T22:21:04.222Z +sunrise : 2018-12-01T13:22:51.106Z +sunset : 2018-11-30T22:21:03.024Z moonrise : 2018-12-01T06:49:33.192Z moonset : 2018-11-30T19:22:02.717Z diff --git a/demo/python/test/culminate_correct.txt b/demo/python/test/culminate_correct.txt index 70ff537d..e76aff3a 100644 --- a/demo/python/test/culminate_correct.txt +++ b/demo/python/test/culminate_correct.txt @@ -1,5 +1,5 @@ search : 2015-02-28T00:00:00.000Z -Sun : 2015-02-28T18:12:32.334Z altitude= 52.13 azimuth= 180.00 +Sun : 2015-02-28T18:12:31.045Z altitude= 52.13 azimuth= 180.00 Moon : 2015-02-28T01:58:15.631Z altitude= 77.88 azimuth= 180.00 Mercury : 2015-02-28T16:31:18.940Z altitude= 42.57 azimuth= 180.00 Venus : 2015-02-28T20:03:56.711Z altitude= 63.24 azimuth= 180.00 diff --git a/demo/python/test/positions_correct.txt b/demo/python/test/positions_correct.txt index b8250b6c..1197cc3e 100644 --- a/demo/python/test/positions_correct.txt +++ b/demo/python/test/positions_correct.txt @@ -1,7 +1,7 @@ UTC date = 2018-11-30T17:55:07.234Z BODY RA DEC AZ ALT -Sun 16.43 -21.68 180.90 22.72 +Sun 16.43 -21.68 180.91 22.72 Moon 11.30 7.81 266.73 14.07 Mercury 15.92 -18.41 189.05 25.50 Venus 13.81 -9.79 224.05 23.85 diff --git a/demo/python/test/riseset_correct.txt b/demo/python/test/riseset_correct.txt index 2247a6ec..a03dadd5 100644 --- a/demo/python/test/riseset_correct.txt +++ b/demo/python/test/riseset_correct.txt @@ -1,5 +1,5 @@ search : 2018-11-30T17:55:07.234Z -sunrise : 2018-12-01T13:22:52.851Z -sunset : 2018-11-30T22:21:04.223Z +sunrise : 2018-12-01T13:22:51.106Z +sunset : 2018-11-30T22:21:03.024Z moonrise : 2018-12-01T06:49:33.193Z moonset : 2018-11-30T19:22:02.718Z diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 9be7b7c9..6881d6c8 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -1885,14 +1885,6 @@ astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time, astro_a vector.z = 0.0; break; - case BODY_SUN: - /* The Sun's heliocentric coordinates are always (0,0,0). No need for light travel correction. */ - vector = CalcEarth(time); - vector.x *= -1.0; - vector.y *= -1.0; - vector.z *= -1.0; - break; - case BODY_MOON: vector = Astronomy_GeoMoon(time); break; @@ -3014,12 +3006,12 @@ astro_angle_result_t Astronomy_LongitudeFromSun(astro_body_t body, astro_time_t if (body == BODY_EARTH) return AngleError(ASTRO_EARTH_NOT_ALLOWED); - sv = Astronomy_GeoVector(BODY_SUN, time, ABERRATION); + sv = Astronomy_GeoVector(BODY_SUN, time, NO_ABERRATION); se = Astronomy_Ecliptic(sv); /* checks for errors in sv */ if (se.status != ASTRO_SUCCESS) return AngleError(se.status); - bv = Astronomy_GeoVector(body, time, ABERRATION); + bv = Astronomy_GeoVector(body, time, NO_ABERRATION); be = Astronomy_Ecliptic(bv); /* checks for errors in bv */ if (be.status != ASTRO_SUCCESS) return AngleError(be.status); @@ -5568,20 +5560,21 @@ static shadow_t PlanetShadow(astro_body_t body, double planet_radius_km, astro_t { astro_vector_t e, p, g; - /* To include light travel time compensation, we use GeoVector instead of HelioVector. */ - + /* Calculate light-travel-corrected vector from Earth to planet. */ g = Astronomy_GeoVector(body, time, ABERRATION); if (g.status != ASTRO_SUCCESS) return ShadowError(g.status); - /* Calculate heliocentric Earth. */ - e = CalcEarth(time); + /* Calculate light-travel-corrected vector from Earth to Sun. */ + e = Astronomy_GeoVector(BODY_SUN, time, ABERRATION); + if (e.status != ASTRO_SUCCESS) + return ShadowError(e.status); - /* Convert light-travel corrected geocentric planet to heliocentric planet. */ + /* Deduce light-travel-corrected vector from Sun to planet. */ p.t = time; - p.x = e.x + g.x; - p.y = e.y + g.y; - p.z = e.z + g.z; + p.x = g.x - e.x; + p.y = g.y - e.y; + p.z = g.z - e.z; /* Calcluate Earth's position from the planet's point of view. */ e.x = -g.x; diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index 534366c8..e902135e 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -2617,11 +2617,6 @@ $ASTRO_IAU_DATA() /* The Earth's geocentric coordinates are always (0,0,0). */ return new AstroVector(0.0, 0.0, 0.0, time); - case Body.Sun: - /* The Sun's heliocentric coordinates are always (0,0,0). No need for light travel correction. */ - vector = CalcEarth(time); - return new AstroVector(-vector.x, -vector.y, -vector.z, time); - case Body.Moon: return GeoMoon(time); @@ -3268,10 +3263,10 @@ $ASTRO_IAU_DATA() if (body == Body.Earth) throw new EarthNotAllowedException(); - AstroVector sv = GeoVector(Body.Sun, time, Aberration.Corrected); + AstroVector sv = GeoVector(Body.Sun, time, Aberration.None); Ecliptic se = EquatorialToEcliptic(sv); - AstroVector bv = GeoVector(body, time, Aberration.Corrected); + AstroVector bv = GeoVector(body, time, Aberration.None); Ecliptic be = EquatorialToEcliptic(bv); return NormalizeLongitude(be.elon - se.elon); diff --git a/generate/template/astronomy.js b/generate/template/astronomy.js index e941428b..bb5a5fb4 100644 --- a/generate/template/astronomy.js +++ b/generate/template/astronomy.js @@ -1685,9 +1685,6 @@ Astronomy.GeoVector = function(body, date, aberration) { } geo = new Vector(h.x-earth.x, h.y-earth.y, h.z-earth.z, time); - if (body === 'Sun') { - return geo; // The Sun's heliocentric coordinates are always (0,0,0). No need to correct. - } let ltime2 = time.AddDays(-geo.Length() / C_AUDAY); dt = Math.abs(ltime2.tt - ltime.tt); if (dt < 1.0e-9) { @@ -1986,10 +1983,10 @@ Astronomy.LongitudeFromSun = function(body, date) { throw 'The Earth does not have a longitude as seen from itself.'; const t = Astronomy.MakeTime(date); - let gb = Astronomy.GeoVector(body, t, true); + let gb = Astronomy.GeoVector(body, t, false); const eb = Astronomy.Ecliptic(gb.x, gb.y, gb.z); - let gs = Astronomy.GeoVector('Sun', t, true); + let gs = Astronomy.GeoVector('Sun', t, false); const es = Astronomy.Ecliptic(gs.x, gs.y, gs.z); return NormalizeLongitude(eb.elon - es.elon); diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index ed396aa9..e5eec68f 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -1497,10 +1497,6 @@ def GeoVector(body, time, aberration): earth = _CalcEarth(ltime) geo = Vector(h.x-earth.x, h.y-earth.y, h.z-earth.z, time) - if body == Body.Sun: - # The Sun's heliocentric coordinates are always (0,0,0). No need to correct. - return geo - ltime2 = time.AddDays(-geo.Length() / _C_AUDAY) dt = abs(ltime2.tt - ltime.tt) if dt < 1.0e-9: @@ -2019,9 +2015,9 @@ def LongitudeFromSun(body, time): """ if body == Body.Earth: raise EarthNotAllowedError() - sv = GeoVector(Body.Sun, time, True) + sv = GeoVector(Body.Sun, time, False) se = Ecliptic(sv) - bv = GeoVector(body, time, True) + bv = GeoVector(body, time, False) be = Ecliptic(bv) return _NormalizeLongitude(be.elon - se.elon) diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 99932a7f..6195ef78 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -3085,14 +3085,6 @@ astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time, astro_a vector.z = 0.0; break; - case BODY_SUN: - /* The Sun's heliocentric coordinates are always (0,0,0). No need for light travel correction. */ - vector = CalcEarth(time); - vector.x *= -1.0; - vector.y *= -1.0; - vector.z *= -1.0; - break; - case BODY_MOON: vector = Astronomy_GeoMoon(time); break; @@ -4214,12 +4206,12 @@ astro_angle_result_t Astronomy_LongitudeFromSun(astro_body_t body, astro_time_t if (body == BODY_EARTH) return AngleError(ASTRO_EARTH_NOT_ALLOWED); - sv = Astronomy_GeoVector(BODY_SUN, time, ABERRATION); + sv = Astronomy_GeoVector(BODY_SUN, time, NO_ABERRATION); se = Astronomy_Ecliptic(sv); /* checks for errors in sv */ if (se.status != ASTRO_SUCCESS) return AngleError(se.status); - bv = Astronomy_GeoVector(body, time, ABERRATION); + bv = Astronomy_GeoVector(body, time, NO_ABERRATION); be = Astronomy_Ecliptic(bv); /* checks for errors in bv */ if (be.status != ASTRO_SUCCESS) return AngleError(be.status); @@ -7223,20 +7215,21 @@ static shadow_t PlanetShadow(astro_body_t body, double planet_radius_km, astro_t { astro_vector_t e, p, g; - /* To include light travel time compensation, we use GeoVector instead of HelioVector. */ - + /* Calculate light-travel-corrected vector from Earth to planet. */ g = Astronomy_GeoVector(body, time, ABERRATION); if (g.status != ASTRO_SUCCESS) return ShadowError(g.status); - /* Calculate heliocentric Earth. */ - e = CalcEarth(time); + /* Calculate light-travel-corrected vector from Earth to Sun. */ + e = Astronomy_GeoVector(BODY_SUN, time, ABERRATION); + if (e.status != ASTRO_SUCCESS) + return ShadowError(e.status); - /* Convert light-travel corrected geocentric planet to heliocentric planet. */ + /* Deduce light-travel-corrected vector from Sun to planet. */ p.t = time; - p.x = e.x + g.x; - p.y = e.y + g.y; - p.z = e.z + g.z; + p.x = g.x - e.x; + p.y = g.y - e.y; + p.z = g.z - e.z; /* Calcluate Earth's position from the planet's point of view. */ e.x = -g.x; diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index 96d87036..59d6e372 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -3821,11 +3821,6 @@ namespace CosineKitty /* The Earth's geocentric coordinates are always (0,0,0). */ return new AstroVector(0.0, 0.0, 0.0, time); - case Body.Sun: - /* The Sun's heliocentric coordinates are always (0,0,0). No need for light travel correction. */ - vector = CalcEarth(time); - return new AstroVector(-vector.x, -vector.y, -vector.z, time); - case Body.Moon: return GeoMoon(time); @@ -4472,10 +4467,10 @@ namespace CosineKitty if (body == Body.Earth) throw new EarthNotAllowedException(); - AstroVector sv = GeoVector(Body.Sun, time, Aberration.Corrected); + AstroVector sv = GeoVector(Body.Sun, time, Aberration.None); Ecliptic se = EquatorialToEcliptic(sv); - AstroVector bv = GeoVector(body, time, Aberration.Corrected); + AstroVector bv = GeoVector(body, time, Aberration.None); Ecliptic be = EquatorialToEcliptic(bv); return NormalizeLongitude(be.elon - se.elon); diff --git a/source/js/astronomy.js b/source/js/astronomy.js index 696a0500..aef60c02 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -2664,9 +2664,6 @@ Astronomy.GeoVector = function(body, date, aberration) { } geo = new Vector(h.x-earth.x, h.y-earth.y, h.z-earth.z, time); - if (body === 'Sun') { - return geo; // The Sun's heliocentric coordinates are always (0,0,0). No need to correct. - } let ltime2 = time.AddDays(-geo.Length() / C_AUDAY); dt = Math.abs(ltime2.tt - ltime.tt); if (dt < 1.0e-9) { @@ -2965,10 +2962,10 @@ Astronomy.LongitudeFromSun = function(body, date) { throw 'The Earth does not have a longitude as seen from itself.'; const t = Astronomy.MakeTime(date); - let gb = Astronomy.GeoVector(body, t, true); + let gb = Astronomy.GeoVector(body, t, false); const eb = Astronomy.Ecliptic(gb.x, gb.y, gb.z); - let gs = Astronomy.GeoVector('Sun', t, true); + let gs = Astronomy.GeoVector('Sun', t, false); const es = Astronomy.Ecliptic(gs.x, gs.y, gs.z); return NormalizeLongitude(eb.elon - es.elon); diff --git a/source/js/astronomy.min.js b/source/js/astronomy.min.js index 91a3f055..cfa4e3d6 100644 --- a/source/js/astronomy.min.js +++ b/source/js/astronomy.min.js @@ -130,10 +130,10 @@ Math.atan2(k[1],k[0])/15,0>c&&(c+=24),24<=c&&(c-=24)):c=0;p=57.29577951308232*Ma c=ka(b,c);a=e.GeoVector(a,b,k);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,ta=Math.cos(V),ua=Math.sin(V));return ma(a,b,c,ta,ua)};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*E(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 va){a:{var c=$jscomp.makeIterator(va[a]);for(a=c.next();!a.done;a=c.next()){a=a.value;var d=a.tt;var k=a.tt+a.ndays;d=(2*b.tt-(k+d))/(k-d);if(-1<=d&&1>=d){var h=[];for(k=0;3>k;++k){var l=1;var g=a.coeff[0][k];var m=d;g+=a.coeff[1][k]*m;for(c=2;cg;++g){k=e.HelioVector(a,l);c&&(d=A(t.Earth,l));h=new r(k.x-d.x,k.y-d.y,k.z-d.z,b);if("Sun"===a)return h;var m=b.AddDays(-h.Length()/173.1446326846693);k=Math.abs(m.tt-l.tt);if(1E-9>k)return h;l=m}throw"Light-travel time solver did not converge: dt="+ +a+'"';};e.HelioDistance=function(a,b){b=e.MakeTime(b);return a in t?F(t[a][2],b.tt/365250):e.HelioVector(a,b).Length()};e.GeoVector=function(a,b,c){b=e.MakeTime(b);if("Moon"===a)return e.GeoMoon(b);if("Earth"===a)return new r(0,0,0,b);var d;c||(d=A(t.Earth,b));for(var k,h,l=b,g=0;10>g;++g){k=e.HelioVector(a,l);c&&(d=A(t.Earth,l));h=new r(k.x-d.x,k.y-d.y,k.z-d.z,b);var m=b.AddDays(-h.Length()/173.1446326846693);k=Math.abs(m.tt-l.tt);if(1E-9>k)return h;l=m}throw"Light-travel time solver did not converge: dt="+ k;};e.Search=function(a,b,c,d){var k=Math.abs((d&&d.dt_tolerance_seconds||1)/86400),h=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){h=q;l=u;b=p;c=r;g=t;f=!1;continue}}}}if(0>h&&0<=g)c=n,l=g;else if(0>g&&0<=l)b=n,h=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 I(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,!0);a=e.Ecliptic(a.x,a.y,a.z);b=e.GeoVector("Sun", -b,!0);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 Ka=function(a,b,c,d,e,h,l,g){this.time=a;this.mag=b;this.phase_angle= +(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){h=q;l=u;b=p;c=r;g=t;f=!1;continue}}}}if(0>h&&0<=g)c=n,l=g;else if(0>g&&0<=l)b=n,h=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 I(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 Ka=function(a,b,c,d,e,h,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=h;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 k=new r(-d.x,-d.y,-d.z,c);b=new r(0,0,0,c);d=0}else"Moon"===a?(k=e.GeoMoon(c),b=new r(d.x+k.x,d.y+k.y,d.z+k.z,c)):(b=e.HelioVector(a,b),k=new r(b.x-d.x,b.y-d.y,b.z-d.z,c)),d=v(k,b);var h=k.Length(),l=b.Length(),g=null;if("Sun"=== a)a=Ga+5*Math.log10(h);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(h/.002573570052980638*l)}else if("Saturn"===a)a=d,g=e.Ecliptic(k.x,k.y,k.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*h),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*h)}return new Ka(c,a,d,l,h,k,b,g)};e.SearchRelativeLongitude=function(a,b,c){function d(c){var d=e.EclipticLongitude(a, diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 8a2b7af1..b561eade 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -3516,10 +3516,6 @@ def GeoVector(body, time, aberration): earth = _CalcEarth(ltime) geo = Vector(h.x-earth.x, h.y-earth.y, h.z-earth.z, time) - if body == Body.Sun: - # The Sun's heliocentric coordinates are always (0,0,0). No need to correct. - return geo - ltime2 = time.AddDays(-geo.Length() / _C_AUDAY) dt = abs(ltime2.tt - ltime.tt) if dt < 1.0e-9: @@ -4038,9 +4034,9 @@ def LongitudeFromSun(body, time): """ if body == Body.Earth: raise EarthNotAllowedError() - sv = GeoVector(Body.Sun, time, True) + sv = GeoVector(Body.Sun, time, False) se = Ecliptic(sv) - bv = GeoVector(body, time, True) + bv = GeoVector(body, time, False) be = Ecliptic(bv) return _NormalizeLongitude(be.elon - se.elon)