jupiters 0.0.3

Jupiter celestial simulation crate for the MilkyWay SolarSystem workspace
Documentation
use jupiters::atmosphere::climate::{
    ClimateState, JUPITERALBEDO, co2doublingforcing, planckspectralradiance, solarfluxatjupiter,
    solarzenith,
};
use jupiters::atmosphere::heatwaves::{
    ThermalAnomaly, dewpointammoniak, saturationvaporpressureammoniana,
};
use jupiters::atmosphere::layers::{
    LAPSERATETROPOSPHERE, ONEBARPRESSURE, ONEBARTEMPK, barometricpressure, effectivetemp,
    meansolarirradiance, olr, scaleheight, stratosphere, thermosphere, troposphere,
};
use jupiters::atmosphere::storms::{JupiterStorm, capeintegrated, capelayer, fujitascale};
use jupiters::atmosphere::winds::{
    OMEGAJUPITER, beauforttoms, coriolisparameter, ekmanspiraldepth, geostrophicwindspeed,
    thermalwindshear, zonalwindprofile,
};
use jupiters::hydrology::glaciers::{icedensity, metallichydrogenmantle, watericecloudcondensate};
use jupiters::hydrology::oceans::{metallichydrogenocean, molecularhydrogenfluid, rockycorelayer};
use jupiters::hydrology::rivers::equatorialjet;

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

    let altitudes = [0.0, 10000.0, 30000.0, 50000.0, 200000.0, 500000.0];
    let mut totaldensity = 0.0;
    for &alt in &altitudes {
        let layer = if alt < 50000.0 {
            &tropo
        } else if alt < 320000.0 {
            &strato
        } else {
            &thermo
        };
        let temp = layer.temperatureat(alt);
        let press = layer.pressureat(alt);
        let dens = layer.densityat(alt);
        totaldensity += dens;
        println!(
            "[{:>14}] alt={:>7.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(30000.0);
    let hscale = scaleheight();
    println!(
        "P(1bar)={:.0}Pa  P(30km)={:.0}Pa  Hscale={:.0}m  Σρ={:.4}",
        psea, phigh, hscale, totaldensity
    );
    println!(
        "constants: P0={} T0={} Γ={}",
        ONEBARPRESSURE, ONEBARTEMPK, LAPSERATETROPOSPHERE
    );

    let state = ClimateState::current();
    let olrval = state.outgoinglongwaveradiation();
    let asr = state.absorbedsolarradiation();
    let totalflux = state.internalplussolarheatflux();
    let imbalance = state.energyimbalance();
    let teff = state.effectiveemissiontemperature();
    let greenhouse = state.ammoniagreenhouseeffectk();
    let nh3rf = state.nh3radiativeforcing();

    println!(
        "NH3: fraction={:.5}  albedo={:.3}  internalHeat={:.2}W/m²",
        state.nh3fraction, state.albedo, state.internalheatwm2
    );
    println!(
        "OLR={:.1}  ASR={:.1}  total={:.1}  imbalance={:.2}W/m²",
        olrval, asr, totalflux, imbalance
    );
    println!(
        "Teff={:.1}K  greenhouse={:.1}K  NH3forcing={:.3}W/m²",
        teff, greenhouse, nh3rf
    );

    let solflux = solarfluxatjupiter();
    let teffglobal = effectivetemp();
    let olrglobal = olr();
    let irradiance = meansolarirradiance();
    let zenith = solarzenith(10.0, 3.13, 0.0);
    let radiance = planckspectralradiance(10e-6, 165.0);
    let doubling = co2doublingforcing();
    println!(
        "Solar: flux={:.1}W/m²  irradiance={:.1}  Teff={:.1}K  OLR={:.1}",
        solflux, irradiance, teffglobal, olrglobal
    );
    println!(
        "zenith(10°)={:.2}°  planck(10µm,165K)={:.6}  2xCO2={:.2}W/m²  albedo={}",
        zenith, radiance, doubling, JUPITERALBEDO
    );

    let grs = JupiterStorm::greatredspot();
    let ob = JupiterStorm::ovalba();
    let ke = grs.kineticenergydensity();
    let vort = grs.vorticityestimate();
    let rossby = grs.rossbydeformationradius();
    let windr = grs.windatradius(20000.0);
    let pi = JupiterStorm::potentialintensityms(5.44, 110.0);
    let fujita = fujitascale(grs.maxwindspeedms);
    let cape1 = capelayer(180.0, 165.0, 1000.0);
    let cape = capeintegrated(170.0, 165.0, 0.002, 500.0, 30);
    println!(
        "GRS: wind={:.0}m/s  R={:.0}km  KE={:.1}J/m³  vort={:.2e}s⁻¹",
        grs.maxwindspeedms, grs.radiuskm, ke, vort
    );
    println!(
        "OvalBA: wind={:.0}m/s  R={:.0}km  Rossby(GRS)={:.0}m",
        ob.maxwindspeedms, ob.radiuskm, rossby
    );
    println!(
        "wind@20Kkm={:.1}m/s  PI={:.1}m/s  Fujita(GRS)={}  CAPE1={:.1}  CAPEint={:.1}J/kg",
        windr, pi, fujita, cape1, cape
    );

    let f30 = coriolisparameter(30.0);
    let vg = geostrophicwindspeed(0.01, 30.0, 0.16);
    let tw = thermalwindshear(5e-6, 165.0, 30.0, 500e2, 750e2);
    let ekman = ekmanspiraldepth(0.05, 30.0);
    let b7 = beauforttoms(7);
    let vzonal = zonalwindprofile(0.0);
    println!(
        "ω={:.5e}  f(30°)={:.5e}  Vg={:.1}  thermalwind={:.4}",
        OMEGAJUPITER, f30, vg, tw
    );
    println!(
        "Ekman={:.0}m  B7={:.1}m/s  zonalwind(eq)={:.0}m/s",
        ekman, b7, vzonal
    );

    let anomaly = ThermalAnomaly {
        temperaturek: 200.0,
        durationdays: 10,
        altitudem: 0.0,
    };
    let excess = anomaly.warmingexcessk();
    let danger = anomaly.dangerlevel();
    let nh3boil = ThermalAnomaly::ammoniaboilingtempk();
    let svp = saturationvaporpressureammoniana(195.0);
    let dp = dewpointammoniak(10000.0);
    println!(
        "Anomaly: T={:.0}K  excess={:.1}K  danger={}  NH3boil={:.1}K",
        anomaly.temperaturek, excess, danger, nh3boil
    );
    println!("SVP(195K)={:.0}Pa  dewpoint(10kPa)={:.1}K", svp, dp);

    let mantle = metallichydrogenmantle();
    let cloud = watericecloudcondensate();
    let rhice = icedensity();
    println!(
        "Ice: ρ={:.0}kg/m³  mantle_thick={:.0}m  cloud_mb={:.4}m/yr",
        rhice,
        mantle.thicknessm,
        cloud.massbalancemyr()
    );

    let h2 = molecularhydrogenfluid();
    let mh = metallichydrogenocean();
    let core = rockycorelayer();
    println!(
        "Interior: H2(ρ={:.0},T={:.0}K)  MH(ρ={:.0},T={:.0}K)  Core(ρ={:.0},T={:.0}K)",
        h2.densitykgm3,
        h2.temperaturek,
        mh.densitykgm3,
        mh.temperaturek,
        core.densitykgm3,
        core.temperaturek
    );
    println!(
        "sound(H2)={:.0}m/s  sound(MH)={:.0}m/s",
        h2.soundspeedms(),
        mh.soundspeedms()
    );

    let jet = equatorialjet();
    let q = jet.dischargem3s();
    let fr = jet.froudenumber();
    println!(
        "EqJet: v={:.0}m/s  Q={:.2e}m³/s  Fr={:.2}",
        jet.velocityms, q, fr
    );
}