diff --git a/configs/test_simple.txt b/configs/test_simple.txt index 04eaa4b..8d7885f 100644 --- a/configs/test_simple.txt +++ b/configs/test_simple.txt @@ -1,22 +1,28 @@ # Simple Test Configuration -# Two planets in circular orbits for easy position verification +# Planets in circular and eccentric orbits # # Orbital periods (for verification): # Earth: ~365 days # Mars: ~687 days +# Comet: ~1444 days # # At t=0: # Earth starts at (1.0 AU, 0, 0) = (1.496e11 m, 0, 0) # Mars starts at (1.5 AU, 0, 0) = (2.244e11 m, 0, 0) -# Both orbit counter-clockwise when viewed from +Z +# Comet starts at perihelion (0.75 AU, 0, 0) = (1.122e11 m, 0, 0) +# All orbit counter-clockwise when viewed from +Z # -# Format: name mass(kg) radius(m) x(m) y(m) z(m) parent_index r g b +# Format: name mass(kg) radius(m) x(m) y(m) z(m) parent_index r g b eccentricity semi_major_axis(m) # Sun at origin (index 0) -Sun 1.989e30 6.96e8 0 0 0 -1 1.0 1.0 0.0 +Sun 1.989e30 6.96e8 0 0 0 -1 1.0 1.0 0.0 0 0 -# Earth at 1 AU on positive X axis (index 1, parent Sun) -Earth 5.972e24 6.371e6 1.496e11 0 0 0 0.0 0.5 1.0 +# Earth at 1 AU - circular orbit (index 1, parent Sun) +Earth 5.972e24 6.371e6 1.496e11 0 0 0 0.0 0.5 1.0 0 1.496e11 -# Mars at 1.5 AU on positive X axis (index 2, parent Sun) -Mars 6.39e23 3.3895e6 2.244e11 0 0 0 0.8 0.3 0.1 +# Mars at 1.5 AU - circular orbit (index 2, parent Sun) +Mars 6.39e23 3.3895e6 2.244e11 0 0 0 0.8 0.3 0.1 0 2.244e11 + +# Comet - eccentric orbit (e=0.7, a=2.5 AU) +# At perihelion: 0.75 AU, aphelion: 4.25 AU +Comet 1e14 5e3 1.122e11 0 0 0 0.5 0.8 1.0 0.7 3.74e11 diff --git a/src/config_loader.cpp b/src/config_loader.cpp index cb0828e..c29abb7 100644 --- a/src/config_loader.cpp +++ b/src/config_loader.cpp @@ -5,20 +5,29 @@ // Parse a single body definition line bool parse_body_line(const char* line, char* name, double* mass, double* radius, - Vec3* pos, int* parent_index, float* r, float* g, float* b) { + Vec3* pos, int* parent_index, float* r, float* g, float* b, + double* eccentricity, double* semi_major_axis) { // Skip empty lines and comments if (line[0] == '\0' || line[0] == '#' || line[0] == '\n') { return false; } - // Parse: name mass(kg) radius(m) x(m) y(m) z(m) parent_index r g b - int result = sscanf(line, "%63s %lf %lf %lf %lf %lf %d %f %f %f", + // Format: name mass radius x y z parent_index r g b eccentricity semi_major_axis + int result = sscanf(line, "%63s %lf %lf %lf %lf %lf %d %f %f %f %lf %lf", name, mass, radius, &pos->x, &pos->y, &pos->z, - parent_index, r, g, b); + parent_index, r, g, b, eccentricity, semi_major_axis); - return result == 10; // All fields must be present + return result == 12; } +struct OrbitParams { + double eccentricity; + double semi_major_axis; +}; + +// Forward declaration +void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params); + // Load system configuration from file bool load_system_config(SimulationState* sim, const char* filepath) { FILE* file = fopen(filepath, "r"); @@ -33,11 +42,18 @@ bool load_system_config(SimulationState* sim, const char* filepath) { Vec3 pos; int parent_index; float r, g, b; + double eccentricity, semi_major_axis; + + OrbitParams orbit_params[100]; while (fgets(line, sizeof(line), file)) { - if (parse_body_line(line, name, &mass, &radius, &pos, &parent_index, &r, &g, &b)) { - Vec3 vel = {0.0, 0.0, 0.0}; // Velocity will be calculated later + if (parse_body_line(line, name, &mass, &radius, &pos, &parent_index, &r, &g, &b, + &eccentricity, &semi_major_axis)) { + Vec3 vel = {0.0, 0.0, 0.0}; add_body(sim, name, mass, radius, pos, vel, parent_index, r, g, b); + + orbit_params[sim->body_count - 1].eccentricity = eccentricity; + orbit_params[sim->body_count - 1].semi_major_axis = semi_major_axis; } } @@ -48,8 +64,8 @@ bool load_system_config(SimulationState* sim, const char* filepath) { return false; } - // Calculate initial velocities for circular orbits - calculate_initial_velocities(sim); + // Calculate initial velocities + calculate_velocities(sim, orbit_params); // Calculate SOI radii calculate_soi_radii(sim); @@ -58,9 +74,8 @@ bool load_system_config(SimulationState* sim, const char* filepath) { return true; } -// Calculate circular orbit velocity: v = sqrt(G * M / r) -// Velocity is perpendicular to position vector -void calculate_initial_velocities(SimulationState* sim) { +// Calculate initial velocities using vis-viva equation +void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params) { // First, handle multiple root bodies (binary stars, etc.) // Find all root bodies and calculate barycentric orbits int root_count = 0; @@ -168,8 +183,16 @@ void calculate_initial_velocities(SimulationState* sim) { continue; } - // Calculate circular orbit speed - double speed = sqrt(G * parent->mass / distance); + double e = orbit_params[i].eccentricity; + double a = orbit_params[i].semi_major_axis; + + // Use vis-viva equation: v = sqrt(G*M*(2/r - 1/a)) + double speed = sqrt(G * parent->mass * (2.0 / distance - 1.0 / a)); + + if (e > 0.01) { + printf(" %s: eccentric orbit (e=%.2f, a=%.3e m), speed at r=%.3e m: %.3f km/s\n", + body->name, e, a, distance, speed / 1000.0); + } // Create velocity perpendicular to position vector // If position is mostly in XY plane, make velocity in XY plane