diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 9f652614..dfb48049 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -78,6 +78,7 @@ _SUN_RADIUS_AU = _SUN_RADIUS_KM / KM_PER_AU _EARTH_FLATTENING = 0.996647180302104 _EARTH_EQUATORIAL_RADIUS_KM = 6378.1366 +_EARTH_POLAR_RADIUS_KM = _EARTH_EQUATORIAL_RADIUS_KM * _EARTH_FLATTENING _EARTH_EQUATORIAL_RADIUS_AU = _EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU _EARTH_MEAN_RADIUS_KM = 6371.0 # mean radius of the Earth's geoid, without atmosphere _EARTH_ATMOSPHERE_KM = 88.0 # effective atmosphere thickness for lunar eclipses @@ -1567,10 +1568,63 @@ def _sidereal_time(time): # return sidereal hours in the half-open range [0, 24). return gst +def _inverse_terra(ovec, st): + # Convert from AU to kilometers + x = ovec[0] * KM_PER_AU + y = ovec[1] * KM_PER_AU + z = ovec[2] * KM_PER_AU + p = math.sqrt(x*x + y*y) + if p < 1.0e-6: + # Special case: within 1 millimeter of a pole! + # Use arbitrary longitude, and latitude determined by polarity of z. + lon_deg = 0.0 + if z > 0.0: + lat_deg = +90.0 + else: + lat_deg = -90.0 + # Elevation is calculated directly from z + height_km = abs(z) - _EARTH_POLAR_RADIUS_KM + else: + stlocl = math.atan2(y, x) + # Calculate exact longitude. + lon_deg = math.degrees(stlocl) - (15.0 * st) + # Normalize longitude to the range (-180, +180]. + while lon_deg <= -180.0: + lon_deg += 360.0 + while lon_deg > +180.0: + lon_deg -= 360.0 + # Numerically solve for exact latitude. + F = _EARTH_FLATTENING ** 2 + def W(context, time): + lat = time.ut + cos = math.cos(lat) + sin = math.sin(lat) + return ((F-1)*_EARTH_EQUATORIAL_RADIUS_KM*sin*cos)/math.sqrt(cos*cos + F*sin*sin) - z*cos + p*sin + # Total hack: use Search(), even though it expects a function of time. + # Pretend that the angle is a time! + # This is not as efficient as I would like, because of all the unnecessary DeltaT calculations. + # But the search is very well tested and converges quickly. + t1 = Time(-math.pi/2) + t2 = Time(+math.pi/2) + dt = 86400.0 * math.radians(1.0e-6) # hack to converge at a millionth of a degree error + tx = Search(W, None, t1, t2, dt) + if tx is None: + raise Error('Cannot solve for geographic latitude') + lat_rad = tx.ut + lat_deg = math.degrees(lat_rad) + # Solve for exact height in meters. + # There are two formulas I can use. Use whichever has the less risky denominator. + sin = math.sin(lat_rad) + cos = math.cos(lat_rad) + adjust = _EARTH_EQUATORIAL_RADIUS_KM/math.sqrt(cos*cos + F*sin*sin) + if abs(sin) > abs(cos): + height_km = z/sin - F*adjust + else: + height_km = p/cos - adjust + return Observer(lat_deg, lon_deg, 1000*height_km) def _terra(observer, st): - df = 1.0 - 0.003352819697896 # flattening of the Earth - df2 = df * df + df2 = _EARTH_FLATTENING ** 2 phi = math.radians(observer.latitude) sinphi = math.sin(phi) cosphi = math.cos(phi) @@ -4238,6 +4292,38 @@ def ObserverVector(time, observer, ofdate): ovec = _precession(ovec, time, _PrecessDir.Into2000) return Vector(ovec[0], ovec[1], ovec[2], time) +def VectorObserver(vector, ofdate): + """Calculates the geographic location corresponding to an equatorial vector. + + This is the inverse function of #ObserverVector. + Instead of converting an #Observer to a #Vector, it converts + a `Vector` to an `Observer`. + + Parameters + ---------- + vector : Vector + The geocentric vector in an equatorial orientation. + The components are expressed in astronomical units (AU). + The time `vector.t` determines the Earth's rotation. + 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 the time `vector.t`. + Or the caller may pass `True` to use the Earth's equator at `vector.t` as the orientation. + + Returns + ------- + Observer + The geographic latitude, longitude, and elevation above sea level + that corresponds to the given equatorial vector. + """ + gast = _sidereal_time(vector.t) + ovec = [vector.x, vector.y, vector.z] + if not ofdate: + ovec = _precession(ovec, vector.t, _PrecessDir.From2000) + ovec = _nutation(ovec, vector.t, _PrecessDir.From2000) + return _inverse_terra(ovec, gast) @enum.unique class Refraction(enum.Enum): diff --git a/demo/python/test/.gitignore b/demo/python/test/.gitignore index 5a1c8f90..200dc141 100644 --- a/demo/python/test/.gitignore +++ b/demo/python/test/.gitignore @@ -10,3 +10,4 @@ lunar_eclipse.txt lunar_angles.txt jupiter_moons.txt galactic.txt +triangulate.txt diff --git a/demo/python/test/triangulate_correct.txt b/demo/python/test/triangulate_correct.txt new file mode 100644 index 00000000..ce34c426 --- /dev/null +++ b/demo/python/test/triangulate_correct.txt @@ -0,0 +1,2 @@ +Solution #1 = (N48.28914299803458, E24.562600140892187, 3862.0102185777796m) +Solution #2 = (N48.28914299803458, E24.562600140892187, 3862.0102185777796m) diff --git a/demo/python/triangulate.py b/demo/python/triangulate.py index 92e472b7..224b6b3e 100755 --- a/demo/python/triangulate.py +++ b/demo/python/triangulate.py @@ -20,6 +20,12 @@ come closest to intersecting. It then prints out the coordinates of that triangulation point, along with the error radius in meters. ''' +Verbose = False + +def Debug(text): + if Verbose: + print(text) + def DotProduct(a, b): return a.x*b.x + a.y*b.y + a.z*b.z @@ -32,20 +38,47 @@ def Intersect(pos1, dir1, pos2, dir2): # 0 = G + Fu - v denom = 1.0 - F*F if denom == 0.0: - print('Cannot solve because directions are parallel.') - else: - u = (F*G - E) / denom - v = G + F*u - print('u={}, v={}'.format(u, v)) + Debug('Cannot solve because directions are parallel.') + return None + + u = (F*G - E) / denom + v = G + F*u + Debug('u={}, v={}'.format(u, v)) + if u < 0 or v < 0: + Debug('Lines of sight do not converge.') + return None + + a = Vector(pos1.x + u*dir1.x, pos1.y + u*dir1.y, pos1.z + u*dir1.z, pos1.t) + b = Vector(pos2.x + v*dir2.x, pos2.y + v*dir2.y, pos2.z + v*dir2.z, pos2.t) + c = Vector((a.x + b.x)/2, (a.y + b.y)/2, (a.z + b.z)/2, a.t) + Debug('c = {}'.format(c)) + # Convert vector back to geographic coordinates + return VectorObserver(c, True) + +def DirectionVector(time, observer, altitude, azimuth): + # Convert horizontal angles to a horizontal unit vector. + hor = Spherical(altitude, azimuth, 1.0) + hvec = VectorFromHorizon(hor, time, Refraction.Airless) + # Find the rotation matrix that converts horizontal vectors to equatorial vectors. + rot = Rotation_HOR_EQD(time, observer) + # Rotate the horizontal (HOR) vector to an equator-of-date (EQD) vector. + evec = RotateVector(rot, hvec) + return evec if __name__ == '__main__': # Validate and parse command line arguments. - if len(sys.argv) != 11: + # The '-v' option enables debug prints. + args = sys.argv[1:] + if len(args) > 0 and args[0] == '-v': + Verbose = True + args = args[1:] + + if len(args) != 10: print(UsageText) sys.exit(1) - lat1, lon1, elv1, az1, alt1, lat2, lon2, elv2, az2, alt2 = [float(x) for x in sys.argv[1:]] + lat1, lon1, elv1, az1, alt1, lat2, lon2, elv2, az2, alt2 = [float(x) for x in args] obs1 = Observer(lat1, lon1, elv1) obs2 = Observer(lat2, lon2, elv2) @@ -56,26 +89,34 @@ if __name__ == '__main__': time = Time(0.0) pos1 = ObserverVector(time, obs1, True) - print('Observer #1 = {}'.format(obs1)) - print('Position #1 = ({:0.6e}, {:0.6e}, {:0.6e})'.format(pos1.x, pos1.y, pos1.z)) + Debug('Observer #1 = {}'.format(obs1)) + Debug('Position #1 = ({:0.6f}, {:0.6f}, {:0.6f})'.format(pos1.x * KM_PER_AU, pos1.y * KM_PER_AU, pos1.z * KM_PER_AU)) pos2 = ObserverVector(time, obs2, True) - print('Observer #2 = {}'.format(obs2)) - print('Position #2 = ({:0.6e}, {:0.6e}, {:0.6e})'.format(pos2.x, pos2.y, pos2.z)) + Debug('Observer #2 = {}'.format(obs2)) + Debug('Position #2 = ({:0.6f}, {:0.6f}, {:0.6f})'.format(pos2.x * KM_PER_AU, pos2.y * KM_PER_AU, pos2.z * KM_PER_AU)) - # Convert horizontal coordinates into unit direction vectors. + # Convert horizontal coordinates into unit direction vectors in EQD orientation. + dir1 = DirectionVector(time, obs1, alt1, az1) + dir2 = DirectionVector(time, obs2, alt2, az2) - hor1 = Spherical(alt1, az1, 1.0) - hor2 = Spherical(alt2, az2, 1.0) - dir1 = VectorFromHorizon(hor1, time, Refraction.Airless) - dir2 = VectorFromHorizon(hor2, time, Refraction.Airless) - print('Direction #1 = ({:0.6e}, {:0.6e}, {:0.6e})'.format(dir1.x, dir1.y, dir1.z)) - print('Direction #2 = ({:0.6e}, {:0.6e}, {:0.6e})'.format(dir2.x, dir2.y, dir2.z)) + Debug('Direction #1 = ({:0.6f}, {:0.6f}, {:0.6f})'.format(dir1.x, dir1.y, dir1.z)) + Debug('Direction #2 = ({:0.6f}, {:0.6f}, {:0.6f})'.format(dir2.x, dir2.y, dir2.z)) # Solve for the target point. - Intersect(pos1, dir1, pos2, dir2) + obs = Intersect(pos1, dir1, pos2, dir2) + if obs is None: + print('ERROR: Could not find intersection.') + sys.exit(1) + + print('Solution #1 = {}'.format(obs)) # Solve again with the inputs reversed, as a sanity check. - Intersect(pos2, dir2, pos1, dir1) + check_obs = Intersect(pos2, dir2, pos1, dir1) + if check_obs is None: + print('INTERNAL ERROR: inconsistent solution.') + sys.exit(1) + + print('Solution #2 = {}'.format(check_obs)) sys.exit(0) diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 6261f034..71f1158b 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -78,6 +78,7 @@ _SUN_RADIUS_AU = _SUN_RADIUS_KM / KM_PER_AU _EARTH_FLATTENING = 0.996647180302104 _EARTH_EQUATORIAL_RADIUS_KM = 6378.1366 +_EARTH_POLAR_RADIUS_KM = _EARTH_EQUATORIAL_RADIUS_KM * _EARTH_FLATTENING _EARTH_EQUATORIAL_RADIUS_AU = _EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU _EARTH_MEAN_RADIUS_KM = 6371.0 # mean radius of the Earth's geoid, without atmosphere _EARTH_ATMOSPHERE_KM = 88.0 # effective atmosphere thickness for lunar eclipses @@ -1031,10 +1032,63 @@ def _sidereal_time(time): # return sidereal hours in the half-open range [0, 24). return gst +def _inverse_terra(ovec, st): + # Convert from AU to kilometers + x = ovec[0] * KM_PER_AU + y = ovec[1] * KM_PER_AU + z = ovec[2] * KM_PER_AU + p = math.sqrt(x*x + y*y) + if p < 1.0e-6: + # Special case: within 1 millimeter of a pole! + # Use arbitrary longitude, and latitude determined by polarity of z. + lon_deg = 0.0 + if z > 0.0: + lat_deg = +90.0 + else: + lat_deg = -90.0 + # Elevation is calculated directly from z + height_km = abs(z) - _EARTH_POLAR_RADIUS_KM + else: + stlocl = math.atan2(y, x) + # Calculate exact longitude. + lon_deg = math.degrees(stlocl) - (15.0 * st) + # Normalize longitude to the range (-180, +180]. + while lon_deg <= -180.0: + lon_deg += 360.0 + while lon_deg > +180.0: + lon_deg -= 360.0 + # Numerically solve for exact latitude. + F = _EARTH_FLATTENING ** 2 + def W(context, time): + lat = time.ut + cos = math.cos(lat) + sin = math.sin(lat) + return ((F-1)*_EARTH_EQUATORIAL_RADIUS_KM*sin*cos)/math.sqrt(cos*cos + F*sin*sin) - z*cos + p*sin + # Total hack: use Search(), even though it expects a function of time. + # Pretend that the angle is a time! + # This is not as efficient as I would like, because of all the unnecessary DeltaT calculations. + # But the search is very well tested and converges quickly. + t1 = Time(-math.pi/2) + t2 = Time(+math.pi/2) + dt = 86400.0 * math.radians(1.0e-6) # hack to converge at a millionth of a degree error + tx = Search(W, None, t1, t2, dt) + if tx is None: + raise Error('Cannot solve for geographic latitude') + lat_rad = tx.ut + lat_deg = math.degrees(lat_rad) + # Solve for exact height in meters. + # There are two formulas I can use. Use whichever has the less risky denominator. + sin = math.sin(lat_rad) + cos = math.cos(lat_rad) + adjust = _EARTH_EQUATORIAL_RADIUS_KM/math.sqrt(cos*cos + F*sin*sin) + if abs(sin) > abs(cos): + height_km = z/sin - F*adjust + else: + height_km = p/cos - adjust + return Observer(lat_deg, lon_deg, 1000*height_km) def _terra(observer, st): - df = 1.0 - 0.003352819697896 # flattening of the Earth - df2 = df * df + df2 = _EARTH_FLATTENING ** 2 phi = math.radians(observer.latitude) sinphi = math.sin(phi) cosphi = math.cos(phi) @@ -2208,6 +2262,38 @@ def ObserverVector(time, observer, ofdate): ovec = _precession(ovec, time, _PrecessDir.Into2000) return Vector(ovec[0], ovec[1], ovec[2], time) +def VectorObserver(vector, ofdate): + """Calculates the geographic location corresponding to an equatorial vector. + + This is the inverse function of #ObserverVector. + Instead of converting an #Observer to a #Vector, it converts + a `Vector` to an `Observer`. + + Parameters + ---------- + vector : Vector + The geocentric vector in an equatorial orientation. + The components are expressed in astronomical units (AU). + The time `vector.t` determines the Earth's rotation. + 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 the time `vector.t`. + Or the caller may pass `True` to use the Earth's equator at `vector.t` as the orientation. + + Returns + ------- + Observer + The geographic latitude, longitude, and elevation above sea level + that corresponds to the given equatorial vector. + """ + gast = _sidereal_time(vector.t) + ovec = [vector.x, vector.y, vector.z] + if not ofdate: + ovec = _precession(ovec, vector.t, _PrecessDir.From2000) + ovec = _nutation(ovec, vector.t, _PrecessDir.From2000) + return _inverse_terra(ovec, gast) @enum.unique class Refraction(enum.Enum): diff --git a/generate/test.py b/generate/test.py index 00e5187b..5d39cd1b 100755 --- a/generate/test.py +++ b/generate/test.py @@ -1692,9 +1692,10 @@ def GeoidTestCase(time, observer, ofdate): dy = astronomy.KM_PER_AU * v((geo_moon.y - surface.y) - topo_moon.vec.y) dz = astronomy.KM_PER_AU * v((geo_moon.z - surface.z) - topo_moon.vec.z) diff = sqrt(dx*dx + dy*dy + dz*dz) - Debug('PY GeoidTestCase: ofdate={}, time={}, lat={}, ht={}, surface=({}, {}, {}), diff = {} km'.format( - ofdate, time, - observer.latitude, observer.longitude, observer.height, + Debug('PY GeoidTestCase: ofdate={}, time={}, obs={}, surface=({}, {}, {}), diff = {} km'.format( + ofdate, + time, + observer, astronomy.KM_PER_AU * surface.x, astronomy.KM_PER_AU * surface.y, astronomy.KM_PER_AU * surface.z, @@ -1703,7 +1704,32 @@ def GeoidTestCase(time, observer, ofdate): # Require 1 millimeter accuracy! (one millionth of a kilometer). if diff > 1.0e-6: - print('JS GeoidTestCase: EXCESSIVE POSITION ERROR.') + print('PY GeoidTestCase: EXCESSIVE POSITION ERROR.') + return 1 + + # Verify that we can convert the surface vector back to an observer. + vobs = astronomy.VectorObserver(surface, ofdate) + lat_diff = vabs(vobs.latitude - observer.latitude) + + # Longitude is meaningless at the poles, so don't bother checking it there. + if -89.99 <= observer.latitude <= +89.99: + lon_diff = vabs(vobs.longitude - observer.longitude) + if lon_diff > 180.0: + lon_diff = 360.0 - lon_diff + lon_diff = vabs(lon_diff * math.cos(math.degrees(observer.latitude))) + if lon_diff > 1.0e-6: + print('PY GeoidTestCase: EXCESSIVE longitude check error = {}'.format(lon_diff)) + return 1 + else: + lon_diff = 0.0 + + h_diff = vabs(vobs.height - observer.height) + Debug('PY GeoidTestCase: vobs={}, lat_diff={}, lon_diff={}, h_diff={}'.format(vobs, lat_diff, lon_diff, h_diff)) + if lat_diff > 1.0e-6: + print('PY GeoidTestCase: EXCESSIVE latitude check error = {}'.format(lat_diff)) + return 1 + if h_diff > 0.001: + print('PY GeoidTestCase: EXCESSIVE height check error = {}'.format(h_diff)) return 1 return 0 diff --git a/source/python/README.md b/source/python/README.md index 50e424bd..61c341a2 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -2566,3 +2566,23 @@ includes the time, as required by all `Time` objects. ### Returns: [`Vector`](#Vector) The vector form of the supplied spherical coordinates. +--- + + +### VectorObserver(vector, ofdate) + +**Calculates the geographic location corresponding to an equatorial vector.** + +This is the inverse function of [`ObserverVector`](#ObserverVector). +Instead of converting an [`Observer`](#Observer) to a [`Vector`](#Vector), it converts +a `Vector` to an `Observer`. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Vector`](#Vector) | `vector` | The geocentric vector in an equatorial orientation. The components are expressed in astronomical units (AU). The time `vector.t` determines the Earth's rotation. | +| `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 the time `vector.t`. Or the caller may pass `True` to use the Earth's equator at `vector.t` as the orientation. | + +### Returns: [`Observer`](#Observer) +The geographic latitude, longitude, and elevation above sea level +that corresponds to the given equatorial vector. + diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 9f652614..dfb48049 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -78,6 +78,7 @@ _SUN_RADIUS_AU = _SUN_RADIUS_KM / KM_PER_AU _EARTH_FLATTENING = 0.996647180302104 _EARTH_EQUATORIAL_RADIUS_KM = 6378.1366 +_EARTH_POLAR_RADIUS_KM = _EARTH_EQUATORIAL_RADIUS_KM * _EARTH_FLATTENING _EARTH_EQUATORIAL_RADIUS_AU = _EARTH_EQUATORIAL_RADIUS_KM / KM_PER_AU _EARTH_MEAN_RADIUS_KM = 6371.0 # mean radius of the Earth's geoid, without atmosphere _EARTH_ATMOSPHERE_KM = 88.0 # effective atmosphere thickness for lunar eclipses @@ -1567,10 +1568,63 @@ def _sidereal_time(time): # return sidereal hours in the half-open range [0, 24). return gst +def _inverse_terra(ovec, st): + # Convert from AU to kilometers + x = ovec[0] * KM_PER_AU + y = ovec[1] * KM_PER_AU + z = ovec[2] * KM_PER_AU + p = math.sqrt(x*x + y*y) + if p < 1.0e-6: + # Special case: within 1 millimeter of a pole! + # Use arbitrary longitude, and latitude determined by polarity of z. + lon_deg = 0.0 + if z > 0.0: + lat_deg = +90.0 + else: + lat_deg = -90.0 + # Elevation is calculated directly from z + height_km = abs(z) - _EARTH_POLAR_RADIUS_KM + else: + stlocl = math.atan2(y, x) + # Calculate exact longitude. + lon_deg = math.degrees(stlocl) - (15.0 * st) + # Normalize longitude to the range (-180, +180]. + while lon_deg <= -180.0: + lon_deg += 360.0 + while lon_deg > +180.0: + lon_deg -= 360.0 + # Numerically solve for exact latitude. + F = _EARTH_FLATTENING ** 2 + def W(context, time): + lat = time.ut + cos = math.cos(lat) + sin = math.sin(lat) + return ((F-1)*_EARTH_EQUATORIAL_RADIUS_KM*sin*cos)/math.sqrt(cos*cos + F*sin*sin) - z*cos + p*sin + # Total hack: use Search(), even though it expects a function of time. + # Pretend that the angle is a time! + # This is not as efficient as I would like, because of all the unnecessary DeltaT calculations. + # But the search is very well tested and converges quickly. + t1 = Time(-math.pi/2) + t2 = Time(+math.pi/2) + dt = 86400.0 * math.radians(1.0e-6) # hack to converge at a millionth of a degree error + tx = Search(W, None, t1, t2, dt) + if tx is None: + raise Error('Cannot solve for geographic latitude') + lat_rad = tx.ut + lat_deg = math.degrees(lat_rad) + # Solve for exact height in meters. + # There are two formulas I can use. Use whichever has the less risky denominator. + sin = math.sin(lat_rad) + cos = math.cos(lat_rad) + adjust = _EARTH_EQUATORIAL_RADIUS_KM/math.sqrt(cos*cos + F*sin*sin) + if abs(sin) > abs(cos): + height_km = z/sin - F*adjust + else: + height_km = p/cos - adjust + return Observer(lat_deg, lon_deg, 1000*height_km) def _terra(observer, st): - df = 1.0 - 0.003352819697896 # flattening of the Earth - df2 = df * df + df2 = _EARTH_FLATTENING ** 2 phi = math.radians(observer.latitude) sinphi = math.sin(phi) cosphi = math.cos(phi) @@ -4238,6 +4292,38 @@ def ObserverVector(time, observer, ofdate): ovec = _precession(ovec, time, _PrecessDir.Into2000) return Vector(ovec[0], ovec[1], ovec[2], time) +def VectorObserver(vector, ofdate): + """Calculates the geographic location corresponding to an equatorial vector. + + This is the inverse function of #ObserverVector. + Instead of converting an #Observer to a #Vector, it converts + a `Vector` to an `Observer`. + + Parameters + ---------- + vector : Vector + The geocentric vector in an equatorial orientation. + The components are expressed in astronomical units (AU). + The time `vector.t` determines the Earth's rotation. + 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 the time `vector.t`. + Or the caller may pass `True` to use the Earth's equator at `vector.t` as the orientation. + + Returns + ------- + Observer + The geographic latitude, longitude, and elevation above sea level + that corresponds to the given equatorial vector. + """ + gast = _sidereal_time(vector.t) + ovec = [vector.x, vector.y, vector.z] + if not ofdate: + ovec = _precession(ovec, vector.t, _PrecessDir.From2000) + ovec = _nutation(ovec, vector.t, _PrecessDir.From2000) + return _inverse_terra(ovec, gast) @enum.unique class Refraction(enum.Enum):