-
-
Notifications
You must be signed in to change notification settings - Fork 210
/
earthlib.py
161 lines (115 loc) · 5.27 KB
/
earthlib.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
"""Formulae for specific earth behaviors and effects."""
from numpy import (abs, arcsin, arccos, arctan2, array, clip, cos,
minimum, pi, sin, sqrt, tan, where, zeros_like)
from .constants import (AU_M, ANGVEL, DAY_S, DEG2RAD, ERAD,
IERS_2010_INVERSE_EARTH_FLATTENING, RAD2DEG, T0, tau)
from .functions import dots
earth_radius_au = ERAD / AU_M
one_minus_flattening = 1.0 - 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
one_minus_flattening_squared = one_minus_flattening * one_minus_flattening
def terra(latitude, longitude, elevation, gast):
"""Deprecated conversion from lat,lon,t -> GCRS; neglects polar motion."""
zero = zeros_like(gast)
sinphi = sin(latitude)
cosphi = cos(latitude)
c = 1.0 / sqrt(cosphi * cosphi +
sinphi * sinphi * one_minus_flattening_squared)
s = one_minus_flattening_squared * c
ach = earth_radius_au * c + elevation
ash = earth_radius_au * s + elevation
# Compute local sidereal time factors at the observer's longitude.
stlocl = 15.0 * DEG2RAD * gast + longitude
sinst = sin(stlocl)
cosst = cos(stlocl)
# Compute position vector components in au.
ac = ach * cosphi
acsst = ac * sinst
accst = ac * cosst
pos = array((accst, acsst, zero + ash * sinphi))
# Compute velocity vector components in au/day.
vel = ANGVEL * DAY_S * array((-acsst, accst, zero))
return pos, vel
def reverse_terra(xyz_au, gast, iterations=3):
"""Deprecated conversion from GCRS -> lat,lon,t; neglects polar motion."""
x, y, z = xyz_au
R = sqrt(x*x + y*y)
lon = (arctan2(y, x) - 15 * DEG2RAD * gast - pi) % tau - pi
lat = arctan2(z, R)
a = ERAD / AU_M
f = 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
e2 = 2.0*f - f*f
i = 0
C = 1.0
while i < iterations:
i += 1
C = 1.0 / sqrt(1.0 - e2 * (sin(lat) ** 2.0))
lat = arctan2(z + a * C * e2 * sin(lat), R)
elevation_m = ((R / cos(lat)) - a * C) * AU_M
return lat, lon, elevation_m
def compute_limb_angle(position_au, observer_au):
"""Determine the angle of an object above or below the Earth's limb.
Given an object's GCRS `position_au` |xyz| vector and the position
of an `observer_au` as a vector in the same coordinate system,
return a tuple that provides `(limb_ang, nadir_ang)`:
limb_angle
Angle of observed object above (+) or below (-) limb in degrees.
nadir_angle
Nadir angle of observed object as a fraction of apparent radius
of limb: <1.0 means below the limb, =1.0 means on the limb, and
>1.0 means above the limb.
"""
# Compute the distance to the object and the distance to the observer.
disobj = sqrt(dots(position_au, position_au))
disobs = sqrt(dots(observer_au, observer_au))
# Compute apparent angular radius of Earth's limb.
aprad = arcsin(minimum(earth_radius_au / disobs, 1.0))
# Compute zenith distance of Earth's limb.
zdlim = pi - aprad
# Compute zenith distance of observed object.
coszd = dots(position_au, observer_au) / (disobj * disobs)
coszd = clip(coszd, -1.0, 1.0)
zdobj = arccos(coszd)
# Angle of object wrt limb is difference in zenith distances.
limb_angle = (zdlim - zdobj) * RAD2DEG
# Nadir angle of object as a fraction of angular radius of limb.
nadir_angle = (pi - zdobj) / aprad
return limb_angle, nadir_angle
def sidereal_time(t):
"""Compute Greenwich Mean Sidereal Time (GMST) in hours at time ``t``."""
theta = earth_rotation_angle(t.whole, t.ut1_fraction)
# The equinox method. See Circular 179, Section 2.6.2.
# Precession-in-RA terms in mean sidereal time taken from third
# reference, eq. (42), with coefficients in arcseconds.
t = (t.whole - T0 + t.tdb_fraction) / 36525.0
st = ( 0.014506 +
(((( - 0.0000000368 * t
- 0.000029956 ) * t
- 0.00000044 ) * t
+ 1.3915817 ) * t
+ 4612.156534 ) * t)
return (st / 54000.0 + theta * 24.0) % 24.0
def earth_rotation_angle(jd_ut1, fraction_ut1=0.0):
"""Return the value of the Earth Rotation Angle (theta) for a UT1 date.
Uses the expression from the note to IAU Resolution B1.8 of 2000.
Returns a fraction between 0.0 and 1.0 whole rotations.
"""
th = 0.7790572732640 + 0.00273781191135448 * (jd_ut1 - T0 + fraction_ut1)
return (th % 1.0 + jd_ut1 % 1.0 + fraction_ut1) % 1.0
def refraction(alt_degrees, temperature_C, pressure_mbar):
"""Given an observed altitude, return how much the image is refracted.
Zero refraction is returned both for objects very near the zenith,
as well as for objects more than one degree below the horizon.
"""
r = 0.016667 / tan((alt_degrees + 7.31 / (alt_degrees + 4.4)) * DEG2RAD)
d = r * (0.28 * pressure_mbar / (temperature_C + 273.0))
return where((-1.0 <= alt_degrees) & (alt_degrees <= 89.9), d, 0.0)
def refract(alt_degrees, temperature_C, pressure_mbar):
"""Given an unrefracted `alt` determine where it will appear in the sky."""
alt = alt_degrees
while True:
alt1 = alt
alt = alt_degrees + refraction(alt, temperature_C, pressure_mbar)
converged = abs(alt - alt1).max() <= 3.0e-5
if converged:
break
return alt