use neptunes::atmosphere::climate::{
    ClimateState, NEPTUNEALBEDO, NEPTUNEEFFECTIVETEMPK, NEPTUNEEQUILIBRIUMTEMPK,
    co2doublingforcing, planck_spectral_radiance, solarfluxatneptune, solarzenith,
};
use neptunes::atmosphere::heatwaves::{
    ThermalAnomaly, dewpointmethanek, saturationvaporpressuremethanepa,
};
use neptunes::atmosphere::layers::{
    AtmosphericLayer, LAPSERATETROPOSPHERE, ONEBARPRESSURE, ONEBARTEMPK, atmospheremolarmass,
    barometricpressure, effectivetemperaturek, meansolarirradiance, scaleheight, stratosphere,
    thermosphere, troposphere,
};
use neptunes::atmosphere::storms::{NeptuneStorm, capeintegrated, capelayer, fujitascale};
use neptunes::atmosphere::winds::{
    OMEGANEPTUNE, beauforttoms, coriolisparameter, ekmanspiraldepth, geostrophicwindspeed,
    gradientwindspeed, thermalwindshear, zonalwindprofile,
};
use neptunes::hydrology::glaciers::{
    GLENFLOWLAWA, ICECREEPEXPONENT, glenstrainrate, icedensity, icerheologyviscosity,
    methaneicecap, shallowicevelocity, watericemantle, watermolarmassfromelements,
};
use neptunes::hydrology::lakes::hypotheticalsubsurfaceocean;
use neptunes::hydrology::oceans::{
    INTERIOROCEANDEPTHM, IONICOCEANDENSITY, MEANSALINITYPSU, dissolvedgashenry, h2omolarmass,
    interiorheatcontentj, ioniclayer, molecularlayer, superioniclayer, surfacelongwaveradiationwm2,
    vaporpressureovermethane,
};
use neptunes::hydrology::rivers::{equatorialjet, neptunejetnetwork, polarjet};

fn main() {
    let tropo = troposphere();
    let strato = stratosphere();
    let thermo = thermosphere();

    let altitudes = [0.0, 10000.0, 50000.0, 100000.0, 500000.0];
    let mut totaldensity = 0.0;
    for &alt in &altitudes {
        let layer: &AtmosphericLayer = if alt < 50000.0 {
            &tropo
        } else if alt < 200000.0 {
            &strato
        } else {
            &thermo
        };
        let temp = layer.temperatureat(alt);
        let press = layer.pressureat(alt);
        let dens = layer.densityat(alt);
        totaldensity += dens;
        println!(
            "[{:>12}] alt={:>8.0}m  T={:.1}K  P={:.0}Pa  ρ={:.4}kg/m³",
            layer.name, alt, temp, press, dens
        );
    }

    let psea = barometricpressure(0.0);
    let phigh = barometricpressure(100000.0);
    let hscale = scaleheight();
    println!(
        "P(1bar)={:.0}Pa  P(100km)={:.0}Pa  Hscale={:.0}m  Σρ={:.4}",
        psea, phigh, hscale, totaldensity
    );
    println!(
        "constants: P0={} T0={} Mair={} Γ={}",
        ONEBARPRESSURE,
        ONEBARTEMPK,
        atmospheremolarmass(),
        LAPSERATETROPOSPHERE
    );

    let current = ClimateState::current();
    let olr = current.outgoinglongwaveradiation();
    let asr = current.absorbedsolarradiation();
    let imbalance = current.energyimbalance();
    let teff = current.effectiveemissiontemperature();
    let greenhouse = current.methanegreenhouseeffectk();
    let ch4forcing = current.ch4radiativeforcing();
    let internalheat = current.internalplussolarheatflux();

    println!(
        "CH4 fraction={:.4}  globalT={:.1}K  albedo={:.2}  internalH={:.1}W/m²",
        current.ch4fraction, current.globalmeantempk, current.albedo, current.internalheatwm2
    );
    println!(
        "OLR={:.1}  ASR={:.1}  imbalance={:.2}  Teff={:.1}K  CH4greenhouse={:.1}K",
        olr, asr, imbalance, teff, greenhouse
    );
    println!(
        "CH4forcing={:.3}W/m²  InternalFlux={:.2}W/m²",
        ch4forcing, internalheat
    );

    let zenith = solarzenith(45.0, 28.32, 0.0);
    let radiance = planck_spectral_radiance(10e-6, 72.0);
    let doubling = co2doublingforcing();
    let irr = meansolarirradiance();
    let teffglob = effectivetemperaturek();
    let sflux = solarfluxatneptune();
    println!(
        "zenith={:.2}°  planck(10µm,72K)={:.6}  2xCO2={:.2}W/m²",
        zenith, radiance, doubling
    );
    println!(
        "solarIrr={:.2}W/m²  Teff={:.1}K  solarFlux={:.3}W/m²",
        irr, teffglob, sflux
    );
    println!(
        "constants: Teq={} Teff={} albedo={}",
        NEPTUNEEQUILIBRIUMTEMPK, NEPTUNEEFFECTIVETEMPK, NEPTUNEALBEDO
    );

    let gds = NeptuneStorm::greatdarkspot();
    let scooter = NeptuneStorm::scooter();
    let sds = NeptuneStorm::smalldarkspot();
    let pi = NeptuneStorm::potentialintensityms(0.4, 52.0);
    println!(
        "GDS: KE={:.2e}J/m³  vort={:.4e}  Rossby={:.0}m  wind@500km={:.1}m/s",
        gds.kineticenergydensity(),
        gds.vorticityestimate(),
        gds.rossbydeformationradius(),
        gds.windatradius(500.0)
    );
    println!(
        "Scooter: KE={:.2e}  vort={:.4e}",
        scooter.kineticenergydensity(),
        scooter.vorticityestimate()
    );
    println!(
        "SDS: KE={:.2e}  Rossby={:.0}m",
        sds.kineticenergydensity(),
        sds.rossbydeformationradius()
    );
    println!("PI(0.4W/m²,52K)={:.1}m/s", pi);

    let fujita = fujitascale(90.0);
    let capesingle = capelayer(80.0, 72.0, 1000.0);
    let cape = capeintegrated(80.0, 72.0, 0.001, 500.0, 40);
    println!(
        "Fujita(90m/s)={}  CAPElayer={:.1}  CAPEintegrated={:.1}J/kg",
        fujita, capesingle, cape
    );

    let f45 = coriolisparameter(45.0);
    let vg = geostrophicwindspeed(0.01, 45.0, 0.45);
    let vgr = gradientwindspeed(0.01, 500000.0, 45.0, 0.45);
    let tw = thermalwindshear(5e-6, 72.0, 45.0, 500e2, 750e2);
    let ekman = ekmanspiraldepth(0.05, 45.0);
    let b7 = beauforttoms(7);
    let zonal30 = zonalwindprofile(30.0);
    let zonal60 = zonalwindprofile(-60.0);
    println!(
        "ω={:.5e}  f(45°)={:.5e}  Vg={:.1}  Vgr={:.1}  thermalwind={:.4}",
        OMEGANEPTUNE, f45, vg, vgr, tw
    );
    println!(
        "Ekman={:.0}m  B7={:.1}m/s  zonal(30°)={:.1}m/s  zonal(-60°)={:.1}m/s",
        ekman, b7, zonal30, zonal60
    );

    let ta = ThermalAnomaly {
        temperaturek: 80.0,
        durationdays: 10,
        altitudem: 50000.0,
    };
    let warming = ta.warmingexcessk();
    let cold = ta.coldexcessk();
    let danger = ta.dangerlevel();
    let methtemp = ThermalAnomaly::methevaportempk();
    let methpress = ThermalAnomaly::methanecondensationpressurepa(72.0);
    let svp = saturationvaporpressuremethanepa(72.0);
    let dp = dewpointmethanek(1e5);
    println!(
        "Anomaly: warm={:.1}K  cold={:.1}K  danger={}  methEvapT={:.1}K",
        warming, cold, danger, methtemp
    );
    println!(
        "CH4condP(72K)={:.0}Pa  SVP(72K)={:.0}Pa  dewpoint(1e5Pa)={:.1}K",
        methpress, svp, dp
    );

    let sup = superioniclayer();
    let ion = ioniclayer();
    let mol = molecularlayer();
    let rhosup = sup.densitykgm3();
    let rhoion = ion.densitykgm3();
    let rhomol = mol.densitykgm3();
    let speedsup = sup.soundspeedms();
    let pdeep = sup.pressureatdepth(3000000.0);
    let ohc = interiorheatcontentj(0.5, 700000.0);
    let lw = surfacelongwaveradiationwm2(72.0, 0.97);
    let vpmethane = vaporpressureovermethane(72.0);
    let co2diss = dissolvedgashenry(40.0, 72.0, 3.4e-4);
    let h2omm = h2omolarmass();
    println!(
        "Ocean: ρsup={:.0} ρion={:.0} ρmol={:.0}  sound={:.0}m/s  P@3Mm={:.2e}Pa",
        rhosup, rhoion, rhomol, speedsup, pdeep
    );
    println!(
        "OHC={:.2e}J  LW={:.1}W/m²  VP_CH4={:.0}Pa  CO2diss={:.6}mol/m³",
        ohc, lw, vpmethane, co2diss
    );
    println!(
        "constants: D={}m  ρion={}  S={}psu  H2O_Mm={}",
        INTERIOROCEANDEPTHM, IONICOCEANDENSITY, MEANSALINITYPSU, h2omm
    );

    let subsea = hypotheticalsubsurfaceocean();
    let depth = subsea.meandepthm();
    let res = subsea.residencetimeyears(1e6);
    let lwlake = subsea.longwaveradiationwm2(0.95);
    println!(
        "SubsurfaceOcean: depth={:.0}m  residence={:.1}yr  LW={:.1}W/m²",
        depth, res, lwlake
    );

    let eqjet = equatorialjet();
    let pjet = polarjet();
    let vmann = eqjet.manningvelocity(3.0, 0.035);
    let fr = eqjet.froudenumber(100.0, 10000.0);
    let re = pjet.reynoldsnumber(50.0, 5000.0);
    let power = eqjet.streampowerwperm();
    let sdischarge = pjet.specificdischargems();
    let net = neptunejetnetwork();
    let outlet = net.outletdischarge();
    let routed = net.routedischarge();
    let routedsum: f64 = routed.iter().sum();
    println!(
        "EqJet: manning={:.1}m/s  Fr={:.4}  power={:.2e}W/m  discharge={:.2e}m/s",
        vmann, fr, power, sdischarge
    );
    println!(
        "PolarJet: Re={:.2e}  Network outlet={:.2e}m³/s  routedΣ={:.2e}",
        re, outlet, routedsum
    );

    let wim = watericemantle();
    let mic = methaneicecap();
    let rhice = icedensity();
    let visc = icerheologyviscosity(70.0, 1e5);
    let glen = glenstrainrate(70.0, 1e5);
    let shallow = shallowicevelocity(5000.0, 0.01, 70.0);
    let wmorph = watermolarmassfromelements();
    println!(
        "WaterIceMantle: massbal={:.4}m/yr  basalτ={:.0}Pa  deform={:.4}m/yr",
        wim.massbalancemyr(),
        wim.basalshearstresspa(),
        wim.deformationvelocitymyr()
    );
    println!(
        "CH4IceCap: massbal={:.4}m/yr  vol={:.2e}m³",
        mic.massbalancemyr(),
        mic.volumem3(1e6)
    );
    println!(
        "ρice={:.0}  η(70K,1e5Pa)={:.2e}Pa·s  glen={:.2e}  shallow={:.4}m/yr",
        rhice, visc, glen, shallow
    );
    println!(
        "constants: n={} A={:.2e} H2Omm={}",
        ICECREEPEXPONENT, GLENFLOWLAWA, wmorph
    );
}
