From e3255c7401aca66d87714cf600a88e073c668cf0 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 23 May 2020 13:08:25 -0400 Subject: [PATCH] Cleaned up and unified Earth and Moon radius constants. In all 4 supported languages, use consistent constant names for Earth and Moon radii. Use Moon's equatorial radius for rise/set timing. Use Moon's mean radius for calculating Moon's umbra radius for detecting solar eclipses. Also use Moon's mean radius for determining whether the Earth's shadow touches the Moon, for finding lunar eclipses. Use the Moon's polar radius for distinguishing between total and annular eclipses, with a 14 meter bias (instead of 1420 meters!) to match Espenak data. Use consistent unit test error threshold of 0.57 minutes for rise/set. Updated demo test data for slight changes to rise/set prediction times. Updated doxygen options to issue an error on any warnings. Fixed the incorrect function name link that doxygen was warning me about. --- demo/csharp/test/riseset_correct.txt | 4 +- demo/nodejs/test/riseset_correct.txt | 4 +- demo/python/test/riseset_correct.txt | 4 +- generate/rise_set_test.js | 2 +- generate/template/astronomy.c | 59 +++++++++++++++------------- generate/template/astronomy.cs | 46 ++++++++++++---------- generate/template/astronomy.js | 31 +++++++++------ generate/template/astronomy.py | 33 ++++++++++------ generate/test.py | 2 +- source/c/Doxyfile | 2 +- source/c/README.md | 2 +- source/c/astronomy.c | 59 +++++++++++++++------------- source/csharp/astronomy.cs | 46 ++++++++++++---------- source/js/astronomy.js | 31 +++++++++------ source/js/astronomy.min.js | 2 +- source/python/astronomy.py | 33 ++++++++++------ 16 files changed, 207 insertions(+), 153 deletions(-) diff --git a/demo/csharp/test/riseset_correct.txt b/demo/csharp/test/riseset_correct.txt index 882d08f5..9977e966 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 -moonrise : 2019-06-14T23:02:43.383Z -moonset : 2019-06-15T09:12:25.377Z +moonrise : 2019-06-14T23:02:43.343Z +moonset : 2019-06-15T09:12:25.416Z diff --git a/demo/nodejs/test/riseset_correct.txt b/demo/nodejs/test/riseset_correct.txt index fa6eb1e2..5337b85c 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 -moonrise : 2018-12-01T06:49:33.231Z -moonset : 2018-11-30T19:22:02.679Z +moonrise : 2018-12-01T06:49:33.192Z +moonset : 2018-11-30T19:22:02.717Z diff --git a/demo/python/test/riseset_correct.txt b/demo/python/test/riseset_correct.txt index 58ff11fc..2247a6ec 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 -moonrise : 2018-12-01T06:49:33.231Z -moonset : 2018-11-30T19:22:02.680Z +moonrise : 2018-12-01T06:49:33.193Z +moonset : 2018-11-30T19:22:02.718Z diff --git a/generate/rise_set_test.js b/generate/rise_set_test.js index d99e4010..14023205 100644 --- a/generate/rise_set_test.js +++ b/generate/rise_set_test.js @@ -91,7 +91,7 @@ function Test() { max_minutes = error_minutes; console.log(`Line ${evt.lnum} : error = ${error_minutes.toFixed(4)}`); } - if (error_minutes > 2) { + if (error_minutes > 0.57) { console.log(`Expected ${evt.date.toISOString()}`); console.log(`Found ${a_date.toString()}`); Fail("Excessive prediction time error."); diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 8b50ceb6..cb80ad7f 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -48,7 +48,6 @@ static const double ASEC2RAD = 4.848136811095359935899141e-6; static const double PI2 = 2.0 * PI; static const double ARC = 3600.0 * 180.0 / PI; /* arcseconds per radian */ static const double C_AUDAY = 173.1446326846693; /* speed of light in AU/day */ -static const double AU = 1.4959787069098932e+11; /* astronomical unit in meters */ static const double KM_PER_AU = 1.4959787069098932e+8; static const double ANGVEL = 7.2921150e-5; static const double SECONDS_PER_DAY = 24.0 * 3600.0; @@ -57,19 +56,24 @@ static const double MEAN_SYNODIC_MONTH = 29.530588; /* average number of day static const double EARTH_ORBITAL_PERIOD = 365.256; static const double NEPTUNE_ORBITAL_PERIOD = 60189.0; static const double REFRACTION_NEAR_HORIZON = 34.0 / 60.0; /* degrees of refractive "lift" seen for objects near horizon */ + static const double SUN_RADIUS_KM = 695700.0; #define SUN_RADIUS_AU (SUN_RADIUS_KM / KM_PER_AU) #define EARTH_FLATTENING 0.996647180302104 -#define EARTH_EQUATORIAL_RADIUS_M 6378136.6 -#define EARTH_POLAR_RADIUS_M (EARTH_FLATTENING * EARTH_EQUATORIAL_RADIUS_M) -#define EARTH_RADIUS_KM 6371.0 /* mean radius of the Earth's geoid, without atmosphere */ +#define EARTH_EQUATORIAL_RADIUS_KM 6378.1366 +#define EARTH_EQUATORIAL_RADIUS_AU (EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU) +#define EARTH_MEAN_RADIUS_KM 6371.0 /* mean radius of the Earth's geoid, without atmosphere */ #define EARTH_ATMOSPHERE_KM 88.0 /* effective atmosphere thickness for lunar eclipses */ -#define EARTH_ECLIPSE_RADIUS_KM (EARTH_RADIUS_KM + EARTH_ATMOSPHERE_KM) +#define EARTH_ECLIPSE_RADIUS_KM (EARTH_MEAN_RADIUS_KM + EARTH_ATMOSPHERE_KM) +/* Note: if we ever need Earth's polar radius, it is (EARTH_FLATTENING * EARTH_EQUATORIAL_RADIUS_KM) */ -static const double MOON_RADIUS_KM = 1737.4; -#define MOON_RADIUS_AU (MOON_RADIUS_KM / KM_PER_AU) -static const double ASEC180 = 180.0 * 60.0 * 60.0; /* arcseconds per 180 degrees (or pi radians) */ +#define MOON_EQUATORIAL_RADIUS_KM 1738.1 +#define MOON_MEAN_RADIUS_KM 1737.4 +#define MOON_POLAR_RADIUS_KM 1736.0 +#define MOON_EQUATORIAL_RADIUS_AU (MOON_EQUATORIAL_RADIUS_KM / KM_PER_AU) + +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. */ @@ -1119,7 +1123,6 @@ static double sidereal_time(astro_time_t *time) static void terra(astro_observer_t observer, double st, double pos[3], double vel[3]) { - double erad_km = EARTH_EQUATORIAL_RADIUS_M / 1000.0; double df2 = EARTH_FLATTENING * EARTH_FLATTENING; double phi = observer.latitude * DEG2RAD; double sinphi = sin(phi); @@ -1127,8 +1130,8 @@ static void terra(astro_observer_t observer, double st, double pos[3], double ve double c = 1.0 / sqrt(cosphi*cosphi + df2*sinphi*sinphi); double s = df2 * c; double ht_km = observer.height / 1000.0; - double ach = erad_km*c + ht_km; - double ash = erad_km*s + ht_km; + double ach = EARTH_EQUATORIAL_RADIUS_KM*c + ht_km; + double ash = EARTH_EQUATORIAL_RADIUS_KM*s + ht_km; double stlocl = (15.0*st + observer.longitude) * DEG2RAD; double sinst = sin(stlocl); double cosst = cos(stlocl); @@ -1382,7 +1385,7 @@ $ASTRO_ADDSOL() *geo_eclip_lon = PI2 * Frac((L0+DLAM/ARC) / PI2); *geo_eclip_lat = lat_seconds * (DEG2RAD / 3600.0); - *distance_au = (ARC * (EARTH_EQUATORIAL_RADIUS_M / AU)) / (0.999953253 * SINPI); + *distance_au = (ARC * EARTH_EQUATORIAL_RADIUS_AU) / (0.999953253 * SINPI); ++_CalcMoonCount; } @@ -3582,9 +3585,9 @@ astro_search_result_t Astronomy_SearchRiseSet( context.observer = observer; switch (body) { - case BODY_SUN: context.body_radius_au = SUN_RADIUS_AU; break; - case BODY_MOON: context.body_radius_au = MOON_RADIUS_AU; break; - default: context.body_radius_au = 0.0; break; + case BODY_SUN: context.body_radius_au = SUN_RADIUS_AU; break; + case BODY_MOON: context.body_radius_au = MOON_EQUATORIAL_RADIUS_AU; break; + default: context.body_radius_au = 0.0; break; } /* @@ -5578,7 +5581,7 @@ static shadow_t MoonShadow(astro_time_t time) m.y += h.y; m.z += h.z; - return CalcShadow(MOON_RADIUS_KM, time, e, m); + return CalcShadow(MOON_MEAN_RADIUS_KM, time, e, m); } @@ -5747,7 +5750,7 @@ astro_lunar_eclipse_t Astronomy_SearchLunarEclipse(astro_time_t startTime) if (shadow.status != ASTRO_SUCCESS) return LunarEclipseError(shadow.status); - if (shadow.r < shadow.p + MOON_RADIUS_KM) + if (shadow.r < shadow.p + MOON_MEAN_RADIUS_KM) { /* This is at least a penumbral eclipse. We will return a result. */ eclipse.status = ASTRO_SUCCESS; @@ -5755,23 +5758,23 @@ astro_lunar_eclipse_t Astronomy_SearchLunarEclipse(astro_time_t startTime) eclipse.center = shadow.time; eclipse.sd_total = 0.0; eclipse.sd_partial = 0.0; - eclipse.sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_RADIUS_KM, 200.0); + eclipse.sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_MEAN_RADIUS_KM, 200.0); if (eclipse.sd_penum <= 0.0) return LunarEclipseError(ASTRO_SEARCH_FAILURE); - if (shadow.r < shadow.k + MOON_RADIUS_KM) + if (shadow.r < shadow.k + MOON_MEAN_RADIUS_KM) { /* This is at least a partial eclipse. */ eclipse.kind = ECLIPSE_PARTIAL; - eclipse.sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM, eclipse.sd_penum); + eclipse.sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_MEAN_RADIUS_KM, eclipse.sd_penum); if (eclipse.sd_partial <= 0.0) return LunarEclipseError(ASTRO_SEARCH_FAILURE); - if (shadow.r + MOON_RADIUS_KM < shadow.k) + if (shadow.r + MOON_MEAN_RADIUS_KM < shadow.k) { /* This is a total eclipse. */ eclipse.kind = ECLIPSE_TOTAL; - eclipse.sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM, eclipse.sd_partial); + eclipse.sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_MEAN_RADIUS_KM, eclipse.sd_partial); if (eclipse.sd_total <= 0.0) return LunarEclipseError(ASTRO_SEARCH_FAILURE); } @@ -5877,7 +5880,7 @@ static astro_global_solar_eclipse_t GeoidIntersect(shadow_t shadow) Solve the quadratic equation that finds whether and where the shadow axis intersects with the Earth in the dilated coordinate system. */ - R = EARTH_EQUATORIAL_RADIUS_M / 1000.0; /* Earth equatorial radius in km */ + R = EARTH_EQUATORIAL_RADIUS_KM; A = v.x*v.x + v.y*v.y + v.z*v.z; B = -2.0 * (v.x*e.x + v.y*e.y + v.z*e.z); C = (e.x*e.x + e.y*e.y + e.z*e.z) - R*R; @@ -5933,7 +5936,7 @@ static astro_global_solar_eclipse_t GeoidIntersect(shadow_t shadow) o.z += shadow.target.z; /* Recalculate the shadow using a vector from the Moon's center toward the observer. */ - surface = CalcShadow(MOON_RADIUS_KM, shadow.time, o, shadow.dir); + surface = CalcShadow(MOON_POLAR_RADIUS_KM, shadow.time, o, shadow.dir); /* If we did everything right, the shadow distance should be very close to zero. */ /* That's because we already determined the observer 'o' is on the shadow axis! */ @@ -5942,8 +5945,8 @@ static astro_global_solar_eclipse_t GeoidIntersect(shadow_t shadow) /* The umbra radius tells us what kind of eclipse the observer sees. */ /* If the umbra radius is positive, this is a total eclipse. Otherwise, it's annular. */ - /* FIXFIXFIX: I added a little bias to match test data. */ - eclipse.kind = (surface.k > 1.42) ? ECLIPSE_TOTAL : ECLIPSE_ANNULAR; + /* HACK: I added a tiny bias (14 meters) to match Espenak test data. */ + eclipse.kind = (surface.k > 0.014) ? ECLIPSE_TOTAL : ECLIPSE_ANNULAR; } return eclipse; @@ -5957,7 +5960,7 @@ static astro_global_solar_eclipse_t GeoidIntersect(shadow_t shadow) * A solar eclipse found may be partial, annular, or total. * See #astro_global_solar_eclipse_t for more information. * To find a series of solar eclipses, call this function once, - * then keep calling #Astronomy_NextSolarEclipse as many times as desired, + * then keep calling #Astronomy_NextGlobalSolarEclipse as many times as desired, * passing in the `peak` value returned from the previous call. * * @param startTime @@ -6000,7 +6003,7 @@ astro_global_solar_eclipse_t Astronomy_SearchGlobalSolarEclipse(astro_time_t sta if (shadow.status != ASTRO_SUCCESS) return GlobalSolarEclipseError(shadow.status); - if (shadow.r < shadow.p + EARTH_RADIUS_KM) + if (shadow.r < shadow.p + EARTH_MEAN_RADIUS_KM) { /* This is at least a partial solar eclipse visible somewhere on Earth. */ /* Try to find an intersection between the shadow axis and the Earth's oblate geoid. */ diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index 59a8c64a..d0728986 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -1068,7 +1068,7 @@ namespace CosineKitty break; case Body.Moon: - this.body_radius_au = Astronomy.MOON_RADIUS_AU; + this.body_radius_au = Astronomy.MOON_EQUATORIAL_RADIUS_AU; break; default: @@ -1378,7 +1378,7 @@ $ASTRO_ADDSOL() return new MoonResult( Astronomy.PI2 * Frac((L0+DLAM/Astronomy.ARC) / Astronomy.PI2), lat_seconds * (Astronomy.DEG2RAD / 3600.0), - (Astronomy.ARC * (Astronomy.ERAD / Astronomy.METERS_PER_AU)) / (0.999953253 * SINPI) + (Astronomy.ARC * Astronomy.EARTH_EQUATORIAL_RADIUS_AU) / (0.999953253 * SINPI) ); } } @@ -1451,9 +1451,23 @@ $ASTRO_ADDSOL() internal const double PI2 = 2.0 * Math.PI; internal const double ARC = 3600.0 * 180.0 / Math.PI; /* arcseconds per radian */ private const double C_AUDAY = 173.1446326846693; /* speed of light in AU/day */ - internal const double ERAD = 6378136.6; /* mean earth radius in meters */ internal const double KM_PER_AU = 1.4959787069098932e+8; - internal const double METERS_PER_AU = KM_PER_AU * 1000.0; /* astronomical unit in meters */ + + internal const double SUN_RADIUS_KM = 695700.0; + internal const double SUN_RADIUS_AU = SUN_RADIUS_KM / KM_PER_AU; + + internal const double EARTH_FLATTENING = 0.996647180302104; + internal const double EARTH_EQUATORIAL_RADIUS_KM = 6378.1366; + internal const double EARTH_EQUATORIAL_RADIUS_AU = EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU; + internal const double EARTH_MEAN_RADIUS_KM = 6371.0; /* mean radius of the Earth's geoid, without atmosphere */ + internal const double EARTH_ATMOSPHERE_KM = 88.0; /* effective atmosphere thickness for lunar eclipses */ + internal const double EARTH_ECLIPSE_RADIUS_KM = EARTH_MEAN_RADIUS_KM + EARTH_ATMOSPHERE_KM; + + internal const double MOON_EQUATORIAL_RADIUS_KM = 1738.1; + internal const double MOON_MEAN_RADIUS_KM = 1737.4; + internal const double MOON_POLAR_RADIUS_KM = 1736.0; + internal const double MOON_EQUATORIAL_RADIUS_AU = (MOON_EQUATORIAL_RADIUS_KM / KM_PER_AU); + private const double ANGVEL = 7.2921150e-5; private const double SECONDS_PER_DAY = 24.0 * 3600.0; private const double SOLAR_DAYS_PER_SIDEREAL_DAY = 0.9972695717592592; @@ -1461,13 +1475,6 @@ $ASTRO_ADDSOL() private const double EARTH_ORBITAL_PERIOD = 365.256; private const double NEPTUNE_ORBITAL_PERIOD = 60189.0; internal const double REFRACTION_NEAR_HORIZON = 34.0 / 60.0; /* degrees of refractive "lift" seen for objects near horizon */ - internal const double SUN_RADIUS_KM = 695700.0; - internal const double SUN_RADIUS_AU = SUN_RADIUS_KM / KM_PER_AU; - internal const double EARTH_RADIUS_KM = 6371.0; /* mean radius of the Earth's geoid, without atmosphere */ - internal const double EARTH_ATMOSPHERE_KM = 88.0; /* effective atmosphere thickness for lunar eclipses */ - internal const double EARTH_ECLIPSE_RADIUS_KM = EARTH_RADIUS_KM + EARTH_ATMOSPHERE_KM; - internal const double MOON_RADIUS_KM = 1737.4; - internal const double MOON_RADIUS_AU = MOON_RADIUS_KM / KM_PER_AU; 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; @@ -2050,7 +2057,6 @@ $ASTRO_IAU_DATA() private static AstroVector terra(Observer observer, double st) { - double erad_km = ERAD / 1000.0; double df = 1.0 - 0.003352819697896; /* flattening of the Earth */ double df2 = df * df; double phi = observer.latitude * DEG2RAD; @@ -2059,8 +2065,8 @@ $ASTRO_IAU_DATA() double c = 1.0 / Math.Sqrt(cosphi*cosphi + df2*sinphi*sinphi); double s = df2 * c; double ht_km = observer.height / 1000.0; - double ach = erad_km*c + ht_km; - double ash = erad_km*s + ht_km; + double ach = EARTH_EQUATORIAL_RADIUS_KM*c + ht_km; + double ash = EARTH_EQUATORIAL_RADIUS_KM*s + ht_km; double stlocl = (15.0*st + observer.longitude) * DEG2RAD; double sinst = Math.Sin(stlocl); double cosst = Math.Cos(stlocl); @@ -4310,25 +4316,25 @@ $ASTRO_IAU_DATA() // is closest to the line passing through the centers of the Sun and Earth. EarthShadowInfo shadow = PeakEarthShadow(fullmoon); - if (shadow.r < shadow.p + MOON_RADIUS_KM) + if (shadow.r < shadow.p + MOON_MEAN_RADIUS_KM) { // This is at least a penumbral eclipse. We will return a result. EclipseKind kind = EclipseKind.Penumbral; double sd_total = 0.0; double sd_partial = 0.0; - double sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_RADIUS_KM, 200.0); + double sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_MEAN_RADIUS_KM, 200.0); - if (shadow.r < shadow.k + MOON_RADIUS_KM) + if (shadow.r < shadow.k + MOON_MEAN_RADIUS_KM) { // This is at least a partial eclipse. kind = EclipseKind.Partial; - sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM, sd_penum); + sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_MEAN_RADIUS_KM, sd_penum); - if (shadow.r + MOON_RADIUS_KM < shadow.k) + if (shadow.r + MOON_MEAN_RADIUS_KM < shadow.k) { // This is a total eclipse. kind = EclipseKind.Total; - sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM, sd_partial); + sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_MEAN_RADIUS_KM, sd_partial); } } return new LunarEclipseInfo(kind, shadow.time, sd_penum, sd_partial, sd_total); diff --git a/generate/template/astronomy.js b/generate/template/astronomy.js index 78a01d4b..3a78688e 100644 --- a/generate/template/astronomy.js +++ b/generate/template/astronomy.js @@ -61,13 +61,22 @@ const MEAN_SYNODIC_MONTH = 29.530588; // average number of days for Moon t const SECONDS_PER_DAY = 24 * 3600; const MILLIS_PER_DAY = SECONDS_PER_DAY * 1000; const SOLAR_DAYS_PER_SIDEREAL_DAY = 0.9972695717592592; + const SUN_RADIUS_KM = 695700.0; const SUN_RADIUS_AU = SUN_RADIUS_KM / KM_PER_AU; -const EARTH_RADIUS_KM = 6371.0; /* mean radius of the Earth's geoid, without atmosphere */ + +const EARTH_FLATTENING = 0.996647180302104; +const EARTH_EQUATORIAL_RADIUS_KM = 6378.1366; +const EARTH_EQUATORIAL_RADIUS_AU = EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU; +const EARTH_MEAN_RADIUS_KM = 6371.0; /* mean radius of the Earth's geoid, without atmosphere */ const EARTH_ATMOSPHERE_KM = 88.0; /* effective atmosphere thickness for lunar eclipses */ -const EARTH_ECLIPSE_RADIUS_KM = EARTH_RADIUS_KM + EARTH_ATMOSPHERE_KM; -const MOON_RADIUS_KM = 1737.4; -const MOON_RADIUS_AU = MOON_RADIUS_KM / KM_PER_AU; +const EARTH_ECLIPSE_RADIUS_KM = EARTH_MEAN_RADIUS_KM + EARTH_ATMOSPHERE_KM; + +const MOON_EQUATORIAL_RADIUS_KM = 1738.1; +const MOON_MEAN_RADIUS_KM = 1737.4; +const MOON_POLAR_RADIUS_KM = 1736.0; +const MOON_EQUATORIAL_RADIUS_AU = (MOON_EQUATORIAL_RADIUS_KM / KM_PER_AU); + 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. */ @@ -2516,7 +2525,7 @@ Astronomy.NextMoonQuarter = function(mq) { Astronomy.SearchRiseSet = function(body, observer, direction, dateStart, limitDays) { // We calculate the apparent angular radius of the Sun and Moon, // but treat all other bodies as points. - let body_radius_au = { Sun:SUN_RADIUS_AU, Moon:MOON_RADIUS_AU }[body] || 0; + let body_radius_au = { Sun:SUN_RADIUS_AU, Moon:MOON_EQUATORIAL_RADIUS_AU }[body] || 0; function peak_altitude(t) { // Return the angular altitude above or below the horizon @@ -4353,22 +4362,22 @@ Astronomy.SearchLunarEclipse = function(date) { /* Search near the full moon for the time when the center of the Moon */ /* is closest to the line passing through the centers of the Sun and Earth. */ const shadow = PeakEarthShadow(fullmoon); - if (shadow.r < shadow.p + MOON_RADIUS_KM) { + if (shadow.r < shadow.p + MOON_MEAN_RADIUS_KM) { /* This is at least a penumbral eclipse. We will return a result. */ let kind = 'penumbral'; let sd_total = 0.0; let sd_partial = 0.0; - let sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_RADIUS_KM, 200.0); + let sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_MEAN_RADIUS_KM, 200.0); - if (shadow.r < shadow.k + MOON_RADIUS_KM) { + if (shadow.r < shadow.k + MOON_MEAN_RADIUS_KM) { /* This is at least a partial eclipse. */ kind = 'partial'; - sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM, sd_penum); + sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_MEAN_RADIUS_KM, sd_penum); - if (shadow.r + MOON_RADIUS_KM < shadow.k) { + if (shadow.r + MOON_MEAN_RADIUS_KM < shadow.k) { /* This is a total eclipse. */ kind = 'total'; - sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM, sd_partial); + sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_MEAN_RADIUS_KM, sd_partial); } } return new LunarEclipseInfo(kind, shadow.time, sd_penum, sd_partial, sd_total); diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 1ac32939..136af4fe 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -55,13 +55,22 @@ _MEAN_SYNODIC_MONTH = 29.530588 _EARTH_ORBITAL_PERIOD = 365.256 _NEPTUNE_ORBITAL_PERIOD = 60189.0 _REFRACTION_NEAR_HORIZON = 34.0 / 60.0 + _SUN_RADIUS_KM = 695700.0 _SUN_RADIUS_AU = _SUN_RADIUS_KM / _KM_PER_AU -_EARTH_RADIUS_KM = 6371.0 # mean radius of the Earth's geoid, without atmosphere -_EARTH_ATMOSPHERE_KM = 88.0 # effective atmosphere thickness for lunar eclipses -_EARTH_ECLIPSE_RADIUS_KM = _EARTH_RADIUS_KM + _EARTH_ATMOSPHERE_KM -_MOON_RADIUS_KM = 1737.4 -_MOON_RADIUS_AU = _MOON_RADIUS_KM / _KM_PER_AU + +_EARTH_FLATTENING = 0.996647180302104 +_EARTH_EQUATORIAL_RADIUS_KM = 6378.1366 +_EARTH_EQUATORIAL_RADIUS_AU = _EARTH_EQUATORIAL_RADIUS_KM / _KM_PER_AU +_EARTH_MEAN_RADIUS_KM = 6371.0 # mean radius of the Earth's geoid, without atmosphere +_EARTH_ATMOSPHERE_KM = 88.0 # effective atmosphere thickness for lunar eclipses +_EARTH_ECLIPSE_RADIUS_KM = _EARTH_MEAN_RADIUS_KM + _EARTH_ATMOSPHERE_KM + +_MOON_EQUATORIAL_RADIUS_KM = 1738.1 +_MOON_MEAN_RADIUS_KM = 1737.4 +_MOON_POLAR_RADIUS_KM = 1736.0 +_MOON_EQUATORIAL_RADIUS_AU = (_MOON_EQUATORIAL_RADIUS_KM / _KM_PER_AU) + _ASEC180 = 180.0 * 60.0 * 60.0 _AU_PER_PARSEC = _ASEC180 / math.pi _EARTH_MOON_MASS_RATIO = 81.30056 @@ -2962,7 +2971,7 @@ def SearchRiseSet(body, observer, direction, startTime, limitDays): elif body == Body.Sun: body_radius = _SUN_RADIUS_AU elif body == Body.Moon: - body_radius = _MOON_RADIUS_AU + body_radius = _MOON_EQUATORIAL_RADIUS_AU else: body_radius = 0.0 @@ -4256,22 +4265,22 @@ def SearchLunarEclipse(startTime): # Search near the full moon for the time when the center of the Moon # is closest to the line passing through the centers of the Sun and Earth. shadow = _PeakEarthShadow(fullmoon) - if shadow.r < shadow.p + _MOON_RADIUS_KM: + if shadow.r < shadow.p + _MOON_MEAN_RADIUS_KM: # This is at least a penumbral eclipse. We will return a result. kind = EclipseKind.Penumbral sd_total = 0.0 sd_partial = 0.0 - sd_penum = _ShadowSemiDurationMinutes(shadow.time, shadow.p + _MOON_RADIUS_KM, 200.0) + sd_penum = _ShadowSemiDurationMinutes(shadow.time, shadow.p + _MOON_MEAN_RADIUS_KM, 200.0) - if shadow.r < shadow.k + _MOON_RADIUS_KM: + if shadow.r < shadow.k + _MOON_MEAN_RADIUS_KM: # This is at least a partial eclipse. kind = EclipseKind.Partial - sd_partial = _ShadowSemiDurationMinutes(shadow.time, shadow.k + _MOON_RADIUS_KM, sd_penum) + sd_partial = _ShadowSemiDurationMinutes(shadow.time, shadow.k + _MOON_MEAN_RADIUS_KM, sd_penum) - if shadow.r + _MOON_RADIUS_KM < shadow.k: + if shadow.r + _MOON_MEAN_RADIUS_KM < shadow.k: # This is a total eclipse. kind = EclipseKind.Total - sd_total = _ShadowSemiDurationMinutes(shadow.time, shadow.k - _MOON_RADIUS_KM, sd_partial) + sd_total = _ShadowSemiDurationMinutes(shadow.time, shadow.k - _MOON_MEAN_RADIUS_KM, sd_partial) return LunarEclipseInfo(kind, shadow.time, sd_penum, sd_partial, sd_total) diff --git a/generate/test.py b/generate/test.py index 0a165085..3e94a8d6 100755 --- a/generate/test.py +++ b/generate/test.py @@ -690,7 +690,7 @@ def Test_RiseSet(filename): error_minutes = (24.0 * 60.0) * abs(a_evt.tt - correct_time.tt) sum_minutes += error_minutes ** 2 max_minutes = max(max_minutes, error_minutes) - if error_minutes > 0.565: + if error_minutes > 0.57: print('PY Test_RiseSet({} line {}): excessive prediction time error = {} minutes.'.format(filename, lnum, error_minutes)) print(' correct = {}, calculated = {}'.format(correct_time, a_evt)) return 1 diff --git a/source/c/Doxyfile b/source/c/Doxyfile index 23ec5805..251c5b33 100644 --- a/source/c/Doxyfile +++ b/source/c/Doxyfile @@ -762,7 +762,7 @@ WARN_NO_PARAMDOC = NO # a warning is encountered. # The default value is: NO. -WARN_AS_ERROR = NO +WARN_AS_ERROR = YES # The WARN_FORMAT tag determines the format of the warning messages that doxygen # can produce. The string should contain the $file, $line, and $text tags, which diff --git a/source/c/README.md b/source/c/README.md index 15a25dff..b689531f 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -1387,7 +1387,7 @@ If the search does not converge within 20 iterations, it will fail with status c -This function finds the first solar eclipse that occurs after `startTime`. A solar eclipse found may be partial, annular, or total. See [`astro_global_solar_eclipse_t`](#astro_global_solar_eclipse_t) for more information. To find a series of solar eclipses, call this function once, then keep calling #Astronomy_NextSolarEclipse as many times as desired, passing in the `peak` value returned from the previous call. +This function finds the first solar eclipse that occurs after `startTime`. A solar eclipse found may be partial, annular, or total. See [`astro_global_solar_eclipse_t`](#astro_global_solar_eclipse_t) for more information. To find a series of solar eclipses, call this function once, then keep calling [`Astronomy_NextGlobalSolarEclipse`](#Astronomy_NextGlobalSolarEclipse) as many times as desired, passing in the `peak` value returned from the previous call. diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 53f32776..0a314b1d 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -48,7 +48,6 @@ static const double ASEC2RAD = 4.848136811095359935899141e-6; static const double PI2 = 2.0 * PI; static const double ARC = 3600.0 * 180.0 / PI; /* arcseconds per radian */ static const double C_AUDAY = 173.1446326846693; /* speed of light in AU/day */ -static const double AU = 1.4959787069098932e+11; /* astronomical unit in meters */ static const double KM_PER_AU = 1.4959787069098932e+8; static const double ANGVEL = 7.2921150e-5; static const double SECONDS_PER_DAY = 24.0 * 3600.0; @@ -57,19 +56,24 @@ static const double MEAN_SYNODIC_MONTH = 29.530588; /* average number of day static const double EARTH_ORBITAL_PERIOD = 365.256; static const double NEPTUNE_ORBITAL_PERIOD = 60189.0; static const double REFRACTION_NEAR_HORIZON = 34.0 / 60.0; /* degrees of refractive "lift" seen for objects near horizon */ + static const double SUN_RADIUS_KM = 695700.0; #define SUN_RADIUS_AU (SUN_RADIUS_KM / KM_PER_AU) #define EARTH_FLATTENING 0.996647180302104 -#define EARTH_EQUATORIAL_RADIUS_M 6378136.6 -#define EARTH_POLAR_RADIUS_M (EARTH_FLATTENING * EARTH_EQUATORIAL_RADIUS_M) -#define EARTH_RADIUS_KM 6371.0 /* mean radius of the Earth's geoid, without atmosphere */ +#define EARTH_EQUATORIAL_RADIUS_KM 6378.1366 +#define EARTH_EQUATORIAL_RADIUS_AU (EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU) +#define EARTH_MEAN_RADIUS_KM 6371.0 /* mean radius of the Earth's geoid, without atmosphere */ #define EARTH_ATMOSPHERE_KM 88.0 /* effective atmosphere thickness for lunar eclipses */ -#define EARTH_ECLIPSE_RADIUS_KM (EARTH_RADIUS_KM + EARTH_ATMOSPHERE_KM) +#define EARTH_ECLIPSE_RADIUS_KM (EARTH_MEAN_RADIUS_KM + EARTH_ATMOSPHERE_KM) +/* Note: if we ever need Earth's polar radius, it is (EARTH_FLATTENING * EARTH_EQUATORIAL_RADIUS_KM) */ -static const double MOON_RADIUS_KM = 1737.4; -#define MOON_RADIUS_AU (MOON_RADIUS_KM / KM_PER_AU) -static const double ASEC180 = 180.0 * 60.0 * 60.0; /* arcseconds per 180 degrees (or pi radians) */ +#define MOON_EQUATORIAL_RADIUS_KM 1738.1 +#define MOON_MEAN_RADIUS_KM 1737.4 +#define MOON_POLAR_RADIUS_KM 1736.0 +#define MOON_EQUATORIAL_RADIUS_AU (MOON_EQUATORIAL_RADIUS_KM / KM_PER_AU) + +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. */ @@ -1197,7 +1201,6 @@ static double sidereal_time(astro_time_t *time) static void terra(astro_observer_t observer, double st, double pos[3], double vel[3]) { - double erad_km = EARTH_EQUATORIAL_RADIUS_M / 1000.0; double df2 = EARTH_FLATTENING * EARTH_FLATTENING; double phi = observer.latitude * DEG2RAD; double sinphi = sin(phi); @@ -1205,8 +1208,8 @@ static void terra(astro_observer_t observer, double st, double pos[3], double ve double c = 1.0 / sqrt(cosphi*cosphi + df2*sinphi*sinphi); double s = df2 * c; double ht_km = observer.height / 1000.0; - double ach = erad_km*c + ht_km; - double ash = erad_km*s + ht_km; + double ach = EARTH_EQUATORIAL_RADIUS_KM*c + ht_km; + double ash = EARTH_EQUATORIAL_RADIUS_KM*s + ht_km; double stlocl = (15.0*st + observer.longitude) * DEG2RAD; double sinst = sin(stlocl); double cosst = cos(stlocl); @@ -1565,7 +1568,7 @@ static void CalcMoon( *geo_eclip_lon = PI2 * Frac((L0+DLAM/ARC) / PI2); *geo_eclip_lat = lat_seconds * (DEG2RAD / 3600.0); - *distance_au = (ARC * (EARTH_EQUATORIAL_RADIUS_M / AU)) / (0.999953253 * SINPI); + *distance_au = (ARC * EARTH_EQUATORIAL_RADIUS_AU) / (0.999953253 * SINPI); ++_CalcMoonCount; } @@ -4782,9 +4785,9 @@ astro_search_result_t Astronomy_SearchRiseSet( context.observer = observer; switch (body) { - case BODY_SUN: context.body_radius_au = SUN_RADIUS_AU; break; - case BODY_MOON: context.body_radius_au = MOON_RADIUS_AU; break; - default: context.body_radius_au = 0.0; break; + case BODY_SUN: context.body_radius_au = SUN_RADIUS_AU; break; + case BODY_MOON: context.body_radius_au = MOON_EQUATORIAL_RADIUS_AU; break; + default: context.body_radius_au = 0.0; break; } /* @@ -7233,7 +7236,7 @@ static shadow_t MoonShadow(astro_time_t time) m.y += h.y; m.z += h.z; - return CalcShadow(MOON_RADIUS_KM, time, e, m); + return CalcShadow(MOON_MEAN_RADIUS_KM, time, e, m); } @@ -7402,7 +7405,7 @@ astro_lunar_eclipse_t Astronomy_SearchLunarEclipse(astro_time_t startTime) if (shadow.status != ASTRO_SUCCESS) return LunarEclipseError(shadow.status); - if (shadow.r < shadow.p + MOON_RADIUS_KM) + if (shadow.r < shadow.p + MOON_MEAN_RADIUS_KM) { /* This is at least a penumbral eclipse. We will return a result. */ eclipse.status = ASTRO_SUCCESS; @@ -7410,23 +7413,23 @@ astro_lunar_eclipse_t Astronomy_SearchLunarEclipse(astro_time_t startTime) eclipse.center = shadow.time; eclipse.sd_total = 0.0; eclipse.sd_partial = 0.0; - eclipse.sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_RADIUS_KM, 200.0); + eclipse.sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_MEAN_RADIUS_KM, 200.0); if (eclipse.sd_penum <= 0.0) return LunarEclipseError(ASTRO_SEARCH_FAILURE); - if (shadow.r < shadow.k + MOON_RADIUS_KM) + if (shadow.r < shadow.k + MOON_MEAN_RADIUS_KM) { /* This is at least a partial eclipse. */ eclipse.kind = ECLIPSE_PARTIAL; - eclipse.sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM, eclipse.sd_penum); + eclipse.sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_MEAN_RADIUS_KM, eclipse.sd_penum); if (eclipse.sd_partial <= 0.0) return LunarEclipseError(ASTRO_SEARCH_FAILURE); - if (shadow.r + MOON_RADIUS_KM < shadow.k) + if (shadow.r + MOON_MEAN_RADIUS_KM < shadow.k) { /* This is a total eclipse. */ eclipse.kind = ECLIPSE_TOTAL; - eclipse.sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM, eclipse.sd_partial); + eclipse.sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_MEAN_RADIUS_KM, eclipse.sd_partial); if (eclipse.sd_total <= 0.0) return LunarEclipseError(ASTRO_SEARCH_FAILURE); } @@ -7532,7 +7535,7 @@ static astro_global_solar_eclipse_t GeoidIntersect(shadow_t shadow) Solve the quadratic equation that finds whether and where the shadow axis intersects with the Earth in the dilated coordinate system. */ - R = EARTH_EQUATORIAL_RADIUS_M / 1000.0; /* Earth equatorial radius in km */ + R = EARTH_EQUATORIAL_RADIUS_KM; A = v.x*v.x + v.y*v.y + v.z*v.z; B = -2.0 * (v.x*e.x + v.y*e.y + v.z*e.z); C = (e.x*e.x + e.y*e.y + e.z*e.z) - R*R; @@ -7588,7 +7591,7 @@ static astro_global_solar_eclipse_t GeoidIntersect(shadow_t shadow) o.z += shadow.target.z; /* Recalculate the shadow using a vector from the Moon's center toward the observer. */ - surface = CalcShadow(MOON_RADIUS_KM, shadow.time, o, shadow.dir); + surface = CalcShadow(MOON_POLAR_RADIUS_KM, shadow.time, o, shadow.dir); /* If we did everything right, the shadow distance should be very close to zero. */ /* That's because we already determined the observer 'o' is on the shadow axis! */ @@ -7597,8 +7600,8 @@ static astro_global_solar_eclipse_t GeoidIntersect(shadow_t shadow) /* The umbra radius tells us what kind of eclipse the observer sees. */ /* If the umbra radius is positive, this is a total eclipse. Otherwise, it's annular. */ - /* FIXFIXFIX: I added a little bias to match test data. */ - eclipse.kind = (surface.k > 1.42) ? ECLIPSE_TOTAL : ECLIPSE_ANNULAR; + /* HACK: I added a tiny bias (14 meters) to match Espenak test data. */ + eclipse.kind = (surface.k > 0.014) ? ECLIPSE_TOTAL : ECLIPSE_ANNULAR; } return eclipse; @@ -7612,7 +7615,7 @@ static astro_global_solar_eclipse_t GeoidIntersect(shadow_t shadow) * A solar eclipse found may be partial, annular, or total. * See #astro_global_solar_eclipse_t for more information. * To find a series of solar eclipses, call this function once, - * then keep calling #Astronomy_NextSolarEclipse as many times as desired, + * then keep calling #Astronomy_NextGlobalSolarEclipse as many times as desired, * passing in the `peak` value returned from the previous call. * * @param startTime @@ -7655,7 +7658,7 @@ astro_global_solar_eclipse_t Astronomy_SearchGlobalSolarEclipse(astro_time_t sta if (shadow.status != ASTRO_SUCCESS) return GlobalSolarEclipseError(shadow.status); - if (shadow.r < shadow.p + EARTH_RADIUS_KM) + if (shadow.r < shadow.p + EARTH_MEAN_RADIUS_KM) { /* This is at least a partial solar eclipse visible somewhere on Earth. */ /* Try to find an intersection between the shadow axis and the Earth's oblate geoid. */ diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index 61582b2c..311686ad 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -1068,7 +1068,7 @@ namespace CosineKitty break; case Body.Moon: - this.body_radius_au = Astronomy.MOON_RADIUS_AU; + this.body_radius_au = Astronomy.MOON_EQUATORIAL_RADIUS_AU; break; default: @@ -1483,7 +1483,7 @@ namespace CosineKitty return new MoonResult( Astronomy.PI2 * Frac((L0+DLAM/Astronomy.ARC) / Astronomy.PI2), lat_seconds * (Astronomy.DEG2RAD / 3600.0), - (Astronomy.ARC * (Astronomy.ERAD / Astronomy.METERS_PER_AU)) / (0.999953253 * SINPI) + (Astronomy.ARC * Astronomy.EARTH_EQUATORIAL_RADIUS_AU) / (0.999953253 * SINPI) ); } } @@ -1556,9 +1556,23 @@ namespace CosineKitty internal const double PI2 = 2.0 * Math.PI; internal const double ARC = 3600.0 * 180.0 / Math.PI; /* arcseconds per radian */ private const double C_AUDAY = 173.1446326846693; /* speed of light in AU/day */ - internal const double ERAD = 6378136.6; /* mean earth radius in meters */ internal const double KM_PER_AU = 1.4959787069098932e+8; - internal const double METERS_PER_AU = KM_PER_AU * 1000.0; /* astronomical unit in meters */ + + internal const double SUN_RADIUS_KM = 695700.0; + internal const double SUN_RADIUS_AU = SUN_RADIUS_KM / KM_PER_AU; + + internal const double EARTH_FLATTENING = 0.996647180302104; + internal const double EARTH_EQUATORIAL_RADIUS_KM = 6378.1366; + internal const double EARTH_EQUATORIAL_RADIUS_AU = EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU; + internal const double EARTH_MEAN_RADIUS_KM = 6371.0; /* mean radius of the Earth's geoid, without atmosphere */ + internal const double EARTH_ATMOSPHERE_KM = 88.0; /* effective atmosphere thickness for lunar eclipses */ + internal const double EARTH_ECLIPSE_RADIUS_KM = EARTH_MEAN_RADIUS_KM + EARTH_ATMOSPHERE_KM; + + internal const double MOON_EQUATORIAL_RADIUS_KM = 1738.1; + internal const double MOON_MEAN_RADIUS_KM = 1737.4; + internal const double MOON_POLAR_RADIUS_KM = 1736.0; + internal const double MOON_EQUATORIAL_RADIUS_AU = (MOON_EQUATORIAL_RADIUS_KM / KM_PER_AU); + private const double ANGVEL = 7.2921150e-5; private const double SECONDS_PER_DAY = 24.0 * 3600.0; private const double SOLAR_DAYS_PER_SIDEREAL_DAY = 0.9972695717592592; @@ -1566,13 +1580,6 @@ namespace CosineKitty private const double EARTH_ORBITAL_PERIOD = 365.256; private const double NEPTUNE_ORBITAL_PERIOD = 60189.0; internal const double REFRACTION_NEAR_HORIZON = 34.0 / 60.0; /* degrees of refractive "lift" seen for objects near horizon */ - internal const double SUN_RADIUS_KM = 695700.0; - internal const double SUN_RADIUS_AU = SUN_RADIUS_KM / KM_PER_AU; - internal const double EARTH_RADIUS_KM = 6371.0; /* mean radius of the Earth's geoid, without atmosphere */ - internal const double EARTH_ATMOSPHERE_KM = 88.0; /* effective atmosphere thickness for lunar eclipses */ - internal const double EARTH_ECLIPSE_RADIUS_KM = EARTH_RADIUS_KM + EARTH_ATMOSPHERE_KM; - internal const double MOON_RADIUS_KM = 1737.4; - internal const double MOON_RADIUS_AU = MOON_RADIUS_KM / KM_PER_AU; 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; @@ -3254,7 +3261,6 @@ namespace CosineKitty private static AstroVector terra(Observer observer, double st) { - double erad_km = ERAD / 1000.0; double df = 1.0 - 0.003352819697896; /* flattening of the Earth */ double df2 = df * df; double phi = observer.latitude * DEG2RAD; @@ -3263,8 +3269,8 @@ namespace CosineKitty double c = 1.0 / Math.Sqrt(cosphi*cosphi + df2*sinphi*sinphi); double s = df2 * c; double ht_km = observer.height / 1000.0; - double ach = erad_km*c + ht_km; - double ash = erad_km*s + ht_km; + double ach = EARTH_EQUATORIAL_RADIUS_KM*c + ht_km; + double ash = EARTH_EQUATORIAL_RADIUS_KM*s + ht_km; double stlocl = (15.0*st + observer.longitude) * DEG2RAD; double sinst = Math.Sin(stlocl); double cosst = Math.Cos(stlocl); @@ -5514,25 +5520,25 @@ namespace CosineKitty // is closest to the line passing through the centers of the Sun and Earth. EarthShadowInfo shadow = PeakEarthShadow(fullmoon); - if (shadow.r < shadow.p + MOON_RADIUS_KM) + if (shadow.r < shadow.p + MOON_MEAN_RADIUS_KM) { // This is at least a penumbral eclipse. We will return a result. EclipseKind kind = EclipseKind.Penumbral; double sd_total = 0.0; double sd_partial = 0.0; - double sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_RADIUS_KM, 200.0); + double sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_MEAN_RADIUS_KM, 200.0); - if (shadow.r < shadow.k + MOON_RADIUS_KM) + if (shadow.r < shadow.k + MOON_MEAN_RADIUS_KM) { // This is at least a partial eclipse. kind = EclipseKind.Partial; - sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM, sd_penum); + sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_MEAN_RADIUS_KM, sd_penum); - if (shadow.r + MOON_RADIUS_KM < shadow.k) + if (shadow.r + MOON_MEAN_RADIUS_KM < shadow.k) { // This is a total eclipse. kind = EclipseKind.Total; - sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM, sd_partial); + sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_MEAN_RADIUS_KM, sd_partial); } } return new LunarEclipseInfo(kind, shadow.time, sd_penum, sd_partial, sd_total); diff --git a/source/js/astronomy.js b/source/js/astronomy.js index 16a3df84..96e10d15 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -61,13 +61,22 @@ const MEAN_SYNODIC_MONTH = 29.530588; // average number of days for Moon t const SECONDS_PER_DAY = 24 * 3600; const MILLIS_PER_DAY = SECONDS_PER_DAY * 1000; const SOLAR_DAYS_PER_SIDEREAL_DAY = 0.9972695717592592; + const SUN_RADIUS_KM = 695700.0; const SUN_RADIUS_AU = SUN_RADIUS_KM / KM_PER_AU; -const EARTH_RADIUS_KM = 6371.0; /* mean radius of the Earth's geoid, without atmosphere */ + +const EARTH_FLATTENING = 0.996647180302104; +const EARTH_EQUATORIAL_RADIUS_KM = 6378.1366; +const EARTH_EQUATORIAL_RADIUS_AU = EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU; +const EARTH_MEAN_RADIUS_KM = 6371.0; /* mean radius of the Earth's geoid, without atmosphere */ const EARTH_ATMOSPHERE_KM = 88.0; /* effective atmosphere thickness for lunar eclipses */ -const EARTH_ECLIPSE_RADIUS_KM = EARTH_RADIUS_KM + EARTH_ATMOSPHERE_KM; -const MOON_RADIUS_KM = 1737.4; -const MOON_RADIUS_AU = MOON_RADIUS_KM / KM_PER_AU; +const EARTH_ECLIPSE_RADIUS_KM = EARTH_MEAN_RADIUS_KM + EARTH_ATMOSPHERE_KM; + +const MOON_EQUATORIAL_RADIUS_KM = 1738.1; +const MOON_MEAN_RADIUS_KM = 1737.4; +const MOON_POLAR_RADIUS_KM = 1736.0; +const MOON_EQUATORIAL_RADIUS_AU = (MOON_EQUATORIAL_RADIUS_KM / KM_PER_AU); + 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. */ @@ -3495,7 +3504,7 @@ Astronomy.NextMoonQuarter = function(mq) { Astronomy.SearchRiseSet = function(body, observer, direction, dateStart, limitDays) { // We calculate the apparent angular radius of the Sun and Moon, // but treat all other bodies as points. - let body_radius_au = { Sun:SUN_RADIUS_AU, Moon:MOON_RADIUS_AU }[body] || 0; + let body_radius_au = { Sun:SUN_RADIUS_AU, Moon:MOON_EQUATORIAL_RADIUS_AU }[body] || 0; function peak_altitude(t) { // Return the angular altitude above or below the horizon @@ -5783,22 +5792,22 @@ Astronomy.SearchLunarEclipse = function(date) { /* Search near the full moon for the time when the center of the Moon */ /* is closest to the line passing through the centers of the Sun and Earth. */ const shadow = PeakEarthShadow(fullmoon); - if (shadow.r < shadow.p + MOON_RADIUS_KM) { + if (shadow.r < shadow.p + MOON_MEAN_RADIUS_KM) { /* This is at least a penumbral eclipse. We will return a result. */ let kind = 'penumbral'; let sd_total = 0.0; let sd_partial = 0.0; - let sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_RADIUS_KM, 200.0); + let sd_penum = ShadowSemiDurationMinutes(shadow.time, shadow.p + MOON_MEAN_RADIUS_KM, 200.0); - if (shadow.r < shadow.k + MOON_RADIUS_KM) { + if (shadow.r < shadow.k + MOON_MEAN_RADIUS_KM) { /* This is at least a partial eclipse. */ kind = 'partial'; - sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM, sd_penum); + sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_MEAN_RADIUS_KM, sd_penum); - if (shadow.r + MOON_RADIUS_KM < shadow.k) { + if (shadow.r + MOON_MEAN_RADIUS_KM < shadow.k) { /* This is a total eclipse. */ kind = 'total'; - sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM, sd_partial); + sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_MEAN_RADIUS_KM, sd_partial); } } return new LunarEclipseInfo(kind, shadow.time, sd_penum, sd_partial, sd_total); diff --git a/source/js/astronomy.min.js b/source/js/astronomy.min.js index 416db01c..c9f2725b 100644 --- a/source/js/astronomy.min.js +++ b/source/js/astronomy.min.js @@ -137,7 +137,7 @@ c))/2;this.helio_dist=e;this.geo_dist=d;this.gc=h;this.hc=k;this.ring_tilt=l};d. 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(k*h)}return new ya(c,a,e,k,h,g,b,l)};d.SearchRelativeLongitude=function(a,b,c){function e(c){var e=d.EclipticLongitude(a,c);c=d.EclipticLongitude("Earth", c);return F(h*(c-e)-b)}var g=z[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 h=g.OrbitalPeriod>z.Earth.OrbitalPeriod?1:-1;g=I(a);c=d.MakeTime(c);var k=e(c);0l;++l){var m=-k/360*g;c=c.AddDays(m);if(1>86400*Math.abs(m))return c;m=k;k=e(c);30>Math.abs(m)&&m!==k&&(m/=m-k,.5m&&(g*=m))}throw"Relative longitude search failed to converge for "+a+ " near "+c.toString()+" (error_angle = "+k+").";};d.MoonPhase=function(a){return d.LongitudeFromSun("Moon",a)};d.SearchMoonPhase=function(a,b,c){function e(b){b=d.MoonPhase(b);return F(b-a)}b=d.MakeTime(b);var g=e(b);0c)return null;c=Math.min(c,h+.9);g=b.AddDays(g);b=b.AddDays(c);return d.Search(e,g,b)};var za=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 za(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 h(e){var f=d.Equator(a,e,b,!0,!0);e=d.Horizon(e,b,f.ra,f.dec).altitude+k/f.dist*57.29577951308232+va;return c*e}var k={Sun:.0046504672612422675,Moon:1.1613801666928728E-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 "+ +b,a,10))&&new za(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 h(e){var f=d.Equator(a,e,b,!0,!0);e=d.Horizon(e,b,f.ra,f.dec).altitude+k/f.dist*57.29577951308232+va;return c*e}var k={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 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=h(e);var n;if(0=f&&0=e.ut+g)return null;p=f.time;f=h(f.time);n=h(q.time)}};var Aa=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 h=Q(e),k=d.Equator(a,e,b,!0,!0);h=(c+k.ra-b.longitude/15-h)%24;1===g?0>h&&(h+=24):-12>h?h+=24:123600*Math.abs(h))return a=d.Horizon(e,b,k.ra,k.dec,"normal"),new Aa(e,a);e=e.AddDays(h/24*.9972695717592592)}};var Ba=function(a,b,c,e){this.mar_equinox=a;this.jun_solstice=b;this.sep_equinox=c;this.dec_solstice=e};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),h=b(270,12,20);return new Ba(c,e,g,h)};var Ca=function(a,b,c,e){this.time=a;this.visibility=b;this.elongation=c;this.ecliptic_separation=e};d.Elongation=function(a,b){b=d.MakeTime(b);var c=d.LongitudeFromSun(a,b);if(180