Struct OceanEndpoint 

pub struct OceanEndpoint {
    pub wind_speed_m_s: f64,
    pub wind_direction: [f64; 2],
    pub fetch_km: f64,
    pub mean_depth_m: f64,
    pub salinity_psu: f64,
    pub surface_temperature_k: f64,
    pub surface_gravity_m_s2: f64,
    pub ior: f64,
    pub absorption_coefficients_rgb_m: [f64; 3],
    pub foam_threshold_wind_m_s: f64,
    pub grid_size: u32,
    pub patch_size_m: f64,
}

Fields

wind_speed_m_s: f64

wind_direction: [f64; 2]

fetch_km: f64

mean_depth_m: f64

salinity_psu: f64

surface_temperature_k: f64

surface_gravity_m_s2: f64

ior: f64

absorption_coefficients_rgb_m: [f64; 3]

foam_threshold_wind_m_s: f64

grid_size: u32

patch_size_m: f64

Implementations

impl OceanEndpoint

pub fn earth_atlantic() -> Self

Examples found in repository

examples/tidal_sim.rs (line 554)

fn main() {
    let orbit = EarthOrbit::new();
    let perioddays = orbit.orbitalperioddays();
    let vperihelion = orbit.velocityatdistance(orbit.perihelionm());
    let vaphelion = orbit.velocityatdistance(orbit.aphelionm());
    let energy = orbit.specific_orbital_energy();
    let angmom = orbit.specific_angular_momentum();
    let vesc = EarthOrbit::escape_velocity_at_surface();
    let fsun = orbit.gravitational_force_sun();
    let r = orbit.current_radius();
    let vmean = orbit.mean_orbital_velocity();
    println!(
        "Orbit: P={:.2}d  vperi={:.0}m/s  vaph={:.0}m/s  vmean={:.0}m/s",
        perioddays, vperihelion, vaphelion, vmean
    );
    println!(
        "  E={:.2e}J/kg  L={:.2e}m²/s  vesc={:.0}m/s  Fsun={:.2e}N  r={:.3e}m",
        energy, angmom, vesc, fsun, r
    );
    println!(
        "  a={:.3e}m  e={}  i={}°  Ω={}°  ω={}°",
        SEMIMAJORAXIS,
        ECCENTRICITY,
        INCLINATIONDEG,
        LONGITUDEASCENDINGNODEDEG,
        ARGUMENTPERIHELIONDEG
    );

    let rot = EarthRotation::new();
    let vequator = rot.surfacevelocityatlatitude(0.0);
    let v45 = rot.surfacevelocityatlatitude(45.0);
    let accel = rot.centripetalaccelerationatlatitude(0.0);
    let coriolis = rot.coriolisparameter(45.0);
    let moi = rot.momentofinertia();
    let ke = rot.rotationalkineticenergy();
    let lrot = rot.angular_momentum();
    let prec = rot.precession_rate_rad_per_year();
    let daylensummer = rot.day_length_variation_due_to_tilt(172, 60.0);
    let daylenwinter = rot.day_length_variation_due_to_tilt(355, 60.0);
    println!(
        "Rotation: veq={:.0}m/s  v45={:.0}m/s  ac={:.4}m/s²  f45={:.5e}",
        vequator, v45, accel, coriolis
    );
    println!(
        "  I={:.3e}kg·m²  KE={:.3e}J  L={:.3e}kg·m²/s  prec={:.2e}rad/yr",
        moi, ke, lrot, prec
    );
    println!(
        "  day@60°N: summer={:.1}h  winter={:.1}h",
        daylensummer, daylenwinter
    );
    println!(
        "  sidereal={:.4}s  solar={}s  tilt={}°={:.4}rad  precperiod={}yr",
        SIDEREALDAYS, SOLARDAYS, AXIALTILTDEG, AXIALTILTRAD, PRECESSIONPERIODYEARS
    );

    let moontide = TidalForce::frommoon();
    let suntide = TidalForce::fromsun();
    let amoon = moontide.tidalacceleration();
    let asun = suntide.tidalacceleration();
    let pot = moontide.tidalpotential(0.0);
    let bulgemoon = moontide.tidalbulgeheight();
    let bulgesun = suntide.tidalbulgeheight();
    let grav = moontide.gravitationalattraction();
    let spring = springtideamplitude();
    let neap = neaptideamplitude();
    let ratio = lunartosolartideratio();
    println!(
        "Tides: amoon={:.2e}  asun={:.2e}  ratio={:.2}",
        amoon, asun, ratio
    );
    println!(
        "  bulgemoon={:.3}m  bulgesun={:.3}m  spring={:.3}m  neap={:.3}m",
        bulgemoon, bulgesun, spring, neap
    );
    println!("  potential(0)={:.2e}  Fgrav={:.2e}N", pot, grav);

    let chix = chicxulubequivalent();
    let tung = tunguskaequivalent();
    let kechix = chix.kineticenergymt();
    let crater = chix.craterdiameterm(2700.0);
    let fireball = chix.fireballradiusm();
    let ejecta = chix.ejectavolumem3(2700.0);
    let vimp = chix.impactvelocity();
    let ketung = tung.kineticenergymt();
    println!(
        "Chicxulub: {:.2e}Mt  crater={:.0}m  fireball={:.0}m  ejecta={:.2e}m³  vimp={:.0}m/s",
        kechix, crater, fireball, ejecta, vimp
    );
    println!("Tunguska: {:.2e}Mt", ketung);

    let mut moon = MoonState::new();
    let source = match moon.source.get() {
        MoonSource::Binary => "binary",
        MoonSource::Simulation => "simulated",
    };
    let pos0 = moon.position();
    moon.step(3600.0);
    let pos1 = moon.position();
    let gmoon = moon.gravityat(EARTHMOONDISTANCE);
    println!(
        "Moon({}): pos0=({:.0},{:.0})  pos1h=({:.0},{:.0})  g@dist={:.4e}m/s²",
        source, pos0.0, pos0.1, pos1.0, pos1.1, gmoon
    );
    println!(
        "  M={:.3e}kg  R={:.4e}m  d={:.3e}m  P={:.0}s",
        LUNARMASS, LUNARRADIUS, EARTHMOONDISTANCE, LUNARORBITALPERIOD
    );

    let iss = ArtificialSatellite::leo("ISS", 420000.0, 408000.0);
    let geo = ArtificialSatellite::geo("GEO-Sat", 300.0);
    let custom = ArtificialSatellite::new("Molniya", 1500.0, 500000.0, 0.7, 1.1);
    let piss = iss.orbitalperiods();
    let viss = iss.orbitalvelocityms();
    let pgeo = geo.orbitalperiods();
    let gsurf = iss.gravityatsurface();
    let posiss = iss.position();
    let posgeo = geo.position();

    let orbittype = |sat: &ArtificialSatellite| -> &str {
        match sat.orbittype {
            OrbitType::LEO => "LEO",
            OrbitType::MEO => "MEO",
            OrbitType::GEO => "GEO",
            OrbitType::HEO => "HEO",
            OrbitType::Custom => "Custom",
        }
    };

    println!(
        "ISS({}): P={:.0}s  v={:.0}m/s  gsurf={:.2}  pos=({:.0},{:.0},{:.0})",
        orbittype(&iss),
        piss,
        viss,
        gsurf,
        posiss.0,
        posiss.1,
        posiss.2
    );
    println!(
        "GEO({}): P={:.0}s  pos=({:.0},{:.0},{:.0})",
        orbittype(&geo),
        pgeo,
        posgeo.0,
        posgeo.1,
        posgeo.2
    );
    println!(
        "Custom({}): e={:.1}  i={:.1}rad",
        orbittype(&custom),
        custom.eccentricity,
        custom.inclinationrad
    );

    let mut constellation = Constellation::new("Starlink-test");
    constellation.add(ArtificialSatellite::leo("S1", 260.0, 550000.0));
    constellation.add(ArtificialSatellite::leo("S2", 260.0, 550000.0));
    let mut satstepped = iss;
    satstepped.step(60.0);
    constellation.stepall(60.0);
    let positions = constellation.positions();
    println!(
        "Constellation: {} sats  positions={}",
        positions.len(),
        positions.len()
    );

    let dt = DateTime::new(2024, 6, 21, 12, 0, 0.0);
    let jd = dt.tojuliandate();
    let dtback = DateTime::fromjuliandate(jd);
    let unix = dt.tounixtimestamp();
    let dtunix = DateTime::fromunixtimestamp(unix);
    println!(
        "DateTime: {}-{}-{} {}:{}:{:.0}  JD={:.4}  unix={:.0}",
        dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second, jd, unix
    );
    println!(
        "  roundtrip JD: {}-{}-{}  roundtrip unix: {}-{}-{}",
        dtback.year, dtback.month, dtback.day, dtunix.year, dtunix.month, dtunix.day
    );
    println!(
        "  J2000JD={} UNIXJD={} SPD={}",
        J2000JD, UNIXEPOCHJD, SECONDSPERDAY
    );

    let mut epoch = Epoch::j2000();
    let epochmjd = Epoch::frommjd(51544.5);
    let centuries = epoch.centuriessincej2000();
    let days = epoch.dayssincej2000();
    let gmst = epoch.gmstdegrees();
    let mjd = epoch.tomjd();
    epoch.advancedays(365.25);
    let centuries1yr = epoch.centuriessincej2000();
    let mut epoch2 = Epoch::fromjd(J2000JD);
    epoch2.advanceseconds(86400.0);
    println!(
        "Epoch: T0={:.4}  days={:.1}  GMST={:.2}°  MJD={:.1}",
        centuries, days, gmst, mjd
    );
    println!(
        "  +1yr: T={:.6}  frommjd.jd={:.1}  epoch2.days={:.1}",
        centuries1yr,
        epochmjd.juliandate,
        epoch2.dayssincej2000()
    );
    println!(
        "  J2000={} J1950={} MJDOFF={}",
        J2000EPOCH, J1950EPOCH, MJDOFFSET
    );

    let mut ts = TimeScale::realtime();
    ts.step(1.0);
    let simdt = ts.simulationdt(1.0);
    let hours = ts.simhours();
    let daysts = ts.simdays();
    let years = ts.simyears();
    ts.pause();
    ts.resume();
    ts.togglepause();
    ts.setspeed(100.0);
    let mut ff = TimeScale::fastforward(10.0);
    ff.step(1.0);
    let sm = TimeScale::slowmotion(2.0);
    println!(
        "TimeScale: dt={:.1}  h={:.6}  d={:.8}  yr={:.10}  ffsim={:.1}s  smspeed={:.1}",
        simdt, hours, daysts, years, ff.simulationtimes, sm.speedmultiplier
    );

    let sunpos = SolarPosition::compute(jd, 48.8, 2.3);
    let above = sunpos.isabovehorizon();
    let distsun = sunpos.distancem();
    println!(
        "Solar: el={:.1}°  az={:.1}°  above={}  dist={:.3e}m  tilt={}°",
        sunpos.elevationdeg, sunpos.azimuthdeg, above, distsun, EARTHAXIALTILTDEG
    );

    let dnc = DayNightCycle::new(jd);
    let stateparis = dnc.stateat(48.8, 2.3);
    let statetext = match stateparis {
        DaylightState::Day => "day",
        DaylightState::Night => "night",
        DaylightState::CivilTwilight => "civiltwilight",
        DaylightState::NauticalTwilight => "nauticaltwilight",
        DaylightState::AstronomicalTwilight => "astrotwilight",
    };
    let terminator = dnc.terminatorpoints(36);
    let ambient = dnc.ambientlight(48.8, 2.3);
    println!(
        "DayNight: Paris={}  ambient={:.2}  terminatorpts={}",
        statetext,
        ambient,
        terminator.len()
    );

    let seasonstate = seasonat(jd, 48.8);
    let seasonname = match seasonstate.seasonnorth {
        Season::Spring => "spring",
        Season::Summer => "summer",
        Season::Autumn => "autumn",
        Season::Winter => "winter",
    };
    let southname = match seasonstate.seasonsouth {
        Season::Spring => "spring",
        Season::Summer => "summer",
        Season::Autumn => "autumn",
        Season::Winter => "winter",
    };
    let (sublat, sublon) = subsolarpoint(jd);
    println!(
        "Season: N={}  S={}  decl={:.2}°  dayh={:.2}  subsolar=({:.1},{:.1})",
        seasonname,
        southname,
        seasonstate.solardeclinationdeg,
        seasonstate.daylengthhours,
        sublat,
        sublon
    );
    println!(
        "  tilt={}°  tropyear={}d  vernaljd={}",
        SEASONSTILT, TROPICALYEARDAYS, VERNALEQUINOXJD
    );

    let paris = LatLon::new(48.8566, 2.3522, 35.0);
    let tokyo = LatLon::new(35.6762, 139.6503, 40.0);
    let ecefparis = paris.toecef();
    let cart = paris.tocartesiansimple();
    let distpt = paris.distanceto(&tokyo);
    let latlonback = ecefparis.tolatlon();
    println!(
        "Paris->ECEF: ({:.0},{:.0},{:.0})  cart=({:.0},{:.0},{:.0})",
        ecefparis.x, ecefparis.y, ecefparis.z, cart[0], cart[1], cart[2]
    );
    println!(
        "  dist(Paris-Tokyo)={:.0}m  roundtriplat={:.4}°",
        distpt, latlonback.latdeg
    );
    println!(
        "  f={:.9}  a={:.0}m  b={:.3}m",
        EARTHFLATTENING, EARTHSEMIMAJORM, EARTHSEMIMINORM
    );

    let regions = RegionDatabase::continents();
    let region = regions.pointinregion(48.8, 2.3);
    if let Some(r) = region {
        let rtype = match r.regiontype {
            RegionType::Continent => "continent",
            RegionType::Ocean => "ocean",
            RegionType::Sea => "sea",
            RegionType::Country => "country",
            RegionType::Island => "island",
        };
        println!(
            "Region(48.8,2.3): {} ({})  area={:.0}km²",
            r.name, rtype, r.areakm2
        );
    }

    let elevprov = ElevationProvider::global(30.0);
    let eleveverest = elevprov.sample(27.988, 86.925);
    let notables = ElevationProvider::notableelevations();
    println!(
        "Elevation(Everest)={:.0}m  notables={}",
        eleveverest,
        notables.len()
    );

    let bathy = BathymetryData::global(60.0);
    let depthmariana = bathy.sample(11.35, 142.2);
    let isocean = bathy.isocean(0.0, -30.0);
    let stats = BathymetryData::oceanstats();
    let basins = OceanBasin::majorbasins();
    println!(
        "Bathy(Mariana)={:.0}m  isocean(0,-30)={}  basins={}  avgdepth={:.0}m",
        depthmariana,
        isocean,
        basins.len(),
        stats.avgdepthm
    );

    let mut hmap = Heightmap::new(64);
    hmap.generateprocedural(42, 6, 0.5, 2.0);
    let hsample = hmap.sample(45.0, 10.0);
    let rat = hmap.radiusat(45.0, 10.0);
    println!(
        "Heightmap(64): sample(45,10)={:.1}m  radius={:.0}m",
        hsample, rat
    );

    let config = LodConfig::default();
    let mut lod = LodTerrain::new(config);
    lod.update([6371000.0 + 1000.0, 0.0, 0.0]);

    let mesh = TerrainMesh::fromregion(0.0, 10.0, 0.0, 10.0, 8, &|lat, lon| hmap.sample(lat, lon));
    println!(
        "Mesh: {} vertices  {} triangles",
        mesh.vertexcount(),
        mesh.trianglecount()
    );

    let classifier = BiomeClassifier::default();
    let biome = classifier.classify(2000.0, 45.0, 0.6);
    let biomename = match biome {
        Biome::Ocean => "ocean",
        Biome::Beach => "beach",
        Biome::Desert => "desert",
        Biome::Grassland => "grassland",
        Biome::Forest => "forest",
        Biome::Tundra => "tundra",
        Biome::Snow => "snow",
        Biome::Mountain => "mountain",
        Biome::Volcanic => "volcanic",
        Biome::Taiga => "taiga",
    };
    let splat = classifier.splat(2000.0, 45.0, 0.6);
    println!(
        "Biome(2000m,45°,0.6)={}  splattop={:.2}",
        biomename, splat.weights[0].1
    );

    let face = Face::PosZ;
    println!("LOD face={:?}", face);

    let all_mats = PbrMaterial::all_earth();
    let total_subsurface: f32 = all_mats.iter().map(|m| m.subsurface).sum();
    let avg_ior: f32 = all_mats.iter().map(|m| m.ior).sum::<f32>() / all_mats.len() as f32;
    let avg_normal: f32 =
        all_mats.iter().map(|m| m.normal_strength).sum::<f32>() / all_mats.len() as f32;
    let emissive_count = all_mats
        .iter()
        .filter(|m| m.emissive[0] > 0.0 || m.emissive[1] > 0.0 || m.emissive[2] > 0.0)
        .count();
    println!(
        "Materials({}) avg_ior={:.3} avg_normal={:.2} total_sss={:.2} emissive={}",
        all_mats.len(),
        avg_ior,
        avg_normal,
        total_subsurface,
        emissive_count
    );
    for m in &all_mats {
        let luminance = m.albedo[0] * 0.3 + m.albedo[1] * 0.59 + m.albedo[2] * 0.11;
        let fresnel_r0 = ((m.ior - 1.0) / (m.ior + 1.0)).powi(2);
        println!(
            "  {}: L={:.3} R={:.2} M={:.2} F₀={:.4} IOR={:.3} SSS={:.2} N={:.2} α={:.2} E=({:.2},{:.2},{:.2})",
            m.name,
            luminance,
            m.roughness,
            m.metallic,
            fresnel_r0,
            m.ior,
            m.subsurface,
            m.normal_strength,
            m.albedo[3],
            m.emissive[0],
            m.emissive[1],
            m.emissive[2]
        );
    }

    let shterrain = ShaderEndpoint::terrain();
    let shatm = ShaderEndpoint::atmosphere();
    let shocean = ShaderEndpoint::ocean();
    let shnights = ShaderEndpoint::night_lights();
    let total_uniform_bytes: usize = [&shterrain, &shatm, &shocean, &shnights]
        .iter()
        .map(|sh| {
            sh.name.len()
                + sh.uniforms
                    .iter()
                    .map(|u| match &u.value {
                        UniformValue::Float(_) => 4,
                        UniformValue::Vec3(_) => 12,
                        UniformValue::Vec4(_) => 16,
                        UniformValue::Int(_) => 4,
                    })
                    .sum::<usize>()
        })
        .sum();
    println!(
        "Shaders: terrain={} atm={} ocean={} night={} total_bytes={}",
        shterrain.uniforms.len(),
        shatm.uniforms.len(),
        shocean.uniforms.len(),
        shnights.uniforms.len(),
        total_uniform_bytes
    );
    for u in &shterrain.uniforms {
        let val = match &u.value {
            UniformValue::Float(f) => format!("f{:.2}", f),
            UniformValue::Vec3(v) => format!("v3({:.1},{:.1},{:.1})", v[0], v[1], v[2]),
            UniformValue::Vec4(v) => format!("v4({:.1},{:.1},{:.1},{:.1})", v[0], v[1], v[2], v[3]),
            UniformValue::Int(i) => format!("i{}", i),
        };
        println!("  {}={}", u.name, val);
    }

    let atm = AtmosphereEndpoint::earth();
    let scale_height_ratio = atm.rayleigh_scale_height_m / atm.mie_scale_height_m;
    let atmo_shell_vol = 4.0 / 3.0
        * std::f64::consts::PI
        * ((atm.planet_radius_m + atm.atmosphere_height_m).powi(3) - atm.planet_radius_m.powi(3));
    let air_density_sl = atm.sea_level_pressure_pa
        / (8.314_462_618 / atm.mean_molar_mass_kg_mol * atm.sea_level_temperature_k);
    let total_atmo_mass_kg = air_density_sl * atmo_shell_vol * 0.5;
    let total_rayleigh_scatter: f64 = atm
        .species
        .iter()
        .map(|s| {
            let n_s = atm.sea_level_number_density_m3 * s.volume_fraction;
            let king =
                (6.0 + 3.0 * s.depolarization_factor) / (6.0 - 7.0 * s.depolarization_factor);
            let dn = s.refractive_index_stp - 1.0;
            8.0 * std::f64::consts::PI.powi(3) * (2.0 * dn).powi(2) * king
                / (3.0 * n_s * (550e-9_f64).powi(4))
        })
        .sum();
    let ozone_optical_depth = atm.ozone_column_du
        * 1e-5
        * (-((25000.0 - atm.ozone_peak_altitude_m) / 5000.0).powi(2)).exp();
    let g = atm.mie_asymmetry_g;
    let hg_forward = (1.0 - g * g) / (1.0 + g * g - 2.0 * g).powf(1.5);
    println!(
        "Atmosphere: M_atm={:.3e}kg  H_ray/H_mie={:.2}  σ_ray550={:.4e}  O₃τ={:.4e}  HG_fwd={:.3}  S₀={:.0}W/m²",
        total_atmo_mass_kg,
        scale_height_ratio,
        total_rayleigh_scatter,
        ozone_optical_depth,
        hg_forward,
        atm.sun_irradiance_w_m2
    );
    for s in &atm.species {
        let king = (6.0 + 3.0 * s.depolarization_factor) / (6.0 - 7.0 * s.depolarization_factor);
        println!(
            "  {} ({}): f={:.6}  M={:.4e}kg/mol  n={:.6}  δ={:.3}  K={:.4}",
            s.name,
            s.symbol,
            s.volume_fraction,
            s.molar_mass_kg_mol,
            s.refractive_index_stp,
            s.depolarization_factor,
            king
        );
    }
    println!(
        "  β_rgb=({:.4e},{:.4e},{:.4e})  mie_coeff={:.2e}",
        atm.rayleigh_coefficients_rgb[0],
        atm.rayleigh_coefficients_rgb[1],
        atm.rayleigh_coefficients_rgb[2],
        atm.mie_coefficient
    );

    let atlantic = OceanEndpoint::earth_atlantic();
    let pacific = OceanEndpoint::earth_pacific();
    let arctic = OceanEndpoint::earth_arctic();
    for ocean in [&atlantic, &pacific, &arctic] {
        let fresnel_r0 = ((ocean.ior - 1.0) / (ocean.ior + 1.0)).powi(2);
        let fetch_m = ocean.fetch_km * 1000.0;
        let grav = ocean.surface_gravity_m_s2;
        let ws = ocean.wind_speed_m_s.max(0.1);
        let omega_p = grav / ws * (0.74 * (grav * fetch_m / (ws * ws)).powf(-0.33));
        let sig_wave_h = 4.0 * (8.1e-3 * grav / omega_p.powi(2)).sqrt();
        let foam_fraction = if ws > ocean.foam_threshold_wind_m_s {
            3.84e-6 * (ws - ocean.foam_threshold_wind_m_s).powf(3.41)
        } else {
            0.0
        };
        let dx = ocean.patch_size_m / ocean.grid_size as f64;
        let penetration_rgb = [
            1.0 / ocean.absorption_coefficients_rgb_m[0],
            1.0 / ocean.absorption_coefficients_rgb_m[1],
            1.0 / ocean.absorption_coefficients_rgb_m[2],
        ];
        let wind_bearing = ocean.wind_direction[1]
            .atan2(ocean.wind_direction[0])
            .to_degrees();
        println!(
            "Ocean(d={:.0}m,S={:.1}psu,T={:.1}K): R₀={:.4}  Hs={:.2}m  foam={:.4}  dx={:.1}m  pen_g={:.1}m  wind={:.0}°",
            ocean.mean_depth_m,
            ocean.salinity_psu,
            ocean.surface_temperature_k,
            fresnel_r0,
            sig_wave_h,
            foam_fraction,
            dx,
            penetration_rgb[1],
            wind_bearing
        );
    }

    let cumulus = CloudLayer::cumulus();
    let stratus = CloudLayer::stratus();
    let cirrus = CloudLayer::cirrus();
    let cb = CloudLayer::cumulonimbus();
    let cloudtype = |ct: &CloudType| match ct {
        CloudType::Cumulus => "Cu",
        CloudType::Stratus => "St",
        CloudType::Cirrus => "Ci",
        CloudType::Cumulonimbus => "Cb",
        CloudType::Altocumulus => "Ac",
        CloudType::Stratocumulus => "Sc",
        CloudType::Nimbostratus => "Ns",
    };
    for l in [&cumulus, &stratus, &cirrus, &cb] {
        let r_eff = l.droplet_radius_um * 1e-6;
        let lwc = 0.3 * l.density;
        let tau = 1.5 * 2.0 * lwc * l.thickness_m * l.coverage / (1000.0 * r_eff);
        let ice_corr = 1.0 - 0.15 * l.ice_fraction;
        let transmittance = (-(tau * ice_corr)).exp();
        let wind_dir = l.wind_direction[1].atan2(l.wind_direction[0]).to_degrees();
        println!(
            "  {} base={:.0}m thick={:.0}m cov={:.2} ρ={:.2} τ={:.2} T={:.4} r={:.0}µm ice={:.1} wind={:.1}m/s@{:.0}°",
            cloudtype(&l.cloud_type),
            l.base_altitude_m,
            l.thickness_m,
            l.coverage,
            l.density,
            tau * ice_corr,
            transmittance,
            l.droplet_radius_um,
            l.ice_fraction,
            l.wind_speed_m_s,
            wind_dir
        );
    }

    let csys = CloudSystemEndpoint::earth_default();
    let stormy = CloudSystemEndpoint::earth_stormy();
    let total_od_default: f64 = csys
        .layers
        .iter()
        .map(|l| {
            let r = l.droplet_radius_um * 1e-6;
            1.5 * 2.0 * 0.3 * l.density * l.thickness_m * l.coverage / (1000.0 * r)
                * (1.0 - 0.15 * l.ice_fraction)
        })
        .sum();
    let total_od_stormy: f64 = stormy
        .layers
        .iter()
        .map(|l| {
            let r = l.droplet_radius_um * 1e-6;
            1.5 * 2.0 * 0.3 * l.density * l.thickness_m * l.coverage / (1000.0 * r)
                * (1.0 - 0.15 * l.ice_fraction)
        })
        .sum();
    println!(
        "CloudSystem: default {} layers τ_total={:.2}  stormy {} layers τ_total={:.2}",
        csys.layers.len(),
        total_od_default,
        stormy.layers.len(),
        total_od_stormy
    );

    let rainforest = tropicalrainforest();
    let boreal = borealforest();
    let reef = coralreef();
    println!("Ecosystems:");
    for eco in [&rainforest, &boreal, &reef] {
        let sp = eco.totalspecies();
        let ind = eco.totalindividuals();
        let shannon = eco.shannonindex();
        let simpson = eco.simpsondiversity();
        let areasp = eco.expectedspeciesfromarea(0.25, 100.0);
        let npp = eco.totalnppgcyr();
        let turnover = eco.biomassturnovertimeyr(15000.0);
        println!(
            "  {} species  {} ind  H'={:.2}  D={:.4}  SAR={:.0}  NPP={:.0}  τ={:.1}yr",
            sp, ind, shannon, simpson, areasp, npp, turnover
        );
    }

    let tb = tropicalbroadleaf();
    let tg = temperategrassland();
    let cf = coniferousforest();
    let photo = tb.photosynthesisrate(400.0, 1500.0, 28.0);
    let canopy = tb.canopyphotosynthesis(photo);
    let transp = tg.transpirationmmday(1.5, 0.3);
    let nppveg = cf.nppkgcm2yr(15.0);
    let crt = cf.carbonresidencetimeyr(nppveg);
    let light = beerlambertlightextinction(2000.0, 4.0, 0.5);
    println!(
        "Vegetation: photo={:.2}  canopy={:.2}  transp={:.2}mm/d  NPP={:.3}  CRT={:.1}yr  light={:.1}",
        photo, canopy, transp, nppveg, crt, light
    );

    let elephant = africanelephant();
    let whale = bluewhale();
    let grel = elephant.growthrate();
    let proj = elephant.projectforward(10.0);
    let metab = whale.metabolicratew();
    let range = elephant.homerangekm2();
    let gentime = whale.generationtimeyears();
    let lifespan = elephant.maxlifespanyears();
    println!(
        "Elephant: r={:.3}  N(+10yr)={:.0}  range={:.1}km²  lifespan={:.0}yr",
        grel, proj, range, lifespan
    );
    println!("Whale: metab={:.0}W  gen={:.1}yr", metab, gentime);

    let mut pp = PredatorPrey {
        prey: africanelephant(),
        predator: bluewhale(),
        attackrate: 0.01,
        conversionefficiency: 0.1,
        predatordeathrate: 0.05,
    };
    let preyr = pp.preygrowthrate();
    let predr = pp.predatorgrowthrate();
    pp.step(0.1);
    println!(
        "Lotka-Volterra: preyr={:.1}  predr={:.1}  -> prey={:.1}  pred={:.2}",
        preyr, predr, pp.prey.count, pp.predator.count
    );
}

pub fn earth_pacific() -> Self

Examples found in repository

examples/tidal_sim.rs (line 555)

fn main() {
    let orbit = EarthOrbit::new();
    let perioddays = orbit.orbitalperioddays();
    let vperihelion = orbit.velocityatdistance(orbit.perihelionm());
    let vaphelion = orbit.velocityatdistance(orbit.aphelionm());
    let energy = orbit.specific_orbital_energy();
    let angmom = orbit.specific_angular_momentum();
    let vesc = EarthOrbit::escape_velocity_at_surface();
    let fsun = orbit.gravitational_force_sun();
    let r = orbit.current_radius();
    let vmean = orbit.mean_orbital_velocity();
    println!(
        "Orbit: P={:.2}d  vperi={:.0}m/s  vaph={:.0}m/s  vmean={:.0}m/s",
        perioddays, vperihelion, vaphelion, vmean
    );
    println!(
        "  E={:.2e}J/kg  L={:.2e}m²/s  vesc={:.0}m/s  Fsun={:.2e}N  r={:.3e}m",
        energy, angmom, vesc, fsun, r
    );
    println!(
        "  a={:.3e}m  e={}  i={}°  Ω={}°  ω={}°",
        SEMIMAJORAXIS,
        ECCENTRICITY,
        INCLINATIONDEG,
        LONGITUDEASCENDINGNODEDEG,
        ARGUMENTPERIHELIONDEG
    );

    let rot = EarthRotation::new();
    let vequator = rot.surfacevelocityatlatitude(0.0);
    let v45 = rot.surfacevelocityatlatitude(45.0);
    let accel = rot.centripetalaccelerationatlatitude(0.0);
    let coriolis = rot.coriolisparameter(45.0);
    let moi = rot.momentofinertia();
    let ke = rot.rotationalkineticenergy();
    let lrot = rot.angular_momentum();
    let prec = rot.precession_rate_rad_per_year();
    let daylensummer = rot.day_length_variation_due_to_tilt(172, 60.0);
    let daylenwinter = rot.day_length_variation_due_to_tilt(355, 60.0);
    println!(
        "Rotation: veq={:.0}m/s  v45={:.0}m/s  ac={:.4}m/s²  f45={:.5e}",
        vequator, v45, accel, coriolis
    );
    println!(
        "  I={:.3e}kg·m²  KE={:.3e}J  L={:.3e}kg·m²/s  prec={:.2e}rad/yr",
        moi, ke, lrot, prec
    );
    println!(
        "  day@60°N: summer={:.1}h  winter={:.1}h",
        daylensummer, daylenwinter
    );
    println!(
        "  sidereal={:.4}s  solar={}s  tilt={}°={:.4}rad  precperiod={}yr",
        SIDEREALDAYS, SOLARDAYS, AXIALTILTDEG, AXIALTILTRAD, PRECESSIONPERIODYEARS
    );

    let moontide = TidalForce::frommoon();
    let suntide = TidalForce::fromsun();
    let amoon = moontide.tidalacceleration();
    let asun = suntide.tidalacceleration();
    let pot = moontide.tidalpotential(0.0);
    let bulgemoon = moontide.tidalbulgeheight();
    let bulgesun = suntide.tidalbulgeheight();
    let grav = moontide.gravitationalattraction();
    let spring = springtideamplitude();
    let neap = neaptideamplitude();
    let ratio = lunartosolartideratio();
    println!(
        "Tides: amoon={:.2e}  asun={:.2e}  ratio={:.2}",
        amoon, asun, ratio
    );
    println!(
        "  bulgemoon={:.3}m  bulgesun={:.3}m  spring={:.3}m  neap={:.3}m",
        bulgemoon, bulgesun, spring, neap
    );
    println!("  potential(0)={:.2e}  Fgrav={:.2e}N", pot, grav);

    let chix = chicxulubequivalent();
    let tung = tunguskaequivalent();
    let kechix = chix.kineticenergymt();
    let crater = chix.craterdiameterm(2700.0);
    let fireball = chix.fireballradiusm();
    let ejecta = chix.ejectavolumem3(2700.0);
    let vimp = chix.impactvelocity();
    let ketung = tung.kineticenergymt();
    println!(
        "Chicxulub: {:.2e}Mt  crater={:.0}m  fireball={:.0}m  ejecta={:.2e}m³  vimp={:.0}m/s",
        kechix, crater, fireball, ejecta, vimp
    );
    println!("Tunguska: {:.2e}Mt", ketung);

    let mut moon = MoonState::new();
    let source = match moon.source.get() {
        MoonSource::Binary => "binary",
        MoonSource::Simulation => "simulated",
    };
    let pos0 = moon.position();
    moon.step(3600.0);
    let pos1 = moon.position();
    let gmoon = moon.gravityat(EARTHMOONDISTANCE);
    println!(
        "Moon({}): pos0=({:.0},{:.0})  pos1h=({:.0},{:.0})  g@dist={:.4e}m/s²",
        source, pos0.0, pos0.1, pos1.0, pos1.1, gmoon
    );
    println!(
        "  M={:.3e}kg  R={:.4e}m  d={:.3e}m  P={:.0}s",
        LUNARMASS, LUNARRADIUS, EARTHMOONDISTANCE, LUNARORBITALPERIOD
    );

    let iss = ArtificialSatellite::leo("ISS", 420000.0, 408000.0);
    let geo = ArtificialSatellite::geo("GEO-Sat", 300.0);
    let custom = ArtificialSatellite::new("Molniya", 1500.0, 500000.0, 0.7, 1.1);
    let piss = iss.orbitalperiods();
    let viss = iss.orbitalvelocityms();
    let pgeo = geo.orbitalperiods();
    let gsurf = iss.gravityatsurface();
    let posiss = iss.position();
    let posgeo = geo.position();

    let orbittype = |sat: &ArtificialSatellite| -> &str {
        match sat.orbittype {
            OrbitType::LEO => "LEO",
            OrbitType::MEO => "MEO",
            OrbitType::GEO => "GEO",
            OrbitType::HEO => "HEO",
            OrbitType::Custom => "Custom",
        }
    };

    println!(
        "ISS({}): P={:.0}s  v={:.0}m/s  gsurf={:.2}  pos=({:.0},{:.0},{:.0})",
        orbittype(&iss),
        piss,
        viss,
        gsurf,
        posiss.0,
        posiss.1,
        posiss.2
    );
    println!(
        "GEO({}): P={:.0}s  pos=({:.0},{:.0},{:.0})",
        orbittype(&geo),
        pgeo,
        posgeo.0,
        posgeo.1,
        posgeo.2
    );
    println!(
        "Custom({}): e={:.1}  i={:.1}rad",
        orbittype(&custom),
        custom.eccentricity,
        custom.inclinationrad
    );

    let mut constellation = Constellation::new("Starlink-test");
    constellation.add(ArtificialSatellite::leo("S1", 260.0, 550000.0));
    constellation.add(ArtificialSatellite::leo("S2", 260.0, 550000.0));
    let mut satstepped = iss;
    satstepped.step(60.0);
    constellation.stepall(60.0);
    let positions = constellation.positions();
    println!(
        "Constellation: {} sats  positions={}",
        positions.len(),
        positions.len()
    );

    let dt = DateTime::new(2024, 6, 21, 12, 0, 0.0);
    let jd = dt.tojuliandate();
    let dtback = DateTime::fromjuliandate(jd);
    let unix = dt.tounixtimestamp();
    let dtunix = DateTime::fromunixtimestamp(unix);
    println!(
        "DateTime: {}-{}-{} {}:{}:{:.0}  JD={:.4}  unix={:.0}",
        dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second, jd, unix
    );
    println!(
        "  roundtrip JD: {}-{}-{}  roundtrip unix: {}-{}-{}",
        dtback.year, dtback.month, dtback.day, dtunix.year, dtunix.month, dtunix.day
    );
    println!(
        "  J2000JD={} UNIXJD={} SPD={}",
        J2000JD, UNIXEPOCHJD, SECONDSPERDAY
    );

    let mut epoch = Epoch::j2000();
    let epochmjd = Epoch::frommjd(51544.5);
    let centuries = epoch.centuriessincej2000();
    let days = epoch.dayssincej2000();
    let gmst = epoch.gmstdegrees();
    let mjd = epoch.tomjd();
    epoch.advancedays(365.25);
    let centuries1yr = epoch.centuriessincej2000();
    let mut epoch2 = Epoch::fromjd(J2000JD);
    epoch2.advanceseconds(86400.0);
    println!(
        "Epoch: T0={:.4}  days={:.1}  GMST={:.2}°  MJD={:.1}",
        centuries, days, gmst, mjd
    );
    println!(
        "  +1yr: T={:.6}  frommjd.jd={:.1}  epoch2.days={:.1}",
        centuries1yr,
        epochmjd.juliandate,
        epoch2.dayssincej2000()
    );
    println!(
        "  J2000={} J1950={} MJDOFF={}",
        J2000EPOCH, J1950EPOCH, MJDOFFSET
    );

    let mut ts = TimeScale::realtime();
    ts.step(1.0);
    let simdt = ts.simulationdt(1.0);
    let hours = ts.simhours();
    let daysts = ts.simdays();
    let years = ts.simyears();
    ts.pause();
    ts.resume();
    ts.togglepause();
    ts.setspeed(100.0);
    let mut ff = TimeScale::fastforward(10.0);
    ff.step(1.0);
    let sm = TimeScale::slowmotion(2.0);
    println!(
        "TimeScale: dt={:.1}  h={:.6}  d={:.8}  yr={:.10}  ffsim={:.1}s  smspeed={:.1}",
        simdt, hours, daysts, years, ff.simulationtimes, sm.speedmultiplier
    );

    let sunpos = SolarPosition::compute(jd, 48.8, 2.3);
    let above = sunpos.isabovehorizon();
    let distsun = sunpos.distancem();
    println!(
        "Solar: el={:.1}°  az={:.1}°  above={}  dist={:.3e}m  tilt={}°",
        sunpos.elevationdeg, sunpos.azimuthdeg, above, distsun, EARTHAXIALTILTDEG
    );

    let dnc = DayNightCycle::new(jd);
    let stateparis = dnc.stateat(48.8, 2.3);
    let statetext = match stateparis {
        DaylightState::Day => "day",
        DaylightState::Night => "night",
        DaylightState::CivilTwilight => "civiltwilight",
        DaylightState::NauticalTwilight => "nauticaltwilight",
        DaylightState::AstronomicalTwilight => "astrotwilight",
    };
    let terminator = dnc.terminatorpoints(36);
    let ambient = dnc.ambientlight(48.8, 2.3);
    println!(
        "DayNight: Paris={}  ambient={:.2}  terminatorpts={}",
        statetext,
        ambient,
        terminator.len()
    );

    let seasonstate = seasonat(jd, 48.8);
    let seasonname = match seasonstate.seasonnorth {
        Season::Spring => "spring",
        Season::Summer => "summer",
        Season::Autumn => "autumn",
        Season::Winter => "winter",
    };
    let southname = match seasonstate.seasonsouth {
        Season::Spring => "spring",
        Season::Summer => "summer",
        Season::Autumn => "autumn",
        Season::Winter => "winter",
    };
    let (sublat, sublon) = subsolarpoint(jd);
    println!(
        "Season: N={}  S={}  decl={:.2}°  dayh={:.2}  subsolar=({:.1},{:.1})",
        seasonname,
        southname,
        seasonstate.solardeclinationdeg,
        seasonstate.daylengthhours,
        sublat,
        sublon
    );
    println!(
        "  tilt={}°  tropyear={}d  vernaljd={}",
        SEASONSTILT, TROPICALYEARDAYS, VERNALEQUINOXJD
    );

    let paris = LatLon::new(48.8566, 2.3522, 35.0);
    let tokyo = LatLon::new(35.6762, 139.6503, 40.0);
    let ecefparis = paris.toecef();
    let cart = paris.tocartesiansimple();
    let distpt = paris.distanceto(&tokyo);
    let latlonback = ecefparis.tolatlon();
    println!(
        "Paris->ECEF: ({:.0},{:.0},{:.0})  cart=({:.0},{:.0},{:.0})",
        ecefparis.x, ecefparis.y, ecefparis.z, cart[0], cart[1], cart[2]
    );
    println!(
        "  dist(Paris-Tokyo)={:.0}m  roundtriplat={:.4}°",
        distpt, latlonback.latdeg
    );
    println!(
        "  f={:.9}  a={:.0}m  b={:.3}m",
        EARTHFLATTENING, EARTHSEMIMAJORM, EARTHSEMIMINORM
    );

    let regions = RegionDatabase::continents();
    let region = regions.pointinregion(48.8, 2.3);
    if let Some(r) = region {
        let rtype = match r.regiontype {
            RegionType::Continent => "continent",
            RegionType::Ocean => "ocean",
            RegionType::Sea => "sea",
            RegionType::Country => "country",
            RegionType::Island => "island",
        };
        println!(
            "Region(48.8,2.3): {} ({})  area={:.0}km²",
            r.name, rtype, r.areakm2
        );
    }

    let elevprov = ElevationProvider::global(30.0);
    let eleveverest = elevprov.sample(27.988, 86.925);
    let notables = ElevationProvider::notableelevations();
    println!(
        "Elevation(Everest)={:.0}m  notables={}",
        eleveverest,
        notables.len()
    );

    let bathy = BathymetryData::global(60.0);
    let depthmariana = bathy.sample(11.35, 142.2);
    let isocean = bathy.isocean(0.0, -30.0);
    let stats = BathymetryData::oceanstats();
    let basins = OceanBasin::majorbasins();
    println!(
        "Bathy(Mariana)={:.0}m  isocean(0,-30)={}  basins={}  avgdepth={:.0}m",
        depthmariana,
        isocean,
        basins.len(),
        stats.avgdepthm
    );

    let mut hmap = Heightmap::new(64);
    hmap.generateprocedural(42, 6, 0.5, 2.0);
    let hsample = hmap.sample(45.0, 10.0);
    let rat = hmap.radiusat(45.0, 10.0);
    println!(
        "Heightmap(64): sample(45,10)={:.1}m  radius={:.0}m",
        hsample, rat
    );

    let config = LodConfig::default();
    let mut lod = LodTerrain::new(config);
    lod.update([6371000.0 + 1000.0, 0.0, 0.0]);

    let mesh = TerrainMesh::fromregion(0.0, 10.0, 0.0, 10.0, 8, &|lat, lon| hmap.sample(lat, lon));
    println!(
        "Mesh: {} vertices  {} triangles",
        mesh.vertexcount(),
        mesh.trianglecount()
    );

    let classifier = BiomeClassifier::default();
    let biome = classifier.classify(2000.0, 45.0, 0.6);
    let biomename = match biome {
        Biome::Ocean => "ocean",
        Biome::Beach => "beach",
        Biome::Desert => "desert",
        Biome::Grassland => "grassland",
        Biome::Forest => "forest",
        Biome::Tundra => "tundra",
        Biome::Snow => "snow",
        Biome::Mountain => "mountain",
        Biome::Volcanic => "volcanic",
        Biome::Taiga => "taiga",
    };
    let splat = classifier.splat(2000.0, 45.0, 0.6);
    println!(
        "Biome(2000m,45°,0.6)={}  splattop={:.2}",
        biomename, splat.weights[0].1
    );

    let face = Face::PosZ;
    println!("LOD face={:?}", face);

    let all_mats = PbrMaterial::all_earth();
    let total_subsurface: f32 = all_mats.iter().map(|m| m.subsurface).sum();
    let avg_ior: f32 = all_mats.iter().map(|m| m.ior).sum::<f32>() / all_mats.len() as f32;
    let avg_normal: f32 =
        all_mats.iter().map(|m| m.normal_strength).sum::<f32>() / all_mats.len() as f32;
    let emissive_count = all_mats
        .iter()
        .filter(|m| m.emissive[0] > 0.0 || m.emissive[1] > 0.0 || m.emissive[2] > 0.0)
        .count();
    println!(
        "Materials({}) avg_ior={:.3} avg_normal={:.2} total_sss={:.2} emissive={}",
        all_mats.len(),
        avg_ior,
        avg_normal,
        total_subsurface,
        emissive_count
    );
    for m in &all_mats {
        let luminance = m.albedo[0] * 0.3 + m.albedo[1] * 0.59 + m.albedo[2] * 0.11;
        let fresnel_r0 = ((m.ior - 1.0) / (m.ior + 1.0)).powi(2);
        println!(
            "  {}: L={:.3} R={:.2} M={:.2} F₀={:.4} IOR={:.3} SSS={:.2} N={:.2} α={:.2} E=({:.2},{:.2},{:.2})",
            m.name,
            luminance,
            m.roughness,
            m.metallic,
            fresnel_r0,
            m.ior,
            m.subsurface,
            m.normal_strength,
            m.albedo[3],
            m.emissive[0],
            m.emissive[1],
            m.emissive[2]
        );
    }

    let shterrain = ShaderEndpoint::terrain();
    let shatm = ShaderEndpoint::atmosphere();
    let shocean = ShaderEndpoint::ocean();
    let shnights = ShaderEndpoint::night_lights();
    let total_uniform_bytes: usize = [&shterrain, &shatm, &shocean, &shnights]
        .iter()
        .map(|sh| {
            sh.name.len()
                + sh.uniforms
                    .iter()
                    .map(|u| match &u.value {
                        UniformValue::Float(_) => 4,
                        UniformValue::Vec3(_) => 12,
                        UniformValue::Vec4(_) => 16,
                        UniformValue::Int(_) => 4,
                    })
                    .sum::<usize>()
        })
        .sum();
    println!(
        "Shaders: terrain={} atm={} ocean={} night={} total_bytes={}",
        shterrain.uniforms.len(),
        shatm.uniforms.len(),
        shocean.uniforms.len(),
        shnights.uniforms.len(),
        total_uniform_bytes
    );
    for u in &shterrain.uniforms {
        let val = match &u.value {
            UniformValue::Float(f) => format!("f{:.2}", f),
            UniformValue::Vec3(v) => format!("v3({:.1},{:.1},{:.1})", v[0], v[1], v[2]),
            UniformValue::Vec4(v) => format!("v4({:.1},{:.1},{:.1},{:.1})", v[0], v[1], v[2], v[3]),
            UniformValue::Int(i) => format!("i{}", i),
        };
        println!("  {}={}", u.name, val);
    }

    let atm = AtmosphereEndpoint::earth();
    let scale_height_ratio = atm.rayleigh_scale_height_m / atm.mie_scale_height_m;
    let atmo_shell_vol = 4.0 / 3.0
        * std::f64::consts::PI
        * ((atm.planet_radius_m + atm.atmosphere_height_m).powi(3) - atm.planet_radius_m.powi(3));
    let air_density_sl = atm.sea_level_pressure_pa
        / (8.314_462_618 / atm.mean_molar_mass_kg_mol * atm.sea_level_temperature_k);
    let total_atmo_mass_kg = air_density_sl * atmo_shell_vol * 0.5;
    let total_rayleigh_scatter: f64 = atm
        .species
        .iter()
        .map(|s| {
            let n_s = atm.sea_level_number_density_m3 * s.volume_fraction;
            let king =
                (6.0 + 3.0 * s.depolarization_factor) / (6.0 - 7.0 * s.depolarization_factor);
            let dn = s.refractive_index_stp - 1.0;
            8.0 * std::f64::consts::PI.powi(3) * (2.0 * dn).powi(2) * king
                / (3.0 * n_s * (550e-9_f64).powi(4))
        })
        .sum();
    let ozone_optical_depth = atm.ozone_column_du
        * 1e-5
        * (-((25000.0 - atm.ozone_peak_altitude_m) / 5000.0).powi(2)).exp();
    let g = atm.mie_asymmetry_g;
    let hg_forward = (1.0 - g * g) / (1.0 + g * g - 2.0 * g).powf(1.5);
    println!(
        "Atmosphere: M_atm={:.3e}kg  H_ray/H_mie={:.2}  σ_ray550={:.4e}  O₃τ={:.4e}  HG_fwd={:.3}  S₀={:.0}W/m²",
        total_atmo_mass_kg,
        scale_height_ratio,
        total_rayleigh_scatter,
        ozone_optical_depth,
        hg_forward,
        atm.sun_irradiance_w_m2
    );
    for s in &atm.species {
        let king = (6.0 + 3.0 * s.depolarization_factor) / (6.0 - 7.0 * s.depolarization_factor);
        println!(
            "  {} ({}): f={:.6}  M={:.4e}kg/mol  n={:.6}  δ={:.3}  K={:.4}",
            s.name,
            s.symbol,
            s.volume_fraction,
            s.molar_mass_kg_mol,
            s.refractive_index_stp,
            s.depolarization_factor,
            king
        );
    }
    println!(
        "  β_rgb=({:.4e},{:.4e},{:.4e})  mie_coeff={:.2e}",
        atm.rayleigh_coefficients_rgb[0],
        atm.rayleigh_coefficients_rgb[1],
        atm.rayleigh_coefficients_rgb[2],
        atm.mie_coefficient
    );

    let atlantic = OceanEndpoint::earth_atlantic();
    let pacific = OceanEndpoint::earth_pacific();
    let arctic = OceanEndpoint::earth_arctic();
    for ocean in [&atlantic, &pacific, &arctic] {
        let fresnel_r0 = ((ocean.ior - 1.0) / (ocean.ior + 1.0)).powi(2);
        let fetch_m = ocean.fetch_km * 1000.0;
        let grav = ocean.surface_gravity_m_s2;
        let ws = ocean.wind_speed_m_s.max(0.1);
        let omega_p = grav / ws * (0.74 * (grav * fetch_m / (ws * ws)).powf(-0.33));
        let sig_wave_h = 4.0 * (8.1e-3 * grav / omega_p.powi(2)).sqrt();
        let foam_fraction = if ws > ocean.foam_threshold_wind_m_s {
            3.84e-6 * (ws - ocean.foam_threshold_wind_m_s).powf(3.41)
        } else {
            0.0
        };
        let dx = ocean.patch_size_m / ocean.grid_size as f64;
        let penetration_rgb = [
            1.0 / ocean.absorption_coefficients_rgb_m[0],
            1.0 / ocean.absorption_coefficients_rgb_m[1],
            1.0 / ocean.absorption_coefficients_rgb_m[2],
        ];
        let wind_bearing = ocean.wind_direction[1]
            .atan2(ocean.wind_direction[0])
            .to_degrees();
        println!(
            "Ocean(d={:.0}m,S={:.1}psu,T={:.1}K): R₀={:.4}  Hs={:.2}m  foam={:.4}  dx={:.1}m  pen_g={:.1}m  wind={:.0}°",
            ocean.mean_depth_m,
            ocean.salinity_psu,
            ocean.surface_temperature_k,
            fresnel_r0,
            sig_wave_h,
            foam_fraction,
            dx,
            penetration_rgb[1],
            wind_bearing
        );
    }

    let cumulus = CloudLayer::cumulus();
    let stratus = CloudLayer::stratus();
    let cirrus = CloudLayer::cirrus();
    let cb = CloudLayer::cumulonimbus();
    let cloudtype = |ct: &CloudType| match ct {
        CloudType::Cumulus => "Cu",
        CloudType::Stratus => "St",
        CloudType::Cirrus => "Ci",
        CloudType::Cumulonimbus => "Cb",
        CloudType::Altocumulus => "Ac",
        CloudType::Stratocumulus => "Sc",
        CloudType::Nimbostratus => "Ns",
    };
    for l in [&cumulus, &stratus, &cirrus, &cb] {
        let r_eff = l.droplet_radius_um * 1e-6;
        let lwc = 0.3 * l.density;
        let tau = 1.5 * 2.0 * lwc * l.thickness_m * l.coverage / (1000.0 * r_eff);
        let ice_corr = 1.0 - 0.15 * l.ice_fraction;
        let transmittance = (-(tau * ice_corr)).exp();
        let wind_dir = l.wind_direction[1].atan2(l.wind_direction[0]).to_degrees();
        println!(
            "  {} base={:.0}m thick={:.0}m cov={:.2} ρ={:.2} τ={:.2} T={:.4} r={:.0}µm ice={:.1} wind={:.1}m/s@{:.0}°",
            cloudtype(&l.cloud_type),
            l.base_altitude_m,
            l.thickness_m,
            l.coverage,
            l.density,
            tau * ice_corr,
            transmittance,
            l.droplet_radius_um,
            l.ice_fraction,
            l.wind_speed_m_s,
            wind_dir
        );
    }

    let csys = CloudSystemEndpoint::earth_default();
    let stormy = CloudSystemEndpoint::earth_stormy();
    let total_od_default: f64 = csys
        .layers
        .iter()
        .map(|l| {
            let r = l.droplet_radius_um * 1e-6;
            1.5 * 2.0 * 0.3 * l.density * l.thickness_m * l.coverage / (1000.0 * r)
                * (1.0 - 0.15 * l.ice_fraction)
        })
        .sum();
    let total_od_stormy: f64 = stormy
        .layers
        .iter()
        .map(|l| {
            let r = l.droplet_radius_um * 1e-6;
            1.5 * 2.0 * 0.3 * l.density * l.thickness_m * l.coverage / (1000.0 * r)
                * (1.0 - 0.15 * l.ice_fraction)
        })
        .sum();
    println!(
        "CloudSystem: default {} layers τ_total={:.2}  stormy {} layers τ_total={:.2}",
        csys.layers.len(),
        total_od_default,
        stormy.layers.len(),
        total_od_stormy
    );

    let rainforest = tropicalrainforest();
    let boreal = borealforest();
    let reef = coralreef();
    println!("Ecosystems:");
    for eco in [&rainforest, &boreal, &reef] {
        let sp = eco.totalspecies();
        let ind = eco.totalindividuals();
        let shannon = eco.shannonindex();
        let simpson = eco.simpsondiversity();
        let areasp = eco.expectedspeciesfromarea(0.25, 100.0);
        let npp = eco.totalnppgcyr();
        let turnover = eco.biomassturnovertimeyr(15000.0);
        println!(
            "  {} species  {} ind  H'={:.2}  D={:.4}  SAR={:.0}  NPP={:.0}  τ={:.1}yr",
            sp, ind, shannon, simpson, areasp, npp, turnover
        );
    }

    let tb = tropicalbroadleaf();
    let tg = temperategrassland();
    let cf = coniferousforest();
    let photo = tb.photosynthesisrate(400.0, 1500.0, 28.0);
    let canopy = tb.canopyphotosynthesis(photo);
    let transp = tg.transpirationmmday(1.5, 0.3);
    let nppveg = cf.nppkgcm2yr(15.0);
    let crt = cf.carbonresidencetimeyr(nppveg);
    let light = beerlambertlightextinction(2000.0, 4.0, 0.5);
    println!(
        "Vegetation: photo={:.2}  canopy={:.2}  transp={:.2}mm/d  NPP={:.3}  CRT={:.1}yr  light={:.1}",
        photo, canopy, transp, nppveg, crt, light
    );

    let elephant = africanelephant();
    let whale = bluewhale();
    let grel = elephant.growthrate();
    let proj = elephant.projectforward(10.0);
    let metab = whale.metabolicratew();
    let range = elephant.homerangekm2();
    let gentime = whale.generationtimeyears();
    let lifespan = elephant.maxlifespanyears();
    println!(
        "Elephant: r={:.3}  N(+10yr)={:.0}  range={:.1}km²  lifespan={:.0}yr",
        grel, proj, range, lifespan
    );
    println!("Whale: metab={:.0}W  gen={:.1}yr", metab, gentime);

    let mut pp = PredatorPrey {
        prey: africanelephant(),
        predator: bluewhale(),
        attackrate: 0.01,
        conversionefficiency: 0.1,
        predatordeathrate: 0.05,
    };
    let preyr = pp.preygrowthrate();
    let predr = pp.predatorgrowthrate();
    pp.step(0.1);
    println!(
        "Lotka-Volterra: preyr={:.1}  predr={:.1}  -> prey={:.1}  pred={:.2}",
        preyr, predr, pp.prey.count, pp.predator.count
    );
}

pub fn earth_arctic() -> Self

Examples found in repository

examples/tidal_sim.rs (line 556)

fn main() {
    let orbit = EarthOrbit::new();
    let perioddays = orbit.orbitalperioddays();
    let vperihelion = orbit.velocityatdistance(orbit.perihelionm());
    let vaphelion = orbit.velocityatdistance(orbit.aphelionm());
    let energy = orbit.specific_orbital_energy();
    let angmom = orbit.specific_angular_momentum();
    let vesc = EarthOrbit::escape_velocity_at_surface();
    let fsun = orbit.gravitational_force_sun();
    let r = orbit.current_radius();
    let vmean = orbit.mean_orbital_velocity();
    println!(
        "Orbit: P={:.2}d  vperi={:.0}m/s  vaph={:.0}m/s  vmean={:.0}m/s",
        perioddays, vperihelion, vaphelion, vmean
    );
    println!(
        "  E={:.2e}J/kg  L={:.2e}m²/s  vesc={:.0}m/s  Fsun={:.2e}N  r={:.3e}m",
        energy, angmom, vesc, fsun, r
    );
    println!(
        "  a={:.3e}m  e={}  i={}°  Ω={}°  ω={}°",
        SEMIMAJORAXIS,
        ECCENTRICITY,
        INCLINATIONDEG,
        LONGITUDEASCENDINGNODEDEG,
        ARGUMENTPERIHELIONDEG
    );

    let rot = EarthRotation::new();
    let vequator = rot.surfacevelocityatlatitude(0.0);
    let v45 = rot.surfacevelocityatlatitude(45.0);
    let accel = rot.centripetalaccelerationatlatitude(0.0);
    let coriolis = rot.coriolisparameter(45.0);
    let moi = rot.momentofinertia();
    let ke = rot.rotationalkineticenergy();
    let lrot = rot.angular_momentum();
    let prec = rot.precession_rate_rad_per_year();
    let daylensummer = rot.day_length_variation_due_to_tilt(172, 60.0);
    let daylenwinter = rot.day_length_variation_due_to_tilt(355, 60.0);
    println!(
        "Rotation: veq={:.0}m/s  v45={:.0}m/s  ac={:.4}m/s²  f45={:.5e}",
        vequator, v45, accel, coriolis
    );
    println!(
        "  I={:.3e}kg·m²  KE={:.3e}J  L={:.3e}kg·m²/s  prec={:.2e}rad/yr",
        moi, ke, lrot, prec
    );
    println!(
        "  day@60°N: summer={:.1}h  winter={:.1}h",
        daylensummer, daylenwinter
    );
    println!(
        "  sidereal={:.4}s  solar={}s  tilt={}°={:.4}rad  precperiod={}yr",
        SIDEREALDAYS, SOLARDAYS, AXIALTILTDEG, AXIALTILTRAD, PRECESSIONPERIODYEARS
    );

    let moontide = TidalForce::frommoon();
    let suntide = TidalForce::fromsun();
    let amoon = moontide.tidalacceleration();
    let asun = suntide.tidalacceleration();
    let pot = moontide.tidalpotential(0.0);
    let bulgemoon = moontide.tidalbulgeheight();
    let bulgesun = suntide.tidalbulgeheight();
    let grav = moontide.gravitationalattraction();
    let spring = springtideamplitude();
    let neap = neaptideamplitude();
    let ratio = lunartosolartideratio();
    println!(
        "Tides: amoon={:.2e}  asun={:.2e}  ratio={:.2}",
        amoon, asun, ratio
    );
    println!(
        "  bulgemoon={:.3}m  bulgesun={:.3}m  spring={:.3}m  neap={:.3}m",
        bulgemoon, bulgesun, spring, neap
    );
    println!("  potential(0)={:.2e}  Fgrav={:.2e}N", pot, grav);

    let chix = chicxulubequivalent();
    let tung = tunguskaequivalent();
    let kechix = chix.kineticenergymt();
    let crater = chix.craterdiameterm(2700.0);
    let fireball = chix.fireballradiusm();
    let ejecta = chix.ejectavolumem3(2700.0);
    let vimp = chix.impactvelocity();
    let ketung = tung.kineticenergymt();
    println!(
        "Chicxulub: {:.2e}Mt  crater={:.0}m  fireball={:.0}m  ejecta={:.2e}m³  vimp={:.0}m/s",
        kechix, crater, fireball, ejecta, vimp
    );
    println!("Tunguska: {:.2e}Mt", ketung);

    let mut moon = MoonState::new();
    let source = match moon.source.get() {
        MoonSource::Binary => "binary",
        MoonSource::Simulation => "simulated",
    };
    let pos0 = moon.position();
    moon.step(3600.0);
    let pos1 = moon.position();
    let gmoon = moon.gravityat(EARTHMOONDISTANCE);
    println!(
        "Moon({}): pos0=({:.0},{:.0})  pos1h=({:.0},{:.0})  g@dist={:.4e}m/s²",
        source, pos0.0, pos0.1, pos1.0, pos1.1, gmoon
    );
    println!(
        "  M={:.3e}kg  R={:.4e}m  d={:.3e}m  P={:.0}s",
        LUNARMASS, LUNARRADIUS, EARTHMOONDISTANCE, LUNARORBITALPERIOD
    );

    let iss = ArtificialSatellite::leo("ISS", 420000.0, 408000.0);
    let geo = ArtificialSatellite::geo("GEO-Sat", 300.0);
    let custom = ArtificialSatellite::new("Molniya", 1500.0, 500000.0, 0.7, 1.1);
    let piss = iss.orbitalperiods();
    let viss = iss.orbitalvelocityms();
    let pgeo = geo.orbitalperiods();
    let gsurf = iss.gravityatsurface();
    let posiss = iss.position();
    let posgeo = geo.position();

    let orbittype = |sat: &ArtificialSatellite| -> &str {
        match sat.orbittype {
            OrbitType::LEO => "LEO",
            OrbitType::MEO => "MEO",
            OrbitType::GEO => "GEO",
            OrbitType::HEO => "HEO",
            OrbitType::Custom => "Custom",
        }
    };

    println!(
        "ISS({}): P={:.0}s  v={:.0}m/s  gsurf={:.2}  pos=({:.0},{:.0},{:.0})",
        orbittype(&iss),
        piss,
        viss,
        gsurf,
        posiss.0,
        posiss.1,
        posiss.2
    );
    println!(
        "GEO({}): P={:.0}s  pos=({:.0},{:.0},{:.0})",
        orbittype(&geo),
        pgeo,
        posgeo.0,
        posgeo.1,
        posgeo.2
    );
    println!(
        "Custom({}): e={:.1}  i={:.1}rad",
        orbittype(&custom),
        custom.eccentricity,
        custom.inclinationrad
    );

    let mut constellation = Constellation::new("Starlink-test");
    constellation.add(ArtificialSatellite::leo("S1", 260.0, 550000.0));
    constellation.add(ArtificialSatellite::leo("S2", 260.0, 550000.0));
    let mut satstepped = iss;
    satstepped.step(60.0);
    constellation.stepall(60.0);
    let positions = constellation.positions();
    println!(
        "Constellation: {} sats  positions={}",
        positions.len(),
        positions.len()
    );

    let dt = DateTime::new(2024, 6, 21, 12, 0, 0.0);
    let jd = dt.tojuliandate();
    let dtback = DateTime::fromjuliandate(jd);
    let unix = dt.tounixtimestamp();
    let dtunix = DateTime::fromunixtimestamp(unix);
    println!(
        "DateTime: {}-{}-{} {}:{}:{:.0}  JD={:.4}  unix={:.0}",
        dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second, jd, unix
    );
    println!(
        "  roundtrip JD: {}-{}-{}  roundtrip unix: {}-{}-{}",
        dtback.year, dtback.month, dtback.day, dtunix.year, dtunix.month, dtunix.day
    );
    println!(
        "  J2000JD={} UNIXJD={} SPD={}",
        J2000JD, UNIXEPOCHJD, SECONDSPERDAY
    );

    let mut epoch = Epoch::j2000();
    let epochmjd = Epoch::frommjd(51544.5);
    let centuries = epoch.centuriessincej2000();
    let days = epoch.dayssincej2000();
    let gmst = epoch.gmstdegrees();
    let mjd = epoch.tomjd();
    epoch.advancedays(365.25);
    let centuries1yr = epoch.centuriessincej2000();
    let mut epoch2 = Epoch::fromjd(J2000JD);
    epoch2.advanceseconds(86400.0);
    println!(
        "Epoch: T0={:.4}  days={:.1}  GMST={:.2}°  MJD={:.1}",
        centuries, days, gmst, mjd
    );
    println!(
        "  +1yr: T={:.6}  frommjd.jd={:.1}  epoch2.days={:.1}",
        centuries1yr,
        epochmjd.juliandate,
        epoch2.dayssincej2000()
    );
    println!(
        "  J2000={} J1950={} MJDOFF={}",
        J2000EPOCH, J1950EPOCH, MJDOFFSET
    );

    let mut ts = TimeScale::realtime();
    ts.step(1.0);
    let simdt = ts.simulationdt(1.0);
    let hours = ts.simhours();
    let daysts = ts.simdays();
    let years = ts.simyears();
    ts.pause();
    ts.resume();
    ts.togglepause();
    ts.setspeed(100.0);
    let mut ff = TimeScale::fastforward(10.0);
    ff.step(1.0);
    let sm = TimeScale::slowmotion(2.0);
    println!(
        "TimeScale: dt={:.1}  h={:.6}  d={:.8}  yr={:.10}  ffsim={:.1}s  smspeed={:.1}",
        simdt, hours, daysts, years, ff.simulationtimes, sm.speedmultiplier
    );

    let sunpos = SolarPosition::compute(jd, 48.8, 2.3);
    let above = sunpos.isabovehorizon();
    let distsun = sunpos.distancem();
    println!(
        "Solar: el={:.1}°  az={:.1}°  above={}  dist={:.3e}m  tilt={}°",
        sunpos.elevationdeg, sunpos.azimuthdeg, above, distsun, EARTHAXIALTILTDEG
    );

    let dnc = DayNightCycle::new(jd);
    let stateparis = dnc.stateat(48.8, 2.3);
    let statetext = match stateparis {
        DaylightState::Day => "day",
        DaylightState::Night => "night",
        DaylightState::CivilTwilight => "civiltwilight",
        DaylightState::NauticalTwilight => "nauticaltwilight",
        DaylightState::AstronomicalTwilight => "astrotwilight",
    };
    let terminator = dnc.terminatorpoints(36);
    let ambient = dnc.ambientlight(48.8, 2.3);
    println!(
        "DayNight: Paris={}  ambient={:.2}  terminatorpts={}",
        statetext,
        ambient,
        terminator.len()
    );

    let seasonstate = seasonat(jd, 48.8);
    let seasonname = match seasonstate.seasonnorth {
        Season::Spring => "spring",
        Season::Summer => "summer",
        Season::Autumn => "autumn",
        Season::Winter => "winter",
    };
    let southname = match seasonstate.seasonsouth {
        Season::Spring => "spring",
        Season::Summer => "summer",
        Season::Autumn => "autumn",
        Season::Winter => "winter",
    };
    let (sublat, sublon) = subsolarpoint(jd);
    println!(
        "Season: N={}  S={}  decl={:.2}°  dayh={:.2}  subsolar=({:.1},{:.1})",
        seasonname,
        southname,
        seasonstate.solardeclinationdeg,
        seasonstate.daylengthhours,
        sublat,
        sublon
    );
    println!(
        "  tilt={}°  tropyear={}d  vernaljd={}",
        SEASONSTILT, TROPICALYEARDAYS, VERNALEQUINOXJD
    );

    let paris = LatLon::new(48.8566, 2.3522, 35.0);
    let tokyo = LatLon::new(35.6762, 139.6503, 40.0);
    let ecefparis = paris.toecef();
    let cart = paris.tocartesiansimple();
    let distpt = paris.distanceto(&tokyo);
    let latlonback = ecefparis.tolatlon();
    println!(
        "Paris->ECEF: ({:.0},{:.0},{:.0})  cart=({:.0},{:.0},{:.0})",
        ecefparis.x, ecefparis.y, ecefparis.z, cart[0], cart[1], cart[2]
    );
    println!(
        "  dist(Paris-Tokyo)={:.0}m  roundtriplat={:.4}°",
        distpt, latlonback.latdeg
    );
    println!(
        "  f={:.9}  a={:.0}m  b={:.3}m",
        EARTHFLATTENING, EARTHSEMIMAJORM, EARTHSEMIMINORM
    );

    let regions = RegionDatabase::continents();
    let region = regions.pointinregion(48.8, 2.3);
    if let Some(r) = region {
        let rtype = match r.regiontype {
            RegionType::Continent => "continent",
            RegionType::Ocean => "ocean",
            RegionType::Sea => "sea",
            RegionType::Country => "country",
            RegionType::Island => "island",
        };
        println!(
            "Region(48.8,2.3): {} ({})  area={:.0}km²",
            r.name, rtype, r.areakm2
        );
    }

    let elevprov = ElevationProvider::global(30.0);
    let eleveverest = elevprov.sample(27.988, 86.925);
    let notables = ElevationProvider::notableelevations();
    println!(
        "Elevation(Everest)={:.0}m  notables={}",
        eleveverest,
        notables.len()
    );

    let bathy = BathymetryData::global(60.0);
    let depthmariana = bathy.sample(11.35, 142.2);
    let isocean = bathy.isocean(0.0, -30.0);
    let stats = BathymetryData::oceanstats();
    let basins = OceanBasin::majorbasins();
    println!(
        "Bathy(Mariana)={:.0}m  isocean(0,-30)={}  basins={}  avgdepth={:.0}m",
        depthmariana,
        isocean,
        basins.len(),
        stats.avgdepthm
    );

    let mut hmap = Heightmap::new(64);
    hmap.generateprocedural(42, 6, 0.5, 2.0);
    let hsample = hmap.sample(45.0, 10.0);
    let rat = hmap.radiusat(45.0, 10.0);
    println!(
        "Heightmap(64): sample(45,10)={:.1}m  radius={:.0}m",
        hsample, rat
    );

    let config = LodConfig::default();
    let mut lod = LodTerrain::new(config);
    lod.update([6371000.0 + 1000.0, 0.0, 0.0]);

    let mesh = TerrainMesh::fromregion(0.0, 10.0, 0.0, 10.0, 8, &|lat, lon| hmap.sample(lat, lon));
    println!(
        "Mesh: {} vertices  {} triangles",
        mesh.vertexcount(),
        mesh.trianglecount()
    );

    let classifier = BiomeClassifier::default();
    let biome = classifier.classify(2000.0, 45.0, 0.6);
    let biomename = match biome {
        Biome::Ocean => "ocean",
        Biome::Beach => "beach",
        Biome::Desert => "desert",
        Biome::Grassland => "grassland",
        Biome::Forest => "forest",
        Biome::Tundra => "tundra",
        Biome::Snow => "snow",
        Biome::Mountain => "mountain",
        Biome::Volcanic => "volcanic",
        Biome::Taiga => "taiga",
    };
    let splat = classifier.splat(2000.0, 45.0, 0.6);
    println!(
        "Biome(2000m,45°,0.6)={}  splattop={:.2}",
        biomename, splat.weights[0].1
    );

    let face = Face::PosZ;
    println!("LOD face={:?}", face);

    let all_mats = PbrMaterial::all_earth();
    let total_subsurface: f32 = all_mats.iter().map(|m| m.subsurface).sum();
    let avg_ior: f32 = all_mats.iter().map(|m| m.ior).sum::<f32>() / all_mats.len() as f32;
    let avg_normal: f32 =
        all_mats.iter().map(|m| m.normal_strength).sum::<f32>() / all_mats.len() as f32;
    let emissive_count = all_mats
        .iter()
        .filter(|m| m.emissive[0] > 0.0 || m.emissive[1] > 0.0 || m.emissive[2] > 0.0)
        .count();
    println!(
        "Materials({}) avg_ior={:.3} avg_normal={:.2} total_sss={:.2} emissive={}",
        all_mats.len(),
        avg_ior,
        avg_normal,
        total_subsurface,
        emissive_count
    );
    for m in &all_mats {
        let luminance = m.albedo[0] * 0.3 + m.albedo[1] * 0.59 + m.albedo[2] * 0.11;
        let fresnel_r0 = ((m.ior - 1.0) / (m.ior + 1.0)).powi(2);
        println!(
            "  {}: L={:.3} R={:.2} M={:.2} F₀={:.4} IOR={:.3} SSS={:.2} N={:.2} α={:.2} E=({:.2},{:.2},{:.2})",
            m.name,
            luminance,
            m.roughness,
            m.metallic,
            fresnel_r0,
            m.ior,
            m.subsurface,
            m.normal_strength,
            m.albedo[3],
            m.emissive[0],
            m.emissive[1],
            m.emissive[2]
        );
    }

    let shterrain = ShaderEndpoint::terrain();
    let shatm = ShaderEndpoint::atmosphere();
    let shocean = ShaderEndpoint::ocean();
    let shnights = ShaderEndpoint::night_lights();
    let total_uniform_bytes: usize = [&shterrain, &shatm, &shocean, &shnights]
        .iter()
        .map(|sh| {
            sh.name.len()
                + sh.uniforms
                    .iter()
                    .map(|u| match &u.value {
                        UniformValue::Float(_) => 4,
                        UniformValue::Vec3(_) => 12,
                        UniformValue::Vec4(_) => 16,
                        UniformValue::Int(_) => 4,
                    })
                    .sum::<usize>()
        })
        .sum();
    println!(
        "Shaders: terrain={} atm={} ocean={} night={} total_bytes={}",
        shterrain.uniforms.len(),
        shatm.uniforms.len(),
        shocean.uniforms.len(),
        shnights.uniforms.len(),
        total_uniform_bytes
    );
    for u in &shterrain.uniforms {
        let val = match &u.value {
            UniformValue::Float(f) => format!("f{:.2}", f),
            UniformValue::Vec3(v) => format!("v3({:.1},{:.1},{:.1})", v[0], v[1], v[2]),
            UniformValue::Vec4(v) => format!("v4({:.1},{:.1},{:.1},{:.1})", v[0], v[1], v[2], v[3]),
            UniformValue::Int(i) => format!("i{}", i),
        };
        println!("  {}={}", u.name, val);
    }

    let atm = AtmosphereEndpoint::earth();
    let scale_height_ratio = atm.rayleigh_scale_height_m / atm.mie_scale_height_m;
    let atmo_shell_vol = 4.0 / 3.0
        * std::f64::consts::PI
        * ((atm.planet_radius_m + atm.atmosphere_height_m).powi(3) - atm.planet_radius_m.powi(3));
    let air_density_sl = atm.sea_level_pressure_pa
        / (8.314_462_618 / atm.mean_molar_mass_kg_mol * atm.sea_level_temperature_k);
    let total_atmo_mass_kg = air_density_sl * atmo_shell_vol * 0.5;
    let total_rayleigh_scatter: f64 = atm
        .species
        .iter()
        .map(|s| {
            let n_s = atm.sea_level_number_density_m3 * s.volume_fraction;
            let king =
                (6.0 + 3.0 * s.depolarization_factor) / (6.0 - 7.0 * s.depolarization_factor);
            let dn = s.refractive_index_stp - 1.0;
            8.0 * std::f64::consts::PI.powi(3) * (2.0 * dn).powi(2) * king
                / (3.0 * n_s * (550e-9_f64).powi(4))
        })
        .sum();
    let ozone_optical_depth = atm.ozone_column_du
        * 1e-5
        * (-((25000.0 - atm.ozone_peak_altitude_m) / 5000.0).powi(2)).exp();
    let g = atm.mie_asymmetry_g;
    let hg_forward = (1.0 - g * g) / (1.0 + g * g - 2.0 * g).powf(1.5);
    println!(
        "Atmosphere: M_atm={:.3e}kg  H_ray/H_mie={:.2}  σ_ray550={:.4e}  O₃τ={:.4e}  HG_fwd={:.3}  S₀={:.0}W/m²",
        total_atmo_mass_kg,
        scale_height_ratio,
        total_rayleigh_scatter,
        ozone_optical_depth,
        hg_forward,
        atm.sun_irradiance_w_m2
    );
    for s in &atm.species {
        let king = (6.0 + 3.0 * s.depolarization_factor) / (6.0 - 7.0 * s.depolarization_factor);
        println!(
            "  {} ({}): f={:.6}  M={:.4e}kg/mol  n={:.6}  δ={:.3}  K={:.4}",
            s.name,
            s.symbol,
            s.volume_fraction,
            s.molar_mass_kg_mol,
            s.refractive_index_stp,
            s.depolarization_factor,
            king
        );
    }
    println!(
        "  β_rgb=({:.4e},{:.4e},{:.4e})  mie_coeff={:.2e}",
        atm.rayleigh_coefficients_rgb[0],
        atm.rayleigh_coefficients_rgb[1],
        atm.rayleigh_coefficients_rgb[2],
        atm.mie_coefficient
    );

    let atlantic = OceanEndpoint::earth_atlantic();
    let pacific = OceanEndpoint::earth_pacific();
    let arctic = OceanEndpoint::earth_arctic();
    for ocean in [&atlantic, &pacific, &arctic] {
        let fresnel_r0 = ((ocean.ior - 1.0) / (ocean.ior + 1.0)).powi(2);
        let fetch_m = ocean.fetch_km * 1000.0;
        let grav = ocean.surface_gravity_m_s2;
        let ws = ocean.wind_speed_m_s.max(0.1);
        let omega_p = grav / ws * (0.74 * (grav * fetch_m / (ws * ws)).powf(-0.33));
        let sig_wave_h = 4.0 * (8.1e-3 * grav / omega_p.powi(2)).sqrt();
        let foam_fraction = if ws > ocean.foam_threshold_wind_m_s {
            3.84e-6 * (ws - ocean.foam_threshold_wind_m_s).powf(3.41)
        } else {
            0.0
        };
        let dx = ocean.patch_size_m / ocean.grid_size as f64;
        let penetration_rgb = [
            1.0 / ocean.absorption_coefficients_rgb_m[0],
            1.0 / ocean.absorption_coefficients_rgb_m[1],
            1.0 / ocean.absorption_coefficients_rgb_m[2],
        ];
        let wind_bearing = ocean.wind_direction[1]
            .atan2(ocean.wind_direction[0])
            .to_degrees();
        println!(
            "Ocean(d={:.0}m,S={:.1}psu,T={:.1}K): R₀={:.4}  Hs={:.2}m  foam={:.4}  dx={:.1}m  pen_g={:.1}m  wind={:.0}°",
            ocean.mean_depth_m,
            ocean.salinity_psu,
            ocean.surface_temperature_k,
            fresnel_r0,
            sig_wave_h,
            foam_fraction,
            dx,
            penetration_rgb[1],
            wind_bearing
        );
    }

    let cumulus = CloudLayer::cumulus();
    let stratus = CloudLayer::stratus();
    let cirrus = CloudLayer::cirrus();
    let cb = CloudLayer::cumulonimbus();
    let cloudtype = |ct: &CloudType| match ct {
        CloudType::Cumulus => "Cu",
        CloudType::Stratus => "St",
        CloudType::Cirrus => "Ci",
        CloudType::Cumulonimbus => "Cb",
        CloudType::Altocumulus => "Ac",
        CloudType::Stratocumulus => "Sc",
        CloudType::Nimbostratus => "Ns",
    };
    for l in [&cumulus, &stratus, &cirrus, &cb] {
        let r_eff = l.droplet_radius_um * 1e-6;
        let lwc = 0.3 * l.density;
        let tau = 1.5 * 2.0 * lwc * l.thickness_m * l.coverage / (1000.0 * r_eff);
        let ice_corr = 1.0 - 0.15 * l.ice_fraction;
        let transmittance = (-(tau * ice_corr)).exp();
        let wind_dir = l.wind_direction[1].atan2(l.wind_direction[0]).to_degrees();
        println!(
            "  {} base={:.0}m thick={:.0}m cov={:.2} ρ={:.2} τ={:.2} T={:.4} r={:.0}µm ice={:.1} wind={:.1}m/s@{:.0}°",
            cloudtype(&l.cloud_type),
            l.base_altitude_m,
            l.thickness_m,
            l.coverage,
            l.density,
            tau * ice_corr,
            transmittance,
            l.droplet_radius_um,
            l.ice_fraction,
            l.wind_speed_m_s,
            wind_dir
        );
    }

    let csys = CloudSystemEndpoint::earth_default();
    let stormy = CloudSystemEndpoint::earth_stormy();
    let total_od_default: f64 = csys
        .layers
        .iter()
        .map(|l| {
            let r = l.droplet_radius_um * 1e-6;
            1.5 * 2.0 * 0.3 * l.density * l.thickness_m * l.coverage / (1000.0 * r)
                * (1.0 - 0.15 * l.ice_fraction)
        })
        .sum();
    let total_od_stormy: f64 = stormy
        .layers
        .iter()
        .map(|l| {
            let r = l.droplet_radius_um * 1e-6;
            1.5 * 2.0 * 0.3 * l.density * l.thickness_m * l.coverage / (1000.0 * r)
                * (1.0 - 0.15 * l.ice_fraction)
        })
        .sum();
    println!(
        "CloudSystem: default {} layers τ_total={:.2}  stormy {} layers τ_total={:.2}",
        csys.layers.len(),
        total_od_default,
        stormy.layers.len(),
        total_od_stormy
    );

    let rainforest = tropicalrainforest();
    let boreal = borealforest();
    let reef = coralreef();
    println!("Ecosystems:");
    for eco in [&rainforest, &boreal, &reef] {
        let sp = eco.totalspecies();
        let ind = eco.totalindividuals();
        let shannon = eco.shannonindex();
        let simpson = eco.simpsondiversity();
        let areasp = eco.expectedspeciesfromarea(0.25, 100.0);
        let npp = eco.totalnppgcyr();
        let turnover = eco.biomassturnovertimeyr(15000.0);
        println!(
            "  {} species  {} ind  H'={:.2}  D={:.4}  SAR={:.0}  NPP={:.0}  τ={:.1}yr",
            sp, ind, shannon, simpson, areasp, npp, turnover
        );
    }

    let tb = tropicalbroadleaf();
    let tg = temperategrassland();
    let cf = coniferousforest();
    let photo = tb.photosynthesisrate(400.0, 1500.0, 28.0);
    let canopy = tb.canopyphotosynthesis(photo);
    let transp = tg.transpirationmmday(1.5, 0.3);
    let nppveg = cf.nppkgcm2yr(15.0);
    let crt = cf.carbonresidencetimeyr(nppveg);
    let light = beerlambertlightextinction(2000.0, 4.0, 0.5);
    println!(
        "Vegetation: photo={:.2}  canopy={:.2}  transp={:.2}mm/d  NPP={:.3}  CRT={:.1}yr  light={:.1}",
        photo, canopy, transp, nppveg, crt, light
    );

    let elephant = africanelephant();
    let whale = bluewhale();
    let grel = elephant.growthrate();
    let proj = elephant.projectforward(10.0);
    let metab = whale.metabolicratew();
    let range = elephant.homerangekm2();
    let gentime = whale.generationtimeyears();
    let lifespan = elephant.maxlifespanyears();
    println!(
        "Elephant: r={:.3}  N(+10yr)={:.0}  range={:.1}km²  lifespan={:.0}yr",
        grel, proj, range, lifespan
    );
    println!("Whale: metab={:.0}W  gen={:.1}yr", metab, gentime);

    let mut pp = PredatorPrey {
        prey: africanelephant(),
        predator: bluewhale(),
        attackrate: 0.01,
        conversionefficiency: 0.1,
        predatordeathrate: 0.05,
    };
    let preyr = pp.preygrowthrate();
    let predr = pp.predatorgrowthrate();
    pp.step(0.1);
    println!(
        "Lotka-Volterra: preyr={:.1}  predr={:.1}  -> prey={:.1}  pred={:.2}",
        preyr, predr, pp.prey.count, pp.predator.count
    );
}

pub fn fresnel_r0(&self) -> f64

pub fn wave_amplitude(&self) -> f64

pub fn depth_color(&self, depth: f64) -> [f64; 3]