From f493707e6dabe25fa41373df852f623696bccd62 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sun, 8 Dec 2019 15:13:55 -0500 Subject: [PATCH] C: Implemented Astronomy_RotateVector. Astronomy_RotateVector translates a vector in one orientation to another orientation, as specified by a rotation matrix. I will use this to implement all the coordinate transforms among EQJ, EQD, ECL, HOR. --- generate/ctest.c | 109 ++++++++++++++++++++++++++++++++++ generate/template/astronomy.c | 31 ++++++++++ source/c/README.md | 25 ++++++++ source/c/astronomy.c | 31 ++++++++++ source/c/astronomy.h | 1 + 5 files changed, 197 insertions(+) diff --git a/generate/ctest.c b/generate/ctest.c index 8c3cd6e0..5a15da85 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -11,6 +11,7 @@ #include "astronomy.h" #define PI 3.14159265358979323846 +static const double DEG2RAD = 0.017453292519943296; #define CHECK(x) do{if(0 != (error = (x))) goto fail;}while(0) @@ -1917,6 +1918,106 @@ static int TestAnglesFromVector(double lat, double lon, double x, double y, doub return 0; } +static astro_rotation_t RotateX(double rot) +{ + astro_rotation_t m; + double ca, sa; + + /* https://en.wikipedia.org/wiki/Rotation_matrix */ + + ca = cos(rot * DEG2RAD); + sa = sin(rot * DEG2RAD); + m.status = ASTRO_SUCCESS; + m.rot[0][0] = 1.0; m.rot[1][0] = 0.0; m.rot[2][0] = 0.0; + m.rot[0][1] = 0.0; m.rot[1][1] = +ca; m.rot[2][1] = -sa; + m.rot[0][2] = 0.0; m.rot[1][2] = +sa; m.rot[2][2] = +ca; + + return m; +} + +static astro_rotation_t RotateY(double rot) +{ + astro_rotation_t m; + double ca, sa; + + /* https://en.wikipedia.org/wiki/Rotation_matrix */ + + ca = cos(rot * DEG2RAD); + sa = sin(rot * DEG2RAD); + m.status = ASTRO_SUCCESS; + m.rot[0][0] = +ca; m.rot[1][0] = 0.0; m.rot[2][0] = +sa; + m.rot[0][1] = 0.0; m.rot[1][1] = 1.0; m.rot[2][1] = 0.0; + m.rot[0][2] = -sa; m.rot[1][2] = 0.0; m.rot[2][2] = +ca; + + return m; +} + +static astro_rotation_t RotateZ(double rot) +{ + astro_rotation_t m; + double ca, sa; + + /* https://en.wikipedia.org/wiki/Rotation_matrix */ + + ca = cos(rot * DEG2RAD); + sa = sin(rot * DEG2RAD); + m.status = ASTRO_SUCCESS; + m.rot[0][0] = +ca; m.rot[1][0] = -sa; m.rot[2][0] = 0.0; + m.rot[0][1] = +sa; m.rot[1][1] = +ca; m.rot[2][1] = 0.0; + m.rot[0][2] = 0.0; m.rot[1][2] = 0.0; m.rot[2][2] = 1.0; + + return m; +} + +static int TestSpin( + double xrot, double yrot, double zrot, + double sx, double sy, double sz, + double tx, double ty, double tz) +{ + astro_rotation_t mx, my, mz, m; + astro_vector_t sv, tv; + double diff, dx, dy, dz; + + /* https://en.wikipedia.org/wiki/Rotation_matrix */ + + mx = RotateX(xrot); + my = RotateY(yrot); + mz = RotateZ(zrot); + m = Astronomy_CombineRotation(mx, my); + m = Astronomy_CombineRotation(m, mz); + + sv.status = ASTRO_SUCCESS; + sv.x = sx; + sv.y = sy; + sv.z = sz; + sv.t = Astronomy_MakeTime(2019, 5, 5, 12, 30, 0.0); + + tv = Astronomy_RotateVector(m, sv); + if (tv.status != ASTRO_SUCCESS) + { + fprintf(stderr, "ERROR(TestSpin): tv.status = %d\n", tv.status); + return 1; + } + + if (tv.t.ut != sv.t.ut) + { + fprintf(stderr, "ERROR(TestSpin): tv time != sv time\n"); + return 1; + } + + dx = tx - tv.x; + dy = ty - tv.y; + dz = tz - tv.z; + diff = sqrt(dx*dx + dy*dy + dz*dz); + printf("TestSpin(xrot=%0.0lf, yrot=%0.0lf, zrot=%0.0lf, sx=%0.0lf, sy=%0.0lf, sz=%0.0lf): diff = %lg\n", xrot, yrot, zrot, sx, sy, sz, diff); + if (diff > 1.0e-15) + { + fprintf(stderr, "TestSpin: EXCESSIVE ERROR\n"); + return 1; + } + return 0; +} + static int RotationTest(void) { int error; @@ -1941,6 +2042,14 @@ static int RotationTest(void) CHECK(TestAnglesFromVector(-90.0, 0.0, 0.0, 0.0, -1.0)); CHECK(TestAnglesFromVector(-30.0, +60.0, 0.43301270189221946, 0.75, -0.5)); + /* Verify rotation of vectors */ + CHECK(TestSpin(0.0, 0.0, 90.0, +1, +2, +3, -2, +1, +3)); + CHECK(TestSpin(0.0, 0.0, 0.0, 1.0, 2.0, -3.0, 1.0, 2.0, -3.0)); + CHECK(TestSpin(90.0, 0.0, 0.0, +1, +2, +3, +1, -3, +2)); + CHECK(TestSpin(0.0, 90.0, 0.0, +1, +2, +3, +3, +2, -1)); + CHECK(TestSpin(180.0, 180.0, 180.0, +1, +2, +3, +1, +2, +3)); + CHECK(TestSpin(0.0, 0.0, -45.0, +1, 0, 0, +0.7071067811865476, -0.7071067811865476, 0)); + error = 0; fail: return error; diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index e1f07179..34173801 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -4239,6 +4239,37 @@ astro_spherical_t Astronomy_SphereFromVector(astro_vector_t vector) return sphere; } +/** + * @brief + * Applies a rotation to a vector, yielding a rotated vector. + * + * This function transforms a vector in one orientation to a vector + * in another orientation. + * + * @param rotation + * A rotation matrix that specifies how the orientation of the vector is to be changed. + * + * @param vector + * The vector whose orientation is to be changed. + * + * @return + * A vector in the orientation specified by `rotation`. + */ +astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t vector) +{ + astro_vector_t target; + + if (rotation.status != ASTRO_SUCCESS || vector.status != ASTRO_SUCCESS) + return VecError(ASTRO_INVALID_PARAMETER, vector.t); + + target.status = ASTRO_SUCCESS; + target.t = vector.t; + target.x = rotation.rot[0][0]*vector.x + rotation.rot[1][0]*vector.y + rotation.rot[2][0]*vector.z; + target.y = rotation.rot[0][1]*vector.x + rotation.rot[1][1]*vector.y + rotation.rot[2][1]*vector.z; + target.z = rotation.rot[0][2]*vector.x + rotation.rot[1][2]*vector.y + rotation.rot[2][2]*vector.z; + + return target; +} #if 0 diff --git a/source/c/README.md b/source/c/README.md index a231ebea..e054b9dd 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -698,6 +698,31 @@ After calling [`Astronomy_SearchMoonQuarter`](#Astronomy_SearchMoonQuarter), thi +--- + + +### Astronomy_RotateVector(rotation, vector) ⇒ [`astro_vector_t`](#astro_vector_t) + +**Applies a rotation to a vector, yielding a rotated vector.** + + + +This function transforms a vector in one orientation to a vector in another orientation. + + + +**Returns:** A vector in the orientation specified by `rotation`. + + + +| Type | Parameter | Description | +| --- | --- | --- | +| [`astro_rotation_t`](#astro_rotation_t) | `rotation` | A rotation matrix that specifies how the orientation of the vector is to be changed. | +| [`astro_vector_t`](#astro_vector_t) | `vector` | The vector whose orientation is to be changed. | + + + + --- diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 1b3ef35b..25cde6ce 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -5478,6 +5478,37 @@ astro_spherical_t Astronomy_SphereFromVector(astro_vector_t vector) return sphere; } +/** + * @brief + * Applies a rotation to a vector, yielding a rotated vector. + * + * This function transforms a vector in one orientation to a vector + * in another orientation. + * + * @param rotation + * A rotation matrix that specifies how the orientation of the vector is to be changed. + * + * @param vector + * The vector whose orientation is to be changed. + * + * @return + * A vector in the orientation specified by `rotation`. + */ +astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t vector) +{ + astro_vector_t target; + + if (rotation.status != ASTRO_SUCCESS || vector.status != ASTRO_SUCCESS) + return VecError(ASTRO_INVALID_PARAMETER, vector.t); + + target.status = ASTRO_SUCCESS; + target.t = vector.t; + target.x = rotation.rot[0][0]*vector.x + rotation.rot[1][0]*vector.y + rotation.rot[2][0]*vector.z; + target.y = rotation.rot[0][1]*vector.x + rotation.rot[1][1]*vector.y + rotation.rot[2][1]*vector.z; + target.z = rotation.rot[0][2]*vector.x + rotation.rot[1][2]*vector.y + rotation.rot[2][2]*vector.z; + + return target; +} #if 0 diff --git a/source/c/astronomy.h b/source/c/astronomy.h index c64b7ab2..546e58f7 100644 --- a/source/c/astronomy.h +++ b/source/c/astronomy.h @@ -616,6 +616,7 @@ 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_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t vector); #if 0 astro_rotation_t Astronomy_Rotation_EQD_EQJ(astro_time_t time);