@ -458,52 +458,3 @@ Vec3 calculate_eccentricity_vector(Vec3 r, Vec3 v, Vec3 h, double mu) {
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 ;
}