Initial Commit

.gitignore

@@ -0,0 +1 @@
+/target

Cargo.toml

@@ -0,0 +1,13 @@
+[package]
+name = "orbital_simulator"
+version = "0.1.0"
+edition = "2021"
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]
+bevy = {version = "0.13.2", default-features = false, features = ["bevy_winit", "bevy_render", "bevy_core_pipeline", "bevy_scene", "bevy_asset", "bevy_gizmos", "bevy_sprite"]}
+bevy_gizmos = "0.13.2"
+nalgebra = "0.32.5"
+ode_solvers = "0.4.0"
+winit = "0.29.15"

src/main.rs

@@ -0,0 +1,32 @@
+use bevy::{
+    ecs::system::Commands,
+    render::camera::ScalingMode,
+    app::{App, Update, Startup},
+    core_pipeline::core_2d::Camera2dBundle,
+    DefaultPlugins,
+};
+
+use std::time::Duration;
+
+mod render;
+mod planets;
+
+fn setup(mut commands: Commands) {
+    let mut camera_bundle = Camera2dBundle::default();
+    camera_bundle.projection.scaling_mode = ScalingMode::FixedVertical(800000000.0);
+    commands.spawn(camera_bundle);
+}
+
+fn main() {
+    App::new()
+        .add_plugins(DefaultPlugins)
+        .insert_resource(render::DrawPathTimer::new(Duration::from_millis(10)))
+        .add_systems(Startup, (setup, planets::add_planets))
+        .add_systems(Update, (
+                planets::gravity::update_gravity,
+                render::add_line,
+                render::draw_planet_positions,
+                render::draw_path,
+            ))
+        .run();
+}

src/planets.rs

@@ -0,0 +1,118 @@
+use bevy::{
+    ecs::{
+        component::Component,
+        system::{Commands, Query, Res},
+        query::With,
+    },
+    time::Time,
+    math::{DVec2, Vec2},
+};
+
+pub(crate) mod gravity;
+
+#[derive(Component)]
+pub(crate) struct Planet;
+
+#[derive(Component, Clone, Copy, Debug)]
+pub(crate) struct Position(pub(crate) DVec2);
+
+impl From<Position> for Vec2 {
+    fn from(position: Position) -> Self {
+        position.0.as_vec2()
+    }
+}
+
+#[derive(Component, Clone, Copy, Default, Debug)]
+pub(crate) struct Velocity(DVec2);
+
+#[derive(Component, Default)]
+pub(crate) struct Acceleration(DVec2);
+
+#[derive(Component)]
+pub(crate) struct Mass(f64);
+
+#[derive(Component)]
+pub(crate) struct Radius(pub(crate) f64);
+
+#[derive(Component, Debug)]
+pub(crate) struct Line {
+    pub(crate) start: Position, 
+    pub(crate) end: Option<Position>,
+}
+
+#[derive(Component, Default, Debug)]
+pub(crate) struct Path {
+    pub(crate) lines: Vec<Line>,
+}
+
+impl<'a> IntoIterator for &'a Path {
+    type Item = Vec2;
+    type IntoIter = PathIterator<'a>;
+    fn into_iter(self) -> Self::IntoIter {
+        PathIterator::new(&self.lines)
+    }
+}
+
+pub(crate) struct PathIterator<'a> {
+    content: &'a Vec<Line>,
+    index: usize,
+    end: usize,
+}
+
+impl<'a> PathIterator<'a> {
+    fn new(content: &'a Vec<Line>) -> Self {
+        Self {content, index: 0, end: content.len()}
+    }
+}
+
+impl Iterator for PathIterator<'_> {
+    type Item = Vec2;
+
+    fn next(&mut self) -> Option<Self::Item> {
+        if self.index < self.end {
+            let item = Some(self.content[self.index].start.0.as_vec2());
+            self.index += 1;
+            item
+        } else {
+            None
+        }
+    }
+}
+
+pub(super) fn add_planets(mut commands: Commands) {
+    /*
+    let moon = (
+            Planet,
+            Position(DVec2::new(384399000.0, 0.0)),
+            Velocity(DVec2::new(0.0, 1022.0)),
+            Acceleration::default(),
+            Mass(7.342e22),
+            Radius(1.7374e6),
+            Path::default(),
+    );
+    */
+
+    let earth = (
+            Planet,
+            Position(DVec2::new(0.0, 0.0)),
+            Velocity(DVec2::new(100.0, -800.5641)),
+            Acceleration::default(),
+            Mass(5.9722e24),
+            Radius(6.371e6),
+            Path::default(),
+    );
+
+    let custom_moon = (
+            Planet,
+            Position(DVec2::new(384399000.0, 0.0)),
+            Velocity(DVec2::new(-1000.0, 400.0)),
+            Acceleration::default(),
+            Mass(7.342e24),
+            Radius(1.7374e6),
+            Path::default(),
+    );
+
+    commands.spawn(earth);
+    commands.spawn(custom_moon);
+
+}

src/planets/gravity.rs

@@ -0,0 +1,114 @@
+use bevy::ecs::system::Query;
+use bevy::ecs::query::With;
+use bevy::math::DVec2;
+
+use crate::planets::{Planet, Position, Velocity, Mass};
+
+use ode_solvers::Rk4;
+
+use nalgebra::{Vector4, Vector2};
+
+use std::error::Error;
+
+const G: f64 = 6.674e-11;
+
+struct GravityODE {
+    alpha: f64,
+    r0: Vector2<f64>,
+    v0: Vector2<f64>,
+}
+
+impl GravityODE {
+    fn new<P, V>(m_1: f64, m_2: f64, p_1: P, p_2: P, v_1: V, v_2: V) -> Self
+    where
+        P: Into<Vector2<f64>>,
+        V: Into<Vector2<f64>>,
+    {
+        let alpha = (m_1 + m_2) * G;
+        let r0: Vector2<f64> = p_2.into() - p_1.into();
+        let v0: Vector2<f64> = v_2.into() - v_1.into();
+
+        GravityODE { alpha, r0, v0 }
+    }
+
+    fn solve(&self, dt: f64, dy: f64) -> Result<(Vector2<f64>, Vector2<f64>), Box<dyn Error>> {
+        let y0 = State::new(self.r0[0], self.r0[1], self.v0[0], self.v0[1]);
+
+        let mut stepper = Rk4::new(self, 0.0, y0, dt, dy);
+
+        stepper.integrate()?;
+        let result = stepper.results().get().1.last().unwrap();
+        let (x_f, y_f) = (result[0], result[1]);
+        let (v_x_f, v_y_f) = (result[2], result[3]);
+
+        Ok((Vector2::new(x_f, y_f), Vector2::new(v_x_f, v_y_f)))
+    } 
+}
+
+type State = Vector4<f64>;
+type Time = f64;
+
+impl ode_solvers::System<Time, Vector4<f64>> for &GravityODE {
+    fn system(&self, t: Time, y: &State, dy: &mut State) {
+        let r = (y[0] * y[0] + y[1] * y[1]).sqrt();
+
+        dy[0] = y[2];
+        dy[1] = y[3];
+        dy[2] = - (self.alpha * y[0]) / r.powi(3);
+        dy[3] = - (self.alpha * y[1]) / r.powi(3);
+    }
+}
+
+impl From<Position> for Vector2<f64> {
+    fn from(position: Position) -> Self {
+        Vector2::new(position.0.x, position.0.y)
+    }
+}
+
+impl From<Velocity> for Vector2<f64> {
+    fn from(position: Velocity) -> Self {
+        Vector2::new(position.0.x, position.0.y)
+    }
+}
+
+impl From<Vector2<f64>> for Position {
+    fn from(v: Vector2<f64>) -> Self {
+        Self(DVec2::new(v[0], v[1]))
+    }
+}
+
+impl From<Vector2<f64>> for Velocity {
+    fn from(v: Vector2<f64>) -> Self {
+        Self(DVec2::new(v[0], v[1]))
+    }
+}
+
+
+pub(crate) fn update_gravity( mut query: Query<(&mut Position, &mut Velocity, &Mass), With<Planet>>,
+) {
+    let mut planets = query.iter_mut();
+
+    let Some((mut p_1, mut v_1, m_1)) = planets.next() else { unreachable!() };
+    let Some((mut p_2, mut v_2, m_2)) = planets.next() else { unreachable!() };
+
+    let system = GravityODE::new(m_1.0, m_2.0, *p_1, *p_2, *v_1, *v_2);
+    let (r, v): (Vector2<f64>, Vector2<f64>) = system.solve(600.0, 60.0).expect("Solving equation");
+
+    let r_1: Vector2<f64> = (*p_1).into();
+    let r_2: Vector2<f64> = (*p_2).into();
+
+    let r_cm: Vector2<f64> = (m_1.0 * r_1 + m_2.0 * r_2) / (m_1.0 + m_2.0);
+
+    let new_p_1 = r_cm - (m_2.0 * r) / (m_1.0 + m_2.0);
+    let new_p_2 = r_cm + (m_1.0 * r) / (m_1.0 + m_2.0);
+
+    let new_v_1 = - (m_2.0 * v) / (m_1.0 + m_2.0);
+    let new_v_2 = (m_1.0 * v) / (m_1.0 + m_2.0);
+
+    *p_1 = new_p_1.into();
+    *p_2 = new_p_2.into();
+
+    *v_1 = new_v_1.into();
+    *v_2 = new_v_2.into();
+
+}

src/render.rs

@@ -0,0 +1,65 @@
+use bevy::{
+    time::{Time, Timer, TimerMode},
+    render::color::Color,
+    gizmos::gizmos::Gizmos,
+    ecs::{
+        system::{Resource, Query, Res, ResMut},
+        query::With,
+    },
+    math::Vec2,
+};
+
+use std::time::Duration;
+
+use crate::planets::{self, Planet, Position, Path, Radius};
+
+#[derive(Resource)]
+pub(crate) struct DrawPathTimer {
+    pub(super) timer: Timer,
+}
+
+impl DrawPathTimer {
+    pub(crate) fn new(time: Duration) -> Self {
+        DrawPathTimer {
+            timer: Timer::new(time, TimerMode::Repeating),
+        }
+    }
+}
+
+pub(crate) fn add_line(
+    time: Res<Time>,
+    mut query: Query<(&planets::Position, &mut planets::Path), With<Planet>>,
+    mut timer: ResMut<DrawPathTimer>,
+) {
+    if timer.timer.tick(time.delta()).just_finished() {
+        for (position, mut path) in query.iter_mut() {
+            if path.lines.len() > 0 {
+                path.lines.last_mut().unwrap().end = Some(*position)
+            }
+
+            path.lines.push(planets::Line {
+                    start: *position,
+                    end: None,
+            });
+        }
+    }
+}
+
+pub(crate) fn draw_path(
+    query: Query<&Path, With<Planet>>,
+    mut gizmos: Gizmos,
+) {
+    for path in query.into_iter() {
+            gizmos.linestrip_2d(path, Color::BLUE);
+    }
+}
+
+pub(crate) fn draw_planet_positions(
+    mut gizmos: Gizmos,
+    query: Query<(&Position, &Radius), With<Planet>>,
+) {
+    gizmos.circle_2d(Vec2::new(0.0, 0.0), 400000.0, Color::BLUE);
+    for (position, radius) in &query {
+        gizmos.circle_2d(position.0.as_vec2(), radius.0 as f32, Color::RED);
+    }
+}