1from numpy import cos
2from skyfield import framelib
3from skyfield.api import Topos, load, wgs84
4from skyfield.constants import AU_M, ERAD
5from skyfield.positionlib import Geocentric
6
7def test_radec_and_altaz_angles_and_rates():
8    # HORIZONS test data in Skyfield repository: authorities/radec-altaz-rates
9    ts = load.timescale()
10    t = ts.utc(2021, 2, 3)
11    top = wgs84.latlon(35.1844866, 248.347300, elevation_m=2106.9128)
12    planets = load('de421.bsp')
13    a = (planets['earth'] + top).at(t).observe(planets['mars']).apparent()
14
15    # First, verify RA and declination.
16
17    frame = framelib.true_equator_and_equinox_of_date
18    dec, ra, distance, dec_rate, ra_rate, range_rate = (
19        a.frame_latlon_and_rates(frame))
20
21    arcseconds = 3600.0
22    assert abs((ra.degrees - 40.75836) * arcseconds) < 0.04
23    assert abs((dec.degrees - 17.16791) * arcseconds) < 0.005
24    assert abs(distance.m - 1.21164331503552 * AU_M) < 120.0
25
26    # Verify RA and declination rates of change.
27
28    assert round(dec_rate.arcseconds.per_hour, 5) == 25.61352
29    assert round(ra_rate.arcseconds.per_hour * cos(dec.radians),
30                 4) == round(75.15571, 4)  # TODO: get last digit to agree?
31    assert abs(range_rate.km_per_s - 16.7926932) < 2e-5
32
33    # Verify altitude and azimuth.
34
35    frame = top
36    alt, az, distance, alt_rate, az_rate, range_rate = (
37        a.frame_latlon_and_rates(frame))
38
39    assert round(alt.degrees, 4) == 65.2758
40    assert round(az.degrees, 4) == 131.8839
41    assert abs(distance.m - 1.21164331503552 * AU_M) < 120.0
42
43    # Verify altitude and azimuth rates of change.
44
45    assert abs(range_rate.km_per_s - 16.7926932) < 2e-5
46    assert round(alt_rate.arcseconds.per_minute, 2) == 548.66
47    assert round(az_rate.arcseconds.per_minute * cos(alt.radians), 2) == 663.55
48
49def test_frame_round_trip():
50    # Does a frame's rotation and twist get applied in the right
51    # directions?  Let's test whether the position and velocity of an
52    # ITRS vector (ERAD,0,0) are restored to the proper orientation.
53    top = Topos(latitude_degrees=0, longitude_degrees=0)
54    ts = load.timescale()
55    t = ts.utc(2020, 11, 27, 15, 34)  # Arbitrary time; LST ~= 20.03.
56    p = top.at(t)
57
58    r = p.frame_xyz(framelib.itrs)
59    assert max(abs(r.m - [ERAD, 0, 0])) < 4e-8 # meters
60
61    r, v = p.frame_xyz_and_velocity(framelib.itrs)
62    assert max(abs(r.m - [ERAD, 0, 0])) < 4e-8 # meters
63    assert max(abs(v.km_per_s)) < 3e-15 # km/s
64
65def test_from_frame_method():
66    ts = load.timescale()
67    t = ts.utc(2020, 11, 27, 15, 34)
68    g1 = Geocentric([1,2,3], [4,5,6], t=t)
69    r, v = g1.frame_xyz_and_velocity(framelib.itrs) # which we trust: see above
70
71    g2 = Geocentric.from_time_and_frame_vectors(t, framelib.itrs, r, v)
72    assert max(abs(g2.position.au - [1,2,3])) < 2e-14
73    assert max(abs(g2.velocity.au_per_d - [4,5,6])) < 3e-14
74
75    # Make sure original vectors were not harmed (for example, by "+=").
76    assert list(g1.position.au) == [1,2,3]
77    assert list(g1.velocity.au_per_d) == [4,5,6]
78
79def test_frame_without_spin():
80    ts = load.timescale()
81    t = ts.utc(2020, 11, 27, 15, 34)
82    g = Geocentric([1,2,3], [4,5,6], t=t)
83
84    # Simply test whether "None" spin raises an exception in either direction.
85    f = framelib.true_equator_and_equinox_of_date
86    r, v = g.frame_xyz_and_velocity(f)
87    Geocentric.from_time_and_frame_vectors(t, f, r, v)
88
89def test_tirs_at_least_runs():
90    # TODO: find an external source for a TIRS vector to test against.
91    # For now, just make sure it doesn't raise an exception.
92    ts = load.timescale()
93    t = ts.utc(2020, 11, 27, 15, 34)
94    g = Geocentric([1,2,3], [4,5,6], t=t)
95    g.frame_xyz_and_velocity(framelib.tirs)
96