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

@@ -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;

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

source/c/README.md

@@ -698,6 +698,31 @@ After calling [`Astronomy_SearchMoonQuarter`](#Astronomy_SearchMoonQuarter), thi
 
 
 
+---
+
+<a name="Astronomy_RotateVector"></a>
+### Astronomy_RotateVector(rotation, vector) &#8658; [`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. | 
+
+
+
+
 ---
 
 <a name="Astronomy_Search"></a>

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

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);