From c91fe513c1558e62dec4b0d17e631344a0473873 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Fri, 19 Nov 2021 21:40:22 -0500 Subject: [PATCH] PY ObserverState Implemented the Python version of the ObserverState function. --- demo/python/astronomy.py | 85 ++++++++++++++++++++++++++++------ generate/template/astronomy.py | 85 ++++++++++++++++++++++++++++------ generate/test.py | 26 +++++++++++ source/python/README.md | 27 +++++++++++ source/python/astronomy.py | 85 ++++++++++++++++++++++++++++------ 5 files changed, 269 insertions(+), 39 deletions(-) diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index e36dad86..dfb9e76e 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -1439,13 +1439,20 @@ def _precession_rot(time, direction): raise Error('Inalid precession direction') +def _rotate(rot, vec): + return [ + rot.rot[0][0]*vec[0] + rot.rot[1][0]*vec[1] + rot.rot[2][0]*vec[2], + rot.rot[0][1]*vec[0] + rot.rot[1][1]*vec[1] + rot.rot[2][1]*vec[2], + rot.rot[0][2]*vec[0] + rot.rot[1][2]*vec[1] + rot.rot[2][2]*vec[2] + ] + def _precession(pos, time, direction): r = _precession_rot(time, direction) - return [ - r.rot[0][0]*pos[0] + r.rot[1][0]*pos[1] + r.rot[2][0]*pos[2], - r.rot[0][1]*pos[0] + r.rot[1][1]*pos[1] + r.rot[2][1]*pos[2], - r.rot[0][2]*pos[0] + r.rot[1][2]*pos[1] + r.rot[2][2]*pos[2] - ] + return _rotate(r, pos) + +def _precession_posvel(state, time, direction): + r = _precession_rot(time, direction) + return RotateState(r, state) class Equatorial: """Equatorial angular coordinates @@ -1537,14 +1544,13 @@ def _nutation_rot(time, direction): raise Error('Invalid nutation direction') - def _nutation(pos, time, direction): r = _nutation_rot(time, direction) - return [ - r.rot[0][0]*pos[0] + r.rot[1][0]*pos[1] + r.rot[2][0]*pos[2], - r.rot[0][1]*pos[0] + r.rot[1][1]*pos[1] + r.rot[2][1]*pos[2], - r.rot[0][2]*pos[0] + r.rot[1][2]*pos[1] + r.rot[2][2]*pos[2] - ] + return _rotate(r, pos) + +def _nutation_posvel(state, time, direction): + r = _nutation_rot(time, direction) + return RotateState(r, state) def _era(time): # Earth Rotation Angle thet1 = 0.7790572732640 + 0.00273781191135448 * time.ut @@ -1628,7 +1634,7 @@ def _inverse_terra(ovec, st): height_km = p/cos - adjust return Observer(lat_deg, lon_deg, 1000*height_km) -def _terra(observer, st): +def _terra_posvel(observer, st): df2 = _EARTH_FLATTENING ** 2 phi = math.radians(observer.latitude) sinphi = math.sin(phi) @@ -1644,9 +1650,15 @@ def _terra(observer, st): return [ ach * cosphi * cosst / KM_PER_AU, ach * cosphi * sinst / KM_PER_AU, - ash * sinphi / KM_PER_AU + ash * sinphi / KM_PER_AU, + -_ANGVEL * ach * cosphi * sinst * 86400 / KM_PER_AU, + +_ANGVEL * ach * cosphi * cosst * 86400 / KM_PER_AU, + 0.0 ] +def _terra(observer, st): + return _terra_posvel(observer, st)[0:3] + def _geo_pos(time, observer): gast = _sidereal_time(time) pos1 = _terra(observer, gast) @@ -4549,6 +4561,53 @@ def ObserverVector(time, observer, ofdate): ovec = _precession(ovec, time, _PrecessDir.Into2000) return Vector(ovec[0], ovec[1], ovec[2], time) +def ObserverState(time, observer, ofdate): + """Calculates geocentric equatorial position and velocity of an observer on the surface of the Earth. + + This function calculates position and velocity vectors of an observer + on or near the surface of the Earth, expressed in equatorial + coordinates. It takes into account the rotation of the Earth at the given + time, along with the given latitude, longitude, and elevation of the observer. + + The caller may pass `ofdate` as `True` to return coordinates relative to the Earth's + equator at the specified time, or `False` to use the J2000 equator. + + The returned position vector has components expressed in astronomical units (AU). + To convert to kilometers, multiply the `x`, `y`, and `z` values by + the constant value #KM_PER_AU. + The returned velocity vector has components expressed in AU/day. + + Parameters + ---------- + time : Time + The date and time for which to calculate the observer's position and velocity vectors. + observer : Observer + The geographic location of a point on or near the surface of the Earth. + ofdate : bool + Selects the date of the Earth's equator in which to express the equatorial coordinates. + The caller may pass `False` to use the orientation of the Earth's equator + at noon UTC on January 1, 2000, in which case this function corrects for precession + and nutation of the Earth as it was at the moment specified by the `time` parameter. + Or the caller may pass `true` to use the Earth's equator at `time` + as the orientation. + + Returns + ------- + StateVector + An equatorial position vector and velocity vector relative to the center of the Earth. + """ + gast = _sidereal_time(time) + ovec = _terra_posvel(observer, gast) + state = StateVector( + ovec[0], ovec[1], ovec[2], + ovec[3], ovec[4], ovec[5], + time + ) + if not ofdate: + state = _nutation_posvel(state, time, _PrecessDir.Into2000) + state = _precession_posvel(state, time, _PrecessDir.Into2000) + return state + def VectorObserver(vector, ofdate): """Calculates the geographic location corresponding to an equatorial vector. diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index ca457084..36dd8d94 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -903,13 +903,20 @@ def _precession_rot(time, direction): raise Error('Inalid precession direction') +def _rotate(rot, vec): + return [ + rot.rot[0][0]*vec[0] + rot.rot[1][0]*vec[1] + rot.rot[2][0]*vec[2], + rot.rot[0][1]*vec[0] + rot.rot[1][1]*vec[1] + rot.rot[2][1]*vec[2], + rot.rot[0][2]*vec[0] + rot.rot[1][2]*vec[1] + rot.rot[2][2]*vec[2] + ] + def _precession(pos, time, direction): r = _precession_rot(time, direction) - return [ - r.rot[0][0]*pos[0] + r.rot[1][0]*pos[1] + r.rot[2][0]*pos[2], - r.rot[0][1]*pos[0] + r.rot[1][1]*pos[1] + r.rot[2][1]*pos[2], - r.rot[0][2]*pos[0] + r.rot[1][2]*pos[1] + r.rot[2][2]*pos[2] - ] + return _rotate(r, pos) + +def _precession_posvel(state, time, direction): + r = _precession_rot(time, direction) + return RotateState(r, state) class Equatorial: """Equatorial angular coordinates @@ -1001,14 +1008,13 @@ def _nutation_rot(time, direction): raise Error('Invalid nutation direction') - def _nutation(pos, time, direction): r = _nutation_rot(time, direction) - return [ - r.rot[0][0]*pos[0] + r.rot[1][0]*pos[1] + r.rot[2][0]*pos[2], - r.rot[0][1]*pos[0] + r.rot[1][1]*pos[1] + r.rot[2][1]*pos[2], - r.rot[0][2]*pos[0] + r.rot[1][2]*pos[1] + r.rot[2][2]*pos[2] - ] + return _rotate(r, pos) + +def _nutation_posvel(state, time, direction): + r = _nutation_rot(time, direction) + return RotateState(r, state) def _era(time): # Earth Rotation Angle thet1 = 0.7790572732640 + 0.00273781191135448 * time.ut @@ -1092,7 +1098,7 @@ def _inverse_terra(ovec, st): height_km = p/cos - adjust return Observer(lat_deg, lon_deg, 1000*height_km) -def _terra(observer, st): +def _terra_posvel(observer, st): df2 = _EARTH_FLATTENING ** 2 phi = math.radians(observer.latitude) sinphi = math.sin(phi) @@ -1108,9 +1114,15 @@ def _terra(observer, st): return [ ach * cosphi * cosst / KM_PER_AU, ach * cosphi * sinst / KM_PER_AU, - ash * sinphi / KM_PER_AU + ash * sinphi / KM_PER_AU, + -_ANGVEL * ach * cosphi * sinst * 86400 / KM_PER_AU, + +_ANGVEL * ach * cosphi * cosst * 86400 / KM_PER_AU, + 0.0 ] +def _terra(observer, st): + return _terra_posvel(observer, st)[0:3] + def _geo_pos(time, observer): gast = _sidereal_time(time) pos1 = _terra(observer, gast) @@ -2507,6 +2519,53 @@ def ObserverVector(time, observer, ofdate): ovec = _precession(ovec, time, _PrecessDir.Into2000) return Vector(ovec[0], ovec[1], ovec[2], time) +def ObserverState(time, observer, ofdate): + """Calculates geocentric equatorial position and velocity of an observer on the surface of the Earth. + + This function calculates position and velocity vectors of an observer + on or near the surface of the Earth, expressed in equatorial + coordinates. It takes into account the rotation of the Earth at the given + time, along with the given latitude, longitude, and elevation of the observer. + + The caller may pass `ofdate` as `True` to return coordinates relative to the Earth's + equator at the specified time, or `False` to use the J2000 equator. + + The returned position vector has components expressed in astronomical units (AU). + To convert to kilometers, multiply the `x`, `y`, and `z` values by + the constant value #KM_PER_AU. + The returned velocity vector has components expressed in AU/day. + + Parameters + ---------- + time : Time + The date and time for which to calculate the observer's position and velocity vectors. + observer : Observer + The geographic location of a point on or near the surface of the Earth. + ofdate : bool + Selects the date of the Earth's equator in which to express the equatorial coordinates. + The caller may pass `False` to use the orientation of the Earth's equator + at noon UTC on January 1, 2000, in which case this function corrects for precession + and nutation of the Earth as it was at the moment specified by the `time` parameter. + Or the caller may pass `true` to use the Earth's equator at `time` + as the orientation. + + Returns + ------- + StateVector + An equatorial position vector and velocity vector relative to the center of the Earth. + """ + gast = _sidereal_time(time) + ovec = _terra_posvel(observer, gast) + state = StateVector( + ovec[0], ovec[1], ovec[2], + ovec[3], ovec[4], ovec[5], + time + ) + if not ofdate: + state = _nutation_posvel(state, time, _PrecessDir.Into2000) + state = _precession_posvel(state, time, _PrecessDir.Into2000) + return state + def VectorObserver(vector, ofdate): """Calculates the geographic location corresponding to an equatorial vector. diff --git a/generate/test.py b/generate/test.py index 44285431..b82699dd 100755 --- a/generate/test.py +++ b/generate/test.py @@ -2016,6 +2016,31 @@ def HelioState(): #----------------------------------------------------------------------------------------------------------- +def TopoStateFunc(body, time): + observer = astronomy.Observer(30.0, -80.0, 1000.0) + observer_state = astronomy.ObserverState(time, observer, False) + if body == _Body_Geo_EMB: + state = astronomy.GeoEmbState(time) + elif body == astronomy.Body.Earth: + state = astronomy.StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time) + else: + raise Exception('PY TopoStateFunction: unsupported body ' + body) + state.x -= observer_state.x + state.y -= observer_state.y + state.z -= observer_state.z + state.vx -= observer_state.vx + state.vy -= observer_state.vy + state.vz -= observer_state.vz + return state + +def TopoState(): + if VerifyStateBody(TopoStateFunc, astronomy.Body.Earth, 'topostate/Earth_N30_W80_1000m.txt', 2.108e-04, 2.430e-04): return 1 + if VerifyStateBody(TopoStateFunc, _Body_Geo_EMB, 'topostate/EMB_N30_W80_1000m.txt', 7.195e-04, 2.497e-04): return 1 + print('PY TopoState: PASS') + return 0 + +#----------------------------------------------------------------------------------------------------------- + def Aberration(): THRESHOLD_SECONDS = 0.4 filename = 'equatorial/Mars_j2000_ofdate_aberration.txt' @@ -2236,6 +2261,7 @@ UnitTests = { 'rotation': Rotation, 'seasons': Seasons, 'time': AstroTime, + 'topostate': TopoState, 'transit': Transit, 'twilight': Twilight, } diff --git a/source/python/README.md b/source/python/README.md index 72f1b22b..0ce19994 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -1771,6 +1771,33 @@ The effective gravitational acceleration expressed in meters per second squared --- + +### ObserverState(time, observer, ofdate) + +**Calculates geocentric equatorial position and velocity of an observer on the surface of the Earth.** + +This function calculates position and velocity vectors of an observer +on or near the surface of the Earth, expressed in equatorial +coordinates. It takes into account the rotation of the Earth at the given +time, along with the given latitude, longitude, and elevation of the observer. +The caller may pass `ofdate` as `True` to return coordinates relative to the Earth's +equator at the specified time, or `False` to use the J2000 equator. +The returned position vector has components expressed in astronomical units (AU). +To convert to kilometers, multiply the `x`, `y`, and `z` values by +the constant value [`KM_PER_AU`](#KM_PER_AU). +The returned velocity vector has components expressed in AU/day. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Time`](#Time) | `time` | The date and time for which to calculate the observer's position and velocity vectors. | +| [`Observer`](#Observer) | `observer` | The geographic location of a point on or near the surface of the Earth. | +| `bool` | `ofdate` | Selects the date of the Earth's equator in which to express the equatorial coordinates. The caller may pass `False` to use the orientation of the Earth's equator at noon UTC on January 1, 2000, in which case this function corrects for precession and nutation of the Earth as it was at the moment specified by the `time` parameter. Or the caller may pass `true` to use the Earth's equator at `time` as the orientation. | + +### Returns: [`StateVector`](#StateVector) +An equatorial position vector and velocity vector relative to the center of the Earth. + +--- + ### ObserverVector(time, observer, ofdate) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index e36dad86..dfb9e76e 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -1439,13 +1439,20 @@ def _precession_rot(time, direction): raise Error('Inalid precession direction') +def _rotate(rot, vec): + return [ + rot.rot[0][0]*vec[0] + rot.rot[1][0]*vec[1] + rot.rot[2][0]*vec[2], + rot.rot[0][1]*vec[0] + rot.rot[1][1]*vec[1] + rot.rot[2][1]*vec[2], + rot.rot[0][2]*vec[0] + rot.rot[1][2]*vec[1] + rot.rot[2][2]*vec[2] + ] + def _precession(pos, time, direction): r = _precession_rot(time, direction) - return [ - r.rot[0][0]*pos[0] + r.rot[1][0]*pos[1] + r.rot[2][0]*pos[2], - r.rot[0][1]*pos[0] + r.rot[1][1]*pos[1] + r.rot[2][1]*pos[2], - r.rot[0][2]*pos[0] + r.rot[1][2]*pos[1] + r.rot[2][2]*pos[2] - ] + return _rotate(r, pos) + +def _precession_posvel(state, time, direction): + r = _precession_rot(time, direction) + return RotateState(r, state) class Equatorial: """Equatorial angular coordinates @@ -1537,14 +1544,13 @@ def _nutation_rot(time, direction): raise Error('Invalid nutation direction') - def _nutation(pos, time, direction): r = _nutation_rot(time, direction) - return [ - r.rot[0][0]*pos[0] + r.rot[1][0]*pos[1] + r.rot[2][0]*pos[2], - r.rot[0][1]*pos[0] + r.rot[1][1]*pos[1] + r.rot[2][1]*pos[2], - r.rot[0][2]*pos[0] + r.rot[1][2]*pos[1] + r.rot[2][2]*pos[2] - ] + return _rotate(r, pos) + +def _nutation_posvel(state, time, direction): + r = _nutation_rot(time, direction) + return RotateState(r, state) def _era(time): # Earth Rotation Angle thet1 = 0.7790572732640 + 0.00273781191135448 * time.ut @@ -1628,7 +1634,7 @@ def _inverse_terra(ovec, st): height_km = p/cos - adjust return Observer(lat_deg, lon_deg, 1000*height_km) -def _terra(observer, st): +def _terra_posvel(observer, st): df2 = _EARTH_FLATTENING ** 2 phi = math.radians(observer.latitude) sinphi = math.sin(phi) @@ -1644,9 +1650,15 @@ def _terra(observer, st): return [ ach * cosphi * cosst / KM_PER_AU, ach * cosphi * sinst / KM_PER_AU, - ash * sinphi / KM_PER_AU + ash * sinphi / KM_PER_AU, + -_ANGVEL * ach * cosphi * sinst * 86400 / KM_PER_AU, + +_ANGVEL * ach * cosphi * cosst * 86400 / KM_PER_AU, + 0.0 ] +def _terra(observer, st): + return _terra_posvel(observer, st)[0:3] + def _geo_pos(time, observer): gast = _sidereal_time(time) pos1 = _terra(observer, gast) @@ -4549,6 +4561,53 @@ def ObserverVector(time, observer, ofdate): ovec = _precession(ovec, time, _PrecessDir.Into2000) return Vector(ovec[0], ovec[1], ovec[2], time) +def ObserverState(time, observer, ofdate): + """Calculates geocentric equatorial position and velocity of an observer on the surface of the Earth. + + This function calculates position and velocity vectors of an observer + on or near the surface of the Earth, expressed in equatorial + coordinates. It takes into account the rotation of the Earth at the given + time, along with the given latitude, longitude, and elevation of the observer. + + The caller may pass `ofdate` as `True` to return coordinates relative to the Earth's + equator at the specified time, or `False` to use the J2000 equator. + + The returned position vector has components expressed in astronomical units (AU). + To convert to kilometers, multiply the `x`, `y`, and `z` values by + the constant value #KM_PER_AU. + The returned velocity vector has components expressed in AU/day. + + Parameters + ---------- + time : Time + The date and time for which to calculate the observer's position and velocity vectors. + observer : Observer + The geographic location of a point on or near the surface of the Earth. + ofdate : bool + Selects the date of the Earth's equator in which to express the equatorial coordinates. + The caller may pass `False` to use the orientation of the Earth's equator + at noon UTC on January 1, 2000, in which case this function corrects for precession + and nutation of the Earth as it was at the moment specified by the `time` parameter. + Or the caller may pass `true` to use the Earth's equator at `time` + as the orientation. + + Returns + ------- + StateVector + An equatorial position vector and velocity vector relative to the center of the Earth. + """ + gast = _sidereal_time(time) + ovec = _terra_posvel(observer, gast) + state = StateVector( + ovec[0], ovec[1], ovec[2], + ovec[3], ovec[4], ovec[5], + time + ) + if not ofdate: + state = _nutation_posvel(state, time, _PrecessDir.Into2000) + state = _precession_posvel(state, time, _PrecessDir.Into2000) + return state + def VectorObserver(vector, ofdate): """Calculates the geographic location corresponding to an equatorial vector.