Skip to content

Instantly share code, notes, and snippets.

@Asmilex
Last active December 3, 2023 09:16
Show Gist options
  • Save Asmilex/da1c68138241b0e5b6c9e1952d1df856 to your computer and use it in GitHub Desktop.
Save Asmilex/da1c68138241b0e5b6c9e1952d1df856 to your computer and use it in GitHub Desktop.
Simple simulation for the two body problem implemented in Rust
use std::{fmt::Debug, io, str::FromStr};
// ─── Helper Functions ────────────────────────────────────────────────────────
/// Calculates the euclidean distance between the two bodies.
/// Given our assumptions, we don't need to be so rigorous with the distance formula
fn distance(a: &CelestialBody, b: &CelestialBody) -> f32 {
let x = a.position.0 - b.position.0;
let y = a.position.1 - b.position.1;
(x.powi(2) + y.powi(2)).sqrt()
}
/// Reads a line from stdin and unwraps it to the desired type.
fn read_input<T: FromStr>() -> T
where
<T as FromStr>::Err: Debug,
{
// We won't be doing any error handling
let mut input = String::new();
io::stdin().read_line(&mut input).unwrap();
input.trim().parse::<T>().unwrap()
}
/// Reads a line from stdin and parses the result the desired types.
fn read_pair<T: FromStr, S: FromStr>() -> (T, S)
where
<T as FromStr>::Err: Debug,
<S as FromStr>::Err: Debug,
{
let mut input = String::new();
io::stdin().read_line(&mut input).unwrap();
let mut iter = input.split_whitespace();
let x = iter.next().unwrap().parse::<T>().unwrap();
let y = iter.next().unwrap().parse::<S>().unwrap();
(x, y)
}
// ─────────────────────────────────────────────────────────────────────────────
#[derive(Debug)]
struct CelestialBody {
// We're not fully utilizing Rust's type system here.
// You should encode dimensions with something like [uom](https://swaits.com/beautiful-experience-with-rust-static-type-system/#hello-uom-crate-so-nice-to-meet-you),
// but this approach works for our toy simulation.
velocity: (f32, f32), // m/s
position: (f32, f32),
mass: f32, // kg
}
#[derive(Debug)]
struct CelestialSystem {
planet: CelestialBody,
meteor: CelestialBody,
gravitational_constant: f32,
}
impl CelestialSystem {
/// Uses Euler's method to run a step in the simulation, which calculates the new state of the meteorite.
/// Based on [Celestial Mechanics on a Graphing Calculator — Euler's method](https://www.ams.org/publicoutreach/feature-column/fcarc-kepler2)
///
/// Key assumptions:
/// 1. The planet is stationary at the origin.
fn simulate(&mut self, delta_time: f32) {
// http://www.physics.drexel.edu/~steve/Courses/Physics-431/two-body.pdf
// Assume that x2 = position from planet = (0, 0). Then, given the meteor position, it follows the next formula:
// dv(t)/dt = -G * M * position / |position|^3
// Discretizing it using Euler's method yields
let distance = distance(&self.planet, &self.meteor);
let new_velocity: (f32, f32) = (
self.meteor.velocity.0
+ delta_time
* (-self.gravitational_constant)
* self.planet.mass
* self.meteor.position.0
/ distance.powi(3),
self.meteor.velocity.1
+ delta_time
* (-self.gravitational_constant)
* self.planet.mass
* self.meteor.position.1
/ distance.powi(3),
);
let new_position: (f32, f32) = (
self.meteor.position.0 + new_velocity.0 * delta_time,
self.meteor.position.1 + new_velocity.1 * delta_time,
);
// I might have screwed this one up — please verify it!
self.meteor.position = new_position;
self.meteor.velocity = new_velocity;
}
fn print_meteor_position(&self) {
println!(
"{:.2} {:.2}",
self.meteor.position.0, self.meteor.position.1
);
}
}
#[derive(Debug)]
struct SimulationStep {
delta_time: f32,
remaining: u32,
}
impl SimulationStep {
/// Consume one of the remaining steps and return the delta time.
fn take(&mut self) -> Option<f32> {
match self.remaining {
0 => None,
_ => {
self.remaining -= 1;
Some(self.delta_time)
}
}
}
}
fn main() {
// ─── Read and parse data ─────────────────────────────────────────────
let planet_mass = read_input::<f32>();
let meteor_position = read_pair::<f32, f32>();
let meteor_velocity = read_pair::<f32, f32>();
let mut world = CelestialSystem {
planet: CelestialBody {
velocity: (0.0, 0.0),
position: (0.0, 0.0),
mass: planet_mass,
},
meteor: CelestialBody {
velocity: meteor_velocity,
position: meteor_position,
mass: 1.0, // we don't really provide a mass for the meteor
},
gravitational_constant: 1.0,
};
// We're hardcoding how many of them we're going to read, since it is a fixed problem
// In reality you'd first ask for how many of them will there be — as an example,
// 3
// 1 2
// 0.5 3
// 0.1 2
let mut simulation_steps: Vec<SimulationStep> = vec![];
for _ in 0..2 {
let (delta_time, times) = read_pair::<f32, u32>();
simulation_steps.push(SimulationStep {
delta_time,
remaining: times,
});
}
// ─── Run the simulation ─────────────────────────────────────────────
for mut step in simulation_steps {
// Until you don't have remaining retries for this delta time...
while step.take().is_some() {
world.simulate(step.delta_time);
world.print_meteor_position();
}
}
}
10
0.0 10.0
0.0 0.0
1.0 2
0.5 2
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment