venuss 0.0.2

Venus celestial simulation crate for the MilkyWay SolarSystem workspace
Documentation
use venuss::{
    GEOTHERMAL_GRADIENT, GLOBAL_STORM_OPTICAL_DEPTH, MEAN_CLOUD_OPTICAL_DEPTH, VENUS_RADIUS,
};
use venuss::{
    atmosphere, biosphere, geodata, geology, hydrology, lighting, physics, rendering, satellites,
    temporal, terrain,
};

const TARGET_FPS: f64 = 60.0;
const REAL_DT_S: f64 = 1.0 / TARGET_FPS;
const FRAME_COUNT: usize = 20_000;

fn biome_material(biome: terrain::texturing::VenusBiome) -> rendering::materials::PbrMaterial {
    match biome {
        terrain::texturing::VenusBiome::Highland | terrain::texturing::VenusBiome::Tessera => {
            rendering::materials::tessera_highland()
        }
        terrain::texturing::VenusBiome::Volcanic => rendering::materials::volcanic_flow(),
        terrain::texturing::VenusBiome::Plains => rendering::materials::weathered_plain(),
    }
}

struct VenusDiagnostics {
    ls_deg: f64,
    local_time_h: f64,
    solar_irradiance_w_m2: f64,
    direct_transmission: f64,
    cloud_optical_depth: f64,
    daylight_state: lighting::day_night::DaylightState,
    season: lighting::seasons::Season,
    sky_color: [f64; 3],
    cloud_top_pressure_pa: f64,
    surface_pressure_pa: f64,
    superrotation_m_s: f64,
    surface_wind_m_s: f64,
    lightning_rate_s: f64,
    solar_tide_bulge_m: f64,
    heat_flow_w_m2: f64,
    rayleigh_number: f64,
    thermal_stress: f64,
    erosion_flux: f64,
    chemical_weathering: f64,
    volcanic_coverage: f64,
    volcanic_volume_km3: f64,
    cloud_habitability: f64,
    biome: terrain::texturing::VenusBiome,
    material_roughness: f32,
    observer_elevation_m: f64,
    orbiters_mean_period_h: f64,
    paleoocean_volume_km3: f64,
    paleolake_volume_km3: f64,
    lava_channel_length_km: f64,
    checksum: f64,
}

impl Default for VenusDiagnostics {
    fn default() -> Self {
        Self {
            ls_deg: 0.0,
            local_time_h: 0.0,
            solar_irradiance_w_m2: 0.0,
            direct_transmission: 0.0,
            cloud_optical_depth: MEAN_CLOUD_OPTICAL_DEPTH,
            daylight_state: lighting::day_night::DaylightState::Night,
            season: lighting::seasons::Season::Spring,
            sky_color: [0.0; 3],
            cloud_top_pressure_pa: 0.0,
            surface_pressure_pa: 0.0,
            superrotation_m_s: 0.0,
            surface_wind_m_s: 0.0,
            lightning_rate_s: 0.0,
            solar_tide_bulge_m: 0.0,
            heat_flow_w_m2: 0.0,
            rayleigh_number: 0.0,
            thermal_stress: 0.0,
            erosion_flux: 0.0,
            chemical_weathering: 0.0,
            volcanic_coverage: 0.0,
            volcanic_volume_km3: 0.0,
            cloud_habitability: 0.0,
            biome: terrain::texturing::VenusBiome::Plains,
            material_roughness: rendering::materials::weathered_plain().roughness,
            observer_elevation_m: 0.0,
            orbiters_mean_period_h: 0.0,
            paleoocean_volume_km3: 0.0,
            paleolake_volume_km3: 0.0,
            lava_channel_length_km: 0.0,
            checksum: 0.0,
        }
    }
}

impl VenusDiagnostics {
    fn recompute_checksum(&mut self) {
        let daylight = match self.daylight_state {
            lighting::day_night::DaylightState::Day => 1.0,
            lighting::day_night::DaylightState::Twilight => 0.5,
            lighting::day_night::DaylightState::Night => 0.0,
        };
        let season = match self.season {
            lighting::seasons::Season::Spring => 0.1,
            lighting::seasons::Season::Summer => 0.2,
            lighting::seasons::Season::Autumn => 0.3,
            lighting::seasons::Season::Winter => 0.4,
        };
        let biome = match self.biome {
            terrain::texturing::VenusBiome::Highland => 1.0,
            terrain::texturing::VenusBiome::Plains => 2.0,
            terrain::texturing::VenusBiome::Volcanic => 3.0,
            terrain::texturing::VenusBiome::Tessera => 4.0,
        };
        self.checksum = self.ls_deg
            + self.local_time_h
            + self.solar_irradiance_w_m2 * 1e-3
            + self.direct_transmission
            + self.cloud_optical_depth * 1e-2
            + daylight
            + season
            + self.sky_color.iter().sum::<f64>()
            + self.cloud_top_pressure_pa * 1e-7
            + self.surface_pressure_pa * 1e-7
            + self.superrotation_m_s * 1e-2
            + self.surface_wind_m_s * 1e-2
            + self.lightning_rate_s * 1e-2
            + self.solar_tide_bulge_m * 1e6
            + self.heat_flow_w_m2
            + self.rayleigh_number * 1e-8
            + self.thermal_stress * 1e-6
            + self.erosion_flux * 1e-4
            + self.chemical_weathering
            + self.volcanic_coverage
            + self.volcanic_volume_km3 * 1e-5
            + self.cloud_habitability
            + biome
            + self.material_roughness as f64
            + self.observer_elevation_m * 1e-4
            + self.orbiters_mean_period_h * 1e-2
            + self.paleoocean_volume_km3 * 1e-9
            + self.paleolake_volume_km3 * 1e-4
            + self.lava_channel_length_km * 1e-3;
    }
}

struct VenusSimulation {
    epoch: temporal::epoch::VenusEpoch,
    time_scale: temporal::time_scale::TimeScale,
    orbit: physics::orbit::VenusOrbit,
    rotation: physics::rotation::VenusRotation,
    rotation_angle_rad: f64,
    climate: atmosphere::climate::VenusClimateState,
    heightmap: terrain::heightmap::Heightmap,
    lod: terrain::lod::LodTerrain,
    atmo_render: rendering::atmosphere_scattering::VenusAtmosphereParams,
    biome_classifier: terrain::texturing::BiomeClassifier,
    terrain_shader: rendering::shaders::ShaderData,
    observer: geodata::coordinates::LatLon,
    current_material: rendering::materials::PbrMaterial,
    orbiters: Vec<satellites::artificial::VenusSatellite>,
    diagnostics: VenusDiagnostics,
}

impl VenusSimulation {
    fn new() -> Self {
        let heightmap = terrain::heightmap::Heightmap::generate(180, 90);
        let observer_alt = heightmap.sample(0.0, 0.0);
        Self {
            epoch: temporal::epoch::VenusEpoch::j2000(),
            time_scale: temporal::time_scale::TimeScale::new(),
            orbit: physics::orbit::VenusOrbit::new(),
            rotation: physics::rotation::VenusRotation::new(),
            rotation_angle_rad: 0.0,
            climate: atmosphere::climate::VenusClimateState::current(),
            heightmap,
            lod: terrain::lod::LodTerrain::new(terrain::lod::LodConfig::default()),
            atmo_render: rendering::atmosphere_scattering::VenusAtmosphereParams::default(),
            biome_classifier: terrain::texturing::BiomeClassifier,
            terrain_shader: rendering::shaders::terrain(),
            observer: geodata::coordinates::LatLon::new(0.0, 0.0, observer_alt),
            current_material: rendering::materials::basaltic_plain(),
            orbiters: vec![
                satellites::artificial::akatsuki(),
                satellites::artificial::venus_express(),
            ],
            diagnostics: VenusDiagnostics::default(),
        }
    }

    fn tick(&mut self, frame: u64, real_dt_s: f64) {
        let sim_dt = self.time_scale.simulation_dt(real_dt_s);
        if sim_dt <= 0.0 {
            return;
        }

        self.epoch.advance_seconds(sim_dt);
        self.orbit.step(sim_dt);
        self.rotation_angle_rad = (self.rotation_angle_rad
            + self.rotation.angular_velocity_rad_s * sim_dt)
            .rem_euclid(2.0 * std::f64::consts::PI);

        let observer_alt = self
            .heightmap
            .sample(self.observer.lat_deg, self.observer.lon_deg);
        self.observer = geodata::coordinates::LatLon::new(
            self.observer.lat_deg,
            self.observer.lon_deg,
            observer_alt,
        );
        let observer_cart = self.observer.to_cartesian();
        let cal = temporal::calendar::VenusDate::from_julian_date(self.epoch.julian_date);
        let ls_deg = cal.ls_deg();
        let season = lighting::seasons::season_at(ls_deg);
        let local_time_h = (((self.rotation_angle_rad / (2.0 * std::f64::consts::PI)) * 24.0)
            + 12.0)
            .rem_euclid(24.0);
        let sun = lighting::solar_position::SolarPosition::compute(
            ls_deg,
            local_time_h,
            self.observer.lat_deg,
            self.observer.lon_deg,
        );
        let day_night = lighting::day_night::DayNightCycle::new(ls_deg);
        let daylight_state = day_night.state_at(self.observer.lat_deg, local_time_h);

        let cloud_top_pressure =
            atmosphere::layers::barometric_pressure(rendering::clouds::cloud_top_altitude_m());
        let surface_pressure = atmosphere::layers::barometric_pressure(observer_alt.max(0.0));
        let mut cloud_optical_depth = self.climate.cloud_optical_depth;
        if matches!(
            season,
            lighting::seasons::Season::Summer | lighting::seasons::Season::Autumn
        ) {
            cloud_optical_depth = (cloud_optical_depth + 0.5 + (frame % 97) as f64 * 0.001)
                .min(GLOBAL_STORM_OPTICAL_DEPTH);
        }
        self.atmo_render.cloud_optical_depth = cloud_optical_depth;

        let solar_irradiance = self.orbit.solar_irradiance();
        let cloud_transmission =
            atmosphere::storms::SulfuricCloudLayer::global().surface_solar_fraction();
        let direct_transmission =
            self.atmo_render.direct_transmission(sun.elevation_deg) * cloud_transmission;
        let sky_color = self.atmo_render.sky_color(sun.elevation_deg);
        let superrotation = atmosphere::winds::superrotation_speed_m_s();
        let surface_wind = atmosphere::winds::mean_surface_wind_speed();
        let lightning = atmosphere::storms::lightning_rate_flashes_per_s();
        let solar_tide_bulge = physics::tides::SolarTide::at_mean_distance().tidal_bulge_height();
        let heat_flow = geology::plate_tectonics::surface_heat_flow_w_m2(3.0, GEOTHERMAL_GRADIENT);
        let interior = geology::plate_tectonics::VenusInterior::new();
        let rayleigh = interior.mantle_rayleigh_number();
        let thermal_stress = geology::erosion::thermal_fatigue_crack_growth_rate(10.0, 8.0e6);
        let erosion_flux = geology::erosion::AeolianErosion::typical().transport_rate(surface_wind);
        let chemical_weathering =
            geology::erosion::chemical_weathering_rate(self.climate.global_mean_temp_k, 0.4);
        let volcanism = geology::volcanism::VenusVolcanism::new();
        let total_volcano_volume = geology::volcanism::VenusVolcanism::major_volcanoes()
            .iter()
            .map(|v| v.volume_km3())
            .sum::<f64>();
        let cloud_hab = biosphere::ecosystems::cloud_layer_microhabitat();
        let biome = self.biome_classifier.classify(
            self.observer.lat_deg,
            observer_alt,
            atmosphere::heatwaves::global_mean_thermal_inertia(),
        );
        self.current_material = biome_material(biome);

        let camera = (observer_cart.x, observer_cart.y, observer_cart.z);
        self.lod.update(camera);
        self.terrain_shader.uniforms[0].1 = rendering::shaders::UniformValue::Float(VENUS_RADIUS);

        let paleoocean = hydrology::oceans::AncientOcean::hypothetical_early_venus();
        let paleolake = hydrology::lakes::hypothetical_early_lake();
        let lava_channel = hydrology::rivers::modeled_lava_channel();
        let orbiters_mean_period_h = self
            .orbiters
            .iter()
            .map(|o| o.orbital_period_s())
            .sum::<f64>()
            / self.orbiters.len() as f64
            / 3600.0;

        self.diagnostics = VenusDiagnostics {
            ls_deg,
            local_time_h,
            solar_irradiance_w_m2: solar_irradiance,
            direct_transmission,
            cloud_optical_depth,
            daylight_state,
            season,
            sky_color,
            cloud_top_pressure_pa: cloud_top_pressure,
            surface_pressure_pa: surface_pressure,
            superrotation_m_s: superrotation,
            surface_wind_m_s: surface_wind,
            lightning_rate_s: lightning,
            solar_tide_bulge_m: solar_tide_bulge,
            heat_flow_w_m2: heat_flow,
            rayleigh_number: rayleigh,
            thermal_stress,
            erosion_flux,
            chemical_weathering,
            volcanic_coverage: volcanism.total_volcanic_coverage(),
            volcanic_volume_km3: total_volcano_volume,
            cloud_habitability: cloud_hab.habitability_score(),
            biome,
            material_roughness: self.current_material.roughness,
            observer_elevation_m: observer_alt,
            orbiters_mean_period_h,
            paleoocean_volume_km3: paleoocean.estimated_volume_km3,
            paleolake_volume_km3: paleolake.volume_km3(),
            lava_channel_length_km: lava_channel.length_km,
            checksum: 0.0,
        };
        self.diagnostics.recompute_checksum();
    }
}

fn main() {
    let mut simulation = VenusSimulation::new();

    for frame in 0..FRAME_COUNT as u64 {
        simulation.tick(frame, REAL_DT_S);
    }
}