diff --git a/src/maneuver.cpp b/src/maneuver.cpp index 22589af..fe6da7c 100644 --- a/src/maneuver.cpp +++ b/src/maneuver.cpp @@ -110,17 +110,7 @@ OrbitalElements preview_burn_result(const Spacecraft* craft, BurnDirection direc return cartesian_to_orbital_elements(pos, new_vel, parent_mass); } -static double normalize_angle(double angle) { - while (angle < 0.0) angle += 2.0 * M_PI; - while (angle >= 2.0 * M_PI) angle -= 2.0 * M_PI; - return angle; -} - -static double angular_distance(double a, double b) { - double diff = fabs(normalize_angle(a) - normalize_angle(b)); - return (diff > M_PI) ? (2.0 * M_PI - diff) : diff; -} - +// Check if target angle lies between current and next angles (accounting for wraparound) static bool angle_between(double current, double next, double target) { double curr_norm = normalize_angle(current); double next_norm = normalize_angle(next); @@ -133,60 +123,17 @@ static bool angle_between(double current, double next, double target) { } } -static double calculate_true_anomaly(Vec3 r, Vec3 v, Vec3 e_vec, double e_mag, double r_mag) { - // For near-circular orbits, eccentricity vector is near-zero - // Compute true anomaly as the angle in the orbital plane - if (e_mag < 1e-10) { - Vec3 h = vec3_cross(r, v); - double h_mag = vec3_magnitude(h); - if (h_mag < 1e-10) return 0.0; - - // Create a coordinate system in the orbital plane - Vec3 z_hat = vec3_scale(h, 1.0 / h_mag); - // Choose x-axis as cross product of Z (world up) and orbit normal - // This gives a consistent reference direction in the orbital plane - Vec3 world_z = {0.0, 0.0, 1.0}; - Vec3 x_hat = vec3_cross(world_z, z_hat); - double x_hat_mag = vec3_magnitude(x_hat); - if (x_hat_mag < 1e-10) { - // Orbit is equatorial, use world X as reference - x_hat = (Vec3){1.0, 0.0, 0.0}; - } else { - x_hat = vec3_scale(x_hat, 1.0 / x_hat_mag); - } - Vec3 y_hat = vec3_cross(z_hat, x_hat); - - // Project position onto this orbital plane coordinate system - double x_proj = vec3_dot(r, x_hat); - double y_proj = vec3_dot(r, y_hat); - - // True anomaly is the angle in the orbital plane - double nu = atan2(y_proj, x_proj); - if (nu < 0) nu += 2.0 * M_PI; - return nu; - } - - // Standard calculation using eccentricity vector - double cos_nu = vec3_dot(e_vec, r) / (e_mag * r_mag); - cos_nu = fmax(-1.0, fmin(1.0, cos_nu)); - double nu = acos(cos_nu); - - // Determine correct quadrant using cross product - Vec3 r_cross_v = vec3_cross(r, v); - double r_cross_v_dot_e = vec3_dot(r_cross_v, e_vec); - if (r_cross_v_dot_e < 0) { - nu = 2.0 * M_PI - nu; - } - - return nu; -} - -static Vec3 calculate_eccentricity_vector(Vec3 r, Vec3 v, Vec3 h, double mu) { - Vec3 v_cross_h = vec3_cross(v, h); - Vec3 v_cross_h_over_mu = vec3_scale(v_cross_h, 1.0 / mu); - double r_mag = vec3_magnitude(r); - Vec3 r_over_mag = vec3_scale(r, 1.0 / r_mag); - return vec3_sub(v_cross_h_over_mu, r_over_mag); +// Propagate orbital state forward by dt and return new position/velocity +static void propagate_state_by_dt(OrbitalElements elements, double parent_mass, double dt, + Vec3* out_r, Vec3* out_v, Vec3* out_h, + Vec3* out_e_vec, double* out_e_mag) { + OrbitalElements future_elements = propagate_orbital_elements(elements, dt, parent_mass); + orbital_elements_to_cartesian(future_elements, parent_mass, out_r, out_v); + + *out_h = vec3_cross(*out_r, *out_v); + double mu = G * parent_mass; + *out_e_vec = calculate_eccentricity_vector(*out_r, *out_v, *out_h, mu); + *out_e_mag = vec3_magnitude(*out_e_vec); } bool check_maneuver_trigger(Maneuver* maneuver, Spacecraft* craft, SimulationState* sim) { @@ -223,18 +170,15 @@ bool check_maneuver_trigger(Maneuver* maneuver, Spacecraft* craft, SimulationSta if (current_diff < 0.01) return true; // Propagate orbit forward by one time step to check if we'll cross trigger - OrbitalElements future_elements = propagate_orbital_elements(craft->orbit, sim->dt, parent->mass); - Vec3 future_r, future_v; - orbital_elements_to_cartesian(future_elements, parent->mass, &future_r, &future_v); + Vec3 future_r, future_v, future_h, future_e_vec; + double future_e_mag; + propagate_state_by_dt(craft->orbit, parent->mass, sim->dt, + &future_r, &future_v, &future_h, + &future_e_vec, &future_e_mag); double future_r_mag = vec3_magnitude(future_r); if (future_r_mag < 1.0) return false; - // Calculate future eccentricity vector for true anomaly calculation - Vec3 future_h = vec3_cross(future_r, future_v); - Vec3 future_e_vec = calculate_eccentricity_vector(future_r, future_v, future_h, mu); - double future_e_mag = vec3_magnitude(future_e_vec); - // Calculate future true anomaly double future_nu = calculate_true_anomaly(future_r, future_v, future_e_vec, future_e_mag, future_r_mag); double future_nu_norm = normalize_angle(future_nu); diff --git a/src/orbital_mechanics.cpp b/src/orbital_mechanics.cpp index 5258eec..51424a6 100644 --- a/src/orbital_mechanics.cpp +++ b/src/orbital_mechanics.cpp @@ -373,3 +373,74 @@ OrbitalElements propagate_orbital_elements(const OrbitalElements& elements, doub return result; } } + +// Normalize angle to [0, 2π) range +double normalize_angle(double angle) { + while (angle < 0.0) angle += 2.0 * M_PI; + while (angle >= 2.0 * M_PI) angle -= 2.0 * M_PI; + return angle; +} + +// Calculate shortest angular distance between two angles (always positive, range [0, π]) +double angular_distance(double a, double b) { + double diff = fabs(normalize_angle(a) - normalize_angle(b)); + return (diff > M_PI) ? (2.0 * M_PI - diff) : diff; +} + +// Calculate eccentricity vector from state vectors +Vec3 calculate_eccentricity_vector(Vec3 r, Vec3 v, Vec3 h, double mu) { + Vec3 v_cross_h = vec3_cross(v, h); + Vec3 v_cross_h_over_mu = vec3_scale(v_cross_h, 1.0 / mu); + double r_mag = vec3_magnitude(r); + Vec3 r_over_mag = vec3_scale(r, 1.0 / r_mag); + return vec3_sub(v_cross_h_over_mu, r_over_mag); +} + +// Calculate true anomaly from position and velocity vectors +double calculate_true_anomaly(Vec3 r, Vec3 v, Vec3 e_vec, double e_mag, double r_mag) { + // For near-circular orbits, eccentricity vector is near-zero + // Compute true anomaly as the angle in the orbital plane + if (e_mag < 1e-10) { + Vec3 h = vec3_cross(r, v); + double h_mag = vec3_magnitude(h); + if (h_mag < 1e-10) return 0.0; + + // Create a coordinate system in the orbital plane + Vec3 z_hat = vec3_scale(h, 1.0 / h_mag); + // Choose x-axis as cross product of Z (world up) and orbit normal + // This gives a consistent reference direction in the orbital plane + Vec3 world_z = {0.0, 0.0, 1.0}; + Vec3 x_hat = vec3_cross(world_z, z_hat); + double x_hat_mag = vec3_magnitude(x_hat); + if (x_hat_mag < 1e-10) { + // Orbit is equatorial, use world X as reference + x_hat = (Vec3){1.0, 0.0, 0.0}; + } else { + x_hat = vec3_scale(x_hat, 1.0 / x_hat_mag); + } + Vec3 y_hat = vec3_cross(z_hat, x_hat); + + // Project position onto this orbital plane coordinate system + double x_proj = vec3_dot(r, x_hat); + double y_proj = vec3_dot(r, y_hat); + + // True anomaly is the angle in the orbital plane + double nu = atan2(y_proj, x_proj); + if (nu < 0) nu += 2.0 * M_PI; + return nu; + } + + // Standard calculation using eccentricity vector + double cos_nu = vec3_dot(e_vec, r) / (e_mag * r_mag); + cos_nu = fmax(-1.0, fmin(1.0, cos_nu)); + double nu = acos(cos_nu); + + // Determine correct quadrant using cross product + Vec3 r_cross_v = vec3_cross(r, v); + double r_cross_v_dot_e = vec3_dot(r_cross_v, e_vec); + if (r_cross_v_dot_e < 0) { + nu = 2.0 * M_PI - nu; + } + + return nu; +} diff --git a/src/orbital_mechanics.h b/src/orbital_mechanics.h index fb191ce..d882a06 100644 --- a/src/orbital_mechanics.h +++ b/src/orbital_mechanics.h @@ -48,4 +48,14 @@ double solve_barker_equation(double mean_anomaly); OrbitalElements propagate_orbital_elements(const OrbitalElements& elements, double dt, double parent_mass); +// Angle utilities +double normalize_angle(double angle); +double angular_distance(double a, double b); + +// Calculate eccentricity vector from state vectors +Vec3 calculate_eccentricity_vector(Vec3 r, Vec3 v, Vec3 h, double mu); + +// Calculate true anomaly from position and velocity vectors +double calculate_true_anomaly(Vec3 r, Vec3 v, Vec3 e_vec, double e_mag, double r_mag); + #endif