use neptunes::geology::earthquakes::{
    Icequake, aftershockrate, annualfrequency, distancefromsptime, icepwavevelocity,
    iceswavevelocity, refractionangle,
};
use neptunes::geology::erosion::{
    ErosionParams, chemicalweatheringrate, sublimationerosionrate, thermalstresserosion,
    winderosionthresholdvelocity,
};
use neptunes::geology::mountains::{
    ICEMANTLEDENSITY, ROCKYCOREDENSITY, interioricepeak, rockycorepeak,
};
use neptunes::geology::plate_tectonics::{
    ICESHELLTHICKNESS, INTERIORVISCOSITY, equatorialcell, iceshellflexuralrigidity,
    iceshellrootdepth, isostaticbalance, polarcell, surfaceheatflow, temperatureatdepth,
};
use neptunes::geology::volcanism::{CryoMagma, cryovolcanicexplosivityindex};

fn main() {
    let eq = Icequake::frommoment(45.0, 120.0, 50.0, 3.9e18);
    let eqmag = Icequake::frommagnitude(45.0, 120.0, 50.0, 5.5);
    let mw = eq.momentmagnitude();
    let energy = eq.energyjoules();
    let vp = icepwavevelocity(920.0, 8.9e9, 3.5e9);
    let vs = iceswavevelocity(920.0, 3.5e9);
    let tp = eq.pwavearrivaltime(200.0, vp);
    let ts = eq.swavearrivaltime(200.0, vs);
    let pga = eq.pgaatdistance(50.0);
    let dist = distancefromsptime(tp, ts, vp, vs);
    let angle = refractionangle(3000.0, 5000.0, 0.5);

    println!(
        "Icequake: Mw={:.1}  energy={:.2e}J  Mw(5.5)={:.1}",
        mw,
        energy,
        eqmag.momentmagnitude()
    );
    println!(
        "Vp={:.0}m/s  Vs={:.0}m/s  tp@200km={:.1}s  ts={:.1}s  PGA@50km={:.4}",
        vp, vs, tp, ts, pga
    );
    println!("distfromSP={:.1}km  refraction={:.3}rad", dist, angle);

    let aftershock = aftershockrate(5.5, 1.0);
    let freqm4 = annualfrequency(4.0);
    println!(
        "aftershocks@1day={:.1}  freq(M4)={:.1}/yr",
        aftershock, freqm4
    );

    let waterammonia = CryoMagma::waterammonia();
    let methanerich = CryoMagma::methanerich();
    let silicatecore = CryoMagma::silicatecore();
    println!(
        "WaterNH3: Mg#={:.2}  alkali={:.2}  Tliq={:.0}°C  η={:.2e}Pa·s  solidΔ={:.1}°C  DI={:.1}",
        waterammonia.mgnumber(),
        waterammonia.totalalkali(),
        waterammonia.liquidustempc(),
        waterammonia.viscositypas(),
        waterammonia.solidusdepressionc(),
        waterammonia.differentiationidx()
    );
    println!(
        "CH4rich:  Mg#={:.2}  alkali={:.2}  Tliq={:.0}°C  η={:.2e}Pa·s",
        methanerich.mgnumber(),
        methanerich.totalalkali(),
        methanerich.liquidustempc(),
        methanerich.viscositypas()
    );
    println!(
        "Silicate: Mg#={:.2}  alkali={:.2}  η={:.2e}Pa·s",
        silicatecore.mgnumber(),
        silicatecore.totalalkali(),
        silicatecore.viscositypas()
    );

    let veismall = cryovolcanicexplosivityindex(0.01);
    let veimed = cryovolcanicexplosivityindex(1.0);
    let veilarge = cryovolcanicexplosivityindex(100.0);
    println!(
        "CryoVEI: 0.01km³={}  1km³={}  100km³={}",
        veismall, veimed, veilarge
    );

    let eqcell = equatorialcell();
    let pcell = polarcell();
    let veq = eqcell.velocitycmperyear(0.0, 90.0);
    let vpol = pcell.velocitycmperyear(80.0, 0.0);
    let vatpoint = eqcell.velocityatpoint(30.0, 60.0);

    println!(
        "EqCell: v@(0,90)={:.3}cm/yr  v@(30,60)={:.4}m/s",
        veq, vatpoint
    );
    println!("PolarCell: v@(80,0)={:.3}cm/yr", vpol);

    let root = iceshellrootdepth(4000.0);
    let flex = iceshellflexuralrigidity(100000.0, 9e9, 0.33);
    let hf = surfaceheatflow(2.0, 0.005);
    let t500km = temperatureatdepth(50.0, 500000.0, 0.005);
    let iso = isostaticbalance(50000.0);
    println!(
        "root(4km)={:.0}m  flex={:.2e}N·m  heatflow={:.4}W/m²",
        root, flex, hf
    );
    println!("T@500km={:.0}K  isostatic={:.0}m", t500km, iso);
    println!(
        "constants: ρice={} ρrock={} shell={}m  η={:.0e}Pa·s",
        ICEMANTLEDENSITY, ROCKYCOREDENSITY, ICESHELLTHICKNESS, INTERIORVISCOSITY
    );

    let ipeak = interioricepeak();
    let rpeak = rockycorepeak();
    println!(
        "IcePeak: h={:.0}m  root={:.0}m  total={:.0}m  iso={:.0}m",
        ipeak.height(),
        ipeak.rootdepthm(),
        ipeak.totalthicknessm(),
        ipeak.isostaticelevation()
    );
    println!(
        "RockPeak: h={:.0}m  root={:.0}m  total={:.0}m  iso={:.0}m",
        rpeak.height(),
        rpeak.rootdepthm(),
        rpeak.totalthicknessm(),
        rpeak.isostaticelevation()
    );

    let icy = ErosionParams {
        windspeedms: 200.0,
        atmosphericdensity: 0.45,
        particlediameterm: 0.001,
        particledensity: 920.0,
    };
    let ewind = icy.winderosionratemmyr();
    let esed = icy.sedimentyieldtkm2yr();
    println!(
        "IcyErosion: wind={:.4}mm/yr  sediment={:.4}t/km²/yr",
        ewind, esed
    );

    let chem = chemicalweatheringrate(72.0, 920.0);
    let thermal = thermalstresserosion(30.0, 5e-5);
    let sublim = sublimationerosionrate(72.0, 100.0);
    let windthresh = winderosionthresholdvelocity(0.001, 920.0);
    println!(
        "chem={:.6}  thermal={:.6}mm/yr  sublim={:.6}mm/yr  windthresh={:.1}m/s",
        chem, thermal, sublim, windthresh
    );
}
