From 4399e72334082b2e330edcc4dbfe93b59d5ffff2 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 9 Dec 2019 11:06:02 -0500 Subject: [PATCH] Added Astronomy_HorizonFromVector, Astronomy_Rotation_EQD_HOR. Astronomy_HorizonFromVector is a specialized variant of Astronomy_SphereFromVector that flips the orientation of the azimuth angle to the more traditional clockwise-from-north direction used in navigation and cartography. It also allows the same optional refraction correction as Astronomy_Horizon. Astronomy_Rotation_EQD_HOR converts a equatorial-of-date vector to a horizontal vector. The horizontal vector has the following components: x = North y = West z = Zenith Removed trailing whitespace in generate.c. --- generate/ctest.c | 60 +++++++++++ generate/generate.c | 74 +++++++------- generate/template/astronomy.c | 182 +++++++++++++++++++++++++++++----- source/c/README.md | 66 +++++++++++- source/c/astronomy.c | 182 +++++++++++++++++++++++++++++----- source/c/astronomy.h | 3 +- 6 files changed, 482 insertions(+), 85 deletions(-) diff --git a/generate/ctest.c b/generate/ctest.c index 8045a431..6a2b0bbb 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -2229,6 +2229,59 @@ fail: return error; } +static int Test_EQD_HOR(astro_body_t body) +{ + astro_time_t time; + astro_observer_t observer; + astro_equatorial_t eqd; + astro_horizon_t hor; + astro_rotation_t rot; + astro_vector_t vec_eqd, vec_hor; + astro_spherical_t sphere; + double diff_az, diff_alt; + int error; + + /* Use existing functions to calculate horizontal coordinates of the body for the time+observer. */ + time = Astronomy_MakeTime(1970, 12, 13, 5, 15, 0.0); + observer = Astronomy_MakeObserver(-37.0, +45.0, 0.0); + CHECK_EQU(eqd, Astronomy_Equator(body, &time, observer, EQUATOR_OF_DATE, ABERRATION)); + printf("Test_EQD_HOR %s: OFDATE ra=%0.3lf, dec=%0.3lf\n", Astronomy_BodyName(body), eqd.ra, eqd.dec); + hor = Astronomy_Horizon(&time, observer, eqd.ra, eqd.dec, REFRACTION_NORMAL); + + /* Calculate the position of the body as an equatorial vector of date. */ + sphere.status = ASTRO_SUCCESS; + sphere.dist = eqd.dist; + sphere.lat = eqd.dec; + sphere.lon = 15.0 * eqd.ra; + CHECK_VECTOR(vec_eqd, Astronomy_VectorFromSphere(sphere, time)); + + /* Calculate rotation matrix to convert equatorial J2000 vector to horizontal vector. */ + rot = Astronomy_Rotation_EQD_HOR(time, observer); + + /* Rotate the equator of date vector to a horizontal vector. */ + CHECK_VECTOR(vec_hor, Astronomy_RotateVector(rot, vec_eqd)); + + /* Convert the horizontal vector to horizontal angular coordinates. */ + sphere = Astronomy_HorizonFromVector(vec_hor, REFRACTION_NORMAL); + CHECK_STATUS(sphere); + + diff_alt = fabs(sphere.lat - hor.altitude); + diff_az = fabs(sphere.lon - hor.azimuth); + + printf("Test_EQD_HOR %s: trusted alt=%0.3lf, az=%0.3lf; test alt=%0.3lf, az=%0.3lf; diff_alt=%lg, diff_az=%lg\n", + Astronomy_BodyName(body), hor.altitude, hor.azimuth, sphere.lat, sphere.lon, diff_alt, diff_az); + + if (diff_alt > 2.0e-14) + { + fprintf(stderr, "Test_EQD_HOR: EXCESSIVE HORIZONTAL ERROR.\n"); + return 1; + } + + error = 0; +fail: + return error; +} + static int RotationTest(void) { int error; @@ -2262,12 +2315,19 @@ static int RotationTest(void) CHECK(TestSpin(0.0, 0.0, -45.0, +1, 0, 0, +0.7071067811865476, -0.7071067811865476, 0)); CHECK(Test_EQJ_ECL()); + CHECK(Test_EQJ_EQD(BODY_MERCURY)); CHECK(Test_EQJ_EQD(BODY_VENUS)); CHECK(Test_EQJ_EQD(BODY_MARS)); CHECK(Test_EQJ_EQD(BODY_JUPITER)); CHECK(Test_EQJ_EQD(BODY_SATURN)); + CHECK(Test_EQD_HOR(BODY_MERCURY)); + CHECK(Test_EQD_HOR(BODY_VENUS)); + CHECK(Test_EQD_HOR(BODY_MARS)); + CHECK(Test_EQD_HOR(BODY_JUPITER)); + CHECK(Test_EQD_HOR(BODY_SATURN)); + error = 0; fail: return error; diff --git a/generate/generate.c b/generate/generate.c index d44d4655..ec591232 100644 --- a/generate/generate.c +++ b/generate/generate.c @@ -97,14 +97,14 @@ static int UnitTestChebyshev(void); #define NEPTUNE_PERIHELION 29.81 #define PLUTO_PERIHELION 29.658 -/* - ErrorRadius[body] is the worst-case (smallest) distance - in AU over which a position error affects Earth-based observations. +/* + ErrorRadius[body] is the worst-case (smallest) distance + in AU over which a position error affects Earth-based observations. For example, an error in the Earth/Moon Barycenter (EMB) affects observations of Venus the most, because that planet comes closest to Earth at 0.277 AU. */ -const double ErrorRadius[] = +const double ErrorRadius[] = { EARTH_PERIHELION - MERCURY_APHELION, /* 0 = Mercury */ EARTH_PERIHELION - VENUS_APHELION, /* 1 = Venus */ @@ -138,7 +138,7 @@ int main(int argc, const char *argv[]) static int PrintUsage(void) { - fprintf(stderr, + fprintf(stderr, "\n" "USAGE:\n" "\n" @@ -188,8 +188,8 @@ static int NovasBodyPos(double jd, int body, double pos[3]) { double factor; - /* - The caller is asking for the Earth's position or the Moon's position. + /* + The caller is asking for the Earth's position or the Moon's position. NOVAS does not directly represent either body. Instead, we have to calculate the Earth or Moon using the Earth/Moon Barycenter (EMB) and the Geocentric Moon (GM). @@ -301,7 +301,7 @@ static int SearchVsop(int body) winner_terms = trunc_terms; winner_arcmin = max_arcmin; winner_threshold = threshold; - + error = SaveVsopFile(&model); if (error) goto fail; } @@ -326,7 +326,7 @@ fail: static int SaveVsopFile(const vsop_model_t *model) { char filename[100]; - + snprintf(filename, sizeof(filename), "output/vsop_%d.txt", (int)model->body); return VsopWriteTrunc(model, filename); } @@ -436,13 +436,13 @@ static int PositionArcminError(int body, double jd, const double a[3], const dou return 0; } - if (body == BODY_MOON) + if (body == BODY_MOON) { /* Measure the position of the moon as seen from the *SURFACE* of the Earth, at the point along a line from the geocenter to the object. This is closer than the geocenter or the barycenter, so the anglular error is larger. - For other bodies, there is no significant difference between + For other bodies, there is no significant difference between topocentric error and geocentric error. */ const double scale = (ERAD / AU) / VectorLength(a); @@ -465,7 +465,7 @@ static int PositionArcminError(int body, double jd, const double a[3], const dou adiff.v[0] = a[0]; adiff.v[1] = a[1]; adiff.v[2] = a[2]; - + bdiff.t = jd; /* doesn't matter, but I don't like leaving undefined memory */ bdiff.v[0] = b[0]; bdiff.v[1] = b[1]; @@ -477,7 +477,7 @@ static int PositionArcminError(int body, double jd, const double a[3], const dou if (body == BODY_EMB || body == BODY_EARTH) { - /* + /* For the Earth/Moon Barycenter (EMB=2), or Earth=11, we use a special estimate of angular error. An error in the Earth's position affects the observed position of all other bodies in the solar system. @@ -528,7 +528,7 @@ calc: adiff.v[0] = a[0] - epos[0]; adiff.v[1] = a[1] - epos[1]; adiff.v[2] = a[2] - epos[2]; - + bdiff.t = jd; /* doesn't matter, but I don't like leaving undefined memory */ bdiff.v[0] = b[0] - epos[0]; bdiff.v[1] = b[1] - epos[1]; @@ -578,7 +578,7 @@ static int MeasureError(const char *inFileName, int nsamples, error_stats_t *sta memset(stats, 0, sizeof(error_stats_t)); error = EphFileOpen(&reader, inFileName); - if (error) + if (error) { fprintf(stderr, "MeasureError: Error %d trying to open file: %s\n", error, inFileName); goto fail; @@ -614,7 +614,7 @@ static int MeasureError(const char *inFileName, int nsamples, error_stats_t *sta /* Calculate "correct" position at the time 'jd' (according to NOVAS). */ error = NovasBodyPos(jd, reader.body, npos); if (error) goto fail; - + /* Calculate our approximation of the position at the time 'jd'. */ x = ChebScale(record.jdStart, record.jdStart + record.jdDelta, jd); ChebApprox(record.numpoly, 3, record.coeff, x, apos); @@ -682,7 +682,7 @@ static int Resample(int body, const char *outFileName, int npoly, int startYear, } error = ChebInit(&encoder, npoly, VECTOR_DIM); - if (error) + if (error) { fprintf(stderr, "ERROR %d returned by ChebInit()\n", error); goto fail; @@ -745,7 +745,7 @@ static int SampleFunc(const void *context, double jd, double pos[CHEB_MAX_DIM]) double jed[2]; double sun_pos[3]; double vel[3]; /* we don't care about velocities... ignored */ - + jed[0] = jd; jed[1] = 0.0; error = state(jed, (short)c->body, pos, vel); @@ -839,9 +839,9 @@ static int CheckEcliptic(const char *filename, int lnum, const char *line, doubl double earth_ecl[3]; double geo[3]; double tt, jd, lon, dist, blon, elon, diff; - + *arcmin = 99999.0; - *outbody = VSOP_INVALID_BODY; + *outbody = VSOP_INVALID_BODY; /* Example line: "e Jupiter opp " */ if (4 != sscanf(line, "e %9[A-Za-z] %9[a-z] %lf %lf", name, event, &tt, &dist)) @@ -883,7 +883,7 @@ static int CheckEcliptic(const char *filename, int lnum, const char *line, doubl { fprintf(stderr, "CheckEcliptic: Invalid distance %0.3lf AU for inferior conjunction; line %d file %s\n", dist, lnum, filename); return 1; - } + } } else if (!strcmp(event, "sup")) { @@ -895,16 +895,16 @@ static int CheckEcliptic(const char *filename, int lnum, const char *line, doubl { fprintf(stderr, "CheckEcliptic: Invalid distance %0.3lf AU for superior conjunction; line %d file %s\n", dist, lnum, filename); return 1; - } + } } else { fprintf(stderr, "CheckEcliptic: Invalid event code '%s' on line %d of file %s\n", event, lnum, filename); return 1; - } + } blon = EclipticLongitude(body_ecl); - elon = EclipticLongitude(earth_ecl); + elon = EclipticLongitude(earth_ecl); diff = LongitudeOffset((elon - blon) - lon); *arcmin = fabs(diff * 60.0); @@ -943,14 +943,14 @@ static int CheckTestVector(const char *filename, int lnum, const char *line, dou } error = NovasBodyPos(jd, body, npos); - if (error) + if (error) { fprintf(stderr, "CheckTestVector: NovasBodyPos returned %d on line %d of file %s\n", error, lnum, filename); return error; } error = PositionArcminError(body, jd, npos, xpos, arcmin); - if (error) + if (error) { fprintf(stderr, "CheckTestVector: PositionArcminError returned %d on line %d of file %s\n", error, lnum, filename); return error; @@ -996,7 +996,7 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const return 1; } - if (body == BODY_EARTH || body == BODY_EMB) + if (body == BODY_EARTH || body == BODY_EMB) { fprintf(stderr, "CheckSkyPos: Cannot calculate sky position of body '%s'\n", name); return 1; @@ -1007,7 +1007,7 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const else if (body == BODY_SUN) bodyIndex = 10; else - bodyIndex = 1 + body; + bodyIndex = 1 + body; error = make_object(0, (short)bodyIndex, name, NULL, &obj); if (error) @@ -1031,12 +1031,12 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const delta_ra = fabs(sky_ast.ra - ra); if (delta_ra > 12.0) { - /* + /* Sometimes the two RA values can straddle the 24-hour mark. For example, one of them is 0.001 and the other 23.999. The actual error is then 0.002 hours, not 23.998. In general, it is never "fair" to call the error greater than - 12 hours or 180 degrees. + 12 hours or 180 degrees. */ delta_ra = 24.0 - delta_ra; } @@ -1065,10 +1065,10 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const return error; } - equ2hor(jd_utc, jd_tt-jd_utc, 1, 0.0, 0.0, - &location->on_surf, sky_dat.ra, sky_dat.dec, 0, + equ2hor(jd_utc, jd_tt-jd_utc, 1, 0.0, 0.0, + &location->on_surf, sky_dat.ra, sky_dat.dec, 0, &check_zenith, &check_azimuth, &check_ra, &check_dec); - + check_altitude = 90.0 - check_zenith; delta_az = fabs(azimuth - check_azimuth); @@ -1077,10 +1077,10 @@ static int CheckSkyPos(observer *location, const char *filename, int lnum, const delta_az *= 60.0; /* Just like RA, diminish the azimuth error as altitude approaches zenith/nadir... */ delta_az *= cos(check_altitude * DEG2RAD); - + delta_alt = 60.0 * (altitude - check_altitude); - *arcmin_hor = sqrt(delta_az*delta_az + delta_alt*delta_alt); + *arcmin_hor = sqrt(delta_az*delta_az + delta_alt*delta_alt); if (*arcmin_hor > 0.9) { @@ -1211,7 +1211,7 @@ static int CheckTestOutput(const char *filename) CHECK(ParseObserver(filename, lnum, line, &location)); break; - case 'v': /* heliocentric cartesian vector represented in J2000 equatorial plane */ + case 'v': /* heliocentric cartesian vector represented in J2000 equatorial plane */ CHECK(CheckTestVector(filename, lnum, line, &arcmin_helio, &body)); UpdateErrorStats(&bundle[body].helio, arcmin_helio); break; @@ -1226,7 +1226,7 @@ static int CheckTestOutput(const char *filename) CHECK(CheckEcliptic(filename, lnum, line, &arcmin_eclip, &body)); UpdateErrorStats(&bundle[body].eclip, arcmin_eclip); break; - + default: fprintf(stderr, "CheckTestOutput: Invalid first character on line %d of file %s\n", lnum, filename); error = 1; diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 8544bc14..88b3286f 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -4191,7 +4191,7 @@ astro_vector_t Astronomy_VectorFromSphere(astro_spherical_t sphere, astro_time_t * Given a Cartesian vector, returns latitude, longitude, and distance. * * @param vector - * Cartesian vector to be converted to spherical coordinates + * Cartesian vector to be converted to spherical coordinates. * * @return * Spherical coordinates that are equivalent to the given vector. @@ -4235,6 +4235,100 @@ astro_spherical_t Astronomy_SphereFromVector(astro_vector_t vector) return sphere; } +/** + * @brief Converts Cartesian coordinates to horizontal coordinates. + * + * Given a horizontal Cartesian vector, returns horizontal azimuth and altitude. + * + * *IMPORTANT:* This function differs from #Astronomy_SphereFromVector in two ways: + * - `Astronomy_SphereFromVector` returns a `lon` value that represents azimuth defined counterclockwise + * from north (e.g., west = +90), but this function represents a clockwise rotation + * (e.g., east = +90). The difference is because `Astronomy_SphereFromVector` is intended + * to preserve the vector "right-hand rule", while this function defines azimuth in a more + * traditional way as used in navigation and cartography. + * - This function optionally corrects for atmospheric refraction, while `Astronomy_SphereFromVector` + * does not. + * + * The returned structure will contain the azimuth in `lon`. + * It is measured in degrees clockwise from north: east = +90 degrees, west = +270 degrees. + * + * The altitude will be stored in `lat`. + * + * The distance to the observed object is stored in `dist`, + * and is expressed in astronomical units (AU). + * + * @param vector + * Cartesian vector to be converted to horizontal coordinates. + * + * @param refraction + * `REFRACTION_NORMAL`: correct altitude for atmospheric refraction (recommended). + * `REFRACTION_NONE`: no atmospheric refraction correction is performed. + * `REFRACTION_JPLHOR`: for JPL Horizons compatibility testing only; not recommended for normal use. + * + * @return + * If successful, `status` hold `ASTRO_SUCCESS` and the other fields are valid as described above. + * Otherwise `status` holds an error code and the other fields are undefined. + */ +astro_spherical_t Astronomy_HorizonFromVector(astro_vector_t vector, astro_refraction_t refraction) +{ + astro_spherical_t sphere; + + sphere = Astronomy_SphereFromVector(vector); + if (sphere.status == ASTRO_SUCCESS) + { + /* Convert azimuth from counterclockwise-from-north to clockwise-from-north. */ + sphere.lon = 360.0 - sphere.lon; + if (sphere.lon >= 360.0) + sphere.lon -= 360.0; + else if (sphere.lon < 0.0) + sphere.lon += 360.0; + + if (refraction == REFRACTION_NONE) + { + /* No atmospheric correction needed; do nothing. */ + } + else if (refraction == REFRACTION_NORMAL || refraction == REFRACTION_JPLHOR) + { + double zd, refr, hd; + + zd = 90.0 - sphere.lat; + + // http://extras.springer.com/1999/978-1-4471-0555-8/chap4/horizons/horizons.pdf + // JPL Horizons says it uses refraction algorithm from + // Meeus "Astronomical Algorithms", 1991, p. 101-102. + // I found the following Go implementation: + // https://github.com/soniakeys/meeus/blob/master/v3/refraction/refract.go + // This is a translation from the function "Saemundsson" there. + // I found experimentally that JPL Horizons clamps the angle to 1 degree below the horizon. + // This is important because the 'refr' formula below goes crazy near hd = -5.11. + hd = sphere.lat; + if (hd < -1.0) + hd = -1.0; + + refr = (1.02 / tan((hd+10.3/(hd+5.11))*DEG2RAD)) / 60.0; + + if (refraction == REFRACTION_NORMAL && zd > 91.0) + { + // In "normal" mode we gradually reduce refraction toward the nadir + // so that we never get an altitude angle less than -90 degrees. + // When horizon angle is -1 degrees, zd = 91, and the factor is exactly 1. + // As zd approaches 180 (the nadir), the fraction approaches 0 linearly. + refr *= (180.0 - zd) / 89.0; + } + + sphere.lat += refr; + } + else + { + /* Invalid refraction option. */ + return SphereError(ASTRO_INVALID_PARAMETER); + } + } + + return sphere; +} + + /** * @brief * Applies a rotation to a vector, yielding a rotated vector. @@ -4268,28 +4362,6 @@ astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t } -#if 0 -/** - * @brief - * Calculates a rotation matrix from equatorial of-date system to equatorial J2000. - * - * This is one of the family of MakeRotation functions that convert - * from one orientation system to another. - * Source: EQD = equatorial system, using equator of date. - * Target: EQJ = equatorial system, using equator at J2000 epoch. - * - * @param time - * The date and time defining the equator of the source orientation. - * - * @return - * A rotation matrix that converts EQD to EQJ. - */ -astro_rotation_t Astronomy_MakeRotation_EQD_EQJ(astro_time_t time) -{ - -} -#endif - /** * @brief * Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL). @@ -4390,6 +4462,70 @@ astro_rotation_t Astronomy_Rotation_EQD_EQJ(astro_time_t time) return Astronomy_CombineRotation(nut, prec); } + +/** + * @brief + * Calculates a rotation matrix from equatorial of-date (EQD) to horizontal (HOR). + * + * This is one of the family of functions that returns a rotation matrix + * for converting from one orientation to another. + * Source: EQD = equatorial system, using equator of the specified date/time. + * Target: HOR = horizontal system. + * + * Use #Astronomy_HorizonFromVector to convert the return value + * to a traditional altitude/azimuth pair. + * + * @param time + * The date and time at which the Earth's equator applies. + * + * @param observer + * A location near the Earth's mean sea level that defines the observer's horizon. + * + * @return + * A rotation matrix that converts EQD to HOR at `time` and for `observer`. + * The components of the horizontal vector are: + * x = north, y = west, z = zenith (straight up from the observer). + * These components are chosen so that the "right-hand rule" works for the vector + * and so that north represents the direction where azimuth = 0. + */ +astro_rotation_t Astronomy_Rotation_EQD_HOR(astro_time_t time, astro_observer_t observer) +{ + astro_rotation_t rot; + double uze[3], une[3], uwe[3]; + double uz[3], un[3], uw[3]; + double spin_angle; + + double sinlat = sin(observer.latitude * DEG2RAD); + double coslat = cos(observer.latitude * DEG2RAD); + double sinlon = sin(observer.longitude * DEG2RAD); + double coslon = cos(observer.longitude * DEG2RAD); + + uze[0] = coslat * coslon; + uze[1] = coslat * sinlon; + uze[2] = sinlat; + + une[0] = -sinlat * coslon; + une[1] = -sinlat * sinlon; + une[2] = coslat; + + uwe[0] = sinlon; + uwe[1] = -coslon; + uwe[2] = 0.0; + + spin_angle = -15.0 * sidereal_time(&time); + spin(spin_angle, uze, uz); + spin(spin_angle, une, un); + spin(spin_angle, uwe, uw); + + rot.rot[0][0] = un[0]; rot.rot[1][0] = un[1]; rot.rot[2][0] = un[2]; + rot.rot[0][1] = uw[0]; rot.rot[1][1] = uw[1]; rot.rot[2][1] = uw[2]; + rot.rot[0][2] = uz[0]; rot.rot[1][2] = uz[1]; rot.rot[2][2] = uz[2]; + + rot.status = ASTRO_SUCCESS; + return rot; +} + + #ifdef __cplusplus } #endif diff --git a/source/c/README.md b/source/c/README.md index d0635aa4..31f82909 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -473,6 +473,43 @@ This function optionally corrects for atmospheric refraction. For most uses, it +--- + + +### Astronomy_HorizonFromVector(vector, refraction) ⇒ [`astro_spherical_t`](#astro_spherical_t) + +**Converts Cartesian coordinates to horizontal coordinates.** + + + +Given a horizontal Cartesian vector, returns horizontal azimuth and altitude. + +*IMPORTANT:* This function differs from [`Astronomy_SphereFromVector`](#Astronomy_SphereFromVector) in two ways: + +- `Astronomy_SphereFromVector` returns a `lon` value that represents azimuth defined counterclockwise from north (e.g., west = +90), but this function represents a clockwise rotation (e.g., east = +90). The difference is because `Astronomy_SphereFromVector` is intended to preserve the vector "right-hand rule", while this function defines azimuth in a more traditional way as used in navigation and cartography. +- This function optionally corrects for atmospheric refraction, while `Astronomy_SphereFromVector` does not. + + +The returned structure will contain the azimuth in `lon`. It is measured in degrees clockwise from north: east = +90 degrees, west = +270 degrees. + +The altitude will be stored in `lat`. + +The distance to the observed object is stored in `dist`, and is expressed in astronomical units (AU). + + + +**Returns:** If successful, `status` hold `ASTRO_SUCCESS` and the other fields are valid as described above. Otherwise `status` holds an error code and the other fields are undefined. + + + +| Type | Parameter | Description | +| --- | --- | --- | +| [`astro_vector_t`](#astro_vector_t) | `vector` | Cartesian vector to be converted to horizontal coordinates. | +| [`astro_refraction_t`](#astro_refraction_t) | `refraction` | `REFRACTION_NORMAL`: correct altitude for atmospheric refraction (recommended). `REFRACTION_NONE`: no atmospheric refraction correction is performed. `REFRACTION_JPLHOR`: for JPL Horizons compatibility testing only; not recommended for normal use. | + + + + --- @@ -764,6 +801,33 @@ This is one of the family of functions that returns a rotation matrix for conver +--- + + +### Astronomy_Rotation_EQD_HOR(time, observer) ⇒ [`astro_rotation_t`](#astro_rotation_t) + +**Calculates a rotation matrix from equatorial of-date (EQD) to horizontal (HOR).** + + + +This is one of the family of functions that returns a rotation matrix for converting from one orientation to another. Source: EQD = equatorial system, using equator of the specified date/time. Target: HOR = horizontal system. + +Use [`Astronomy_HorizonFromVector`](#Astronomy_HorizonFromVector) to convert the return value to a traditional altitude/azimuth pair. + + + +**Returns:** A rotation matrix that converts EQD to HOR at `time` and for `observer`. The components of the horizontal vector are: x = north, y = west, z = zenith (straight up from the observer). These components are chosen so that the "right-hand rule" works for the vector and so that north represents the direction where azimuth = 0. + + + +| Type | Parameter | Description | +| --- | --- | --- | +| [`astro_time_t`](#astro_time_t) | `time` | The date and time at which the Earth's equator applies. | +| [`astro_observer_t`](#astro_observer_t) | `observer` | A location near the Earth's mean sea level that defines the observer's horizon. | + + + + --- @@ -1162,7 +1226,7 @@ Given a Cartesian vector, returns latitude, longitude, and distance. | Type | Parameter | Description | | --- | --- | --- | -| [`astro_vector_t`](#astro_vector_t) | `vector` | Cartesian vector to be converted to spherical coordinates | +| [`astro_vector_t`](#astro_vector_t) | `vector` | Cartesian vector to be converted to spherical coordinates. | diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 86b81460..affe55cb 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -5430,7 +5430,7 @@ astro_vector_t Astronomy_VectorFromSphere(astro_spherical_t sphere, astro_time_t * Given a Cartesian vector, returns latitude, longitude, and distance. * * @param vector - * Cartesian vector to be converted to spherical coordinates + * Cartesian vector to be converted to spherical coordinates. * * @return * Spherical coordinates that are equivalent to the given vector. @@ -5474,6 +5474,100 @@ astro_spherical_t Astronomy_SphereFromVector(astro_vector_t vector) return sphere; } +/** + * @brief Converts Cartesian coordinates to horizontal coordinates. + * + * Given a horizontal Cartesian vector, returns horizontal azimuth and altitude. + * + * *IMPORTANT:* This function differs from #Astronomy_SphereFromVector in two ways: + * - `Astronomy_SphereFromVector` returns a `lon` value that represents azimuth defined counterclockwise + * from north (e.g., west = +90), but this function represents a clockwise rotation + * (e.g., east = +90). The difference is because `Astronomy_SphereFromVector` is intended + * to preserve the vector "right-hand rule", while this function defines azimuth in a more + * traditional way as used in navigation and cartography. + * - This function optionally corrects for atmospheric refraction, while `Astronomy_SphereFromVector` + * does not. + * + * The returned structure will contain the azimuth in `lon`. + * It is measured in degrees clockwise from north: east = +90 degrees, west = +270 degrees. + * + * The altitude will be stored in `lat`. + * + * The distance to the observed object is stored in `dist`, + * and is expressed in astronomical units (AU). + * + * @param vector + * Cartesian vector to be converted to horizontal coordinates. + * + * @param refraction + * `REFRACTION_NORMAL`: correct altitude for atmospheric refraction (recommended). + * `REFRACTION_NONE`: no atmospheric refraction correction is performed. + * `REFRACTION_JPLHOR`: for JPL Horizons compatibility testing only; not recommended for normal use. + * + * @return + * If successful, `status` hold `ASTRO_SUCCESS` and the other fields are valid as described above. + * Otherwise `status` holds an error code and the other fields are undefined. + */ +astro_spherical_t Astronomy_HorizonFromVector(astro_vector_t vector, astro_refraction_t refraction) +{ + astro_spherical_t sphere; + + sphere = Astronomy_SphereFromVector(vector); + if (sphere.status == ASTRO_SUCCESS) + { + /* Convert azimuth from counterclockwise-from-north to clockwise-from-north. */ + sphere.lon = 360.0 - sphere.lon; + if (sphere.lon >= 360.0) + sphere.lon -= 360.0; + else if (sphere.lon < 0.0) + sphere.lon += 360.0; + + if (refraction == REFRACTION_NONE) + { + /* No atmospheric correction needed; do nothing. */ + } + else if (refraction == REFRACTION_NORMAL || refraction == REFRACTION_JPLHOR) + { + double zd, refr, hd; + + zd = 90.0 - sphere.lat; + + // http://extras.springer.com/1999/978-1-4471-0555-8/chap4/horizons/horizons.pdf + // JPL Horizons says it uses refraction algorithm from + // Meeus "Astronomical Algorithms", 1991, p. 101-102. + // I found the following Go implementation: + // https://github.com/soniakeys/meeus/blob/master/v3/refraction/refract.go + // This is a translation from the function "Saemundsson" there. + // I found experimentally that JPL Horizons clamps the angle to 1 degree below the horizon. + // This is important because the 'refr' formula below goes crazy near hd = -5.11. + hd = sphere.lat; + if (hd < -1.0) + hd = -1.0; + + refr = (1.02 / tan((hd+10.3/(hd+5.11))*DEG2RAD)) / 60.0; + + if (refraction == REFRACTION_NORMAL && zd > 91.0) + { + // In "normal" mode we gradually reduce refraction toward the nadir + // so that we never get an altitude angle less than -90 degrees. + // When horizon angle is -1 degrees, zd = 91, and the factor is exactly 1. + // As zd approaches 180 (the nadir), the fraction approaches 0 linearly. + refr *= (180.0 - zd) / 89.0; + } + + sphere.lat += refr; + } + else + { + /* Invalid refraction option. */ + return SphereError(ASTRO_INVALID_PARAMETER); + } + } + + return sphere; +} + + /** * @brief * Applies a rotation to a vector, yielding a rotated vector. @@ -5507,28 +5601,6 @@ astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t } -#if 0 -/** - * @brief - * Calculates a rotation matrix from equatorial of-date system to equatorial J2000. - * - * This is one of the family of MakeRotation functions that convert - * from one orientation system to another. - * Source: EQD = equatorial system, using equator of date. - * Target: EQJ = equatorial system, using equator at J2000 epoch. - * - * @param time - * The date and time defining the equator of the source orientation. - * - * @return - * A rotation matrix that converts EQD to EQJ. - */ -astro_rotation_t Astronomy_MakeRotation_EQD_EQJ(astro_time_t time) -{ - -} -#endif - /** * @brief * Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL). @@ -5629,6 +5701,70 @@ astro_rotation_t Astronomy_Rotation_EQD_EQJ(astro_time_t time) return Astronomy_CombineRotation(nut, prec); } + +/** + * @brief + * Calculates a rotation matrix from equatorial of-date (EQD) to horizontal (HOR). + * + * This is one of the family of functions that returns a rotation matrix + * for converting from one orientation to another. + * Source: EQD = equatorial system, using equator of the specified date/time. + * Target: HOR = horizontal system. + * + * Use #Astronomy_HorizonFromVector to convert the return value + * to a traditional altitude/azimuth pair. + * + * @param time + * The date and time at which the Earth's equator applies. + * + * @param observer + * A location near the Earth's mean sea level that defines the observer's horizon. + * + * @return + * A rotation matrix that converts EQD to HOR at `time` and for `observer`. + * The components of the horizontal vector are: + * x = north, y = west, z = zenith (straight up from the observer). + * These components are chosen so that the "right-hand rule" works for the vector + * and so that north represents the direction where azimuth = 0. + */ +astro_rotation_t Astronomy_Rotation_EQD_HOR(astro_time_t time, astro_observer_t observer) +{ + astro_rotation_t rot; + double uze[3], une[3], uwe[3]; + double uz[3], un[3], uw[3]; + double spin_angle; + + double sinlat = sin(observer.latitude * DEG2RAD); + double coslat = cos(observer.latitude * DEG2RAD); + double sinlon = sin(observer.longitude * DEG2RAD); + double coslon = cos(observer.longitude * DEG2RAD); + + uze[0] = coslat * coslon; + uze[1] = coslat * sinlon; + uze[2] = sinlat; + + une[0] = -sinlat * coslon; + une[1] = -sinlat * sinlon; + une[2] = coslat; + + uwe[0] = sinlon; + uwe[1] = -coslon; + uwe[2] = 0.0; + + spin_angle = -15.0 * sidereal_time(&time); + spin(spin_angle, uze, uz); + spin(spin_angle, une, un); + spin(spin_angle, uwe, uw); + + rot.rot[0][0] = un[0]; rot.rot[1][0] = un[1]; rot.rot[2][0] = un[2]; + rot.rot[0][1] = uw[0]; rot.rot[1][1] = uw[1]; rot.rot[2][1] = uw[2]; + rot.rot[0][2] = uz[0]; rot.rot[1][2] = uz[1]; rot.rot[2][2] = uz[2]; + + rot.status = ASTRO_SUCCESS; + return rot; +} + + #ifdef __cplusplus } #endif diff --git a/source/c/astronomy.h b/source/c/astronomy.h index 4fb7290e..a7d3758c 100644 --- a/source/c/astronomy.h +++ b/source/c/astronomy.h @@ -616,11 +616,12 @@ astro_rotation_t Astronomy_InverseRotation(astro_rotation_t rotation); astro_rotation_t Astronomy_CombineRotation(astro_rotation_t a, astro_rotation_t b); astro_vector_t Astronomy_VectorFromSphere(astro_spherical_t sphere, astro_time_t time); astro_spherical_t Astronomy_SphereFromVector(astro_vector_t vector); +astro_spherical_t Astronomy_HorizonFromVector(astro_vector_t vector, astro_refraction_t refraction); astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t vector); astro_rotation_t Astronomy_Rotation_EQD_EQJ(astro_time_t time); //astro_rotation_t Astronomy_Rotation_EQD_ECL(astro_time_t time); -//astro_rotation_t Astronomy_Rotation_EQD_HOR(astro_time_t time, astro_observer_t observer); +astro_rotation_t Astronomy_Rotation_EQD_HOR(astro_time_t time, astro_observer_t observer); astro_rotation_t Astronomy_Rotation_EQJ_EQD(astro_time_t time); astro_rotation_t Astronomy_Rotation_EQJ_ECL(void); //astro_rotation_t Astronomy_Rotation_EQJ_HOR(astro_time_t time, astro_observer_t observer);