From 7b543249b13da01f8409ea9573c0ba92e4bb9841 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 21 Jun 2021 15:34:56 -0400 Subject: [PATCH] Implemented C version of VectorObserver. --- demo/python/astronomy.py | 2 +- generate/ctest.c | 45 ++++++++++++- generate/template/astronomy.c | 118 +++++++++++++++++++++++++++++++++ generate/template/astronomy.py | 2 +- source/c/README.md | 25 +++++++ source/c/astronomy.c | 118 +++++++++++++++++++++++++++++++++ source/c/astronomy.h | 2 + source/python/README.md | 2 +- source/python/astronomy.py | 2 +- 9 files changed, 311 insertions(+), 5 deletions(-) diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 6845f0e3..e662086c 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -4309,7 +4309,7 @@ def VectorObserver(vector, ofdate): The components are expressed in astronomical units (AU). The time `vector.t` determines the Earth's rotation. ofdate : bool - Selects the date of the Earth's equator in which to express the equatorial coordinates. + Selects the date of the Earth's equator in which `vector` is expressed. The caller may pass `False` to use the orientation of the Earth's equator at noon UTC on January 1, 2000, in which case this function corrects for precession and nutation of the Earth as it was at the moment specified by the the time `vector.t`. diff --git a/generate/ctest.c b/generate/ctest.c index 03858ebb..8b83a59c 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -3913,7 +3913,8 @@ static int GeoidTestCase(astro_time_t time, astro_observer_t observer, astro_equ astro_vector_t surface; astro_vector_t geo_moon; astro_equatorial_t topo_moon; - double dx, dy, dz, diff; + double dx, dy, dz, diff, lat_diff, lon_diff, h_diff; + astro_observer_t vobs; topo_moon = Astronomy_Equator(BODY_MOON, &time, observer, equdate, NO_ABERRATION); CHECK_STATUS(topo_moon); @@ -3951,6 +3952,35 @@ static int GeoidTestCase(astro_time_t time, astro_observer_t observer, astro_equ if (diff > 1.0e-6) FAIL("C GeoidTestCase: EXCESSIVE POSITION ERROR.\n"); + /* Verify that we can convert the surface vector back to an observer. */ + vobs = Astronomy_VectorObserver(surface, equdate); + lat_diff = ABS(vobs.latitude - observer.latitude); + + /* Longitude is meaningless at the poles, so don't bother checking it there. */ + if (-89.99 <= observer.latitude && observer.latitude <= +89.99) + { + lon_diff = ABS(vobs.longitude - observer.longitude); + if (lon_diff > 180.0) + lon_diff = 360.0 - lon_diff; + lon_diff = ABS(lon_diff * cos(DEG2RAD * observer.latitude)); + if (lon_diff > 1.0e-6) + FAIL("C GeoidTestCase: EXCESSIVE longitude check error = %lf\n", lon_diff); + } + else + { + lon_diff = 0.0; + } + + h_diff = ABS(vobs.height - observer.height); + DEBUG("C GeoidTestCase: vobs=(lat=%lf, lon=%lf, height=%lf), lat_diff=%lg, lon_diff=%lg, h_diff=%lg\n", + vobs.latitude, vobs.longitude, vobs.height, lat_diff, lon_diff, h_diff); + + if (lat_diff > 1.0e-6) + FAIL("C GeoidTestCase: EXCESSIVE latitude check error = %lf\n", lat_diff); + + if (h_diff > 0.001) + FAIL("C GeoidTestCase: EXCESSIVE height check error = %lf\n", h_diff); + error = 0; fail: return error; @@ -3963,6 +3993,8 @@ static int GeoidTest(void) int tindex, oindex; astro_time_t time; astro_vector_t vec; + astro_observer_t observer; + int lat, lon; const astro_time_t time_list[] = { @@ -4001,6 +4033,17 @@ static int GeoidTest(void) } } + /* More exhaustive tests for a single time value across many different geographic coordinates. */ + time = Astronomy_MakeTime(2021, 06, 20, 15, 8, 0.0); + for (lat = -90; lat <= +90; lat += 1) + { + for (lon = -175; lon <= +180; lon += 5) + { + observer = Astronomy_MakeObserver(lat, lon, 0.0); + CHECK(GeoidTestCase(time, observer, EQUATOR_OF_DATE)); + } + } + printf("C GeoidTest: PASS\n"); error = 0; fail: diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index c0c1549f..03287605 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -163,6 +163,7 @@ static const double SUN_RADIUS_KM = 695700.0; #define EARTH_FLATTENING 0.996647180302104 #define EARTH_EQUATORIAL_RADIUS_KM 6378.1366 #define EARTH_EQUATORIAL_RADIUS_AU (EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU) +#define EARTH_POLAR_RADIUS_KM (EARTH_EQUATORIAL_RADIUS_KM * EARTH_FLATTENING) #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_MEAN_RADIUS_KM + EARTH_ATMOSPHERE_KM) @@ -1446,6 +1447,77 @@ static double sidereal_time(astro_time_t *time) return gst; /* return sidereal hours in the half-open range [0, 24). */ } +static astro_observer_t inverse_terra(const double ovec[3], double st) +{ + double x, y, z, p, F, W, D, c, s, c2, s2; + double lon_deg, lat_deg, lat, radicand, factor, denom, adjust; + double height_km, stlocl; + astro_observer_t observer; + + /* Convert from AU to kilometers. */ + x = ovec[0] * KM_PER_AU; + y = ovec[1] * KM_PER_AU; + z = ovec[2] * KM_PER_AU; + p = sqrt(x*x + y*y); + if (p < 1.0e-6) + { + /* Special case: within 1 millimeter of a pole! */ + /* Use arbitrary longitude, and latitude determined by polarity of z. */ + lon_deg = 0.0; + lat_deg = (z > 0.0) ? +90.0 : -90.0; + /* Elevation is calculated directly from z */ + height_km = fabs(z) - EARTH_POLAR_RADIUS_KM; + } + else + { + stlocl = atan2(y, x); + /* Calculate exact longitude. */ + lon_deg = RAD2DEG*stlocl - (15.0 * st); + /* Normalize longitude to the range (-180, +180]. */ + while (lon_deg <= -180.0) + lon_deg += 360.0; + while (lon_deg > +180.0) + lon_deg -= 360.0; + /* Numerically solve for exact latitude, using Newton's Method. */ + F = EARTH_FLATTENING * EARTH_FLATTENING; + /* Start with initial latitude estimate, based on a spherical Earth. */ + lat = atan2(z, p); + for(;;) + { + /* Calculate the error function W(lat). */ + /* We try to find the root of W, meaning where the error is 0. */ + c = cos(lat); + s = sin(lat); + factor = (F-1)*EARTH_EQUATORIAL_RADIUS_KM; + c2 = c*c; + s2 = s*s; + radicand = c2 + F*s2; + denom = sqrt(radicand); + W = (factor*s*c)/denom - z*c + p*s; + if (fabs(W) < 1.0e-12) + break; /* The error is now negligible. */ + /* Error is still too large. Find the next estimate. */ + /* Calculate D = the derivative of W with respect to lat. */ + D = factor*((c2 - s2)/denom - s2*c2*(F-1)/(factor*radicand)) + z*s + p*c; + lat -= W/D; + } + /* We now have a solution for the latitude in radians. */ + lat_deg = lat * RAD2DEG; + /* Solve for exact height in meters. */ + /* There are two formulas I can use. Use whichever has the less risky denominator. */ + adjust = EARTH_EQUATORIAL_RADIUS_KM / denom; + if (fabs(s) > fabs(c)) + height_km = z/s - F*adjust; + else + height_km = p/c - adjust; + } + + observer.latitude = lat_deg; + observer.longitude = lon_deg; + observer.height = 1000.0 * height_km; + return observer; +} + static void terra(astro_observer_t observer, double st, double pos[3]) { double df2 = EARTH_FLATTENING * EARTH_FLATTENING; @@ -2920,6 +2992,52 @@ astro_vector_t Astronomy_ObserverVector( } +/** + * @brief Calculates the geographic location corresponding to an equatorial vector. + * + * This is the inverse function of #Astronomy_ObserverVector. + * Instead of converting #astro_observer_t to #astro_vector_t, + * it converts `astro_vector_t` to `astro_observer_t`. + * + * @param vector + * The date and time for which to calculate the observer's position vector. + * The components are expressed in Astronomical Units (AU). + * You can calculate AU by dividing kilometers by the constant #KM_PER_AU. + * The time `vector.t` determines the Earth's rotation. + * + * @param equdate + * Selects the date of the Earth's equator in which `vector` is expressed. + * The caller may select `EQUATOR_J2000` to use the orientation of the Earth's equator + * at noon UTC on January 1, 2000, in which case this function corrects for precession + * and nutation of the Earth as it was at the moment specified by `vector.t`. + * Or the caller may select `EQUATOR_OF_DATE` to use the Earth's equator at `vector.t` + * as the orientation. + * + * @return + * The geographic latitude, longitude, and elevation above sea level + * that corresponds to the given equatorial vector. + */ +astro_observer_t Astronomy_VectorObserver( + astro_vector_t vector, + astro_equator_date_t equdate) +{ + double gast; + double pos1[3]; + double pos2[3]; + + gast = sidereal_time(&vector.t); + pos1[0] = vector.x; + pos1[1] = vector.y; + pos1[2] = vector.z; + if (equdate == EQUATOR_J2000) + { + precession(pos1, vector.t, FROM_2000, pos2); + nutation(pos2, &vector.t, FROM_2000, pos1); + } + return inverse_terra(pos1, gast); +} + + /** * @brief Calculates the apparent location of a body relative to the local horizon of an observer on the Earth. * diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 78af857f..e03a8526 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -2279,7 +2279,7 @@ def VectorObserver(vector, ofdate): The components are expressed in astronomical units (AU). The time `vector.t` determines the Earth's rotation. ofdate : bool - Selects the date of the Earth's equator in which to express the equatorial coordinates. + Selects the date of the Earth's equator in which `vector` is expressed. The caller may pass `False` to use the orientation of the Earth's equator at noon UTC on January 1, 2000, in which case this function corrects for precession and nutation of the Earth as it was at the moment specified by the the time `vector.t`. diff --git a/source/c/README.md b/source/c/README.md index 1fcb6ae7..42f16884 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -2332,6 +2332,31 @@ Calculates the non-negative length of the given vector. The length is expressed + +--- + + +### Astronomy_VectorObserver(vector, equdate) ⇒ [`astro_observer_t`](#astro_observer_t) + +**Calculates the geographic location corresponding to an equatorial vector.** + + + +This is the inverse function of [`Astronomy_ObserverVector`](#Astronomy_ObserverVector). Instead of converting [`astro_observer_t`](#astro_observer_t) to [`astro_vector_t`](#astro_vector_t), it converts `[`astro_vector_t`](#astro_vector_t)` to `[`astro_observer_t`](#astro_observer_t)`. + + + +**Returns:** The geographic latitude, longitude, and elevation above sea level that corresponds to the given equatorial vector. + + + +| Type | Parameter | Description | +| --- | --- | --- | +| [`astro_vector_t`](#astro_vector_t) | `vector` | The date and time for which to calculate the observer's position vector. The components are expressed in Astronomical Units (AU). You can calculate AU by dividing kilometers by the constant [`KM_PER_AU`](#KM_PER_AU). The time `vector.t` determines the Earth's rotation. | +| [`astro_equator_date_t`](#astro_equator_date_t) | `equdate` | Selects the date of the Earth's equator in which `vector` is expressed. The caller may select `EQUATOR_J2000` to use the orientation of the Earth's equator at noon UTC on January 1, 2000, in which case this function corrects for precession and nutation of the Earth as it was at the moment specified by `vector.t`. Or the caller may select `EQUATOR_OF_DATE` to use the Earth's equator at `vector.t` as the orientation. | + + + ## Constants diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 79b8fb42..1d7096b3 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -163,6 +163,7 @@ static const double SUN_RADIUS_KM = 695700.0; #define EARTH_FLATTENING 0.996647180302104 #define EARTH_EQUATORIAL_RADIUS_KM 6378.1366 #define EARTH_EQUATORIAL_RADIUS_AU (EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU) +#define EARTH_POLAR_RADIUS_KM (EARTH_EQUATORIAL_RADIUS_KM * EARTH_FLATTENING) #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_MEAN_RADIUS_KM + EARTH_ATMOSPHERE_KM) @@ -1524,6 +1525,77 @@ static double sidereal_time(astro_time_t *time) return gst; /* return sidereal hours in the half-open range [0, 24). */ } +static astro_observer_t inverse_terra(const double ovec[3], double st) +{ + double x, y, z, p, F, W, D, c, s, c2, s2; + double lon_deg, lat_deg, lat, radicand, factor, denom, adjust; + double height_km, stlocl; + astro_observer_t observer; + + /* Convert from AU to kilometers. */ + x = ovec[0] * KM_PER_AU; + y = ovec[1] * KM_PER_AU; + z = ovec[2] * KM_PER_AU; + p = sqrt(x*x + y*y); + if (p < 1.0e-6) + { + /* Special case: within 1 millimeter of a pole! */ + /* Use arbitrary longitude, and latitude determined by polarity of z. */ + lon_deg = 0.0; + lat_deg = (z > 0.0) ? +90.0 : -90.0; + /* Elevation is calculated directly from z */ + height_km = fabs(z) - EARTH_POLAR_RADIUS_KM; + } + else + { + stlocl = atan2(y, x); + /* Calculate exact longitude. */ + lon_deg = RAD2DEG*stlocl - (15.0 * st); + /* Normalize longitude to the range (-180, +180]. */ + while (lon_deg <= -180.0) + lon_deg += 360.0; + while (lon_deg > +180.0) + lon_deg -= 360.0; + /* Numerically solve for exact latitude, using Newton's Method. */ + F = EARTH_FLATTENING * EARTH_FLATTENING; + /* Start with initial latitude estimate, based on a spherical Earth. */ + lat = atan2(z, p); + for(;;) + { + /* Calculate the error function W(lat). */ + /* We try to find the root of W, meaning where the error is 0. */ + c = cos(lat); + s = sin(lat); + factor = (F-1)*EARTH_EQUATORIAL_RADIUS_KM; + c2 = c*c; + s2 = s*s; + radicand = c2 + F*s2; + denom = sqrt(radicand); + W = (factor*s*c)/denom - z*c + p*s; + if (fabs(W) < 1.0e-12) + break; /* The error is now negligible. */ + /* Error is still too large. Find the next estimate. */ + /* Calculate D = the derivative of W with respect to lat. */ + D = factor*((c2 - s2)/denom - s2*c2*(F-1)/(factor*radicand)) + z*s + p*c; + lat -= W/D; + } + /* We now have a solution for the latitude in radians. */ + lat_deg = lat * RAD2DEG; + /* Solve for exact height in meters. */ + /* There are two formulas I can use. Use whichever has the less risky denominator. */ + adjust = EARTH_EQUATORIAL_RADIUS_KM / denom; + if (fabs(s) > fabs(c)) + height_km = z/s - F*adjust; + else + height_km = p/c - adjust; + } + + observer.latitude = lat_deg; + observer.longitude = lon_deg; + observer.height = 1000.0 * height_km; + return observer; +} + static void terra(astro_observer_t observer, double st, double pos[3]) { double df2 = EARTH_FLATTENING * EARTH_FLATTENING; @@ -4146,6 +4218,52 @@ astro_vector_t Astronomy_ObserverVector( } +/** + * @brief Calculates the geographic location corresponding to an equatorial vector. + * + * This is the inverse function of #Astronomy_ObserverVector. + * Instead of converting #astro_observer_t to #astro_vector_t, + * it converts `astro_vector_t` to `astro_observer_t`. + * + * @param vector + * The date and time for which to calculate the observer's position vector. + * The components are expressed in Astronomical Units (AU). + * You can calculate AU by dividing kilometers by the constant #KM_PER_AU. + * The time `vector.t` determines the Earth's rotation. + * + * @param equdate + * Selects the date of the Earth's equator in which `vector` is expressed. + * The caller may select `EQUATOR_J2000` to use the orientation of the Earth's equator + * at noon UTC on January 1, 2000, in which case this function corrects for precession + * and nutation of the Earth as it was at the moment specified by `vector.t`. + * Or the caller may select `EQUATOR_OF_DATE` to use the Earth's equator at `vector.t` + * as the orientation. + * + * @return + * The geographic latitude, longitude, and elevation above sea level + * that corresponds to the given equatorial vector. + */ +astro_observer_t Astronomy_VectorObserver( + astro_vector_t vector, + astro_equator_date_t equdate) +{ + double gast; + double pos1[3]; + double pos2[3]; + + gast = sidereal_time(&vector.t); + pos1[0] = vector.x; + pos1[1] = vector.y; + pos1[2] = vector.z; + if (equdate == EQUATOR_J2000) + { + precession(pos1, vector.t, FROM_2000, pos2); + nutation(pos2, &vector.t, FROM_2000, pos1); + } + return inverse_terra(pos1, gast); +} + + /** * @brief Calculates the apparent location of a body relative to the local horizon of an observer on the Earth. * diff --git a/source/c/astronomy.h b/source/c/astronomy.h index 305e3d4d..a98e5ee3 100644 --- a/source/c/astronomy.h +++ b/source/c/astronomy.h @@ -978,6 +978,8 @@ astro_vector_t Astronomy_ObserverVector( astro_equator_date_t equdate ); +astro_observer_t Astronomy_VectorObserver(astro_vector_t vector, astro_equator_date_t equdate); + astro_ecliptic_t Astronomy_SunPosition(astro_time_t time); astro_ecliptic_t Astronomy_Ecliptic(astro_vector_t equ); astro_angle_result_t Astronomy_EclipticLongitude(astro_body_t body, astro_time_t time); diff --git a/source/python/README.md b/source/python/README.md index 61c341a2..b222e1e6 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -2580,7 +2580,7 @@ a `Vector` to an `Observer`. | Type | Parameter | Description | | --- | --- | --- | | [`Vector`](#Vector) | `vector` | The geocentric vector in an equatorial orientation. The components are expressed in astronomical units (AU). The time `vector.t` determines the Earth's rotation. | -| `bool` | `ofdate` | Selects the date of the Earth's equator in which to express the equatorial coordinates. The caller may pass `False` to use the orientation of the Earth's equator at noon UTC on January 1, 2000, in which case this function corrects for precession and nutation of the Earth as it was at the moment specified by the the time `vector.t`. Or the caller may pass `True` to use the Earth's equator at `vector.t` as the orientation. | +| `bool` | `ofdate` | Selects the date of the Earth's equator in which `vector` is expressed. The caller may pass `False` to use the orientation of the Earth's equator at noon UTC on January 1, 2000, in which case this function corrects for precession and nutation of the Earth as it was at the moment specified by the the time `vector.t`. Or the caller may pass `True` to use the Earth's equator at `vector.t` as the orientation. | ### Returns: [`Observer`](#Observer) The geographic latitude, longitude, and elevation above sea level diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 6845f0e3..e662086c 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -4309,7 +4309,7 @@ def VectorObserver(vector, ofdate): The components are expressed in astronomical units (AU). The time `vector.t` determines the Earth's rotation. ofdate : bool - Selects the date of the Earth's equator in which to express the equatorial coordinates. + Selects the date of the Earth's equator in which `vector` is expressed. The caller may pass `False` to use the orientation of the Earth's equator at noon UTC on January 1, 2000, in which case this function corrects for precession and nutation of the Earth as it was at the moment specified by the the time `vector.t`.