Browse Source
- Added orbital_mechanics module with orbital_elements_to_cartesian() - Updated initialize_orbital_objects() to use orbital mechanics - Added validate_initial_positions() for post-initialization checking - Fixed test files to use global_position/global_velocity - Updated config loader to support spacecraft altitude parameter - Fixed orbital_mechanics.cpp velocity calculation bug (removed duplicate scaling) - Updated Makefile to include orbital_mechanics.o in test build - Renamed simulation.h OrbitalMetrics to OrbitalAnalysis to avoid conflict - Added docs/parabolic_union_implementation.md for parabolic orbit support plan Note: Test configs still need manual fix for orbit table TOML syntaxmain
13 changed files with 276 additions and 93 deletions
@ -0,0 +1,152 @@
|
||||
# Parabolic Orbit Union Implementation Plan |
||||
|
||||
## Overview |
||||
Add support for parabolic orbits (e≈1.0) using semi-latus rectum parameter `p` instead of the current hacky `semi_major_axis = 1.0e30` infinity approximation. |
||||
|
||||
## Problem |
||||
Current implementation uses `semi_major_axis = 1.0e30` to approximate infinity for parabolic orbits, causing: |
||||
1. Numerical precision issues with extremely large distances (~6.68e18 AU) |
||||
2. Velocities approaching zero (1.6e-08 km/s instead of ~42 km/s escape velocity) |
||||
3. Test failures due to floating-point equality (final_distance ≈ initial_distance) |
||||
|
||||
## Solution |
||||
Use a union in `OrbitalElements` struct to support both `semi_major_axis` (for elliptical/hyperbolic) and `semi_latus_rectum` (for parabolic). |
||||
|
||||
## Mathematical Background |
||||
For parabolic orbits (e=1.0), the semi-major axis is theoretically infinity. Using semi-latus rectum `p` is mathematically correct: |
||||
|
||||
Position: `r = p / (1 + cos(ν))` |
||||
Velocity: `v = √(2μ / r)` |
||||
|
||||
Where: |
||||
- `p` = semi-latus rectum |
||||
- `ν` = true anomaly |
||||
- `μ` = GM (gravitational parameter) |
||||
|
||||
For parabolic orbits: `p = 2q` where `q` is perihelion distance |
||||
|
||||
## Implementation Steps |
||||
|
||||
### Phase 1: Update OrbitalElements Struct |
||||
**File: `src/orbital_mechanics.h`** |
||||
|
||||
```cpp |
||||
struct OrbitalElements { |
||||
union { |
||||
double semi_major_axis; // for elliptical (e<1) and hyperbolic (e>1) |
||||
double semi_latus_rectum; // for parabolic (e≈1) |
||||
}; |
||||
double eccentricity; |
||||
double true_anomaly; |
||||
double inclination; |
||||
double longitude_of_ascending_node; |
||||
double argument_of_periapsis; |
||||
}; |
||||
``` |
||||
|
||||
### Phase 2: Update Config Loader |
||||
**File: `src/config_loader.cpp`** |
||||
|
||||
Add to `parse_toml_body()` and `parse_toml_spacecraft()`: |
||||
|
||||
1. Parse both `semi_major_axis` and `semi_latus_rectum` from orbit table |
||||
2. Initialize union field based on which is specified |
||||
3. Validate exactly one is present per eccentricity range |
||||
|
||||
**Validation Logic:** |
||||
```cpp |
||||
bool has_semi_major = (semi_major.type == TOML_FP64); |
||||
bool has_semi_latus = (semi_latus.type == TOML_FP64); |
||||
|
||||
if (fabs(elements.eccentricity - 1.0) < 0.005) { |
||||
// Parabolic orbit - requires semi_latus_rectum |
||||
if (!has_semi_latus) { |
||||
printf("Error: Parabolic orbit requires 'semi_latus_rectum'\n"); |
||||
return false; |
||||
} |
||||
if (has_semi_major) { |
||||
printf("Error: Parabolic orbit cannot have 'semi_major_axis'\n"); |
||||
return false; |
||||
} |
||||
elements.semi_latus_rectum = semi_latus.u.fp64; |
||||
} else { |
||||
// Elliptical or hyperbolic - requires semi_major_axis |
||||
if (!has_semi_major) { |
||||
printf("Error: Elliptical/hyperbolic orbit requires 'semi_major_axis'\n"); |
||||
return false; |
||||
} |
||||
if (has_semi_latus) { |
||||
printf("Error: Elliptical/hyperbolic orbit cannot have 'semi_latus_rectum'\n"); |
||||
return false; |
||||
} |
||||
elements.semi_major_axis = semi_major.u.fp64; |
||||
} |
||||
``` |
||||
|
||||
### Phase 3: Update orbital_mechanics.cpp |
||||
**File: `src/orbital_mechanics.cpp`** |
||||
|
||||
Update parabolic case (line 21-23): |
||||
|
||||
```cpp |
||||
} else if (fabs(e - 1.0) < 0.005) { |
||||
double p = elements.semi_latus_rectum; |
||||
r = p / (1.0 + cos(nu)); |
||||
v_mag = sqrt(2.0 * mu / r); |
||||
} |
||||
``` |
||||
|
||||
Remove the `2.0 * a` approximation that requires `a=1.0e30`. |
||||
|
||||
### Phase 4: Update Test Configs |
||||
**File: `tests/configs/parabolic_comet.toml`** |
||||
|
||||
Replace `semi_major_axis = 1.0e30` with `semi_latus_rectum = 1.496e11` (p = 1 AU): |
||||
|
||||
```toml |
||||
[[bodies]] |
||||
name = "ParabolicComet" |
||||
mass = 1.0e14 |
||||
radius = 5.0e3 |
||||
parent_index = 0 |
||||
color = { r = 0.7, g = 0.8, b = 0.9 } |
||||
orbit = { |
||||
semi_latus_rectum = 1.496e11, |
||||
eccentricity = 1.0, |
||||
true_anomaly = 0.0 |
||||
} |
||||
``` |
||||
|
||||
### Phase 5: Update Documentation |
||||
**File: `docs/technical_reference.md`** |
||||
|
||||
1. Update `OrbitalElements` struct documentation to show union |
||||
2. Add note about `semi_latus_rectum` being required for parabolic orbits (e≈1.0) |
||||
3. Document `semi_latus_rectum` in config format section |
||||
|
||||
**File: `docs/unified_orbital_elements_plan.md`** |
||||
|
||||
Mark union implementation as complete in Phase 7 status. |
||||
|
||||
## Validation Steps |
||||
|
||||
1. Build: `make clean && make` |
||||
2. Run parabolic test: `./orbit_test '[parabolic]'` |
||||
3. Verify velocity is correct: should be ~42,127 m/s escape velocity at 1 AU |
||||
4. Verify energy is ~0 (parabolic orbits have total energy = 0) |
||||
|
||||
## Decisions Made |
||||
|
||||
### Default Behavior |
||||
No backward compatibility for `semi_major_axis` on parabolic orbits - require explicit `semi_latus_rectum` for all parabolic configs. This is cleaner than trying to auto-convert `p = 2*a`. |
||||
|
||||
### Spacecraft Altitude Parameter |
||||
Spacecraft `altitude` parameter is not supported for parabolic orbits in this implementation. If user specifies `altitude` with `eccentricity ≈ 1.0`, the config loader will require `semi_latus_rectum` instead and reject `altitude` or `semi_major_axis`. Added to future todos for Phase 8+. |
||||
|
||||
### Parabolic Detection Tolerance |
||||
Using `|e - 1.0| < 0.005` as threshold for detecting parabolic orbits. This matches tolerance used elsewhere in the codebase. |
||||
|
||||
## Future Enhancements (TODO) |
||||
|
||||
- Spacecraft `altitude` parameter for parabolic orbits: parse `altitude` and convert to `semi_latus_rectum = parent_radius + altitude` when eccentricity is parabolic |
||||
- Consider adding explicit `perihelion` parameter to config file, then derive `semi_latus_rectum = 2 * perihelion` for parabolic orbits |
||||
Loading…
Reference in new issue