From af7cd4c7ea8bae1598155bb9ecccf7996afacd7a Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 16 May 2020 09:36:31 -0400 Subject: [PATCH] C LunarEclipse: Optimize by pruning full moons with large ecliptic latitudes. When the full moon's ecliptic latitude is larger than 1.8 degrees, even a penumbral eclipse is not possible. Thus there is no need to search for the minimum shadow distance in that case. This decreased unit test CalcMoon() count to 127155. Improvement ratio over original algorithm = 4.55. --- generate/ctest.c | 3 ++ generate/template/astronomy.c | 69 ++++++++++++++++++++--------------- source/c/astronomy.c | 69 ++++++++++++++++++++--------------- 3 files changed, 83 insertions(+), 58 deletions(-) diff --git a/generate/ctest.c b/generate/ctest.c index 46a108ae..32813da4 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -2709,6 +2709,9 @@ static int LunarEclipseTest(void) if (sum_diff_minutes > 0.274) FAIL("C LunarEclipseTest: EXCESSIVE AVERAGE TIME ERROR: %lf\n", sum_diff_minutes); + if (skip_count > 9) + FAIL("C LunarEclipseTest: Skipped too many penumbral eclipses: %d\n", skip_count); + printf("C LunarEclipseTest: PASS (verified %d, skipped %d, moon calcs %d, max_diff_minutes = %0.3lf, avg_diff_minutes = %0.3lf)\n", lnum, skip_count, _CalcMoonCount, max_diff_minutes, sum_diff_minutes); error = 0; fail: diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index e1e678d7..d9eb6ead 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -5640,12 +5640,14 @@ static double ShadowSemiDurationMinutes(astro_time_t center_time, double radius_ */ astro_lunar_eclipse_t Astronomy_SearchLunarEclipse(astro_time_t startTime) { + const double PruneLatitude = 1.8; /* full Moon's ecliptic latitude above which eclipse is impossible */ astro_time_t fmtime; astro_lunar_eclipse_t eclipse; astro_search_result_t fullmoon; earth_shadow_t shadow; int fmcount; double r1, r2; + double eclip_lat, eclip_lon, distance; /* Iterate through consecutive full moons until we find any kind of lunar eclipse. */ fmtime = startTime; @@ -5656,44 +5658,53 @@ astro_lunar_eclipse_t Astronomy_SearchLunarEclipse(astro_time_t startTime) if (fullmoon.status != ASTRO_SUCCESS) return EclipseError(fullmoon.status); - /* 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.time); - if (shadow.status != ASTRO_SUCCESS) - return EclipseError(shadow.status); - - r1 = fabs(shadow.r - MOON_RADIUS_KM); - r2 = fabs(shadow.r + MOON_RADIUS_KM); - if (r1 < shadow.p) + /* + Pruning: if the full Moon's ecliptic latitude is too large, + a lunar eclipse is not possible. Avoid needless work searching for + the minimum moon distance. + */ + CalcMoon(fullmoon.time.tt / 36525.0, &eclip_lon, &eclip_lat, &distance); + if (RAD2DEG * fabs(eclip_lat) < PruneLatitude) { - /* This is at least a penumbral eclipse. We will return a result. */ - eclipse.status = ASTRO_SUCCESS; - eclipse.kind = ECLIPSE_PENUMBRAL; - 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); - if (eclipse.sd_penum <= 0.0) - return EclipseError(ASTRO_SEARCH_FAILURE); + /* 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.time); + if (shadow.status != ASTRO_SUCCESS) + return EclipseError(shadow.status); - if (r1 < shadow.k) + r1 = fabs(shadow.r - MOON_RADIUS_KM); + r2 = fabs(shadow.r + MOON_RADIUS_KM); + if (r1 < shadow.p) { - /* This is at least a partial eclipse. */ - eclipse.kind = ECLIPSE_PARTIAL; - eclipse.sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM); - if (eclipse.sd_partial <= 0.0) + /* This is at least a penumbral eclipse. We will return a result. */ + eclipse.status = ASTRO_SUCCESS; + eclipse.kind = ECLIPSE_PENUMBRAL; + 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); + if (eclipse.sd_penum <= 0.0) return EclipseError(ASTRO_SEARCH_FAILURE); - if (r2 < shadow.k) + if (r1 < shadow.k) { - /* This is a total eclipse. */ - eclipse.kind = ECLIPSE_TOTAL; - eclipse.sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM); - if (eclipse.sd_total <= 0.0) + /* This is at least a partial eclipse. */ + eclipse.kind = ECLIPSE_PARTIAL; + eclipse.sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM); + if (eclipse.sd_partial <= 0.0) return EclipseError(ASTRO_SEARCH_FAILURE); + + if (r2 < shadow.k) + { + /* This is a total eclipse. */ + eclipse.kind = ECLIPSE_TOTAL; + eclipse.sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM); + if (eclipse.sd_total <= 0.0) + return EclipseError(ASTRO_SEARCH_FAILURE); + } } + return eclipse; } - return eclipse; } /* We didn't find an eclipse on this full moon, so search for the next one. */ diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 6473ebf6..42372da3 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -7295,12 +7295,14 @@ static double ShadowSemiDurationMinutes(astro_time_t center_time, double radius_ */ astro_lunar_eclipse_t Astronomy_SearchLunarEclipse(astro_time_t startTime) { + const double PruneLatitude = 1.8; /* full Moon's ecliptic latitude above which eclipse is impossible */ astro_time_t fmtime; astro_lunar_eclipse_t eclipse; astro_search_result_t fullmoon; earth_shadow_t shadow; int fmcount; double r1, r2; + double eclip_lat, eclip_lon, distance; /* Iterate through consecutive full moons until we find any kind of lunar eclipse. */ fmtime = startTime; @@ -7311,44 +7313,53 @@ astro_lunar_eclipse_t Astronomy_SearchLunarEclipse(astro_time_t startTime) if (fullmoon.status != ASTRO_SUCCESS) return EclipseError(fullmoon.status); - /* 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.time); - if (shadow.status != ASTRO_SUCCESS) - return EclipseError(shadow.status); - - r1 = fabs(shadow.r - MOON_RADIUS_KM); - r2 = fabs(shadow.r + MOON_RADIUS_KM); - if (r1 < shadow.p) + /* + Pruning: if the full Moon's ecliptic latitude is too large, + a lunar eclipse is not possible. Avoid needless work searching for + the minimum moon distance. + */ + CalcMoon(fullmoon.time.tt / 36525.0, &eclip_lon, &eclip_lat, &distance); + if (RAD2DEG * fabs(eclip_lat) < PruneLatitude) { - /* This is at least a penumbral eclipse. We will return a result. */ - eclipse.status = ASTRO_SUCCESS; - eclipse.kind = ECLIPSE_PENUMBRAL; - 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); - if (eclipse.sd_penum <= 0.0) - return EclipseError(ASTRO_SEARCH_FAILURE); + /* 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.time); + if (shadow.status != ASTRO_SUCCESS) + return EclipseError(shadow.status); - if (r1 < shadow.k) + r1 = fabs(shadow.r - MOON_RADIUS_KM); + r2 = fabs(shadow.r + MOON_RADIUS_KM); + if (r1 < shadow.p) { - /* This is at least a partial eclipse. */ - eclipse.kind = ECLIPSE_PARTIAL; - eclipse.sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM); - if (eclipse.sd_partial <= 0.0) + /* This is at least a penumbral eclipse. We will return a result. */ + eclipse.status = ASTRO_SUCCESS; + eclipse.kind = ECLIPSE_PENUMBRAL; + 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); + if (eclipse.sd_penum <= 0.0) return EclipseError(ASTRO_SEARCH_FAILURE); - if (r2 < shadow.k) + if (r1 < shadow.k) { - /* This is a total eclipse. */ - eclipse.kind = ECLIPSE_TOTAL; - eclipse.sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM); - if (eclipse.sd_total <= 0.0) + /* This is at least a partial eclipse. */ + eclipse.kind = ECLIPSE_PARTIAL; + eclipse.sd_partial = ShadowSemiDurationMinutes(shadow.time, shadow.k + MOON_RADIUS_KM); + if (eclipse.sd_partial <= 0.0) return EclipseError(ASTRO_SEARCH_FAILURE); + + if (r2 < shadow.k) + { + /* This is a total eclipse. */ + eclipse.kind = ECLIPSE_TOTAL; + eclipse.sd_total = ShadowSemiDurationMinutes(shadow.time, shadow.k - MOON_RADIUS_KM); + if (eclipse.sd_total <= 0.0) + return EclipseError(ASTRO_SEARCH_FAILURE); + } } + return eclipse; } - return eclipse; } /* We didn't find an eclipse on this full moon, so search for the next one. */