cinnaboot
/
claudes_game
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Move orbital calculation helpers to orbital_mechanics module, add propagate_state_by_dt helper

main
cinnaboot 5 months ago
parent
db802f1983
commit
e23a1c874b
3 changed files with 98 additions and 73 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 90
      src/maneuver.cpp
  2. 71
      src/orbital_mechanics.cpp
  3. 10
      src/orbital_mechanics.h

90
src/maneuver.cpp
Unescape View File

@ -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);

71
src/orbital_mechanics.cpp
Unescape View File

@ -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;
}

10
src/orbital_mechanics.h
Unescape View File

@ -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

Write Preview
Loading…
Cancel
Save