jupiters 0.0.3

use jupiters::physics::collisions::{Impactor, asteroidimpactor, cometimpact};
use jupiters::physics::orbit::JupiterOrbit;
use jupiters::physics::rotation::{JupiterRotation, SIDEREALDAYS};
use jupiters::physics::tides::{
    TidalForce, iotosolartideratio, neaptideamplitude, springtideamplitude,
};
use jupiters::satellites::artificial::{ArtificialSatellite, Constellation};
use jupiters::satellites::io::IoState;

fn main() {
    let orbit = JupiterOrbit::new();
    let period = orbit.orbitalperioddays();
    let vperi = orbit.velocityatdistance(orbit.perihelionm());
    let vaph = orbit.velocityatdistance(orbit.aphelionm());
    let vmean = orbit.mean_orbital_velocity();
    let energy = orbit.specific_orbital_energy();
    let angular = orbit.specific_angular_momentum();
    let vesc = JupiterOrbit::escape_velocity_at_surface();
    let fgrav = orbit.gravitational_force_sun();
    println!(
        "Orbit: P={:.1}d  q={:.3e}m  Q={:.3e}m",
        period,
        orbit.perihelionm(),
        orbit.aphelionm()
    );
    println!(
        "v(peri)={:.0}m/s  v(aph)={:.0}m/s  v(mean)={:.0}m/s  v(esc)={:.0}m/s",
        vperi, vaph, vmean, vesc
    );
    println!(
        "E={:.3e}J/kg  L={:.3e}m²/s  F(sun)={:.3e}N",
        energy, angular, fgrav
    );

    let rot = JupiterRotation::new();
    let veq = rot.surfacevelocityatlatitude(0.0);
    let v45 = rot.surfacevelocityatlatitude(45.0);
    let moi = rot.momentofinertia();
    let ke = rot.rotationalkineticenergy();
    let angmom = rot.angular_momentum();
    let prec = rot.precession_rate_rad_per_year();
    println!(
        "Rotation: day={:.0}s  v(eq)={:.0}m/s  v(45°)={:.0}m/s",
        SIDEREALDAYS, veq, v45
    );
    println!(
        "I={:.3e}kg·m²  KE={:.3e}J  L={:.3e}kg·m²/s  prec={:.2e}rad/yr",
        moi, ke, angmom, prec
    );

    let latitudes = [0.0, 15.0, 30.0, 45.0, 60.0, 75.0];
    for &lat in &latitudes {
        let daylength = rot.day_length_variation_due_to_tilt(172, lat);
        let centripetal = rot.centripetalaccelerationatlatitude(lat);
        let coriolis = rot.coriolisparameter(lat);
        println!(
            "lat={:.0}°: daylight={:.1}h  centripetal={:.3}m/s²  f={:.5e}",
            lat, daylength, centripetal, coriolis
        );
    }

    let iotide = TidalForce::fromio();
    let suntide = TidalForce::fromsun();
    let ioheight = iotide.tidalbulgeheight();
    let sunheight = suntide.tidalbulgeheight();
    let ioaccel = iotide.tidalacceleration();
    let sunaccel = suntide.tidalacceleration();
    let pot0 = iotide.tidalpotential(0.0);
    let pot90 = iotide.tidalpotential(std::f64::consts::FRAC_PI_2);
    let spring = springtideamplitude();
    let neap = neaptideamplitude();
    let ratio = iotosolartideratio();
    println!(
        "Io tide: accel={:.3e}m/s²  bulge={:.3e}m  F={:.3e}N",
        ioaccel,
        ioheight,
        iotide.gravitationalattraction()
    );
    println!(
        "Sun tide: accel={:.3e}m/s²  bulge={:.3e}m  F={:.3e}N",
        sunaccel,
        sunheight,
        suntide.gravitationalattraction()
    );
    println!(
        "Io/Sun ratio={:.1}  spring={:.3e}m  neap={:.3e}m",
        ratio, spring, neap
    );
    println!("Io potential: θ=0→{:.3e}  θ=π/2→{:.3e}", pot0, pot90);

    let mut io = IoState::simulated();
    let dt = 3600.0;
    let mut mindist = f64::MAX;
    let mut maxdist = 0.0_f64;
    for _ in 0..200 {
        io.step(dt);
        let d = io.distance();
        if d < mindist {
            mindist = d;
        }
        if d > maxdist {
            maxdist = d;
        }
    }
    let (x, y, z) = io.position();
    let r = (x * x + y * y + z * z).sqrt();
    println!(
        "Io after 200h: r={:.0}m  pos=({:.0},{:.0},{:.0})  min={:.0}  max={:.0}",
        r, x, y, z, mindist, maxdist
    );

    let sat = ArtificialSatellite::loworbit("Juno", 1500.0, 200000.0);
    let speriod = sat.orbitalperiods();
    let svel = sat.orbitalvelocityms();
    let gsurf = sat.gravityatsurface();
    println!(
        "Juno: alt={:.0}m  P={:.0}s  v={:.0}m/s  g={:.2}m/s²",
        sat.altitudem, speriod, svel, gsurf
    );

    let mut constellation = Constellation::new("Jovian Network");
    for i in 0..4 {
        let alt = 500_000.0 + (i as f64) * 200_000.0;
        constellation.add(ArtificialSatellite::loworbit(
            &format!("Sat{}", i + 1),
            100.0,
            alt,
        ));
    }
    constellation.stepall(3600.0);
    let positions = constellation.positions();
    for (i, &(x, y, z)) in positions.iter().enumerate() {
        let r = (x * x + y * y + z * z).sqrt();
        println!("  Sat{}: r={:.0}m", i + 1, r);
    }

    let comet = cometimpact();
    let ast = asteroidimpactor();
    let custom = Impactor::asteroid(2000.0, 25.0);
    println!(
        "Comet(5km): KE={:.2e}MT  v_impact={:.0}m/s  crater(ρ=1000)={:.0}m",
        comet.kineticenergymt(),
        comet.impactvelocity(),
        comet.craterdiameterm(1000.0)
    );
    println!(
        "Asteroid(100m): KE={:.2e}MT  fireball={:.0}m  ejecta={:.2e}m³",
        ast.kineticenergymt(),
        ast.fireballradiusm(),
        ast.ejectavolumem3(1000.0)
    );
    println!(
        "Custom(2km,25km/s): KE={:.2e}MT  crater={:.0}m",
        custom.kineticenergymt(),
        custom.craterdiameterm(1000.0)
    );
}