solarsystems 0.0.1

pub mod asteroids;
pub mod belts;
pub mod comets;
pub mod config;
pub mod darkmatter;
pub mod dwarfplanets;
pub mod dynamics;
pub mod exports;
pub mod gravity;
pub mod orbits;
pub mod orchestrator;
pub mod planets;
pub mod rendering;
pub mod satellites;
pub mod stars;

use std::fmt;
use std::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign};

use sciforge::hub::domain::common::constants::{C, G};

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Vec3 {
    pub x: f64,
    pub y: f64,
    pub z: f64,
}

impl Vec3 {
    pub const ZERO: Self = Self {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };

    pub fn new(x: f64, y: f64, z: f64) -> Self {
        Self { x, y, z }
    }

    pub fn magnitude(&self) -> f64 {
        self.magnitude_squared().sqrt()
    }

    pub fn magnitude_squared(&self) -> f64 {
        self.x * self.x + self.y * self.y + self.z * self.z
    }

    pub fn normalize(&self) -> Self {
        let m = self.magnitude();
        if m == 0.0 {
            return Self::ZERO;
        }
        *self / m
    }

    pub fn dot(&self, o: &Self) -> f64 {
        self.x * o.x + self.y * o.y + self.z * o.z
    }

    pub fn cross(&self, o: &Self) -> Self {
        Self {
            x: self.y * o.z - self.z * o.y,
            y: self.z * o.x - self.x * o.z,
            z: self.x * o.y - self.y * o.x,
        }
    }

    pub fn distance(&self, o: &Self) -> f64 {
        (*self - *o).magnitude()
    }

    pub fn angle_between(&self, o: &Self) -> f64 {
        let d = self.dot(o) / (self.magnitude() * o.magnitude());
        d.clamp(-1.0, 1.0).acos()
    }

    pub fn rotate_x(&self, a: f64) -> Self {
        let (s, c) = a.sin_cos();
        Self {
            x: self.x,
            y: self.y * c - self.z * s,
            z: self.y * s + self.z * c,
        }
    }

    pub fn rotate_y(&self, a: f64) -> Self {
        let (s, c) = a.sin_cos();
        Self {
            x: self.x * c + self.z * s,
            y: self.y,
            z: -self.x * s + self.z * c,
        }
    }

    pub fn rotate_z(&self, a: f64) -> Self {
        let (s, c) = a.sin_cos();
        Self {
            x: self.x * c - self.y * s,
            y: self.x * s + self.y * c,
            z: self.z,
        }
    }

    pub fn lerp(&self, o: &Self, t: f64) -> Self {
        *self * (1.0 - t) + *o * t
    }
}

impl Add for Vec3 {
    type Output = Self;
    fn add(self, r: Self) -> Self {
        Self {
            x: self.x + r.x,
            y: self.y + r.y,
            z: self.z + r.z,
        }
    }
}

impl AddAssign for Vec3 {
    fn add_assign(&mut self, r: Self) {
        self.x += r.x;
        self.y += r.y;
        self.z += r.z;
    }
}

impl Sub for Vec3 {
    type Output = Self;
    fn sub(self, r: Self) -> Self {
        Self {
            x: self.x - r.x,
            y: self.y - r.y,
            z: self.z - r.z,
        }
    }
}

impl SubAssign for Vec3 {
    fn sub_assign(&mut self, r: Self) {
        self.x -= r.x;
        self.y -= r.y;
        self.z -= r.z;
    }
}

impl Mul<f64> for Vec3 {
    type Output = Self;
    fn mul(self, s: f64) -> Self {
        Self {
            x: self.x * s,
            y: self.y * s,
            z: self.z * s,
        }
    }
}

impl Mul<Vec3> for f64 {
    type Output = Vec3;
    fn mul(self, v: Vec3) -> Vec3 {
        Vec3 {
            x: self * v.x,
            y: self * v.y,
            z: self * v.z,
        }
    }
}

impl MulAssign<f64> for Vec3 {
    fn mul_assign(&mut self, s: f64) {
        self.x *= s;
        self.y *= s;
        self.z *= s;
    }
}

impl Div<f64> for Vec3 {
    type Output = Self;
    fn div(self, s: f64) -> Self {
        Self {
            x: self.x / s,
            y: self.y / s,
            z: self.z / s,
        }
    }
}

impl Neg for Vec3 {
    type Output = Self;
    fn neg(self) -> Self {
        Self {
            x: -self.x,
            y: -self.y,
            z: -self.z,
        }
    }
}

impl fmt::Display for Vec3 {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "({:.6e}, {:.6e}, {:.6e})", self.x, self.y, self.z)
    }
}

#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum BodyId {
    Sun,
    Mercury,
    Venus,
    Earth,
    Mars,
    Jupiter,
    Saturn,
    Uranus,
    Neptune,
    Moon,
    Phobos,
    Deimos,
    Io,
    Europa,
    Ganymede,
    Callisto,
    Titan,
    Enceladus,
    Titania,
    Oberon,
    Triton,
    // Dwarf planets
    Pluto,
    Eris,
    Haumea,
    Makemake,
    Sedna,
    // Asteroids
    Ceres,
    Vesta,
    Pallas,
    Hygiea,
    Juno,
    // Comets
    Halley,
    Encke,
    HaleBopp,
    ChuryumovGerasimenko,
    Tempel1,
    // Belt particles (procedurally generated)
    BeltAsteroid(u16),
    KuiperObject(u16),
}

impl fmt::Display for BodyId {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "{:?}", self)
    }
}

impl BodyId {
    pub fn all_planets() -> &'static [BodyId] {
        &[
            BodyId::Mercury,
            BodyId::Venus,
            BodyId::Earth,
            BodyId::Mars,
            BodyId::Jupiter,
            BodyId::Saturn,
            BodyId::Uranus,
            BodyId::Neptune,
        ]
    }

    pub fn all_satellites() -> &'static [BodyId] {
        &[
            BodyId::Moon,
            BodyId::Phobos,
            BodyId::Deimos,
            BodyId::Io,
            BodyId::Europa,
            BodyId::Ganymede,
            BodyId::Callisto,
            BodyId::Titan,
            BodyId::Enceladus,
            BodyId::Titania,
            BodyId::Oberon,
            BodyId::Triton,
        ]
    }

    pub fn all_dwarf_planets() -> &'static [BodyId] {
        &[
            BodyId::Pluto,
            BodyId::Eris,
            BodyId::Haumea,
            BodyId::Makemake,
            BodyId::Sedna,
        ]
    }

    pub fn all_asteroids() -> &'static [BodyId] {
        &[
            BodyId::Ceres,
            BodyId::Vesta,
            BodyId::Pallas,
            BodyId::Hygiea,
            BodyId::Juno,
        ]
    }

    pub fn all_comets() -> &'static [BodyId] {
        &[
            BodyId::Halley,
            BodyId::Encke,
            BodyId::HaleBopp,
            BodyId::ChuryumovGerasimenko,
            BodyId::Tempel1,
        ]
    }

    pub fn all() -> &'static [BodyId] {
        &[
            BodyId::Sun,
            BodyId::Mercury,
            BodyId::Venus,
            BodyId::Earth,
            BodyId::Mars,
            BodyId::Jupiter,
            BodyId::Saturn,
            BodyId::Uranus,
            BodyId::Neptune,
            BodyId::Moon,
            BodyId::Phobos,
            BodyId::Deimos,
            BodyId::Io,
            BodyId::Europa,
            BodyId::Ganymede,
            BodyId::Callisto,
            BodyId::Titan,
            BodyId::Enceladus,
            BodyId::Titania,
            BodyId::Oberon,
            BodyId::Triton,
            BodyId::Pluto,
            BodyId::Eris,
            BodyId::Haumea,
            BodyId::Makemake,
            BodyId::Sedna,
            BodyId::Ceres,
            BodyId::Vesta,
            BodyId::Pallas,
            BodyId::Hygiea,
            BodyId::Juno,
            BodyId::Halley,
            BodyId::Encke,
            BodyId::HaleBopp,
            BodyId::ChuryumovGerasimenko,
            BodyId::Tempel1,
        ]
    }

    pub fn parent(&self) -> Option<BodyId> {
        match self {
            BodyId::Mercury
            | BodyId::Venus
            | BodyId::Earth
            | BodyId::Mars
            | BodyId::Jupiter
            | BodyId::Saturn
            | BodyId::Uranus
            | BodyId::Neptune => Some(BodyId::Sun),
            BodyId::Moon => Some(BodyId::Earth),
            BodyId::Phobos | BodyId::Deimos => Some(BodyId::Mars),
            BodyId::Io | BodyId::Europa | BodyId::Ganymede | BodyId::Callisto => {
                Some(BodyId::Jupiter)
            }
            BodyId::Titan | BodyId::Enceladus => Some(BodyId::Saturn),
            BodyId::Titania | BodyId::Oberon => Some(BodyId::Uranus),
            BodyId::Triton => Some(BodyId::Neptune),
            BodyId::Pluto | BodyId::Eris | BodyId::Haumea | BodyId::Makemake | BodyId::Sedna => {
                Some(BodyId::Sun)
            }
            BodyId::Ceres | BodyId::Vesta | BodyId::Pallas | BodyId::Hygiea | BodyId::Juno => {
                Some(BodyId::Sun)
            }
            BodyId::Halley
            | BodyId::Encke
            | BodyId::HaleBopp
            | BodyId::ChuryumovGerasimenko
            | BodyId::Tempel1 => Some(BodyId::Sun),
            BodyId::BeltAsteroid(_) | BodyId::KuiperObject(_) => Some(BodyId::Sun),
            BodyId::Sun => None,
        }
    }

    pub fn is_planet(&self) -> bool {
        Self::all_planets().contains(self)
    }

    pub fn is_satellite(&self) -> bool {
        Self::all_satellites().contains(self)
    }

    pub fn is_dwarf_planet(&self) -> bool {
        Self::all_dwarf_planets().contains(self)
    }

    pub fn is_asteroid(&self) -> bool {
        matches!(self, BodyId::BeltAsteroid(_)) || Self::all_asteroids().contains(self)
    }

    pub fn is_comet(&self) -> bool {
        Self::all_comets().contains(self)
    }

    pub fn is_belt_asteroid(&self) -> bool {
        matches!(self, BodyId::BeltAsteroid(_))
    }

    pub fn is_kuiper_object(&self) -> bool {
        matches!(self, BodyId::KuiperObject(_))
    }

    pub fn is_belt_particle(&self) -> bool {
        matches!(self, BodyId::BeltAsteroid(_) | BodyId::KuiperObject(_))
    }
}

#[derive(Clone, Debug)]
pub struct CelestialBody {
    pub id: BodyId,
    pub name: &'static str,
    pub mass: f64,
    pub radius: f64,
    pub position: Vec3,
    pub velocity: Vec3,
    pub acceleration: Vec3,
    pub j2: f64,
}

impl CelestialBody {
    pub fn gravitational_parameter(&self) -> f64 {
        G * self.mass
    }

    pub fn escape_velocity(&self) -> f64 {
        (2.0 * G * self.mass / self.radius).sqrt()
    }

    pub fn surface_gravity(&self) -> f64 {
        G * self.mass / (self.radius * self.radius)
    }

    pub fn mean_density(&self) -> f64 {
        self.mass / ((4.0 / 3.0) * std::f64::consts::PI * self.radius.powi(3))
    }

    pub fn hill_sphere(&self, parent_mass: f64, distance: f64) -> f64 {
        distance * (self.mass / (3.0 * parent_mass)).powf(1.0 / 3.0)
    }

    pub fn roche_limit(&self, satellite_density: f64) -> f64 {
        self.radius * 2.456 * (self.mean_density() / satellite_density).powf(1.0 / 3.0)
    }

    pub fn sphere_of_influence(&self, parent_mass: f64, distance: f64) -> f64 {
        distance * (self.mass / parent_mass).powf(2.0 / 5.0)
    }

    pub fn kinetic_energy(&self) -> f64 {
        0.5 * self.mass * self.velocity.magnitude_squared()
    }

    pub fn potential_energy(&self, other: &CelestialBody) -> f64 {
        let r = self.position.distance(&other.position);
        -G * self.mass * other.mass / r
    }

    pub fn angular_momentum(&self) -> Vec3 {
        self.mass * self.position.cross(&self.velocity)
    }

    pub fn specific_angular_momentum(&self) -> Vec3 {
        self.position.cross(&self.velocity)
    }

    pub fn specific_orbital_energy(&self, central_mass: f64) -> f64 {
        let v2 = self.velocity.magnitude_squared();
        let r = self.position.magnitude();
        0.5 * v2 - G * central_mass / r
    }

    pub fn orbital_period(&self, central_mass: f64) -> f64 {
        let a = self.position.magnitude();
        let mu = G * (central_mass + self.mass);
        2.0 * std::f64::consts::PI * (a.powi(3) / mu).sqrt()
    }

    pub fn vis_viva_speed(&self, central_mass: f64, semi_major_axis: f64) -> f64 {
        let r = self.position.magnitude();
        let mu = G * (central_mass + self.mass);
        (mu * (2.0 / r - 1.0 / semi_major_axis)).abs().sqrt()
    }

    pub fn schwarzschild_radius(&self) -> f64 {
        2.0 * G * self.mass / (C * C)
    }
}