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.
118 lines
3.7 KiB
118 lines
3.7 KiB
#!/usr/bin/env python3 |
|
""" |
|
Precalculate expected values for test_parabolic_orbit. |
|
Simulates a parabolic comet orbiting the Sun for 300 days. |
|
""" |
|
|
|
import math |
|
import sys |
|
sys.path.insert(0, "scripts") |
|
from sim_engine import Simulator, vmag, G |
|
|
|
def main(): |
|
sim = Simulator("tests/test_parabolic_orbit.toml", dt=60.0) |
|
|
|
comet = sim.get_body("ParabolicComet") |
|
sun = sim.get_body("Sun") |
|
|
|
# Initial state |
|
r0 = vmag(comet.global_pos) |
|
v0 = vmag(comet.global_vel) |
|
mu = G * sun.mass |
|
escape_v0 = math.sqrt(2.0 * mu / r0) |
|
circular_v0 = math.sqrt(mu / r0) |
|
|
|
print(f"// === Initial Conditions (SI units) ===") |
|
print(f"// Distance: {r0:.6f} m ({r0 / 1.496e11:.6f} AU)") |
|
print(f"// Velocity: {v0:.6f} m/s ({v0 / 1000.0:.6f} km/s)") |
|
print(f"// Escape velocity: {escape_v0:.6f} m/s ({escape_v0 / 1000.0:.6f} km/s)") |
|
print(f"// Circular velocity: {circular_v0:.6f} m/s ({circular_v0 / 1000.0:.6f} km/s)") |
|
print(f"// Velocity error from escape: {(abs(v0 - escape_v0) / escape_v0) * 100.0:.6f}%") |
|
print(f"// Eccentricity: {comet.orbit.e:.6f}") |
|
print() |
|
|
|
# Energy at start (local frame, comet relative to sun) |
|
KE0 = 0.5 * comet.mass * v0**2 |
|
PE0 = -mu * comet.mass / r0 |
|
E0 = KE0 + PE0 |
|
|
|
print(f"// === Energy (Joules) ===") |
|
print(f"// Initial KE: {KE0:.6e}") |
|
print(f"// Initial PE: {PE0:.6e}") |
|
print(f"// Initial total E: {E0:.6e}") |
|
print() |
|
|
|
# Run simulation for 300 days |
|
total_seconds = 300.0 * 86400.0 |
|
steps = int(total_seconds / sim.dt) |
|
print(f"// Total steps: {steps}") |
|
print() |
|
|
|
# Record states every 1000 steps |
|
distances = [] |
|
velocities = [] |
|
energies = [] |
|
|
|
for i in range(steps): |
|
sim._step() |
|
|
|
if i % 1000 == 0: |
|
r = vmag(comet.global_pos) |
|
v = vmag(comet.global_vel) |
|
KE = 0.5 * comet.mass * v**2 |
|
PE = -mu * comet.mass / r |
|
E = KE + PE |
|
distances.append(r) |
|
velocities.append(v) |
|
energies.append(E) |
|
|
|
print(f"// Step {i}: t={sim.time/86400.0:.1f} days, r={r:.6f} m ({r/1.496e11:.4f} AU), v={v:.6f} m/s ({v/1000.0:.4f} km/s), E={E:.6e}") |
|
|
|
# Final state |
|
rf = vmag(comet.global_pos) |
|
vf = vmag(comet.global_vel) |
|
KEf = 0.5 * comet.mass * vf**2 |
|
PEf = -mu * comet.mass / rf |
|
Ef = KEf + PEf |
|
|
|
print() |
|
print(f"// === Final State (t=300 days) ===") |
|
print(f"// Distance: {rf:.6f} m ({rf / 1.496e11:.6f} AU)") |
|
print(f"// Velocity: {vf:.6f} m/s ({vf / 1000.0:.6f} km/s)") |
|
print(f"// Final KE: {KEf:.6e}") |
|
print(f"// Final PE: {PEf:.6e}") |
|
print(f"// Final total E: {Ef:.6e}") |
|
print() |
|
|
|
# Energy drift |
|
avg_KE = (KE0 + KEf) / 2.0 |
|
energy_drift = abs(Ef - E0) |
|
energy_drift_pct = (energy_drift / avg_KE) * 100.0 if avg_KE > 0 else 0.0 |
|
|
|
print(f"// === Energy Drift ===") |
|
print(f"// Absolute drift: {energy_drift:.6e} J") |
|
print(f"// Drift percent: {energy_drift_pct:.6f}%") |
|
print() |
|
|
|
# Velocity trend |
|
vel_decreases = 0 |
|
for i in range(1, len(velocities)): |
|
if velocities[i] < velocities[i-1]: |
|
vel_decreases += 1 |
|
total_checks = len(velocities) - 1 |
|
|
|
print(f"// === Velocity Trend ===") |
|
print(f"// Velocity decreases: {vel_decreases} / {total_checks}") |
|
print(f"// Ratio: {vel_decreases / total_checks:.2%}") |
|
print() |
|
|
|
# Assertions summary |
|
print(f"// === Assertions ===") |
|
print(f"// final_distance ({rf:.2f} m) > initial_distance ({r0:.2f} m): {rf > r0}") |
|
print(f"// final_velocity ({vf:.2f} m/s) < initial_velocity ({v0:.2f} m/s): {vf < v0}") |
|
print(f"// E0 >= -1e25: {E0 >= -1e25}") |
|
print(f"// energy_drift_pct < 1.0: {energy_drift_pct < 1.0}") |
|
print(f"// vel_decreases > total/2: {vel_decreases > total_checks // 2}") |
|
|
|
if __name__ == "__main__": |
|
main()
|
|
|