Fix orbital_elements_to_cartesian: circular orbit velocity and refactor for clarity

src/orbital_mechanics.cpp

@@ -13,16 +13,20 @@ void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass,
     double r, v_mag;
 
     if (fabs(e) < 1e-10) {
+        // Circular orbit
         r = a;
         v_mag = sqrt(mu / a);
     } else if (e < 1.0) {
+        // Elliptical orbit
         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) < 0.005) {
+        // Parabolic orbit
         double p = elements.semi_latus_rectum;
         r = p / (1.0 + cos(nu));
         v_mag = sqrt(2.0 * mu / r);
     } else {
+        // Hyperbolic orbit
         r = a * (1.0 - e * e) / (1.0 + e * cos(nu));
         v_mag = sqrt(mu * (2.0 / r - 1.0 / a));
     }
@@ -35,21 +39,25 @@ void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass,
     double sin_nu = sin(nu);
     double cos_nu = cos(nu);
 
+    double p;
+    if (fabs(e - 1.0) < 0.005) {
+        p = elements.semi_latus_rectum;
+    } else {
+        p = a * (1.0 - e * e);
+    }
+
     double vx_orbital, vy_orbital;
 
     if (fabs(e) < 1e-10) {
-        vx_orbital = 0.0;
+        // Circular orbit: velocity rotates with position
-        vy_orbital = v_mag;
+        vx_orbital = -v_mag * sin_nu;
-    } else if (e < 1.0) {
+        vy_orbital = v_mag * cos_nu;
-        double p = a * (1.0 - e * e);
-        vx_orbital = -sqrt(mu / p) * sin_nu;
-        vy_orbital = sqrt(mu / p) * (e + cos_nu);
     } else if (fabs(e - 1.0) < 0.005) {
-        double p = elements.semi_latus_rectum;
+        // Parabolic orbit: use (1 + cos_nu) term
         vx_orbital = -sqrt(mu / p) * sin_nu;
         vy_orbital = sqrt(mu / p) * (1.0 + cos_nu);
     } else {
-        double p = a * (1.0 - e * e);
+        // Elliptical or hyperbolic orbit: use (e + cos_nu) term
         vx_orbital = -sqrt(mu / p) * sin_nu;
         vy_orbital = sqrt(mu / p) * (e + cos_nu);
     }