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