From 55d4cc3dde4022fff95fb7139a7df9457eec2da2 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 9 Dec 2019 16:53:02 -0500 Subject: [PATCH] C: Added Astronomy_Rotation_EQD_ECL, ECL_EQD, HOR_ECL, ECL_HOR. These are the final 4 rotation functions to complete every possible coordinate transform. Finished unit tests of verifying that all triangular cycles of transitive rotation are consistent. --- generate/ctest.c | 55 +++++++++++++----- generate/template/astronomy.c | 104 ++++++++++++++++++++++++++++++++++ source/c/README.md | 100 ++++++++++++++++++++++++++++++++ source/c/astronomy.c | 104 ++++++++++++++++++++++++++++++++++ source/c/astronomy.h | 8 +-- 5 files changed, 354 insertions(+), 17 deletions(-) diff --git a/generate/ctest.c b/generate/ctest.c index de0eed91..110d113e 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -1775,7 +1775,7 @@ static int CheckUnitVector(int lnum, const char *name, astro_rotation_t r, int i sum += x*x; } x = fabs(sum - 1.0); - if (x > 5e-16) + if (x > 1.0e-15) { fprintf(stderr, "ERROR(%s line %d): unit error = %lg for i0=%d, j0=%d, di=%d, dj=%d\n", name, lnum, x, i0, j0, di, dj); return 1; @@ -2349,6 +2349,24 @@ fail: #define CHECK_INVERSE(a,b) CHECK(CheckInverse(#a, #b, a, b)) +static int CheckCycle( + const char *cyclename, + astro_rotation_t arot, astro_rotation_t brot, astro_rotation_t crot) +{ + int error; + astro_rotation_t xrot; + + xrot = Astronomy_CombineRotation(arot, brot); + CHECK_ROTMAT(xrot); + xrot = Astronomy_InverseRotation(xrot); + CHECK(CompareMatrices(cyclename, crot, xrot, 2.0e-15)); + error = 0; +fail: + return error; +} + +#define CHECK_CYCLE(a,b,c) CHECK(CheckCycle(("CheckCycle: " #a ", " #b ", " #c), a, b, c)) + static int Test_RotRoundTrip(void) { int error; @@ -2358,18 +2376,15 @@ static int Test_RotRoundTrip(void) astro_rotation_t eqj_hor, hor_eqj; astro_rotation_t eqj_eqd, eqd_eqj; astro_rotation_t hor_eqd, eqd_hor; - //astro_rotation_t eqd_ecl, ecl_eqd; - //astro_rotation_t hor_ecl, ecl_hor; + astro_rotation_t eqd_ecl, ecl_eqd; + astro_rotation_t hor_ecl, ecl_hor; time = Astronomy_MakeTime(2067, 5, 30, 14, 45, 0.0); observer = Astronomy_MakeObserver(+28.0, -82.0, 0.0); /* - In each round trip, we calculate a forward rotation and a backward rotation. + In each round trip, calculate a forward rotation and a backward rotation. Verify the two are inverse matrices. - Then verify that combining different sequences of rotations result - in the expected combination. - For example, (EQJ ==> HOR ==> ECL) must be the same matrix as (EQJ ==> ECL). */ /* Round trip #1: EQJ <==> EQD. */ @@ -2393,14 +2408,28 @@ static int Test_RotRoundTrip(void) CHECK_INVERSE(eqd_hor, hor_eqd); /* Round trip #5: EQD <==> ECL. */ - //eqd_ecl = Astronomy_Rotation_EQD_ECL(time, observer); - //ecl_eqd = Astronomy_Rotation_ECL_EQD(time, observer); - //CHECK_INVERSE(eqd_ecl, ecl_eqd); + eqd_ecl = Astronomy_Rotation_EQD_ECL(time); + ecl_eqd = Astronomy_Rotation_ECL_EQD(time); + CHECK_INVERSE(eqd_ecl, ecl_eqd); /* Round trip #6: HOR <==> ECL. */ - //hor_ecl = Astronomy_Rotation_HOR_ECL(time, observer); - //ecl_hor = Astronomy_Rotation_ECL_HOR(time, observer); - //CHECK_INVERSE(hor_ecl, ecl_hor); + hor_ecl = Astronomy_Rotation_HOR_ECL(time, observer); + ecl_hor = Astronomy_Rotation_ECL_HOR(time, observer); + CHECK_INVERSE(hor_ecl, ecl_hor); + + /* + Verify that combining different sequences of rotations result + in the expected combination. + For example, (EQJ ==> HOR ==> ECL) must be the same matrix as (EQJ ==> ECL). + Each of these is a "triangle" of relationships between 3 orientations. + There are 4 possible ways to pick 3 orientations from the 4 to form a triangle. + Because we have just proved that each transformation is reversible, + we only need to verify the triangle in one cyclic direction. + */ + CHECK_CYCLE(eqj_ecl, ecl_eqd, eqd_eqj); /* excluded corner = HOR */ + CHECK_CYCLE(eqj_hor, hor_ecl, ecl_eqj); /* excluded corner = EQD */ + CHECK_CYCLE(eqj_hor, hor_eqd, eqd_eqj); /* excluded corner = ECL */ + CHECK_CYCLE(ecl_eqd, eqd_hor, hor_ecl); /* excluded corner = EQJ */ printf("Test_RotRoundTrip: PASS\n"); error = 0; diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 3c197456..d23d22c8 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -4617,6 +4617,110 @@ astro_rotation_t Astronomy_Rotation_EQJ_HOR(astro_time_t time, astro_observer_t } +/** + * @brief + * Calculates a rotation matrix from equatorial of-date (EQD) to ecliptic J2000 (ECL). + * + * 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 date. + * Target: ECL = ecliptic system, using equator at J2000 epoch. + * + * @param time + * The date and time of the source equator. + * + * @return + * A rotation matrix that converts EQD to ECL. + */ +astro_rotation_t Astronomy_Rotation_EQD_ECL(astro_time_t time) +{ + astro_rotation_t eqd_eqj; + astro_rotation_t eqj_ecl; + + eqd_eqj = Astronomy_Rotation_EQD_EQJ(time); + eqj_ecl = Astronomy_Rotation_EQJ_ECL(); + return Astronomy_CombineRotation(eqd_eqj, eqj_ecl); +} + + +/** + * @brief + * Calculates a rotation matrix from ecliptic J2000 (ECL) to equatorial of-date (EQD). + * + * This is one of the family of functions that returns a rotation matrix + * for converting from one orientation to another. + * Source: ECL = ecliptic system, using equator at J2000 epoch. + * Target: EQD = equatorial system, using equator of date. + * + * @param time + * The date and time of the desired equator. + * + * @return + * A rotation matrix that converts ECL to EQD. + */ +astro_rotation_t Astronomy_Rotation_ECL_EQD(astro_time_t time) +{ + astro_rotation_t rot = Astronomy_Rotation_EQD_ECL(time); + return Astronomy_InverseRotation(rot); +} + +/** + * @brief + * Calculates a rotation matrix from ecliptic J2000 (ECL) to horizontal (HOR). + * + * This is one of the family of functions that returns a rotation matrix + * for converting from one orientation to another. + * Source: ECL = ecliptic system, using equator at J2000 epoch. + * Target: HOR = horizontal system. + * + * Use #Astronomy_HorizonFromVector to convert the return value + * to a traditional altitude/azimuth pair. + * + * @param time + * The date and time of the desired horizontal orientation. + * + * @param observer + * A location near the Earth's mean sea level that defines the observer's horizon. + * + * @return + * A rotation matrix that converts ECL 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_ECL_HOR(astro_time_t time, astro_observer_t observer) +{ + astro_rotation_t ecl_eqd = Astronomy_Rotation_ECL_EQD(time); + astro_rotation_t eqd_hor = Astronomy_Rotation_EQD_HOR(time, observer); + return Astronomy_CombineRotation(ecl_eqd, eqd_hor); +} + +/** + * @brief + * Calculates a rotation matrix from horizontal (HOR) to ecliptic J2000 (ECL). + * + * This is one of the family of functions that returns a rotation matrix + * for converting from one orientation to another. + * Source: HOR = horizontal system. + * Target: ECL = ecliptic system, using equator at J2000 epoch. + * + * @param time + * The date and time of the horizontal observation. + * + * @param observer + * The location of the horizontal observer. + * + * @return + * A rotation matrix that converts HOR to ECL. + */ +astro_rotation_t Astronomy_Rotation_HOR_ECL(astro_time_t time, astro_observer_t observer) +{ + astro_rotation_t rot = Astronomy_Rotation_ECL_HOR(time, observer); + return Astronomy_InverseRotation(rot); +} + + #ifdef __cplusplus } #endif diff --git a/source/c/README.md b/source/c/README.md index e532ae21..ea4535a8 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -760,6 +760,30 @@ This function transforms a vector in one orientation to a vector in another orie +--- + + +### Astronomy_Rotation_ECL_EQD(time) ⇒ [`astro_rotation_t`](#astro_rotation_t) + +**Calculates a rotation matrix from ecliptic J2000 (ECL) to equatorial of-date (EQD).** + + + +This is one of the family of functions that returns a rotation matrix for converting from one orientation to another. Source: ECL = ecliptic system, using equator at J2000 epoch. Target: EQD = equatorial system, using equator of date. + + + +**Returns:** A rotation matrix that converts ECL to EQD. + + + +| Type | Parameter | Description | +| --- | --- | --- | +| [`astro_time_t`](#astro_time_t) | `time` | The date and time of the desired equator. | + + + + --- @@ -777,6 +801,57 @@ This is one of the family of functions that returns a rotation matrix for conver +--- + + +### Astronomy_Rotation_ECL_HOR(time, observer) ⇒ [`astro_rotation_t`](#astro_rotation_t) + +**Calculates a rotation matrix from ecliptic J2000 (ECL) to horizontal (HOR).** + + + +This is one of the family of functions that returns a rotation matrix for converting from one orientation to another. Source: ECL = ecliptic system, using equator at J2000 epoch. 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 ECL 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 of the desired horizontal orientation. | +| [`astro_observer_t`](#astro_observer_t) | `observer` | A location near the Earth's mean sea level that defines the observer's horizon. | + + + + +--- + + +### Astronomy_Rotation_EQD_ECL(time) ⇒ [`astro_rotation_t`](#astro_rotation_t) + +**Calculates a rotation matrix from equatorial of-date (EQD) to ecliptic J2000 (ECL).** + + + +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 date. Target: ECL = ecliptic system, using equator at J2000 epoch. + + + +**Returns:** A rotation matrix that converts EQD to ECL. + + + +| Type | Parameter | Description | +| --- | --- | --- | +| [`astro_time_t`](#astro_time_t) | `time` | The date and time of the source equator. | + + + + --- @@ -896,6 +971,31 @@ Use [`Astronomy_HorizonFromVector`](#Astronomy_HorizonFromVector) to convert the +--- + + +### Astronomy_Rotation_HOR_ECL(time, observer) ⇒ [`astro_rotation_t`](#astro_rotation_t) + +**Calculates a rotation matrix from horizontal (HOR) to ecliptic J2000 (ECL).** + + + +This is one of the family of functions that returns a rotation matrix for converting from one orientation to another. Source: HOR = horizontal system. Target: ECL = ecliptic system, using equator at J2000 epoch. + + + +**Returns:** A rotation matrix that converts HOR to ECL. + + + +| Type | Parameter | Description | +| --- | --- | --- | +| [`astro_time_t`](#astro_time_t) | `time` | The date and time of the horizontal observation. | +| [`astro_observer_t`](#astro_observer_t) | `observer` | The location of the horizontal observer. | + + + + --- diff --git a/source/c/astronomy.c b/source/c/astronomy.c index ecfd14b0..d0c8bbe3 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -5856,6 +5856,110 @@ astro_rotation_t Astronomy_Rotation_EQJ_HOR(astro_time_t time, astro_observer_t } +/** + * @brief + * Calculates a rotation matrix from equatorial of-date (EQD) to ecliptic J2000 (ECL). + * + * 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 date. + * Target: ECL = ecliptic system, using equator at J2000 epoch. + * + * @param time + * The date and time of the source equator. + * + * @return + * A rotation matrix that converts EQD to ECL. + */ +astro_rotation_t Astronomy_Rotation_EQD_ECL(astro_time_t time) +{ + astro_rotation_t eqd_eqj; + astro_rotation_t eqj_ecl; + + eqd_eqj = Astronomy_Rotation_EQD_EQJ(time); + eqj_ecl = Astronomy_Rotation_EQJ_ECL(); + return Astronomy_CombineRotation(eqd_eqj, eqj_ecl); +} + + +/** + * @brief + * Calculates a rotation matrix from ecliptic J2000 (ECL) to equatorial of-date (EQD). + * + * This is one of the family of functions that returns a rotation matrix + * for converting from one orientation to another. + * Source: ECL = ecliptic system, using equator at J2000 epoch. + * Target: EQD = equatorial system, using equator of date. + * + * @param time + * The date and time of the desired equator. + * + * @return + * A rotation matrix that converts ECL to EQD. + */ +astro_rotation_t Astronomy_Rotation_ECL_EQD(astro_time_t time) +{ + astro_rotation_t rot = Astronomy_Rotation_EQD_ECL(time); + return Astronomy_InverseRotation(rot); +} + +/** + * @brief + * Calculates a rotation matrix from ecliptic J2000 (ECL) to horizontal (HOR). + * + * This is one of the family of functions that returns a rotation matrix + * for converting from one orientation to another. + * Source: ECL = ecliptic system, using equator at J2000 epoch. + * Target: HOR = horizontal system. + * + * Use #Astronomy_HorizonFromVector to convert the return value + * to a traditional altitude/azimuth pair. + * + * @param time + * The date and time of the desired horizontal orientation. + * + * @param observer + * A location near the Earth's mean sea level that defines the observer's horizon. + * + * @return + * A rotation matrix that converts ECL 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_ECL_HOR(astro_time_t time, astro_observer_t observer) +{ + astro_rotation_t ecl_eqd = Astronomy_Rotation_ECL_EQD(time); + astro_rotation_t eqd_hor = Astronomy_Rotation_EQD_HOR(time, observer); + return Astronomy_CombineRotation(ecl_eqd, eqd_hor); +} + +/** + * @brief + * Calculates a rotation matrix from horizontal (HOR) to ecliptic J2000 (ECL). + * + * This is one of the family of functions that returns a rotation matrix + * for converting from one orientation to another. + * Source: HOR = horizontal system. + * Target: ECL = ecliptic system, using equator at J2000 epoch. + * + * @param time + * The date and time of the horizontal observation. + * + * @param observer + * The location of the horizontal observer. + * + * @return + * A rotation matrix that converts HOR to ECL. + */ +astro_rotation_t Astronomy_Rotation_HOR_ECL(astro_time_t time, astro_observer_t observer) +{ + astro_rotation_t rot = Astronomy_Rotation_ECL_HOR(time, observer); + return Astronomy_InverseRotation(rot); +} + + #ifdef __cplusplus } #endif diff --git a/source/c/astronomy.h b/source/c/astronomy.h index b7a1f734..b9534caa 100644 --- a/source/c/astronomy.h +++ b/source/c/astronomy.h @@ -620,17 +620,17 @@ astro_spherical_t Astronomy_HorizonFromVector(astro_vector_t vector, astro_refra 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_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_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); -//astro_rotation_t Astronomy_Rotation_ECL_EQD(astro_time_t time); +astro_rotation_t Astronomy_Rotation_ECL_EQD(astro_time_t time); astro_rotation_t Astronomy_Rotation_ECL_EQJ(void); -//astro_rotation_t Astronomy_Rotation_ECL_HOR(astro_time_t time, astro_observer_t observer); +astro_rotation_t Astronomy_Rotation_ECL_HOR(astro_time_t time, astro_observer_t observer); astro_rotation_t Astronomy_Rotation_HOR_EQD(astro_time_t time, astro_observer_t observer); astro_rotation_t Astronomy_Rotation_HOR_EQJ(astro_time_t time, astro_observer_t observer); -//astro_rotation_t Astronomy_Rotation_HOR_ECL(astro_time_t time, astro_observer_t observer); +astro_rotation_t Astronomy_Rotation_HOR_ECL(astro_time_t time, astro_observer_t observer); #ifdef __cplusplus }