From 85d9113f7785e1279243571efe31581205bc50e4 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sun, 15 Dec 2019 21:24:32 -0500 Subject: [PATCH] Python: Added more rotation functions and unit tests. --- generate/template/astronomy.py | 146 +++++++++++++++++++++++++++++++++ generate/test.py | 47 +++++++++++ source/python/README.md | 102 +++++++++++++++++++++++ source/python/astronomy.py | 146 +++++++++++++++++++++++++++++++++ 4 files changed, 441 insertions(+) diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 995044fb..473c1c05 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -3022,6 +3022,107 @@ def NextLunarApsis(apsis): raise InternalError() return next + +def VectorFromSphere(sphere, time): + """Converts spherical coordinates to Cartesian coordinates. + + Given spherical coordinates and a time at which they are valid, + returns a vector of Cartesian coordinates. The returned value + includes the time, as required by all `Time` objects. + + Parameters + ---------- + sphere : Spherical + Spherical coordinates to be converted. + time : Time + The time that should be included in the returned vector. + + Returns + ------- + Vector + The vector form of the supplied spherical coordinates. + """ + radlat = math.radians(sphere.lat) + radlon = math.radians(sphere.lon) + rcoslat = sphere.dist * math.cos(radlat) + return Vector( + rcoslat * math.cos(radlon), + rcoslat * math.sin(radlon), + sphere.dist * math.sin(radlat), + time + ) + + +def VectorFromEquator(equ, time): + """Given angular equatorial coordinates in `equ`, calculates equatorial vector. + + Parameters + ---------- + equ : Equatorial + Angular equatorial coordinates to be converted to a vector. + time : Time + The date and time of the observation. This is needed because the returned + vector object requires a valid time value when passed to certain other functions. + + Returns + ------- + Vector + A vector in the equatorial system. + """ + return VectorFromSphere(Spherical(equ.dec, 15.0 * equ.ra, equ.dist), time) + + +def EquatorFromVector(vec): + """Given an equatorial vector, calculates equatorial angular coordinates. + + Parameters + ---------- + vec : Vector + A vector in an equatorial coordinate system. + + Returns + ------- + Equatorial + Angular coordinates expressed in the same equatorial system as `vec`. + """ + sphere = SphereFromVector(vec) + return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist) + + +def SphereFromVector(vector): + """Converts Cartesian coordinates to spherical coordinates. + + Given a Cartesian vector, returns latitude, longitude, and distance. + + Parameters + ---------- + vector : Vector + Cartesian vector to be converted to spherical coordinates. + + Returns + ------- + Spherical + Spherical coordinates that are equivalent to the given vector. + """ + xyproj = vector.x*vector.x + vector.y*vector.y + dist = math.sqrt(xyproj + vector.z*vector.z) + if xyproj == 0.0: + if vector.z == 0.0: + raise Exception('Zero-length vector not allowed.') + lon = 0.0 + if vector.z < 0.0: + lat = -90.0 + else: + lat = +90.0 + else: + lon = math.degrees(math.atan2(vector.y, vector.x)) + if lon < 0.0: + lon += 360.0 + lat = math.degrees(math.atan2(vector.z, math.sqrt(xyproj))) + return Spherical(lat, lon, dist) + + + def InverseRotation(rotation): """Calculates the inverse of a rotation matrix. @@ -3159,3 +3260,48 @@ def Rotation_ECL_EQJ(): [ 0, +c, +s], [ 0, -s, +c] ]) + +def Rotation_EQJ_EQD(time): + """Calculates a rotation matrix from equatorial J2000 (EQJ) 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: EQJ = equatorial system, using equator at J2000 epoch. + Target: EQD = equatorial system, using equator of the specified date/time. + + Parameters + ---------- + time : Time + The date and time at which the Earth's equator defines the target orientation. + + Returns + ------- + RotationMatrix + A rotation matrix that converts EQJ to EQD at `time`. + """ + prec = _precession_rot(0.0, time.tt) + nut = _nutation_rot(time, 0) + return CombineRotation(prec, nut) + + +def Rotation_EQD_EQJ(time): + """Calculates a rotation matrix from equatorial of-date (EQD) 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: EQD = equatorial system, using equator of the specified date/time. + Target: EQJ = equatorial system, using equator at J2000 epoch. + + Parameters + ---------- + time : Time + The date and time at which the Earth's equator defines the source orientation. + + Returns + ------- + RotationMatrix + A rotation matrix that converts EQD at `time` to EQJ. + """ + nut = _nutation_rot(time, 1) + prec = _precession_rot(time.tt, 0.0) + return CombineRotation(nut, prec) diff --git a/generate/test.py b/generate/test.py index 112f7b3f..667a6a7c 100755 --- a/generate/test.py +++ b/generate/test.py @@ -812,10 +812,57 @@ def Test_EQJ_ECL(): print('Test_EQJ_ECL: EXCESSIVE REVERSE ROTATION ERROR') sys.exit(1) + +def Test_EQJ_EQD(body): + # Verify conversion of equatorial J2000 to equatorial of-date, and back. + # Use established functions to calculate spherical coordinates for the body, in both EQJ and EQD. + time = astronomy.Time.Make(2019, 12, 8, 20, 50, 0) + observer = astronomy.Observer(+35, -85, 0) + eq2000 = astronomy.Equator(body, time, observer, False, True) + eqdate = astronomy.Equator(body, time, observer, True, True) + + # Convert EQJ spherical coordinates to vector. + v2000 = astronomy.VectorFromEquator(eq2000, time) + + # Find rotation matrix. + r = astronomy.Rotation_EQJ_EQD(time) + + # Rotate EQJ vector to EQD vector. + vdate = astronomy.RotateVector(r, v2000) + + # Convert vector back to angular equatorial coordinates. + equcheck = astronomy.EquatorFromVector(vdate) + + # Compare the result with the eqdate. + ra_diff = abs(equcheck.ra - eqdate.ra) + dec_diff = abs(equcheck.dec - eqdate.dec) + dist_diff = abs(equcheck.dist - eqdate.dist) + print('Test_EQJ_EQD: {} ra={}, dec={}, dist={}, ra_diff={}, dec_diff={}, dist_diff={}'.format( + body.name, eqdate.ra, eqdate.dec, eqdate.dist, ra_diff, dec_diff, dist_diff + )) + if ra_diff > 1.0e-14 or dec_diff > 1.0e-14 or dist_diff > 4.0e-15: + print('Test_EQJ_EQD: EXCESSIVE ERROR') + sys.exit(1) + + # Perform the inverse conversion back to equatorial J2000 coordinates. + ir = astronomy.Rotation_EQD_EQJ(time) + t2000 = astronomy.RotateVector(ir, vdate) + diff = VectorDiff(t2000, v2000) + print('Test_EQJ_EQD: {} inverse diff = {}'.format(body.name, diff)) + if diff > 3.0e-15: + print('Test_EQJ_EQD: EXCESSIVE INVERSE ERROR') + sys.exit(1) + + def Test_Rotation(): Rotation_MatrixInverse() Rotation_MatrixMultiply() Test_EQJ_ECL() + Test_EQJ_EQD(astronomy.Body.Mercury) + Test_EQJ_EQD(astronomy.Body.Venus) + Test_EQJ_EQD(astronomy.Body.Mars) + Test_EQJ_EQD(astronomy.Body.Jupiter) + Test_EQJ_EQD(astronomy.Body.Saturn) print('Python Test_Rotation: PASS') return 0 diff --git a/source/python/README.md b/source/python/README.md index ba6e2231..c00783dc 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -626,6 +626,20 @@ Equatorial coordinates in the specified frame of reference. --- + +### EquatorFromVector(vec) + +**Given an equatorial vector, calculates equatorial angular coordinates.** + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Vector`](#Vector) | `vec` | A vector in an equatorial coordinate system. | + +### Returns: Equatorial +Angular coordinates expressed in the same equatorial system as `vec`. + +--- + ### GeoMoon(time) @@ -962,6 +976,25 @@ A rotation matrix that converts ECL to EQJ. --- + +### Rotation_EQD_EQJ(time) + +**Calculates a rotation matrix from equatorial of-date (EQD) 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: EQD = equatorial system, using equator of the specified date/time. +Target: EQJ = equatorial system, using equator at J2000 epoch. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Time`](#Time) | `time` | The date and time at which the Earth's equator defines the source orientation. | + +### Returns: RotationMatrix +A rotation matrix that converts EQD at `time` to EQJ. + +--- + ### Rotation_EQJ_ECL() @@ -977,6 +1010,25 @@ A rotation matrix that converts EQJ to ECL. --- + +### Rotation_EQJ_EQD(time) + +**Calculates a rotation matrix from equatorial J2000 (EQJ) 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: EQJ = equatorial system, using equator at J2000 epoch. +Target: EQD = equatorial system, using equator of the specified date/time. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Time`](#Time) | `time` | The date and time at which the Earth's equator defines the target orientation. | + +### Returns: RotationMatrix +A rotation matrix that converts EQJ to EQD at `time`. + +--- + ### Search(func, context, t1, t2, dt_tolerance_seconds) @@ -1334,6 +1386,22 @@ of winter in the southern hemisphere. --- + +### SphereFromVector(vector) + +**Converts Cartesian coordinates to spherical coordinates.** + +Given a Cartesian vector, returns latitude, longitude, and distance. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Vector`](#Vector) | `vector` | Cartesian vector to be converted to spherical coordinates. | + +### Returns: Spherical +Spherical coordinates that are equivalent to the given vector. + +--- + ### SunPosition(time) @@ -1359,3 +1427,37 @@ In fact, the function #Seasons does use this function for that purpose. ### Returns: #EclipticCoordinates The ecliptic coordinates of the Sun using the Earth's true equator of date. +--- + + +### VectorFromEquator(equ, time) + +**Given angular equatorial coordinates in `equ`, calculates equatorial vector.** + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Equatorial`](#Equatorial) | `equ` | Angular equatorial coordinates to be converted to a vector. | +| [`Time`](#Time) | `time` | The date and time of the observation. This is needed because the returned vector object requires a valid time value when passed to certain other functions. | + +### Returns: Vector +A vector in the equatorial system. + +--- + + +### VectorFromSphere(sphere, time) + +**Converts spherical coordinates to Cartesian coordinates.** + +Given spherical coordinates and a time at which they are valid, +returns a vector of Cartesian coordinates. The returned value +includes the time, as required by all `Time` objects. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Spherical`](#Spherical) | `sphere` | Spherical coordinates to be converted. | +| [`Time`](#Time) | `time` | The time that should be included in the returned vector. | + +### Returns: Vector +The vector form of the supplied spherical coordinates. + diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 43270a20..3827237f 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -5083,6 +5083,107 @@ def NextLunarApsis(apsis): raise InternalError() return next + +def VectorFromSphere(sphere, time): + """Converts spherical coordinates to Cartesian coordinates. + + Given spherical coordinates and a time at which they are valid, + returns a vector of Cartesian coordinates. The returned value + includes the time, as required by all `Time` objects. + + Parameters + ---------- + sphere : Spherical + Spherical coordinates to be converted. + time : Time + The time that should be included in the returned vector. + + Returns + ------- + Vector + The vector form of the supplied spherical coordinates. + """ + radlat = math.radians(sphere.lat) + radlon = math.radians(sphere.lon) + rcoslat = sphere.dist * math.cos(radlat) + return Vector( + rcoslat * math.cos(radlon), + rcoslat * math.sin(radlon), + sphere.dist * math.sin(radlat), + time + ) + + +def VectorFromEquator(equ, time): + """Given angular equatorial coordinates in `equ`, calculates equatorial vector. + + Parameters + ---------- + equ : Equatorial + Angular equatorial coordinates to be converted to a vector. + time : Time + The date and time of the observation. This is needed because the returned + vector object requires a valid time value when passed to certain other functions. + + Returns + ------- + Vector + A vector in the equatorial system. + """ + return VectorFromSphere(Spherical(equ.dec, 15.0 * equ.ra, equ.dist), time) + + +def EquatorFromVector(vec): + """Given an equatorial vector, calculates equatorial angular coordinates. + + Parameters + ---------- + vec : Vector + A vector in an equatorial coordinate system. + + Returns + ------- + Equatorial + Angular coordinates expressed in the same equatorial system as `vec`. + """ + sphere = SphereFromVector(vec) + return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist) + + +def SphereFromVector(vector): + """Converts Cartesian coordinates to spherical coordinates. + + Given a Cartesian vector, returns latitude, longitude, and distance. + + Parameters + ---------- + vector : Vector + Cartesian vector to be converted to spherical coordinates. + + Returns + ------- + Spherical + Spherical coordinates that are equivalent to the given vector. + """ + xyproj = vector.x*vector.x + vector.y*vector.y + dist = math.sqrt(xyproj + vector.z*vector.z) + if xyproj == 0.0: + if vector.z == 0.0: + raise Exception('Zero-length vector not allowed.') + lon = 0.0 + if vector.z < 0.0: + lat = -90.0 + else: + lat = +90.0 + else: + lon = math.degrees(math.atan2(vector.y, vector.x)) + if lon < 0.0: + lon += 360.0 + lat = math.degrees(math.atan2(vector.z, math.sqrt(xyproj))) + return Spherical(lat, lon, dist) + + + def InverseRotation(rotation): """Calculates the inverse of a rotation matrix. @@ -5220,3 +5321,48 @@ def Rotation_ECL_EQJ(): [ 0, +c, +s], [ 0, -s, +c] ]) + +def Rotation_EQJ_EQD(time): + """Calculates a rotation matrix from equatorial J2000 (EQJ) 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: EQJ = equatorial system, using equator at J2000 epoch. + Target: EQD = equatorial system, using equator of the specified date/time. + + Parameters + ---------- + time : Time + The date and time at which the Earth's equator defines the target orientation. + + Returns + ------- + RotationMatrix + A rotation matrix that converts EQJ to EQD at `time`. + """ + prec = _precession_rot(0.0, time.tt) + nut = _nutation_rot(time, 0) + return CombineRotation(prec, nut) + + +def Rotation_EQD_EQJ(time): + """Calculates a rotation matrix from equatorial of-date (EQD) 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: EQD = equatorial system, using equator of the specified date/time. + Target: EQJ = equatorial system, using equator at J2000 epoch. + + Parameters + ---------- + time : Time + The date and time at which the Earth's equator defines the source orientation. + + Returns + ------- + RotationMatrix + A rotation matrix that converts EQD at `time` to EQJ. + """ + nut = _nutation_rot(time, 1) + prec = _precession_rot(time.tt, 0.0) + return CombineRotation(nut, prec)