Simple simulation for the two body problem implemented in Rust

main.rs

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();
        }
    }
}

parameters

10
0.0 10.0
0.0 0.0
1.0 2
0.5 2