You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
250 lines
10 KiB
250 lines
10 KiB
#!/usr/env python3 |
|
""" |
|
Precalculate expected values for test_cartesian_to_elements_advanced.cpp. |
|
|
|
Replicates all test cases: convert elements -> cartesian -> back to elements, |
|
then report round-trip errors for each assertion. |
|
|
|
Usage: |
|
python3 scripts/precalc_cartesian_to_elements_advanced.py |
|
""" |
|
|
|
import sys, math |
|
sys.path.insert(0, 'scripts') |
|
from sim_engine import orbital_to_cartesian, cartesian_to_orbital_elements, vmag, OrbitalElements, G, normalize_angle |
|
|
|
M_sun = 1.989e30 |
|
mu = G * M_sun |
|
|
|
def make_elements(a, e, nu, inc, Omega, omega, semi_latus_rectum=None): |
|
el = OrbitalElements(a=a, e=e, nu=nu, inc=inc, Omega=Omega, omega=omega) |
|
if semi_latus_rectum is not None: |
|
el.p = semi_latus_rectum # use 'p' field for semi_latus_rectum |
|
return el |
|
|
|
def roundtrip(el): |
|
pos, vel = orbital_to_cartesian(el, M_sun) |
|
return cartesian_to_orbital_elements(pos, vel, M_sun) |
|
|
|
def ang_diff(a, b): |
|
"""Shortest angular distance.""" |
|
return abs(normalize_angle(a) - normalize_angle(b)) |
|
|
|
def report(name, original, recovered, fields): |
|
"""Print round-trip errors for specified fields.""" |
|
print(f" # {name}:") |
|
for field in fields: |
|
orig_val = getattr(original, field) |
|
rec_val = getattr(recovered, field) |
|
err = abs(orig_val - rec_val) |
|
print(f" {field:20s} = {orig_val:20.15e} -> {rec_val:20.15e} error = {err:.2e}") |
|
|
|
# ============================================================================= |
|
# SECTION: eccentricity spectrum |
|
# ============================================================================= |
|
print("=" * 70) |
|
print("SECTION: eccentricity spectrum: circular to highly hyperbolic") |
|
print("=" * 70) |
|
|
|
r = 1.496e11 |
|
v_circ = math.sqrt(mu / r) |
|
|
|
# 1. Circular orbit (e=0) |
|
circular = make_elements(r, 0.0, 0.0, 0.0, 0.0, 0.0) |
|
rec_circ = roundtrip(circular) |
|
print("\n [1] Circular orbit (e=0):") |
|
print(f" ecc error = {abs(rec_circ.e - 0.0):.2e} (test tol: 1e-10)") |
|
print(f" a error = {abs(rec_circ.a - r):.2e} (test tol: 1e-2)") |
|
|
|
# 2. Near-circular (e=0.001) |
|
near_circ = make_elements(1.496e11, 0.001, 0.5, 0.0, 0.0, 0.0) |
|
rec_near_circ = roundtrip(near_circ) |
|
print(f"\n [2] Near-circular (e=0.001):") |
|
print(f" ecc error = {abs(rec_near_circ.e - 0.001):.2e} (test tol: 1e-6)") |
|
print(f" a error = {abs(rec_near_circ.a - 1.496e11):.2e} (test tol: 1e-2)") |
|
|
|
# 3. Elliptical (e=0.5) |
|
elliptical = make_elements(1.0e11, 0.5, 0.8, 0.0, 0.0, 0.0) |
|
rec_elliptical = roundtrip(elliptical) |
|
print(f"\n [3] Elliptical (e=0.5):") |
|
print(f" ecc error = {abs(rec_elliptical.e - 0.5):.2e} (test tol: 1e-4)") |
|
print(f" a error = {abs(rec_elliptical.a - 1.0e11):.2e} (test tol: 1e-2)") |
|
|
|
# 4. Highly elliptical (e=0.95) |
|
high_ell = make_elements(1.0e11, 0.95, 0.1, 0.0, 0.0, 0.0) |
|
rec_high_ell = roundtrip(high_ell) |
|
print(f"\n [4] Highly elliptical (e=0.95):") |
|
print(f" ecc error = {abs(rec_high_ell.e - 0.95):.2e} (test tol: 1e-3)") |
|
print(f" a error = {abs(rec_high_ell.a - 1.0e11):.2e} (test tol: 1e-2)") |
|
|
|
# 5. Near-parabolic (e=0.999) |
|
near_par = make_elements(1.0e11, 0.999, 0.05, 0.0, 0.0, 0.0) |
|
rec_near_par = roundtrip(near_par) |
|
print(f"\n [5] Near-parabolic (e=0.999):") |
|
print(f" ecc error = {abs(rec_near_par.e - 0.999):.2e} (test tol: 1e-3)") |
|
|
|
# 6. Parabolic (e=1.0) |
|
parabolic = make_elements(0.0, 1.0, 0.5, 0.0, 0.0, 0.0, semi_latus_rectum=1.0e11) |
|
rec_parabolic = roundtrip(parabolic) |
|
print(f"\n [6] Parabolic (e=1.0):") |
|
print(f" ecc error = {abs(rec_parabolic.e - 1.0):.2e} (test tol: 1e-2)") |
|
# semi_latus_rectum is stored in 'p' field in the script |
|
print(f" p error = {abs(rec_parabolic.p - 1.0e11):.2e} (test tol: 1e-2)") |
|
|
|
# 7. Hyperbolic (e=2.0) |
|
hyper = make_elements(-1.0e11, 2.0, 0.5, 0.0, 0.0, 0.0) |
|
rec_hyper = roundtrip(hyper) |
|
print(f"\n [7] Hyperbolic (e=2.0):") |
|
print(f" ecc error = {abs(rec_hyper.e - 2.0):.2e} (test tol: 1e-3)") |
|
print(f" a error = {abs(rec_hyper.a - (-1.0e11)):.2e} (test tol: 1e-2)") |
|
|
|
# 8. Highly hyperbolic (e=10.0) |
|
high_hyper = make_elements(-1.0e10, 10.0, 0.8, 0.0, 0.0, 0.0) |
|
rec_high_hyper = roundtrip(high_hyper) |
|
print(f"\n [8] Highly hyperbolic (e=10.0):") |
|
print(f" ecc error = {abs(rec_high_hyper.e - 10.0):.2e} (test tol: 1e-3)") |
|
print(f" a error = {abs(rec_high_hyper.a - (-1.0e10)):.2e} (test tol: 1e-2)") |
|
|
|
# ============================================================================= |
|
# SECTION: inclination |
|
# ============================================================================= |
|
print("\n" + "=" * 70) |
|
print("SECTION: inclination: zero, polar, and retrograde") |
|
print("=" * 70) |
|
|
|
# 1. Zero inclination |
|
eq = make_elements(1.0e11, 0.3, 0.5, 0.0, 0.0, 0.0) |
|
rec_eq = roundtrip(eq) |
|
print(f"\n [1] Equatorial (inc=0):") |
|
print(f" inc error = {abs(rec_eq.inc - 0.0):.2e} (test tol: 1e-6)") |
|
print(f" ecc error = {abs(rec_eq.e - 0.3):.2e} (test tol: 1e-4)") |
|
|
|
# 2. Polar (inc=90 deg) |
|
polar = make_elements(1.0e11, 0.2, 0.6, math.pi / 2.0, 0.5, 0.3) |
|
rec_polar = roundtrip(polar) |
|
print(f"\n [2] Polar (inc=90 deg):") |
|
print(f" inc error = {abs(rec_polar.inc - math.pi/2.0):.2e} (test tol: 1e-4)") |
|
print(f" Omega error = {abs(rec_polar.Omega - 0.5):.2e} (test tol: 1e-4)") |
|
print(f" omega error = {abs(rec_polar.omega - 0.3):.2e} (test tol: 1e-4)") |
|
|
|
# 3. Retrograde (inc=180 deg) |
|
retro = make_elements(1.0e11, 0.2, 0.6, math.pi, 0.5, 0.3) |
|
rec_retro = roundtrip(retro) |
|
print(f"\n [3] Retrograde (inc=180 deg):") |
|
print(f" inc error = {abs(rec_retro.inc - math.pi):.2e} (test tol: 1e-4)") |
|
|
|
# ============================================================================= |
|
# SECTION: true anomaly at key orbital positions |
|
# ============================================================================= |
|
print("\n" + "=" * 70) |
|
print("SECTION: true anomaly at key orbital positions") |
|
print("=" * 70) |
|
|
|
nu_tests = [ |
|
(0.0, 0.0, "periapsis"), |
|
(math.pi, math.pi, "apoapsis"), |
|
(math.pi / 2.0, math.pi / 2.0, "quadrature +90"), |
|
(-math.pi / 2.0, 3.0 * math.pi / 2.0, "quadrature -90"), |
|
(3.0 * math.pi / 2.0, 3.0 * math.pi / 2.0, "quadrature +270"), |
|
(-3.0 * math.pi / 2.0, math.pi / 2.0, "quadrature -270"), |
|
] |
|
|
|
for i, (nu_in, nu_exp, label) in enumerate(nu_tests): |
|
el = make_elements(1.0e11, 0.5, nu_in, 0.0, 0.0, 0.0) |
|
rec = roundtrip(el) |
|
nu_err = abs(rec.nu - nu_exp) |
|
e_err = abs(rec.e - 0.5) |
|
print(f"\n [{i+1}] {label} (input nu={nu_in:.6f}):") |
|
print(f" nu error = {nu_err:.2e} (test tol: 1e-6)") |
|
print(f" ecc error = {e_err:.2e} (test tol: 1e-4)") |
|
|
|
# ============================================================================= |
|
# SECTION: quadrature at various eccentricities |
|
# ============================================================================= |
|
print("\n" + "=" * 70) |
|
print("SECTION: quadrature at various eccentricities") |
|
print("=" * 70) |
|
|
|
e_tests = [(0.9, 1e-3, 1e-5), (0.1, 1e-5, 1e-6)] |
|
for i, (e, e_tol, nu_tol) in enumerate(e_tests): |
|
el = make_elements(1.0e11, e, math.pi / 2.0, 0.0, 0.0, 0.0) |
|
rec = roundtrip(el) |
|
e_err = abs(rec.e - e) |
|
a_err = abs(rec.a - 1.0e11) |
|
nu_err = abs(rec.nu - math.pi / 2.0) |
|
print(f"\n [{i+1}] e={e}:") |
|
print(f" ecc error = {e_err:.2e} (test tol: {e_tol:.0e}) {'PASS' if e_err <= e_tol else 'FAIL'}") |
|
print(f" a error = {a_err:.2e} (test tol: 1e-2)") |
|
print(f" nu error = {nu_err:.2e} (test tol: {nu_tol:.0e}) {'PASS' if nu_err <= nu_tol else 'FAIL'}") |
|
|
|
# ============================================================================= |
|
# SECTION: large true anomaly values |
|
# ============================================================================= |
|
print("\n" + "=" * 70) |
|
print("SECTION: large true anomaly values") |
|
print("=" * 70) |
|
|
|
large_nu_tests = [ |
|
(5.0, 5.0, 1e-6, "nu=5.0"), |
|
(-5.0, 1.28318530717958623, 1e-6, "nu=-5.0"), |
|
(10.0, 10.0 - 2.0 * math.pi, 1e-5, "nu=10.0"), |
|
] |
|
|
|
for i, (nu_in, nu_exp, tol, label) in enumerate(large_nu_tests): |
|
el = make_elements(1.0e11, 0.5, nu_in, 0.0, 0.0, 0.0) |
|
rec = roundtrip(el) |
|
nu_err = abs(rec.nu - nu_exp) |
|
e_err = abs(rec.e - 0.5) |
|
a_err = abs(rec.a - 1.0e11) |
|
print(f"\n [{i+1}] {label}:") |
|
print(f" ecc error = {e_err:.2e} (test tol: 1e-4)") |
|
print(f" a error = {a_err:.2e} (test tol: 1e-2)") |
|
print(f" nu error = {nu_err:.2e} (test tol: {tol:.0e}) {'PASS' if nu_err <= tol else 'FAIL'}") |
|
|
|
# ============================================================================= |
|
# SECTION: 3D orientation with quadrature point |
|
# ============================================================================= |
|
print("\n" + "=" * 70) |
|
print("SECTION: 3D orientation with quadrature point") |
|
print("=" * 70) |
|
|
|
el = make_elements(1.0e11, 0.5, math.pi / 2.0, math.pi / 3.0, math.pi / 4.0, math.pi / 6.0) |
|
rec = roundtrip(el) |
|
print(f" ecc error = {abs(rec.e - 0.5):.2e} (test tol: 1e-4)") |
|
print(f" a error = {abs(rec.a - 1.0e11):.2e} (test tol: 1e-2)") |
|
print(f" nu error = {abs(rec.nu - math.pi/2.0):.2e} (test tol: 1e-5)") |
|
print(f" inc error = {abs(rec.inc - math.pi/3.0):.2e} (test tol: 1e-4)") |
|
print(f" Omega error = {abs(rec.Omega - math.pi/4.0):.2e} (test tol: 1e-4)") |
|
print(f" omega error = {abs(rec.omega - math.pi/6.0):.2e} (test tol: 1e-4)") |
|
|
|
# ============================================================================= |
|
# SECTION: multiple true anomaly points in sequence |
|
# ============================================================================= |
|
print("\n" + "=" * 70) |
|
print("SECTION: multiple true anomaly points in sequence") |
|
print("=" * 70) |
|
|
|
nu_seq = [0.0, math.pi / 4.0, math.pi / 2.0, 3.0 * math.pi / 4.0, math.pi] |
|
for i, nu in enumerate(nu_seq): |
|
el = make_elements(1.0e11, 0.5, nu, 0.0, 0.0, 0.0) |
|
rec = roundtrip(el) |
|
nu_err = abs(rec.nu - nu) |
|
e_err = abs(rec.e - 0.5) |
|
a_err = abs(rec.a - 1.0e11) |
|
print(f"\n [{i+1}] nu={nu:.6f}:") |
|
print(f" ecc error = {e_err:.2e} (test tol: 1e-4)") |
|
print(f" a error = {a_err:.2e} (test tol: 1e-2)") |
|
print(f" nu error = {nu_err:.2e} (test tol: 1e-6)") |
|
|
|
# ============================================================================= |
|
# SECTION: hyperbolic orbit at quadrature point |
|
# ============================================================================= |
|
print("\n" + "=" * 70) |
|
print("SECTION: hyperbolic orbit at quadrature point") |
|
print("=" * 70) |
|
|
|
el = make_elements(-1.0e11, 2.0, math.pi / 2.0, 0.0, 0.0, 0.0) |
|
rec = roundtrip(el) |
|
print(f" ecc error = {abs(rec.e - 2.0):.2e} (test tol: 1e-3)") |
|
print(f" a error = {abs(rec.a - (-1.0e11)):.2e} (test tol: 1e-2)") |
|
print(f" nu error = {abs(rec.nu - math.pi/2.0):.2e} (test tol: 1e-5)")
|
|
|