From 3e2c39f40a2122e009d328ad5e8d952023713ff8 Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Tue, 24 Sep 2024 13:35:30 -0600 Subject: [PATCH] Refactor tests to simplify sampling LEO and GEO orbits --- tests/ssapy_test_helpers.py | 43 +++++ tests/test_accel.py | 162 ++----------------- tests/test_orbit.py | 304 ++++-------------------------------- 3 files changed, 85 insertions(+), 424 deletions(-) diff --git a/tests/ssapy_test_helpers.py b/tests/ssapy_test_helpers.py index e63f231..d625699 100644 --- a/tests/ssapy_test_helpers.py +++ b/tests/ssapy_test_helpers.py @@ -1,4 +1,5 @@ import numpy as np +import ssapy def timer(f): @@ -34,3 +35,45 @@ def checkSphere(lon1, lat1, lon2, lat2, atol=1e-14, verbose=False): if verbose: print(f"max angle difference {np.rad2deg(np.max(da))*3600} arcsec") np.testing.assert_allclose(da, 0, rtol=0, atol=atol) + + +def sample_orbit(t, r_low, r_high, v_low, v_high): + """ + Sample a random orbit by drawing a position vector with magnitude in + (r_low, r_high), a velocity vector with magnitude in (v_low, v_high), and + a direction, at time t. + """ + def sample_vector(lower, upper): + """ + Helper function for sampling a position and velocity vector. + """ + # Sample magnitude of vector + mag = np.random.uniform(lower, upper) + # Sample a random direction (not uniform on sphere) + theta = np.random.uniform(0, 2*np.pi) + phi = np.random.uniform(0, np.pi) + vec_x = mag * np.cos(theta) * np.sin(phi) + vec_y = mag * np.sin(theta) * np.sin(phi) + vec_z = mag * np.cos(phi) + vec = np.array([vec_x, vec_y, vec_z]) + return vec + + # Sample position and velocity + r = sample_vector(r_low, r_high) + v = sample_vector(v_low, v_high) + + return ssapy.Orbit(r, v, t) + + +def sample_GEO_orbit(t, r_low=4e7, r_high=5e7, v_low=2.7e3, v_high=3.3e3): + """ + Returns a ssapy.Orbit object with a distance near RGEO and a velocity near VGEO. + """ + return sample_orbit(t, r_low, r_high, v_low, v_high) + + +def sample_LEO_orbit(t, r_low=7e6, r_high=1e7, v_low=7.7e3, v_high=7.9e3): + """ + Returns a ssapy.Orbit object with a position near LEO and a velocity near VLEO. + """ + return sample_orbit(t, r_low, r_high, v_low, v_high) diff --git a/tests/test_accel.py b/tests/test_accel.py index 8dc1b49..1a1f00d 100644 --- a/tests/test_accel.py +++ b/tests/test_accel.py @@ -4,7 +4,7 @@ import ssapy from ssapy.utils import norm, iers_interp -from ssapy_test_helpers import timer +from ssapy_test_helpers import timer, sample_LEO_orbit, sample_GEO_orbit iers_interp(0.0) # Prime the IERS interpolant cache earth = ssapy.get_body("earth") @@ -715,25 +715,7 @@ def test_RK4_vs_analytic(): for _ in range(10): while True: - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - orbit = ssapy.Orbit(r, v, t0) + orbit = sample_GEO_orbit(t0) if norm(orbit.periapsis) > 1e7: break @@ -747,25 +729,7 @@ def test_RK4_vs_analytic(): times = t0 + np.linspace(0, 5, 1000)*u.h for _ in range(10): while True: - # Repeat near LEO - r = np.random.uniform(7e6, 1e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VLEO - v = np.random.uniform(7.7e3, 7.9e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - orbit = ssapy.Orbit(r, v, t0.gps) + orbit = sample_LEO_orbit(t0) if norm(orbit.periapsis) > 6475e3: break @@ -787,25 +751,7 @@ def test_scipy_propagator(): for _ in range(10): while True: - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - orbit = ssapy.Orbit(r, v, t0) + orbit = sample_GEO_orbit(t0) if norm(orbit.periapsis) > 1e7: break @@ -819,25 +765,7 @@ def test_scipy_propagator(): times = t0 + np.linspace(0, 5, 1000)*u.h for _ in range(10): while True: - # Repeat near LEO - r = np.random.uniform(7e6, 1e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VLEO - v = np.random.uniform(7.7e3, 7.9e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - orbit = ssapy.Orbit(r, v, t0.gps) + orbit = sample_LEO_orbit(t0) if norm(orbit.periapsis) > 6475e3: break @@ -860,25 +788,7 @@ def test_RK8(): for _ in range(10): while True: - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - orbit = ssapy.Orbit(r, v, t0) + orbit = sample_GEO_orbit(t0) if norm(orbit.periapsis) > 1e7: break @@ -895,25 +805,7 @@ def test_RK8(): times = t0 + np.arange(0, 500)*30*u.s for _ in range(10): while True: - # Repeat near LEO - r = np.random.uniform(7e6, 1e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VLEO - v = np.random.uniform(7.7e3, 7.9e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - orbit = ssapy.Orbit(r, v, t0.gps) + orbit = sample_LEO_orbit(t0) if norm(orbit.periapsis) > 6475e3: break @@ -939,25 +831,7 @@ def test_RK78(): for _ in range(10): while True: - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - orbit = ssapy.Orbit(r, v, t0) + orbit = sample_GEO_orbit(t0) if norm(orbit.periapsis) > 1e7: break @@ -975,25 +849,7 @@ def test_RK78(): times = t0 + np.arange(0, 500)*30*u.s for _ in range(10): while True: - # Repeat near LEO - r = np.random.uniform(7e6, 1e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VLEO - v = np.random.uniform(7.7e3, 7.9e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - orbit = ssapy.Orbit(r, v, t0.gps) + orbit = sample_LEO_orbit(t0) if norm(orbit.periapsis) > 6475e3: break diff --git a/tests/test_orbit.py b/tests/test_orbit.py index 944a666..96f6f23 100755 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -15,7 +15,7 @@ from ssapy.orbit import _ellipticalEccentricToMeanLongitude, _ellipticalMeanToEccentricLongitude from ssapy.orbit import _hyperbolicEccentricToMeanLongitude, _hyperbolicMeanToEccentricLongitude from ssapy.utils import normed, norm, teme_to_gcrf -from ssapy_test_helpers import timer, checkAngle, checkSphere +from ssapy_test_helpers import timer, checkAngle, checkSphere, sample_orbit, sample_LEO_orbit, sample_GEO_orbit try: @@ -142,28 +142,10 @@ def test_orbit_ctor(): eqElts = [] kElts = [] for _ in range(1000): - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - rs.append(r) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - vs.append(v) + orbit = sample_GEO_orbit(t=0.0) + rs.append(orbit.r) + vs.append(orbit.v) - orbit = ssapy.Orbit(r, v, 0.0) eqElts.append(orbit.equinoctialElements) kElts.append(orbit.keplerianElements) orbit2 = ssapy.Orbit.fromEquinoctialElements(*orbit.equinoctialElements, orbit.t) @@ -182,32 +164,14 @@ def test_orbit_ctor(): np.testing.assert_allclose(orb.hx, np.tan(orb.i/2)*np.cos(orb.raan)) np.testing.assert_allclose(orb.hy, np.tan(orb.i/2)*np.sin(orb.raan)) checkAngle(orb.lv, orb.trueAnomaly + lonPa, atol=1e-13) - np.testing.assert_allclose(r, orb.r, atol=1e-4, rtol=0) - np.testing.assert_allclose(v, orb.v, atol=1e-8, rtol=0) + np.testing.assert_allclose(orbit.r, orb.r, atol=1e-4, rtol=0) + np.testing.assert_allclose(orbit.v, orb.v, atol=1e-8, rtol=0) np.testing.assert_equal(0, orb.t) - # Repeat near LEO - r = np.random.uniform(7e6, 1e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - rs.append(r) - # Pick a velocity near VLEO - v = np.random.uniform(7.7e3, 7.9e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - vs.append(v) + orbit = sample_LEO_orbit(t=0.0) + rs.append(orbit.r) + vs.append(orbit.v) - orbit = ssapy.Orbit(r, v, 0.0) eqElts.append(orbit.equinoctialElements) kElts.append(orbit.keplerianElements) orbit2 = ssapy.Orbit.fromEquinoctialElements(*orbit.equinoctialElements, orbit.t) @@ -226,8 +190,8 @@ def test_orbit_ctor(): np.testing.assert_allclose(orb.hx, np.tan(orb.i/2)*np.cos(orb.raan)) np.testing.assert_allclose(orb.hy, np.tan(orb.i/2)*np.sin(orb.raan)) checkAngle(orb.lv, orb.trueAnomaly + lonPa, atol=1e-13) - np.testing.assert_allclose(r, orb.r, atol=1e-4, rtol=0) - np.testing.assert_allclose(v, orb.v, atol=1e-8, rtol=0) + np.testing.assert_allclose(orbit.r, orb.r, atol=1e-4, rtol=0) + np.testing.assert_allclose(orbit.v, orb.v, atol=1e-8, rtol=0) np.testing.assert_equal(0, orb.t) # Construct an "Orbit" with array arguments @@ -367,27 +331,8 @@ def test_orbit_rv(): np.random.seed(57721) for _ in range(1000): - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbit = ssapy.Orbit(r, v, 0.0) - + # Test near GEO + orbit = sample_GEO_orbit(t=0.0) true_anomalies = np.random.uniform(-2*np.pi, 2*np.pi, size=1) true_longitudes = true_anomalies + orbit.pa + orbit.raan r1, v1 = orbit._rvFromKeplerian(np.atleast_2d(true_anomalies)) @@ -396,31 +341,11 @@ def test_orbit_rv(): np.testing.assert_allclose(v1, v2, atol=1e-6, rtol=1e-11) # Repeat near LEO - r = np.random.uniform(7e6, 1e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VLEO - v = np.random.uniform(7.7e3, 7.9e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbit = ssapy.Orbit(r, v, 0.0) - + orbit = sample_LEO_orbit(t=0.0) true_anomalies = np.random.uniform(-2*np.pi, 2*np.pi, size=1) true_longitudes = true_anomalies + orbit.pa + orbit.raan r1, v1 = orbit._rvFromKeplerian(np.atleast_2d(true_anomalies)) r2, v2 = orbit._rvFromEquinoctial(lv=np.atleast_2d(true_longitudes)) - np.testing.assert_allclose(r1, r2, atol=1e-4, rtol=1e-13) np.testing.assert_allclose(v1, v2, atol=1e-6, rtol=1e-11) @@ -448,27 +373,8 @@ def vec3ToArr(vec3): t0 = Time('J2000.0') for i in range(100): - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbit = ssapy.Orbit(r, v, t0) + orbit = sample_GEO_orbit(t0) true_anomalies = np.random.uniform(-2*np.pi, 2*np.pi, size=1) # Now construct an orekit orbit and compare @@ -648,26 +554,8 @@ def test_rv(): orbits = [] for _ in range(NORBIT): - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbits.append(ssapy.Orbit(r, v, Time("J2000"))) + orbit = sample_GEO_orbit(t=Time("J2000")) + orbits.append(orbit) for prop in [ ssapy.KeplerianPropagator(), ssapy.SeriesPropagator(0), @@ -742,26 +630,8 @@ def test_groundTrack(): orbits = [] for _ in range(NORBIT): - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbits.append(ssapy.Orbit(r, v, Time("J2000"))) + orbit = sample_GEO_orbit(t=Time("J2000")) + orbits.append(orbit) for prop in [ ssapy.KeplerianPropagator(), ssapy.SeriesPropagator(0), @@ -834,26 +704,8 @@ def test_dircos(): # Just checking broadcasting behavior for the moment orbits = [] for _ in range(NORBIT): - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbits.append(ssapy.Orbit(r, v, Time("J2000"))) + orbit = sample_GEO_orbit(t=Time("J2000")) + orbits.append(orbit) times = Time("J2000") + np.linspace(-2, 2, NTIME)*u.year @@ -977,26 +829,8 @@ def test_radec(): # Just checking broadcasting behavior for the moment orbits = [] for _ in range(30): - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbits.append(ssapy.Orbit(r, v, Time("J2000"))) + orbit = sample_GEO_orbit(t=Time("J2000")) + orbits.append(orbit) times = Time("J2000") + np.linspace(-2, 2, NTIME)*u.year @@ -1151,33 +985,16 @@ def test_radecRate(): orbits = [] for _ in range(NORBIT): while True: - # Pick a distance between LEO and GEO - r = np.random.uniform(7e6, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near between LEO and GEO - v = np.random.uniform(2.7e3, 8.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - orbit = ssapy.Orbit(r, v, Time("J2000")) - # exclude if intersects earth or is unbound; + # For sampling orbits, pick a distance between LEO and GEO + orbit = sample_orbit(t=Time("J2000"), r_low=7e6, r_high=5e7, v_low=2.7e3, v_high=8.3e3) + # Exclude orbit if it intersects earth or is unbound if norm(orbit.periapsis) < 6.4e6: continue if orbit.energy > 0: continue break - orbits.append(ssapy.Orbit(r, v, Time("J2000"))) + orbits.append(orbit) times = Time("J2000") + np.linspace(-2, 2, NTIME)*u.year @@ -1368,26 +1185,8 @@ def test_altaz(): # Just checking broadcasting behavior for the moment orbits = [] for _ in range(NORBIT): - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbits.append(ssapy.Orbit(r, v, Time("J2000"))) + orbit = sample_GEO_orbit(t=Time("J2000")) + orbits.append(orbit) times = Time("J2000") + np.linspace(-2, 2, NTIME)*u.year @@ -1454,26 +1253,8 @@ def test_multiprocessing(): orbits = [] for _ in range(NORBIT): - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbits.append(ssapy.Orbit(r, v, Time("J2000"))) + orbit = sample_GEO_orbit(t=Time("J2000")) + orbits.append(orbit) time = Time("J2000") + np.linspace(-2, 2, NTIME)*u.year ff = f(time) @@ -1604,27 +1385,8 @@ def test_sgp4_vector(): np.random.seed(57721566490153 % 2**32) for _ in range(10): while True: - # Pick a distance near GEO - r = np.random.uniform(4e7, 5e7) - # Pick a random direction (not uniform on sphere) - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - x = r * np.cos(theta) * np.sin(phi) - y = r * np.sin(theta) * np.sin(phi) - z = r * np.cos(phi) - r = np.array([x, y, z]) - # Pick a velocity near VGEO - v = np.random.uniform(2.7e3, 3.3e3) - # and a randomish direction - theta = np.random.uniform(0, 2*np.pi) - phi = np.random.uniform(0, np.pi) - vx = v * np.cos(theta) * np.sin(phi) - vy = v * np.sin(theta) * np.sin(phi) - vz = v * np.cos(phi) - v = np.array([vx, vy, vz]) - - orbit = ssapy.Orbit(r, v, 0.0) - # Skip orbits with perigees close to Earth's surface + orbit = sample_GEO_orbit(t=0.0) + # Skip orbit if perigee is close to Earth's surface if np.sqrt(np.sum(orbit.periapsis**2)) > 1e7: break