solarsystems 0.0.1

use std::f64::consts::PI;

use dwarfplanetsfactory::types::cold_classical::ColdClassical;
use dwarfplanetsfactory::types::plutino::Plutino;
use dwarfplanetsfactory::types::scattered_disk::ScatteredDisk;

use crate::config::parameters::MU_SUN;
use crate::orbits::keplerian::OrbitalElements;
use crate::{BodyId, CelestialBody, Vec3};

const AU: f64 = 1.495978707e11;

/// Simple deterministic PRNG (xorshift64) for reproducible belt generation.
struct Rng(u64);

impl Rng {
    fn next_u64(&mut self) -> u64 {
        self.0 ^= self.0 << 13;
        self.0 ^= self.0 >> 7;
        self.0 ^= self.0 << 17;
        self.0
    }

    fn f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 / ((1u64 << 53) as f64)
    }

    fn range(&mut self, lo: f64, hi: f64) -> f64 {
        lo + self.f64() * (hi - lo)
    }
}

/// Generate `count` Kuiper belt objects using `dwarfplanetsfactory` types.
///
/// Three dynamical populations (matching observations):
/// - **Cold Classical** (40%): a = 42–47 AU, e = 0–0.1, i = 0–5°
/// - **Plutinos / 3:2 resonance** (20%): a ≈ 39.4 AU, e = 0.1–0.25, i = 0–15°
/// - **Scattered disk** (40%): a = 35–100 AU, e = 0.2–0.6, i = 0–30°
pub fn generate_kuiper_belt(count: u16) -> Vec<CelestialBody> {
    let mut rng = Rng(0xCE1E_5714_4B01_BE17);
    let mut bodies = Vec::with_capacity(count as usize);

    let n_classical = (count as f64 * 0.4) as u16;
    let n_plutino = (count as f64 * 0.2) as u16;
    let n_scattered = count - n_classical - n_plutino;

    let mut idx: u16 = 0;

    // ── Cold Classical KBOs ──
    for _ in 0..n_classical {
        let radius = 10.0_f64.powf(rng.range(5.3, 6.1)); // 200 km – 1260 km
        let a_au = rng.range(42.0, 47.0);
        let ecc = rng.range(0.0, 0.1);
        let inc_deg = rng.range(0.0, 5.0);

        let kbo = ColdClassical::new(radius, a_au)
            .with_eccentricity(ecc)
            .with_inclination(inc_deg)
            .with_density(rng.range(800.0, 1500.0));

        let elements = make_elements(
            kbo.semi_major_axis,
            kbo.eccentricity,
            kbo.inclination.to_radians(), // factory stores degrees
            &mut rng,
        );
        let (position, velocity) = elements.to_state_vectors(MU_SUN);

        bodies.push(CelestialBody {
            id: BodyId::KuiperObject(idx),
            name: "Cold Classical KBO",
            mass: kbo.mass(),
            radius: kbo.radius,
            position,
            velocity,
            acceleration: Vec3::ZERO,
            j2: 0.0,
        });
        idx += 1;
    }

    // ── Plutinos (3:2 Neptune resonance) ──
    for _ in 0..n_plutino {
        let radius = 10.0_f64.powf(rng.range(5.3, 6.1));
        let ecc = rng.range(0.1, 0.25);
        let inc_deg = rng.range(0.0, 15.0);

        let kbo = Plutino::new(radius, ecc)
            .with_inclination(inc_deg)
            .with_density(rng.range(1200.0, 2200.0));

        let elements = make_elements(
            kbo.semi_major_axis, // fixed ~39.4 AU
            kbo.eccentricity,
            kbo.inclination.to_radians(),
            &mut rng,
        );
        let (position, velocity) = elements.to_state_vectors(MU_SUN);

        bodies.push(CelestialBody {
            id: BodyId::KuiperObject(idx),
            name: "Plutino",
            mass: kbo.mass(),
            radius: kbo.radius,
            position,
            velocity,
            acceleration: Vec3::ZERO,
            j2: 0.0,
        });
        idx += 1;
    }

    // ── Scattered Disk Objects ──
    for _ in 0..n_scattered {
        let radius = 10.0_f64.powf(rng.range(5.3, 6.1));
        let a_au = rng.range(35.0, 100.0);
        let ecc = rng.range(0.2, 0.6);
        let inc_deg = rng.range(0.0, 30.0);

        let kbo = ScatteredDisk::new(radius, a_au, ecc)
            .with_inclination(inc_deg)
            .with_density(rng.range(1000.0, 2000.0));

        let elements = make_elements(
            kbo.semi_major_axis,
            kbo.eccentricity,
            kbo.inclination.to_radians(),
            &mut rng,
        );
        let (position, velocity) = elements.to_state_vectors(MU_SUN);

        bodies.push(CelestialBody {
            id: BodyId::KuiperObject(idx),
            name: "Scattered Disk Object",
            mass: kbo.mass(),
            radius: kbo.radius,
            position,
            velocity,
            acceleration: Vec3::ZERO,
            j2: 0.0,
        });
        idx += 1;
    }

    bodies
}

fn make_elements(a_au: f64, ecc: f64, inc_rad: f64, rng: &mut Rng) -> OrbitalElements {
    OrbitalElements {
        semi_major_axis: a_au * AU,
        eccentricity: ecc,
        inclination: inc_rad,
        longitude_ascending_node: rng.range(0.0, 2.0 * PI),
        argument_periapsis: rng.range(0.0, 2.0 * PI),
        true_anomaly: rng.range(0.0, 2.0 * PI),
    }
}