From 74dd133391f5ef8118ca5ea54d79abcbe3a5ed0b Mon Sep 17 00:00:00 2001 From: Don Cross Date: Fri, 24 May 2019 16:58:03 -0400 Subject: [PATCH] Implemented C function SearchLunarApsis, but not yet tested. --- generate/template/astronomy.c | 110 ++++++++++++++++++++++++++++++++- generate/template/astronomy.js | 2 +- source/c/astronomy.c | 110 ++++++++++++++++++++++++++++++++- source/c/astronomy.h | 18 ++++++ source/js/astronomy.js | 2 +- 5 files changed, 238 insertions(+), 4 deletions(-) diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 33188da2..66bd58fa 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -264,6 +264,18 @@ static astro_illum_t IllumError(astro_status_t status) return result; } +static astro_apsis_t ApsisError(astro_status_t status) +{ + astro_apsis_t result; + + result.status = status; + result.time.tt = result.time.ut = NAN; + result.kind = APSIS_INVALID; + result.dist_km = result.dist_au = NAN; + + return result; +} + static astro_func_result_t SynodicPeriod(astro_body_t body) { static const double Te = 365.256; /* Earth's orbital period in days */ @@ -3038,6 +3050,102 @@ astro_illum_t Astronomy_SearchPeakMagnitude(astro_body_t body, astro_time_t star return IllumError(ASTRO_SEARCH_FAILURE); } +static double MoonDistance(astro_time_t t) +{ + double lon, lat, dist; + CalcMoon(t.tt / 36525.0, &lon, &lat, &dist); + return dist; +} + +static astro_func_result_t distance_slope(void *context, astro_time_t time) +{ + static const double dt = 0.001; + astro_time_t t1 = Astronomy_AddDays(time, -dt/2.0); + astro_time_t t2 = Astronomy_AddDays(time, +dt/2.0); + double dist1, dist2; + int direction = *((int *)context); + astro_func_result_t result; + + dist1 = MoonDistance(t1); + dist2 = MoonDistance(t2); + result.value = direction * (dist2 - dist1) / dt; + result.status = ASTRO_SUCCESS; + return result; +} + +astro_apsis_t Astronomy_SearchLunarApsis(astro_time_t startTime) +{ + astro_time_t t1, t2; + astro_search_result_t search; + astro_func_result_t m1, m2; + int positive_direction = +1; + int negative_direction = -1; + const double increment = 5.0; /* number of days to skip in each iteration */ + astro_apsis_t result; + + /* + Check the rate of change of the distance dr/dt at the start time. + If it is positive, the Moon is currently getting farther away, + so start looking for apogee. + Conversely, if dr/dt < 0, start looking for perigee. + Either way, the polarity of the slope will change, so the product will be negative. + Handle the crazy corner case of exactly touching zero by checking for m1*m2 <= 0. + */ + + t1 = startTime; + m1 = distance_slope(&positive_direction, t1); + if (m1.status != ASTRO_SUCCESS) + return ApsisError(m1.status); + + for(;;) + { + t2 = Astronomy_AddDays(t1, increment); + m2 = distance_slope(&positive_direction, t2); + if (m2.status != ASTRO_SUCCESS) + return ApsisError(m2.status); + + if (m1.value * m2.value <= 0.0) + { + /* There is a change of slope polarity within the time range [t1, t2]. */ + /* Therefore this time range contains an apsis. */ + /* Figure out whether it is perigee or apogee. */ + + if (m1.value < 0.0 || m2.value > 0.0) + { + /* We found a minimum-distance event: perigee. */ + /* Search the time range for the time when the slope goes from negative to positive. */ + search = Astronomy_Search(distance_slope, &positive_direction, t1, t2, 1.0); + result.kind = APSIS_PERICENTER; + } + else if (m1.value > 0.0 || m2.value < 0.0) + { + /* We found a maximum-distance event: apogee. */ + /* Search the time range for the time when the slope goes from positive to negative. */ + search = Astronomy_Search(distance_slope, &negative_direction, t1, t2, 1.0); + result.kind = APSIS_APOCENTER; + } + else + { + /* This should never happen. It should not be possible for both slopes to be zero. */ + return ApsisError(ASTRO_INTERNAL_ERROR); + } + + if (search.status != ASTRO_SUCCESS) + return ApsisError(search.status); + + result.status = ASTRO_SUCCESS; + result.time = search.time; + result.dist_au = MoonDistance(search.time); + result.dist_km = result.dist_au * KM_PER_AU; + return result; + } + + /* We have not yet found a slope polarity change. Keep searching. */ + t1 = t2; + m1 = m2; + } +} + #ifdef __cplusplus } #endif @@ -3049,5 +3157,5 @@ astro_illum_t Astronomy_SearchPeakMagnitude(astro_body_t body, astro_time_t star ------------------------------------------- X NextLunarApsis - X SearchLunarApsis + - SearchLunarApsis */ diff --git a/generate/template/astronomy.js b/generate/template/astronomy.js index 8d97f4e8..44c99100 100644 --- a/generate/template/astronomy.js +++ b/generate/template/astronomy.js @@ -3302,7 +3302,7 @@ Astronomy.SearchLunarApsis = function(startDate) { // Check the rate of change of the distance dr/dt at the start time. // If it is positive, the Moon is currently getting farther away, // so start looking for apogee. - // Conversely, if dr/dt < 0, start looking for apogee. + // Conversely, if dr/dt < 0, start looking for perigee. // Either way, the polarity of the slope will change, so the product will be negative. // Handle the crazy corner case of exactly touching zero by checking for m1*m2 <= 0. diff --git a/source/c/astronomy.c b/source/c/astronomy.c index ad17bebe..f591c629 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -264,6 +264,18 @@ static astro_illum_t IllumError(astro_status_t status) return result; } +static astro_apsis_t ApsisError(astro_status_t status) +{ + astro_apsis_t result; + + result.status = status; + result.time.tt = result.time.ut = NAN; + result.kind = APSIS_INVALID; + result.dist_km = result.dist_au = NAN; + + return result; +} + static astro_func_result_t SynodicPeriod(astro_body_t body) { static const double Te = 365.256; /* Earth's orbital period in days */ @@ -4094,6 +4106,102 @@ astro_illum_t Astronomy_SearchPeakMagnitude(astro_body_t body, astro_time_t star return IllumError(ASTRO_SEARCH_FAILURE); } +static double MoonDistance(astro_time_t t) +{ + double lon, lat, dist; + CalcMoon(t.tt / 36525.0, &lon, &lat, &dist); + return dist; +} + +static astro_func_result_t distance_slope(void *context, astro_time_t time) +{ + static const double dt = 0.001; + astro_time_t t1 = Astronomy_AddDays(time, -dt/2.0); + astro_time_t t2 = Astronomy_AddDays(time, +dt/2.0); + double dist1, dist2; + int direction = *((int *)context); + astro_func_result_t result; + + dist1 = MoonDistance(t1); + dist2 = MoonDistance(t2); + result.value = direction * (dist2 - dist1) / dt; + result.status = ASTRO_SUCCESS; + return result; +} + +astro_apsis_t Astronomy_SearchLunarApsis(astro_time_t startTime) +{ + astro_time_t t1, t2; + astro_search_result_t search; + astro_func_result_t m1, m2; + int positive_direction = +1; + int negative_direction = -1; + const double increment = 5.0; /* number of days to skip in each iteration */ + astro_apsis_t result; + + /* + Check the rate of change of the distance dr/dt at the start time. + If it is positive, the Moon is currently getting farther away, + so start looking for apogee. + Conversely, if dr/dt < 0, start looking for perigee. + Either way, the polarity of the slope will change, so the product will be negative. + Handle the crazy corner case of exactly touching zero by checking for m1*m2 <= 0. + */ + + t1 = startTime; + m1 = distance_slope(&positive_direction, t1); + if (m1.status != ASTRO_SUCCESS) + return ApsisError(m1.status); + + for(;;) + { + t2 = Astronomy_AddDays(t1, increment); + m2 = distance_slope(&positive_direction, t2); + if (m2.status != ASTRO_SUCCESS) + return ApsisError(m2.status); + + if (m1.value * m2.value <= 0.0) + { + /* There is a change of slope polarity within the time range [t1, t2]. */ + /* Therefore this time range contains an apsis. */ + /* Figure out whether it is perigee or apogee. */ + + if (m1.value < 0.0 || m2.value > 0.0) + { + /* We found a minimum-distance event: perigee. */ + /* Search the time range for the time when the slope goes from negative to positive. */ + search = Astronomy_Search(distance_slope, &positive_direction, t1, t2, 1.0); + result.kind = APSIS_PERICENTER; + } + else if (m1.value > 0.0 || m2.value < 0.0) + { + /* We found a maximum-distance event: apogee. */ + /* Search the time range for the time when the slope goes from positive to negative. */ + search = Astronomy_Search(distance_slope, &negative_direction, t1, t2, 1.0); + result.kind = APSIS_APOCENTER; + } + else + { + /* This should never happen. It should not be possible for both slopes to be zero. */ + return ApsisError(ASTRO_INTERNAL_ERROR); + } + + if (search.status != ASTRO_SUCCESS) + return ApsisError(search.status); + + result.status = ASTRO_SUCCESS; + result.time = search.time; + result.dist_au = MoonDistance(search.time); + result.dist_km = result.dist_au * KM_PER_AU; + return result; + } + + /* We have not yet found a slope polarity change. Keep searching. */ + t1 = t2; + m1 = m2; + } +} + #ifdef __cplusplus } #endif @@ -4105,5 +4213,5 @@ astro_illum_t Astronomy_SearchPeakMagnitude(astro_body_t body, astro_time_t star ------------------------------------------- X NextLunarApsis - X SearchLunarApsis + - SearchLunarApsis */ diff --git a/source/c/astronomy.h b/source/c/astronomy.h index 8dc35e41..e02e75f0 100644 --- a/source/c/astronomy.h +++ b/source/c/astronomy.h @@ -212,6 +212,23 @@ typedef struct } astro_illum_t; +typedef enum +{ + APSIS_PERICENTER, + APSIS_APOCENTER, + APSIS_INVALID +} +astro_apsis_kind_t; + +typedef struct +{ + astro_status_t status; + astro_time_t time; + astro_apsis_kind_t kind; + double dist_au; + double dist_km; +} +astro_apsis_t; /*---------- functions ----------*/ @@ -282,6 +299,7 @@ astro_search_result_t Astronomy_SearchRiseSet( astro_seasons_t Astronomy_Seasons(int calendar_year); astro_illum_t Astronomy_Illumination(astro_body_t body, astro_time_t time); astro_illum_t Astronomy_SearchPeakMagnitude(astro_body_t body, astro_time_t startDate); +astro_apsis_t Astronomy_SearchLunarApsis(astro_time_t startTime); #ifdef __cplusplus } diff --git a/source/js/astronomy.js b/source/js/astronomy.js index 6ee6c8ae..82dcf401 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -4140,7 +4140,7 @@ Astronomy.SearchLunarApsis = function(startDate) { // Check the rate of change of the distance dr/dt at the start time. // If it is positive, the Moon is currently getting farther away, // so start looking for apogee. - // Conversely, if dr/dt < 0, start looking for apogee. + // Conversely, if dr/dt < 0, start looking for perigee. // Either way, the polarity of the slope will change, so the product will be negative. // Handle the crazy corner case of exactly touching zero by checking for m1*m2 <= 0.