From d572edb10fe4f3c3779c2644233ae01de85abeda Mon Sep 17 00:00:00 2001
From: Don Cross <cosinekitty@gmail.com>
Date: Sun, 15 Dec 2019 20:41:03 -0500
Subject: [PATCH] Python: Implemented Rotation_ECL_EQJ.

---
 generate/template/astronomy.py | 23 +++++++++++++++++++
 generate/test.py               | 41 ++++++++++++++++++++++++++++++++++
 source/python/README.md        | 15 +++++++++++++
 source/python/astronomy.py     | 23 +++++++++++++++++++
 4 files changed, 102 insertions(+)

diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py
index 75edefa3..995044fb 100644
--- a/generate/template/astronomy.py
+++ b/generate/template/astronomy.py
@@ -3136,3 +3136,26 @@ def Rotation_EQJ_ECL():
         [ 0, +c, -s],
         [ 0, +s, +c]
     ])
+
+
+def Rotation_ECL_EQJ():
+    """Calculates a rotation matrix from ecliptic J2000 (ECL) to equatorial J2000 (EQJ).
+
+    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: EQJ = equatorial system, using equator at J2000 epoch.
+
+    Returns
+    -------
+    RotationMatrix
+        A rotation matrix that converts ECL to EQJ.
+    """
+    # ob = mean obliquity of the J2000 ecliptic = 0.40909260059599012 radians.
+    c = 0.9174821430670688    # cos(ob)
+    s = 0.3977769691083922    # sin(ob)
+    return RotationMatrix([
+        [ 1,  0,  0],
+        [ 0, +c, +s],
+        [ 0, -s, +c]
+    ])
diff --git a/generate/test.py b/generate/test.py
index c4274be1..112f7b3f 100755
--- a/generate/test.py
+++ b/generate/test.py
@@ -736,6 +736,7 @@ def CompareMatrices(caller, a, b, tolerance):
                 print('ERROR({}): matrix[{}][{}] = {}, expected {}, diff {}'.format(caller, i, j, a.rot[i][j], b.rot[i][j], diff))
                 sys.exit(1)
 
+
 def Rotation_MatrixInverse():
     a = astronomy.RotationMatrix([
         [1, 4, 7],
@@ -750,6 +751,7 @@ def Rotation_MatrixInverse():
     b = astronomy.InverseRotation(a)
     CompareMatrices('Rotation_MatrixInverse', b, v, 0)
 
+
 def Rotation_MatrixMultiply():
     a = astronomy.RotationMatrix([
         [1, 4, 7],
@@ -772,9 +774,48 @@ def Rotation_MatrixMultiply():
     c = astronomy.CombineRotation(b, a)
     CompareMatrices('Rotation_MatrixMultiply', c, v, 0)
 
+
+def VectorDiff(a, b):
+    dx = a.x - b.x
+    dy = a.y - b.y
+    dz = a.z - b.z
+    return math.sqrt(dx*dx + dy*dy + dz*dz)
+
+
+def Test_EQJ_ECL():
+    r = astronomy.Rotation_EQJ_ECL()
+    # Calculate heliocentric Earth position at a test time.
+    time = astronomy.Time.Make(2019, 12, 8, 19, 39, 15)
+    ev = astronomy.HelioVector(astronomy.Body.Earth, time)
+
+    # Use the existing astronomy.Ecliptic() to calculate ecliptic vector and angles.
+    ecl = astronomy.Ecliptic(ev)
+    print('Test_EQJ_ECL ecl = ({}, {}, {})'.format(ecl.ex, ecl.ey, ecl.ez))
+
+    # Now compute the same vector via rotation matrix.
+    ee = astronomy.RotateVector(r, ev)
+    dx = ee.x - ecl.ex
+    dy = ee.y - ecl.ey
+    dz = ee.z - ecl.ez
+    diff = math.sqrt(dx*dx + dy*dy + dz*dz)
+    print('Test_EQJ_ECL ee = ({}, {}, {}); diff = {}'.format(ee.x, ee.y, ee.z, diff))
+    if diff > 1.0e-16:
+        print('Test_EQJ_ECL: EXCESSIVE VECTOR ERROR')
+        sys.exit(1)
+
+    # Reverse the test: go from ecliptic back to equatorial.
+    ir = astronomy.Rotation_ECL_EQJ()
+    et = astronomy.RotateVector(ir, ee)
+    idiff = VectorDiff(et, ev)
+    print('Test_EQJ_ECL ev diff = {}'.format(idiff))
+    if idiff > 1.0e-16:
+        print('Test_EQJ_ECL: EXCESSIVE REVERSE ROTATION ERROR')
+        sys.exit(1)
+
 def Test_Rotation():
     Rotation_MatrixInverse()
     Rotation_MatrixMultiply()
+    Test_EQJ_ECL()
     print('Python Test_Rotation: PASS')
     return 0
 
diff --git a/source/python/README.md b/source/python/README.md
index 0a070750..ba6e2231 100644
--- a/source/python/README.md
+++ b/source/python/README.md
@@ -947,6 +947,21 @@ A vector in the orientation specified by `rotation`.
 
 ---
 
+<a name="Rotation_ECL_EQJ"></a>
+### Rotation_ECL_EQJ()
+
+**Calculates a rotation matrix from ecliptic J2000 (ECL) to equatorial J2000 (EQJ).**
+
+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: EQJ = equatorial system, using equator at J2000 epoch.
+
+### Returns: RotationMatrix
+A rotation matrix that converts ECL to EQJ.
+
+---
+
 <a name="Rotation_EQJ_ECL"></a>
 ### Rotation_EQJ_ECL()
 
diff --git a/source/python/astronomy.py b/source/python/astronomy.py
index ebbf02a4..43270a20 100644
--- a/source/python/astronomy.py
+++ b/source/python/astronomy.py
@@ -5197,3 +5197,26 @@ def Rotation_EQJ_ECL():
         [ 0, +c, -s],
         [ 0, +s, +c]
     ])
+
+
+def Rotation_ECL_EQJ():
+    """Calculates a rotation matrix from ecliptic J2000 (ECL) to equatorial J2000 (EQJ).
+
+    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: EQJ = equatorial system, using equator at J2000 epoch.
+
+    Returns
+    -------
+    RotationMatrix
+        A rotation matrix that converts ECL to EQJ.
+    """
+    # ob = mean obliquity of the J2000 ecliptic = 0.40909260059599012 radians.
+    c = 0.9174821430670688    # cos(ob)
+    s = 0.3977769691083922    # sin(ob)
+    return RotationMatrix([
+        [ 1,  0,  0],
+        [ 0, +c, +s],
+        [ 0, -s, +c]
+    ])