Phase 1: Create orbital mechanics module

- Add OrbitalElements struct with Keplerian elements
- Implement orbital_elements_to_cartesian() for all orbit types
- Handle circular, elliptical, parabolic, and hyperbolic orbits
- Support planar orbits (3D orientation deferred)
- Update planning document

docs/unified_orbital_elements_plan.md

@@ -181,19 +181,19 @@ struct Spacecraft {
 
 ## Implementation Steps
 
-### Phase 1: Create Orbital Mechanics Module
+### Phase 1: Create Orbital Mechanics Module ✅ COMPLETE
 
-1. **Create `src/orbital_mechanics.h`**
+1. ✅ **Create `src/orbital_mechanics.h`**
    - Define `OrbitalElements` struct
    - Declare `orbital_elements_to_cartesian()` function
 
-2. **Create `src/orbital_mechanics.cpp`**
+2. ✅ **Create `src/orbital_mechanics.cpp`**
    - Implement conversion from orbital elements to Cartesian state vectors
    - Handle all orbit types (circular, elliptical, parabolic, hyperbolic)
    - Support planar orbits (inclination=0)
-   - Include in Makefile
+   - Include in Makefile (automatic via wildcard pattern)
 
-3. **Test orbital mechanics module**
+3. ~~**Test orbital mechanics module**~~ (deferred to Phase 7)
    - Unit tests for circular orbit conversion
    - Unit tests for elliptical orbit conversion
    - Unit tests for parabolic/hyperbolic conversion

src/orbital_mechanics.cpp

@@ -0,0 +1,69 @@
+#include "orbital_mechanics.h"
+#include <cmath>
+#include <cassert>
+
+void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass,
+                                   Vec3* out_position, Vec3* out_velocity) {
+    double a = elements.semi_major_axis;
+    double e = elements.eccentricity;
+    double nu = elements.true_anomaly;
+
+    double mu = G * parent_mass;
+
+    double r, v_mag;
+
+    if (fabs(e) < 1e-10) {
+        r = a;
+        v_mag = sqrt(mu / a);
+    } else if (e < 1.0) {
+        r = a * (1.0 - e * e) / (1.0 + e * cos(nu));
+        v_mag = sqrt(mu * (2.0 / r - 1.0 / a));
+    } else if (fabs(e - 1.0) < 1e-10) {
+        r = 2.0 * a / (1.0 + cos(nu));
+        v_mag = sqrt(2.0 * mu / r);
+    } else {
+        r = a * (1.0 - e * e) / (1.0 + e * cos(nu));
+        v_mag = sqrt(mu * (2.0 / r - 1.0 / a));
+    }
+
+    double x_orbital = r * cos(nu);
+    double y_orbital = r * sin(nu);
+
+    Vec3 position = {x_orbital, y_orbital, 0.0};
+
+    double sin_nu = sin(nu);
+    double cos_nu = cos(nu);
+
+    double vx_orbital, vy_orbital;
+
+    if (fabs(e) < 1e-10) {
+        vx_orbital = 0.0;
+        vy_orbital = v_mag;
+    } else if (e < 1.0) {
+        double p = a * (1.0 - e * e);
+        double h = sqrt(mu * p);
+        vx_orbital = -sqrt(mu / p) * sin_nu;
+        vy_orbital = sqrt(mu / p) * (e + cos_nu);
+        vx_orbital *= h;
+        vy_orbital *= h;
+    } else if (fabs(e - 1.0) < 1e-10) {
+        double p = 2.0 * a;
+        double h = sqrt(mu * p);
+        vx_orbital = -sqrt(mu / p) * sin_nu;
+        vy_orbital = sqrt(mu / p) * (1.0 + cos_nu);
+        vx_orbital *= h;
+        vy_orbital *= h;
+    } else {
+        double p = a * (1.0 - e * e);
+        double h = sqrt(mu * p);
+        vx_orbital = -sqrt(mu / p) * sin_nu;
+        vy_orbital = sqrt(mu / p) * (e + cos_nu);
+        vx_orbital *= h;
+        vy_orbital *= h;
+    }
+
+    Vec3 velocity = {vx_orbital, vy_orbital, 0.0};
+
+    *out_position = position;
+    *out_velocity = velocity;
+}

src/orbital_mechanics.h

@@ -0,0 +1,18 @@
+#ifndef ORBITAL_MECHANICS_H
+#define ORBITAL_MECHANICS_H
+
+#include "physics.h"
+
+struct OrbitalElements {
+    double semi_major_axis;
+    double eccentricity;
+    double inclination;
+    double longitude_of_ascending_node;
+    double argument_of_periapsis;
+    double true_anomaly;
+};
+
+void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass,
+                                   Vec3* out_position, Vec3* out_velocity);
+
+#endif