Vendor Bevy rendering crates (Phase 1 complete)

Closes #6, #7, #8, #9, #10
Refs #2, #122

Vendored bevy_render, bevy_core_pipeline, and bevy_pbr from Bevy v0.17.2
(commit 566358363126dd69f6e457e47f306c68f8041d2a) into libmarathon.

- ~51K LOC vendored to crates/libmarathon/src/render/
- Merged bevy_render_macros into crates/macros/
- Fixed 773→0 compilation errors
- Updated dependencies (encase 0.10→0.11, added 4 new deps)
- Removed bevy_render/pbr/core_pipeline from app Cargo features

All builds passing, macOS smoke test successful.

Signed-off-by: Sienna Meridian Satterwhite <sienna@r3t.io>

crates/libmarathon/src/render/pbr/atmosphere/aerial_view_lut.wgsl

@@ -0,0 +1,65 @@
+#import bevy_pbr::{
+    mesh_view_types::{Lights, DirectionalLight},
+    atmosphere::{
+        types::{Atmosphere, AtmosphereSettings},
+        bindings::{atmosphere, settings, view, lights, aerial_view_lut_out},
+        functions::{
+            sample_transmittance_lut, sample_atmosphere, rayleigh, henyey_greenstein,
+            sample_multiscattering_lut, AtmosphereSample, sample_local_inscattering,
+            uv_to_ndc, max_atmosphere_distance, uv_to_ray_direction, 
+            MIDPOINT_RATIO, get_view_position
+        },
+    }
+}
+
+
+@group(0) @binding(13) var aerial_view_lut_out: texture_storage_3d<rgba16float, write>;
+
+@compute
+@workgroup_size(16, 16, 1)
+fn main(@builtin(global_invocation_id) idx: vec3<u32>) {
+    if any(idx.xy > settings.aerial_view_lut_size.xy) { return; }
+
+    let uv = (vec2<f32>(idx.xy) + 0.5) / vec2<f32>(settings.aerial_view_lut_size.xy);
+    let ray_dir = uv_to_ray_direction(uv);
+    let world_pos = get_view_position();
+    
+    let r = length(world_pos);
+    let t_max = settings.aerial_view_lut_max_distance;
+
+    var prev_t = 0.0;
+    var total_inscattering = vec3(0.0);
+    var throughput = vec3(1.0);
+
+    for (var slice_i: u32 = 0; slice_i < settings.aerial_view_lut_size.z; slice_i++) {
+        for (var step_i: u32 = 0; step_i < settings.aerial_view_lut_samples; step_i++) {
+            let t_i = t_max * (f32(slice_i) + ((f32(step_i) + MIDPOINT_RATIO) / f32(settings.aerial_view_lut_samples))) / f32(settings.aerial_view_lut_size.z);
+            let dt = (t_i - prev_t);
+            prev_t = t_i;
+
+            let sample_pos = world_pos + ray_dir * t_i;
+            let local_r = length(sample_pos);
+            let local_up = normalize(sample_pos);
+
+            let local_atmosphere = sample_atmosphere(local_r);
+            let sample_optical_depth = local_atmosphere.extinction * dt;
+            let sample_transmittance = exp(-sample_optical_depth);
+
+            // evaluate one segment of the integral
+            var inscattering = sample_local_inscattering(local_atmosphere, ray_dir, sample_pos);
+
+            // Analytical integration of the single scattering term in the radiance transfer equation
+            let s_int = (inscattering - inscattering * sample_transmittance) / local_atmosphere.extinction;
+            total_inscattering += throughput * s_int;
+
+            throughput *= sample_transmittance;
+            if all(throughput < vec3(0.001)) {
+                break;
+            }
+        }
+
+        // Store in log space to allow linear interpolation of exponential values between slices
+        let log_inscattering = log(max(total_inscattering, vec3(1e-6)));
+        textureStore(aerial_view_lut_out, vec3(vec2<u32>(idx.xy), slice_i), vec4(log_inscattering, 0.0));
+    }
+}

crates/libmarathon/src/render/pbr/atmosphere/bindings.wgsl

@@ -0,0 +1,22 @@
+#define_import_path bevy_pbr::atmosphere::bindings
+
+#import bevy_render::view::View;
+
+#import bevy_pbr::{
+    mesh_view_types::Lights,
+    atmosphere::types::{Atmosphere, AtmosphereSettings, AtmosphereTransforms}
+}
+
+@group(0) @binding(0) var<uniform> atmosphere: Atmosphere;
+@group(0) @binding(1) var<uniform> settings: AtmosphereSettings;
+@group(0) @binding(2) var<uniform> atmosphere_transforms: AtmosphereTransforms;
+@group(0) @binding(3) var<uniform> view: View;
+@group(0) @binding(4) var<uniform> lights: Lights;
+@group(0) @binding(5) var transmittance_lut: texture_2d<f32>;
+@group(0) @binding(6) var transmittance_lut_sampler: sampler;
+@group(0) @binding(7) var multiscattering_lut: texture_2d<f32>;
+@group(0) @binding(8) var multiscattering_lut_sampler: sampler;
+@group(0) @binding(9) var sky_view_lut: texture_2d<f32>;
+@group(0) @binding(10) var sky_view_lut_sampler: sampler;
+@group(0) @binding(11) var aerial_view_lut: texture_3d<f32>;
+@group(0) @binding(12) var aerial_view_lut_sampler: sampler;

crates/libmarathon/src/render/pbr/atmosphere/bruneton_functions.wgsl

@@ -0,0 +1,139 @@
+// Copyright (c) 2017 Eric Bruneton
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+// 1. Redistributions of source code must retain the above copyright
+//    notice, this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright
+//    notice, this list of conditions and the following disclaimer in the
+//    documentation and/or other materials provided with the distribution.
+// 3. Neither the name of the copyright holders nor the names of its
+//    contributors may be used to endorse or promote products derived from
+//    this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+// THE POSSIBILITY OF SUCH DAMAGE.
+//
+// Precomputed Atmospheric Scattering
+// Copyright (c) 2008 INRIA
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+// 1. Redistributions of source code must retain the above copyright
+//    notice, this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright
+//    notice, this list of conditions and the following disclaimer in the
+//    documentation and/or other materials provided with the distribution.
+// 3. Neither the name of the copyright holders nor the names of its
+//    contributors may be used to endorse or promote products derived from
+//    this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+// THE POSSIBILITY OF SUCH DAMAGE.
+
+#define_import_path bevy_pbr::atmosphere::bruneton_functions
+
+#import bevy_pbr::atmosphere::{
+    types::Atmosphere,
+    bindings::atmosphere,
+}
+
+// Mapping from view height (r) and zenith cos angle (mu) to UV coordinates in the transmittance LUT
+// Assuming r between ground and top atmosphere boundary, and mu= cos(zenith_angle)
+// Chosen to increase precision near the ground and to work around a discontinuity at the horizon
+// See Bruneton and Neyret 2008, "Precomputed Atmospheric Scattering" section 4
+fn transmittance_lut_r_mu_to_uv(r: f32, mu: f32) -> vec2<f32> {
+  // Distance along a horizontal ray from the ground to the top atmosphere boundary
+    let H = sqrt(atmosphere.top_radius * atmosphere.top_radius - atmosphere.bottom_radius * atmosphere.bottom_radius);
+
+  // Distance from a point at height r to the horizon
+  // ignore the case where r <= atmosphere.bottom_radius
+    let rho = sqrt(max(r * r - atmosphere.bottom_radius * atmosphere.bottom_radius, 0.0));
+
+  // Distance from a point at height r to the top atmosphere boundary at zenith angle mu
+    let d = distance_to_top_atmosphere_boundary(r, mu);
+
+  // Minimum and maximum distance to the top atmosphere boundary from a point at height r
+    let d_min = atmosphere.top_radius - r; // length of the ray straight up to the top atmosphere boundary
+    let d_max = rho + H; // length of the ray to the top atmosphere boundary and grazing the horizon
+
+    let u = (d - d_min) / (d_max - d_min);
+    let v = rho / H;
+    return vec2<f32>(u, v);
+}
+
+// Inverse of the mapping above, mapping from UV coordinates in the transmittance LUT to view height (r) and zenith cos angle (mu)
+fn transmittance_lut_uv_to_r_mu(uv: vec2<f32>) -> vec2<f32> {
+  // Distance to top atmosphere boundary for a horizontal ray at ground level
+    let H = sqrt(atmosphere.top_radius * atmosphere.top_radius - atmosphere.bottom_radius * atmosphere.bottom_radius);
+
+  // Distance to the horizon, from which we can compute r:
+    let rho = H * uv.y;
+    let r = sqrt(rho * rho + atmosphere.bottom_radius * atmosphere.bottom_radius);
+
+  // Distance to the top atmosphere boundary for the ray (r,mu), and its minimum
+  // and maximum values over all mu- obtained for (r,1) and (r,mu_horizon) -
+  // from which we can recover mu:
+    let d_min = atmosphere.top_radius - r;
+    let d_max = rho + H;
+    let d = d_min + uv.x * (d_max - d_min);
+
+    var mu: f32;
+    if d == 0.0 {
+        mu = 1.0;
+    } else {
+        mu = (H * H - rho * rho - d * d) / (2.0 * r * d);
+    }
+
+    mu = clamp(mu, -1.0, 1.0);
+
+    return vec2<f32>(r, mu);
+}
+
+/// Simplified ray-sphere intersection
+/// where:
+/// Ray origin, o = [0,0,r] with r <= atmosphere.top_radius
+/// mu is the cosine of spherical coordinate theta (-1.0 <= mu <= 1.0)
+/// so ray direction in spherical coordinates is [1,acos(mu),0] which needs to be converted to cartesian
+/// Direction of ray, u = [0,sqrt(1-mu*mu),mu]
+/// Center of sphere, c = [0,0,0]
+/// Radius of sphere, r = atmosphere.top_radius
+/// This function solves the quadratic equation for line-sphere intersection simplified under these assumptions
+fn distance_to_top_atmosphere_boundary(r: f32, mu: f32) -> f32 {
+  // ignore the case where r > atmosphere.top_radius
+    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atmosphere.top_radius * atmosphere.top_radius, 0.0);
+    return max(-r * mu + sqrt(positive_discriminant), 0.0);
+}
+
+/// Simplified ray-sphere intersection
+/// as above for intersections with the ground
+fn distance_to_bottom_atmosphere_boundary(r: f32, mu: f32) -> f32 {
+    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atmosphere.bottom_radius * atmosphere.bottom_radius, 0.0);
+    return max(-r * mu - sqrt(positive_discriminant), 0.0);
+}
+
+fn ray_intersects_ground(r: f32, mu: f32) -> bool {
+    return mu < 0.0 && r * r * (mu * mu - 1.0) + atmosphere.bottom_radius * atmosphere.bottom_radius >= 0.0;
+}

crates/libmarathon/src/render/pbr/atmosphere/environment.rs

@@ -0,0 +1,332 @@
+use crate::render::pbr::{
+    resources::{
+        AtmosphereSamplers, AtmosphereTextures, AtmosphereTransform, AtmosphereTransforms,
+        AtmosphereTransformsOffset,
+    },
+    GpuAtmosphereSettings, GpuLights, LightMeta, ViewLightsUniformOffset,
+};
+use bevy_asset::{load_embedded_asset, AssetServer, Assets, Handle, RenderAssetUsages};
+use bevy_ecs::{
+    component::Component,
+    entity::Entity,
+    query::{QueryState, With, Without},
+    resource::Resource,
+    system::{lifetimeless::Read, Commands, Query, Res, ResMut},
+    world::{FromWorld, World},
+};
+use bevy_image::Image;
+use bevy_light::{AtmosphereEnvironmentMapLight, GeneratedEnvironmentMapLight};
+use bevy_math::{Quat, UVec2};
+use crate::render::{
+    extract_component::{ComponentUniforms, DynamicUniformIndex, ExtractComponent},
+    render_asset::RenderAssets,
+    render_graph::{Node, NodeRunError, RenderGraphContext},
+    render_resource::{binding_types::*, *},
+    renderer::{RenderContext, RenderDevice},
+    texture::{CachedTexture, GpuImage},
+    view::{ViewUniform, ViewUniformOffset, ViewUniforms},
+};
+use bevy_utils::default;
+use tracing::warn;
+
+use super::Atmosphere;
+
+// Render world representation of an environment map light for the atmosphere
+#[derive(Component, ExtractComponent, Clone)]
+pub struct AtmosphereEnvironmentMap {
+    pub environment_map: Handle<Image>,
+    pub size: UVec2,
+}
+
+#[derive(Component)]
+pub struct AtmosphereProbeTextures {
+    pub environment: TextureView,
+    pub transmittance_lut: CachedTexture,
+    pub multiscattering_lut: CachedTexture,
+    pub sky_view_lut: CachedTexture,
+    pub aerial_view_lut: CachedTexture,
+}
+
+#[derive(Component)]
+pub(crate) struct AtmosphereProbeBindGroups {
+    pub environment: BindGroup,
+}
+
+#[derive(Resource)]
+pub struct AtmosphereProbeLayouts {
+    pub environment: BindGroupLayout,
+}
+
+#[derive(Resource)]
+pub struct AtmosphereProbePipeline {
+    pub environment: CachedComputePipelineId,
+}
+
+pub fn init_atmosphere_probe_layout(mut commands: Commands, render_device: Res<RenderDevice>) {
+    let environment = render_device.create_bind_group_layout(
+        "environment_bind_group_layout",
+        &BindGroupLayoutEntries::sequential(
+            ShaderStages::COMPUTE,
+            (
+                uniform_buffer::<Atmosphere>(true),
+                uniform_buffer::<GpuAtmosphereSettings>(true),
+                uniform_buffer::<AtmosphereTransform>(true),
+                uniform_buffer::<ViewUniform>(true),
+                uniform_buffer::<GpuLights>(true),
+                texture_2d(TextureSampleType::Float { filterable: true }), //transmittance lut and sampler
+                sampler(SamplerBindingType::Filtering),
+                texture_2d(TextureSampleType::Float { filterable: true }), //multiscattering lut and sampler
+                sampler(SamplerBindingType::Filtering),
+                texture_2d(TextureSampleType::Float { filterable: true }), //sky view lut and sampler
+                sampler(SamplerBindingType::Filtering),
+                texture_3d(TextureSampleType::Float { filterable: true }), //aerial view lut ans sampler
+                sampler(SamplerBindingType::Filtering),
+                texture_storage_2d_array(
+                    // output 2D array texture
+                    TextureFormat::Rgba16Float,
+                    StorageTextureAccess::WriteOnly,
+                ),
+            ),
+        ),
+    );
+
+    commands.insert_resource(AtmosphereProbeLayouts { environment });
+}
+
+pub(super) fn prepare_atmosphere_probe_bind_groups(
+    probes: Query<(Entity, &AtmosphereProbeTextures), With<AtmosphereEnvironmentMap>>,
+    render_device: Res<RenderDevice>,
+    layouts: Res<AtmosphereProbeLayouts>,
+    samplers: Res<AtmosphereSamplers>,
+    view_uniforms: Res<ViewUniforms>,
+    lights_uniforms: Res<LightMeta>,
+    atmosphere_transforms: Res<AtmosphereTransforms>,
+    atmosphere_uniforms: Res<ComponentUniforms<Atmosphere>>,
+    settings_uniforms: Res<ComponentUniforms<GpuAtmosphereSettings>>,
+    mut commands: Commands,
+) {
+    for (entity, textures) in &probes {
+        let environment = render_device.create_bind_group(
+            "environment_bind_group",
+            &layouts.environment,
+            &BindGroupEntries::sequential((
+                atmosphere_uniforms.binding().unwrap(),
+                settings_uniforms.binding().unwrap(),
+                atmosphere_transforms.uniforms().binding().unwrap(),
+                view_uniforms.uniforms.binding().unwrap(),
+                lights_uniforms.view_gpu_lights.binding().unwrap(),
+                &textures.transmittance_lut.default_view,
+                &samplers.transmittance_lut,
+                &textures.multiscattering_lut.default_view,
+                &samplers.multiscattering_lut,
+                &textures.sky_view_lut.default_view,
+                &samplers.sky_view_lut,
+                &textures.aerial_view_lut.default_view,
+                &samplers.aerial_view_lut,
+                &textures.environment,
+            )),
+        );
+
+        commands
+            .entity(entity)
+            .insert(AtmosphereProbeBindGroups { environment });
+    }
+}
+
+pub(super) fn prepare_probe_textures(
+    view_textures: Query<&AtmosphereTextures, With<Atmosphere>>,
+    probes: Query<
+        (Entity, &AtmosphereEnvironmentMap),
+        (
+            With<AtmosphereEnvironmentMap>,
+            Without<AtmosphereProbeTextures>,
+        ),
+    >,
+    gpu_images: Res<RenderAssets<GpuImage>>,
+    mut commands: Commands,
+) {
+    for (probe, render_env_map) in &probes {
+        let environment = gpu_images.get(&render_env_map.environment_map).unwrap();
+        // create a cube view
+        let environment_view = environment.texture.create_view(&TextureViewDescriptor {
+            dimension: Some(TextureViewDimension::D2Array),
+            ..Default::default()
+        });
+        // Get the first view entity's textures to borrow
+        if let Some(view_textures) = view_textures.iter().next() {
+            commands.entity(probe).insert(AtmosphereProbeTextures {
+                environment: environment_view,
+                transmittance_lut: view_textures.transmittance_lut.clone(),
+                multiscattering_lut: view_textures.multiscattering_lut.clone(),
+                sky_view_lut: view_textures.sky_view_lut.clone(),
+                aerial_view_lut: view_textures.aerial_view_lut.clone(),
+            });
+        }
+    }
+}
+
+pub fn init_atmosphere_probe_pipeline(
+    pipeline_cache: Res<PipelineCache>,
+    layouts: Res<AtmosphereProbeLayouts>,
+    asset_server: Res<AssetServer>,
+    mut commands: Commands,
+) {
+    let environment = pipeline_cache.queue_compute_pipeline(ComputePipelineDescriptor {
+        label: Some("environment_pipeline".into()),
+        layout: vec![layouts.environment.clone()],
+        shader: load_embedded_asset!(asset_server.as_ref(), "environment.wgsl"),
+        ..default()
+    });
+    commands.insert_resource(AtmosphereProbePipeline { environment });
+}
+
+// Ensure power-of-two dimensions to avoid edge update issues on cubemap faces
+pub fn validate_environment_map_size(size: UVec2) -> UVec2 {
+    let new_size = UVec2::new(
+        size.x.max(1).next_power_of_two(),
+        size.y.max(1).next_power_of_two(),
+    );
+    if new_size != size {
+        warn!(
+            "Non-power-of-two AtmosphereEnvironmentMapLight size {}, correcting to {new_size}",
+            size
+        );
+    }
+    new_size
+}
+
+pub fn prepare_atmosphere_probe_components(
+    probes: Query<(Entity, &AtmosphereEnvironmentMapLight), (Without<AtmosphereEnvironmentMap>,)>,
+    mut commands: Commands,
+    mut images: ResMut<Assets<Image>>,
+) {
+    for (entity, env_map_light) in &probes {
+        // Create a cubemap image in the main world that we can reference
+        let new_size = validate_environment_map_size(env_map_light.size);
+        let mut environment_image = Image::new_fill(
+            Extent3d {
+                width: new_size.x,
+                height: new_size.y,
+                depth_or_array_layers: 6,
+            },
+            TextureDimension::D2,
+            &[0; 8],
+            TextureFormat::Rgba16Float,
+            RenderAssetUsages::all(),
+        );
+
+        environment_image.texture_view_descriptor = Some(TextureViewDescriptor {
+            dimension: Some(TextureViewDimension::Cube),
+            ..Default::default()
+        });
+
+        environment_image.texture_descriptor.usage = TextureUsages::TEXTURE_BINDING
+            | TextureUsages::STORAGE_BINDING
+            | TextureUsages::COPY_SRC;
+
+        // Add the image to assets to get a handle
+        let environment_handle = images.add(environment_image);
+
+        commands.entity(entity).insert(AtmosphereEnvironmentMap {
+            environment_map: environment_handle.clone(),
+            size: new_size,
+        });
+
+        commands
+            .entity(entity)
+            .insert(GeneratedEnvironmentMapLight {
+                environment_map: environment_handle,
+                intensity: env_map_light.intensity,
+                rotation: Quat::IDENTITY,
+                affects_lightmapped_mesh_diffuse: env_map_light.affects_lightmapped_mesh_diffuse,
+            });
+    }
+}
+
+pub(super) struct EnvironmentNode {
+    main_view_query: QueryState<(
+        Read<DynamicUniformIndex<Atmosphere>>,
+        Read<DynamicUniformIndex<GpuAtmosphereSettings>>,
+        Read<AtmosphereTransformsOffset>,
+        Read<ViewUniformOffset>,
+        Read<ViewLightsUniformOffset>,
+    )>,
+    probe_query: QueryState<(
+        Read<AtmosphereProbeBindGroups>,
+        Read<AtmosphereEnvironmentMap>,
+    )>,
+}
+
+impl FromWorld for EnvironmentNode {
+    fn from_world(world: &mut World) -> Self {
+        Self {
+            main_view_query: QueryState::new(world),
+            probe_query: QueryState::new(world),
+        }
+    }
+}
+
+impl Node for EnvironmentNode {
+    fn update(&mut self, world: &mut World) {
+        self.main_view_query.update_archetypes(world);
+        self.probe_query.update_archetypes(world);
+    }
+
+    fn run(
+        &self,
+        graph: &mut RenderGraphContext,
+        render_context: &mut RenderContext,
+        world: &World,
+    ) -> Result<(), NodeRunError> {
+        let pipeline_cache = world.resource::<PipelineCache>();
+        let pipelines = world.resource::<AtmosphereProbePipeline>();
+        let view_entity = graph.view_entity();
+
+        let Some(environment_pipeline) = pipeline_cache.get_compute_pipeline(pipelines.environment)
+        else {
+            return Ok(());
+        };
+
+        let (Ok((
+            atmosphere_uniforms_offset,
+            settings_uniforms_offset,
+            atmosphere_transforms_offset,
+            view_uniforms_offset,
+            lights_uniforms_offset,
+        )),) = (self.main_view_query.get_manual(world, view_entity),)
+        else {
+            return Ok(());
+        };
+
+        for (bind_groups, env_map_light) in self.probe_query.iter_manual(world) {
+            let mut pass =
+                render_context
+                    .command_encoder()
+                    .begin_compute_pass(&ComputePassDescriptor {
+                        label: Some("environment_pass"),
+                        timestamp_writes: None,
+                    });
+
+            pass.set_pipeline(environment_pipeline);
+            pass.set_bind_group(
+                0,
+                &bind_groups.environment,
+                &[
+                    atmosphere_uniforms_offset.index(),
+                    settings_uniforms_offset.index(),
+                    atmosphere_transforms_offset.index(),
+                    view_uniforms_offset.offset,
+                    lights_uniforms_offset.offset,
+                ],
+            );
+
+            pass.dispatch_workgroups(
+                env_map_light.size.x / 8,
+                env_map_light.size.y / 8,
+                6, // 6 cubemap faces
+            );
+        }
+
+        Ok(())
+    }
+}

crates/libmarathon/src/render/pbr/atmosphere/environment.wgsl

@@ -0,0 +1,39 @@
+#import bevy_pbr::{
+    atmosphere::{
+        functions::{direction_world_to_atmosphere, sample_sky_view_lut, get_view_position},
+    },
+    utils::sample_cube_dir
+}
+
+@group(0) @binding(13) var output: texture_storage_2d_array<rgba16float, write>;
+
+@compute @workgroup_size(8, 8, 1)
+fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
+    let dimensions = textureDimensions(output);
+    let slice_index = global_id.z;
+    
+    if (global_id.x >= dimensions.x || global_id.y >= dimensions.y || slice_index >= 6u) {
+        return;
+    }
+    
+    // Calculate normalized UV coordinates for this pixel
+    let uv = vec2<f32>(
+        (f32(global_id.x) + 0.5) / f32(dimensions.x),
+        (f32(global_id.y) + 0.5) / f32(dimensions.y)
+    );
+
+    var ray_dir_ws = sample_cube_dir(uv, slice_index);
+    
+    // invert the z direction to account for cubemaps being lefthanded
+    ray_dir_ws.z = -ray_dir_ws.z;
+
+    let world_pos = get_view_position();
+    let r = length(world_pos);
+    let up = normalize(world_pos);
+
+    let ray_dir_as = direction_world_to_atmosphere(ray_dir_ws.xyz, up);
+    let inscattering = sample_sky_view_lut(r, ray_dir_as);
+    let color = vec4<f32>(inscattering, 1.0);
+
+    textureStore(output, vec2<i32>(global_id.xy), i32(slice_index), color);
+}

crates/libmarathon/src/render/pbr/atmosphere/functions.wgsl

@@ -0,0 +1,528 @@
+#define_import_path bevy_pbr::atmosphere::functions
+
+#import bevy_render::maths::{PI, HALF_PI, PI_2, fast_acos, fast_acos_4, fast_atan2, ray_sphere_intersect}
+
+#import bevy_pbr::atmosphere::{
+    types::Atmosphere,
+    bindings::{
+        atmosphere, settings, view, lights, transmittance_lut, transmittance_lut_sampler, 
+        multiscattering_lut, multiscattering_lut_sampler, sky_view_lut, sky_view_lut_sampler,
+        aerial_view_lut, aerial_view_lut_sampler, atmosphere_transforms
+    },
+    bruneton_functions::{
+        transmittance_lut_r_mu_to_uv, transmittance_lut_uv_to_r_mu, 
+        ray_intersects_ground, distance_to_top_atmosphere_boundary, 
+        distance_to_bottom_atmosphere_boundary
+    },
+}
+
+// NOTE FOR CONVENTIONS: 
+// r:
+//   radius, or distance from planet center 
+//
+// altitude:
+//   distance from planet **surface**
+//
+// mu:
+//   cosine of the zenith angle of a ray with
+//   respect to the planet normal
+//
+// atmosphere space:
+//   abbreviated as "as" (contrast with vs, cs, ws), this space is similar
+//   to view space, but with the camera positioned horizontally on the planet
+//   surface, so the horizon is a horizontal line centered vertically in the
+//   frame. This enables the non-linear latitude parametrization the paper uses 
+//   to concentrate detail near the horizon 
+
+
+// CONSTANTS
+const FRAC_PI: f32 = 0.3183098862; // 1 / π
+const FRAC_2_PI: f32 = 0.15915494309;  // 1 / (2π)
+const FRAC_3_16_PI: f32 = 0.0596831036594607509; // 3 / (16π)
+const FRAC_4_PI: f32 = 0.07957747154594767; // 1 / (4π)
+const ROOT_2: f32 = 1.41421356; // √2
+const EPSILON: f32 = 1.0; // 1 meter
+
+// During raymarching, each segment is sampled at a single point. This constant determines
+// where in the segment that sample is taken (0.0 = start, 0.5 = middle, 1.0 = end).
+// We use 0.3 to sample closer to the start of each segment, which better approximates
+// the exponential falloff of atmospheric density.
+const MIDPOINT_RATIO: f32 = 0.3;
+
+// LUT UV PARAMETERIZATIONS
+
+fn unit_to_sub_uvs(val: vec2<f32>, resolution: vec2<f32>) -> vec2<f32> {
+    return (val + 0.5f / resolution) * (resolution / (resolution + 1.0f));
+}
+
+fn sub_uvs_to_unit(val: vec2<f32>, resolution: vec2<f32>) -> vec2<f32> {
+    return (val - 0.5f / resolution) * (resolution / (resolution - 1.0f));
+}
+
+fn multiscattering_lut_r_mu_to_uv(r: f32, mu: f32) -> vec2<f32> {
+    let u = 0.5 + 0.5 * mu;
+    let v = saturate((r - atmosphere.bottom_radius) / (atmosphere.top_radius - atmosphere.bottom_radius)); //TODO
+    return unit_to_sub_uvs(vec2(u, v), vec2<f32>(settings.multiscattering_lut_size));
+}
+
+fn multiscattering_lut_uv_to_r_mu(uv: vec2<f32>) -> vec2<f32> {
+    let adj_uv = sub_uvs_to_unit(uv, vec2<f32>(settings.multiscattering_lut_size));
+    let r = mix(atmosphere.bottom_radius, atmosphere.top_radius, adj_uv.y);
+    let mu = adj_uv.x * 2 - 1;
+    return vec2(r, mu);
+}
+
+fn sky_view_lut_r_mu_azimuth_to_uv(r: f32, mu: f32, azimuth: f32) -> vec2<f32> {
+    let u = (azimuth * FRAC_2_PI) + 0.5;
+
+    let v_horizon = sqrt(r * r - atmosphere.bottom_radius * atmosphere.bottom_radius);
+    let cos_beta = v_horizon / r;
+    // Using fast_acos_4 for better precision at small angles
+    // to avoid artifacts at the horizon
+    let beta = fast_acos_4(cos_beta);
+    let horizon_zenith = PI - beta;
+    let view_zenith = fast_acos_4(mu);
+
+    // Apply non-linear transformation to compress more texels 
+    // near the horizon where high-frequency details matter most
+    // l is latitude in [-π/2, π/2] and v is texture coordinate in [0,1]
+    let l = view_zenith - horizon_zenith;
+    let abs_l = abs(l);
+
+    let v = 0.5 + 0.5 * sign(l) * sqrt(abs_l / HALF_PI);
+
+    return unit_to_sub_uvs(vec2(u, v), vec2<f32>(settings.sky_view_lut_size));
+}
+
+fn sky_view_lut_uv_to_zenith_azimuth(r: f32, uv: vec2<f32>) -> vec2<f32> {
+    let adj_uv = sub_uvs_to_unit(vec2(uv.x, 1.0 - uv.y), vec2<f32>(settings.sky_view_lut_size));
+    let azimuth = (adj_uv.x - 0.5) * PI_2;
+
+    // Horizon parameters
+    let v_horizon = sqrt(r * r - atmosphere.bottom_radius * atmosphere.bottom_radius);
+    let cos_beta = v_horizon / r;
+    let beta = fast_acos_4(cos_beta);
+    let horizon_zenith = PI - beta;
+
+    // Inverse of horizon-detail mapping to recover original latitude from texture coordinate
+    let t = abs(2.0 * (adj_uv.y - 0.5));
+    let l = sign(adj_uv.y - 0.5) * HALF_PI * t * t;
+
+    return vec2(horizon_zenith - l, azimuth);
+}
+
+// LUT SAMPLING
+
+fn sample_transmittance_lut(r: f32, mu: f32) -> vec3<f32> {
+    let uv = transmittance_lut_r_mu_to_uv(r, mu);
+    return textureSampleLevel(transmittance_lut, transmittance_lut_sampler, uv, 0.0).rgb;
+}
+
+// NOTICE: This function is copyrighted by Eric Bruneton and INRIA, and falls
+// under the license reproduced in bruneton_functions.wgsl (variant of MIT license)
+//
+// FIXME: this function should be in bruneton_functions.wgsl, but because naga_oil doesn't 
+// support cyclic imports it's stuck here
+fn sample_transmittance_lut_segment(r: f32, mu: f32, t: f32) -> vec3<f32> {
+    let r_t = get_local_r(r, mu, t);
+    let mu_t = clamp((r * mu + t) / r_t, -1.0, 1.0);
+
+    if ray_intersects_ground(r, mu) {
+        return min(
+            sample_transmittance_lut(r_t, -mu_t) / sample_transmittance_lut(r, -mu),
+            vec3(1.0)
+        );
+    } else {
+        return min(
+            sample_transmittance_lut(r, mu) / sample_transmittance_lut(r_t, mu_t), vec3(1.0)
+        );
+    }
+}
+
+fn sample_multiscattering_lut(r: f32, mu: f32) -> vec3<f32> {
+    let uv = multiscattering_lut_r_mu_to_uv(r, mu);
+    return textureSampleLevel(multiscattering_lut, multiscattering_lut_sampler, uv, 0.0).rgb;
+}
+
+fn sample_sky_view_lut(r: f32, ray_dir_as: vec3<f32>) -> vec3<f32> {
+    let mu = ray_dir_as.y;
+    let azimuth = fast_atan2(ray_dir_as.x, -ray_dir_as.z);
+    let uv = sky_view_lut_r_mu_azimuth_to_uv(r, mu, azimuth);
+    return textureSampleLevel(sky_view_lut, sky_view_lut_sampler, uv, 0.0).rgb;
+}
+
+fn ndc_to_camera_dist(ndc: vec3<f32>) -> f32 {
+    let view_pos = view.view_from_clip * vec4(ndc, 1.0);
+    let t = length(view_pos.xyz / view_pos.w) * settings.scene_units_to_m;
+    return t;
+}
+
+// RGB channels: total inscattered light along the camera ray to the current sample.
+// A channel: average transmittance across all wavelengths to the current sample.
+fn sample_aerial_view_lut(uv: vec2<f32>, t: f32) -> vec3<f32> {
+    let t_max = settings.aerial_view_lut_max_distance;
+    let num_slices = f32(settings.aerial_view_lut_size.z);
+    // Each texel stores the value of the scattering integral over the whole slice,
+    // which requires us to offset the w coordinate by half a slice. For
+    // example, if we wanted the value of the integral at the boundary between slices,
+    // we'd need to sample at the center of the previous slice, and vice-versa for
+    // sampling in the center of a slice.
+    let uvw = vec3(uv, saturate(t / t_max - 0.5 / num_slices));
+    let sample = textureSampleLevel(aerial_view_lut, aerial_view_lut_sampler, uvw, 0.0);
+    // Since sampling anywhere between w=0 and w=t_slice will clamp to the first slice,
+    // we need to do a linear step over the first slice towards zero at the camera's
+    // position to recover the correct integral value.
+    let t_slice = t_max / num_slices;
+    let fade = saturate(t / t_slice);
+    // Recover the values from log space
+    return exp(sample.rgb) * fade;
+}
+
+// PHASE FUNCTIONS
+
+// -(L . V) == (L . -V). -V here is our ray direction, which points away from the view 
+// instead of towards it (which would be the *view direction*, V)
+
+// evaluates the rayleigh phase function, which describes the likelihood
+// of a rayleigh scattering event scattering light from the light direction towards the view
+fn rayleigh(neg_LdotV: f32) -> f32 {
+    return FRAC_3_16_PI * (1 + (neg_LdotV * neg_LdotV));
+}
+
+// evaluates the henyey-greenstein phase function, which describes the likelihood
+// of a mie scattering event scattering light from the light direction towards the view
+fn henyey_greenstein(neg_LdotV: f32) -> f32 {
+    let g = atmosphere.mie_asymmetry;
+    let denom = 1.0 + g * g - 2.0 * g * neg_LdotV;
+    return FRAC_4_PI * (1.0 - g * g) / (denom * sqrt(denom));
+}
+
+// ATMOSPHERE SAMPLING
+
+struct AtmosphereSample {
+    /// units: m^-1
+    rayleigh_scattering: vec3<f32>,
+
+    /// units: m^-1
+    mie_scattering: f32,
+
+    /// the sum of scattering and absorption. Since the phase function doesn't
+    /// matter for this, we combine rayleigh and mie extinction to a single 
+    //  value.
+    //
+    /// units: m^-1
+    extinction: vec3<f32>
+}
+
+/// Samples atmosphere optical densities at a given radius
+fn sample_atmosphere(r: f32) -> AtmosphereSample {
+    let altitude = clamp(r, atmosphere.bottom_radius, atmosphere.top_radius) - atmosphere.bottom_radius;
+
+    // atmosphere values at altitude
+    let mie_density = exp(-atmosphere.mie_density_exp_scale * altitude);
+    let rayleigh_density = exp(-atmosphere.rayleigh_density_exp_scale * altitude);
+    var ozone_density: f32 = max(0.0, 1.0 - (abs(altitude - atmosphere.ozone_layer_altitude) / (atmosphere.ozone_layer_width * 0.5)));
+
+    let mie_scattering = mie_density * atmosphere.mie_scattering;
+    let mie_absorption = mie_density * atmosphere.mie_absorption;
+    let mie_extinction = mie_scattering + mie_absorption;
+
+    let rayleigh_scattering = rayleigh_density * atmosphere.rayleigh_scattering;
+    // no rayleigh absorption
+    // rayleigh extinction is the sum of scattering and absorption
+
+    // ozone doesn't contribute to scattering
+    let ozone_absorption = ozone_density * atmosphere.ozone_absorption;
+
+    var sample: AtmosphereSample;
+    sample.rayleigh_scattering = rayleigh_scattering;
+    sample.mie_scattering = mie_scattering;
+    sample.extinction = rayleigh_scattering + mie_extinction + ozone_absorption;
+
+    return sample;
+}
+
+/// evaluates L_scat, equation 3 in the paper, which gives the total single-order scattering towards the view at a single point
+fn sample_local_inscattering(local_atmosphere: AtmosphereSample, ray_dir: vec3<f32>, world_pos: vec3<f32>) -> vec3<f32> {
+    let local_r = length(world_pos);
+    let local_up = normalize(world_pos);
+    var inscattering = vec3(0.0);
+    for (var light_i: u32 = 0u; light_i < lights.n_directional_lights; light_i++) {
+        let light = &lights.directional_lights[light_i];
+
+        let mu_light = dot((*light).direction_to_light, local_up);
+
+        // -(L . V) == (L . -V). -V here is our ray direction, which points away from the view
+        // instead of towards it (as is the convention for V)
+        let neg_LdotV = dot((*light).direction_to_light, ray_dir);
+
+        // Phase functions give the proportion of light
+        // scattered towards the camera for each scattering type
+        let rayleigh_phase = rayleigh(neg_LdotV);
+        let mie_phase = henyey_greenstein(neg_LdotV);
+        let scattering_coeff = local_atmosphere.rayleigh_scattering * rayleigh_phase + local_atmosphere.mie_scattering * mie_phase;
+
+        let transmittance_to_light = sample_transmittance_lut(local_r, mu_light);
+        let shadow_factor = transmittance_to_light * f32(!ray_intersects_ground(local_r, mu_light));
+
+        // Transmittance from scattering event to light source
+        let scattering_factor = shadow_factor * scattering_coeff;
+
+        // Additive factor from the multiscattering LUT
+        let psi_ms = sample_multiscattering_lut(local_r, mu_light);
+        let multiscattering_factor = psi_ms * (local_atmosphere.rayleigh_scattering + local_atmosphere.mie_scattering);
+
+        inscattering += (*light).color.rgb * (scattering_factor + multiscattering_factor);
+    }
+    return inscattering;
+}
+
+fn sample_sun_radiance(ray_dir_ws: vec3<f32>) -> vec3<f32> {
+    let view_pos = get_view_position();
+    let r = length(view_pos);
+    let up = normalize(view_pos);
+    let mu_view = dot(ray_dir_ws, up);
+    let shadow_factor = f32(!ray_intersects_ground(r, mu_view));
+    var sun_radiance = vec3(0.0);
+    for (var light_i: u32 = 0u; light_i < lights.n_directional_lights; light_i++) {
+        let light = &lights.directional_lights[light_i];
+        let neg_LdotV = dot((*light).direction_to_light, ray_dir_ws);
+        let angle_to_sun = fast_acos(clamp(neg_LdotV, -1.0, 1.0));
+        let w = max(0.5 * fwidth(angle_to_sun), 1e-6);
+        let sun_angular_size = (*light).sun_disk_angular_size;
+        let sun_intensity = (*light).sun_disk_intensity;
+        if sun_angular_size > 0.0 && sun_intensity > 0.0 {
+            let factor = 1 - smoothstep(sun_angular_size * 0.5 - w, sun_angular_size * 0.5 + w, angle_to_sun);
+            let sun_solid_angle = (sun_angular_size * sun_angular_size) * 0.25 * PI;
+            sun_radiance += ((*light).color.rgb / sun_solid_angle) * sun_intensity * factor * shadow_factor;
+        }
+    }
+    return sun_radiance;
+}
+
+// TRANSFORM UTILITIES
+
+fn max_atmosphere_distance(r: f32, mu: f32) -> f32 {
+    let t_top = distance_to_top_atmosphere_boundary(r, mu);
+    let t_bottom = distance_to_bottom_atmosphere_boundary(r, mu);
+    let hits = ray_intersects_ground(r, mu);
+    return mix(t_top, t_bottom, f32(hits));
+}
+
+/// Returns the observer's position in the atmosphere
+fn get_view_position() -> vec3<f32> {
+    var world_pos = view.world_position * settings.scene_units_to_m + vec3(0.0, atmosphere.bottom_radius, 0.0);
+    
+    // If the camera is underground, clamp it to the ground surface along the local up.
+    let r = length(world_pos);
+    // Nudge r above ground to avoid sqrt cancellation, zero-length segments where 
+    // r is equal to bottom_radius, which show up as black pixels
+    let min_radius = atmosphere.bottom_radius + EPSILON;
+    if r < min_radius {
+        let up = normalize(world_pos);
+        world_pos = up * min_radius;
+    }
+
+    return world_pos;
+}
+
+// We assume the `up` vector at the view position is the y axis, since the world is locally flat/level.
+// t = distance along view ray in atmosphere space
+// NOTE: this means that if your world is actually spherical, this will be wrong.
+fn get_local_up(r: f32, t: f32, ray_dir: vec3<f32>) -> vec3<f32> {
+    return normalize(vec3(0.0, r, 0.0) + t * ray_dir);
+}
+
+// Given a ray starting at radius r, with mu = cos(zenith angle),
+// and a t = distance along the ray, gives the new radius at point t
+fn get_local_r(r: f32, mu: f32, t: f32) -> f32 {
+    return sqrt(t * t + 2.0 * r * mu * t + r * r);
+}
+
+// Convert uv [0.0 .. 1.0] coordinate to ndc space xy [-1.0 .. 1.0]
+fn uv_to_ndc(uv: vec2<f32>) -> vec2<f32> {
+    return uv * vec2(2.0, -2.0) + vec2(-1.0, 1.0);
+}
+
+/// Convert ndc space xy coordinate [-1.0 .. 1.0] to uv [0.0 .. 1.0]
+fn ndc_to_uv(ndc: vec2<f32>) -> vec2<f32> {
+    return ndc * vec2(0.5, -0.5) + vec2(0.5);
+}
+
+/// Converts a direction in world space to atmosphere space
+fn direction_world_to_atmosphere(dir_ws: vec3<f32>, up: vec3<f32>) -> vec3<f32> {
+    // Camera forward in world space (-Z in view to world transform)
+    let forward_ws = (view.world_from_view * vec4(0.0, 0.0, -1.0, 0.0)).xyz;
+    let tangent_z = normalize(up * dot(forward_ws, up) - forward_ws);
+    let tangent_x = cross(up, tangent_z);
+    return vec3(
+        dot(dir_ws, tangent_x),
+        dot(dir_ws, up),
+        dot(dir_ws, tangent_z),
+    );
+}
+
+/// Converts a direction in atmosphere space to world space
+fn direction_atmosphere_to_world(dir_as: vec3<f32>) -> vec3<f32> {
+    let dir_ws = atmosphere_transforms.world_from_atmosphere * vec4(dir_as, 0.0);
+    return dir_ws.xyz;
+}
+
+// Modified from skybox.wgsl. For this pass we don't need to apply a separate sky transform or consider camera viewport.
+// Returns a normalized ray direction in world space.
+fn uv_to_ray_direction(uv: vec2<f32>) -> vec3<f32> {
+    // Using world positions of the fragment and camera to calculate a ray direction
+    // breaks down at large translations. This code only needs to know the ray direction.
+    // The ray direction is along the direction from the camera to the fragment position.
+    // In view space, the camera is at the origin, so the view space ray direction is
+    // along the direction of the fragment position - (0,0,0) which is just the
+    // fragment position.
+    // Use the position on the near clipping plane to avoid -inf world position
+    // because the far plane of an infinite reverse projection is at infinity.
+    let view_position_homogeneous = view.view_from_clip * vec4(
+        uv_to_ndc(uv),
+        1.0,
+        1.0,
+    );
+
+    let view_ray_direction = view_position_homogeneous.xyz / view_position_homogeneous.w;
+    // Transforming the view space ray direction by the inverse view matrix, transforms the
+    // direction to world space. Note that the w element is set to 0.0, as this is a
+    // vector direction, not a position, That causes the matrix multiplication to ignore
+    // the translations from the view matrix.
+    let ray_direction = (view.world_from_view * vec4(view_ray_direction, 0.0)).xyz;
+
+    return normalize(ray_direction);
+}
+
+fn zenith_azimuth_to_ray_dir(zenith: f32, azimuth: f32) -> vec3<f32> {
+    let sin_zenith = sin(zenith);
+    let mu = cos(zenith);
+    let sin_azimuth = sin(azimuth);
+    let cos_azimuth = cos(azimuth);
+    return vec3(sin_azimuth * sin_zenith, mu, -cos_azimuth * sin_zenith);
+}
+
+struct RaymarchSegment {
+    start: f32,
+    end: f32,
+}
+
+fn get_raymarch_segment(r: f32, mu: f32) -> RaymarchSegment {
+    // Get both intersection points with atmosphere
+    let atmosphere_intersections = ray_sphere_intersect(r, mu, atmosphere.top_radius);
+    let ground_intersections = ray_sphere_intersect(r, mu, atmosphere.bottom_radius);
+    
+    var segment: RaymarchSegment;
+    
+    if r < atmosphere.bottom_radius {
+        // Inside planet - start from bottom of atmosphere
+        segment.start = ground_intersections.y; // Use second intersection point with ground
+        segment.end = atmosphere_intersections.y;
+    } else if r < atmosphere.top_radius {
+        // Inside atmosphere
+        segment.start = 0.0;
+        segment.end = select(atmosphere_intersections.y, ground_intersections.x, ray_intersects_ground(r, mu));
+    } else {
+        // Outside atmosphere
+        if atmosphere_intersections.x < 0.0 {
+            // No intersection with atmosphere
+            return segment;
+        }
+        // Start at atmosphere entry, end at exit or ground
+        segment.start = atmosphere_intersections.x;
+        segment.end = select(atmosphere_intersections.y, ground_intersections.x, ray_intersects_ground(r, mu));
+    }
+
+    return segment;
+}
+
+struct RaymarchResult {
+    inscattering: vec3<f32>,
+    transmittance: vec3<f32>,
+}
+
+fn raymarch_atmosphere(
+    pos: vec3<f32>,
+    ray_dir: vec3<f32>,
+    t_max: f32,
+    max_samples: u32,
+    uv: vec2<f32>,
+    ground: bool
+) -> RaymarchResult {
+    let r = length(pos);
+    let up = normalize(pos);
+    let mu = dot(ray_dir, up);
+    
+    // Optimization: Reduce sample count at close proximity to the scene
+    let sample_count = mix(1.0, f32(max_samples), saturate(t_max * 0.01));
+    
+    let segment = get_raymarch_segment(r, mu);
+    let t_start = segment.start;
+    var t_end = segment.end;
+    
+    t_end = min(t_end, t_max);
+    let t_total = t_end - t_start;
+    
+    var result: RaymarchResult;
+    result.inscattering = vec3(0.0);
+    result.transmittance = vec3(1.0);
+    
+    // Skip if invalid segment
+    if t_total <= 0.0 {
+        return result;
+    }
+    
+    var prev_t = t_start;
+    var optical_depth = vec3(0.0);
+    for (var s = 0.0; s < sample_count; s += 1.0) {
+        // Linear distribution from atmosphere entry to exit/ground
+        let t_i = t_start + t_total * (s + MIDPOINT_RATIO) / sample_count;
+        let dt_i = (t_i - prev_t);
+        prev_t = t_i;
+
+        let sample_pos = pos + ray_dir * t_i;
+        let local_r = length(sample_pos);
+        let local_up = normalize(sample_pos);
+        let local_atmosphere = sample_atmosphere(local_r);
+
+        let sample_optical_depth = local_atmosphere.extinction * dt_i;
+        optical_depth += sample_optical_depth;
+        let sample_transmittance = exp(-sample_optical_depth);
+
+        let inscattering = sample_local_inscattering(
+            local_atmosphere,
+            ray_dir,
+            sample_pos
+        );
+
+        let s_int = (inscattering - inscattering * sample_transmittance) / local_atmosphere.extinction;
+        result.inscattering += result.transmittance * s_int;
+
+        result.transmittance *= sample_transmittance;
+        if all(result.transmittance < vec3(0.001)) {
+            break;
+        }
+    }
+
+    // include reflected luminance from planet ground 
+    if ground && ray_intersects_ground(r, mu) {
+        for (var light_i: u32 = 0u; light_i < lights.n_directional_lights; light_i++) {
+            let light = &lights.directional_lights[light_i];
+            let light_dir = (*light).direction_to_light;
+            let light_color = (*light).color.rgb;
+            let transmittance_to_ground = exp(-optical_depth);
+            // position on the sphere and get the sphere normal (up)
+            let sphere_point = pos + ray_dir * t_end;
+            let sphere_normal = normalize(sphere_point);
+            let mu_light = dot(light_dir, sphere_normal);
+            let transmittance_to_light = sample_transmittance_lut(0.0, mu_light);
+            let light_luminance = transmittance_to_light * max(mu_light, 0.0) * light_color;
+            // Normalized Lambert BRDF
+            let ground_luminance = transmittance_to_ground * atmosphere.ground_albedo / PI;
+            result.inscattering += ground_luminance * light_luminance;
+        }
+    }
+
+    return result;
+}

crates/libmarathon/src/render/pbr/atmosphere/mod.rs

@@ -0,0 +1,510 @@
+//! Procedural Atmospheric Scattering.
+//!
+//! This plugin implements [Hillaire's 2020 paper](https://sebh.github.io/publications/egsr2020.pdf)
+//! on real-time atmospheric scattering. While it *will* work simply as a
+//! procedural skybox, it also does much more. It supports dynamic time-of-
+//! -day, multiple directional lights, and since it's applied as a post-processing
+//! effect *on top* of the existing skybox, a starry skybox would automatically
+//! show based on the time of day. Scattering in front of terrain (similar
+//! to distance fog, but more complex) is handled as well, and takes into
+//! account the directional light color and direction.
+//!
+//! Adding the [`Atmosphere`] component to a 3d camera will enable the effect,
+//! which by default is set to look similar to Earth's atmosphere. See the
+//! documentation on the component itself for information regarding its fields.
+//!
+//! Performance-wise, the effect should be fairly cheap since the LUTs (Look
+//! Up Tables) that encode most of the data are small, and take advantage of the
+//! fact that the atmosphere is symmetric. Performance is also proportional to
+//! the number of directional lights in the scene. In order to tune
+//! performance more finely, the [`AtmosphereSettings`] camera component
+//! manages the size of each LUT and the sample count for each ray.
+//!
+//! Given how similar it is to [`crate::volumetric_fog`], it might be expected
+//! that these two modules would work together well. However for now using both
+//! at once is untested, and might not be physically accurate. These may be
+//! integrated into a single module in the future.
+//!
+//! On web platforms, atmosphere rendering will look slightly different. Specifically, when calculating how light travels
+//! through the atmosphere, we use a simpler averaging technique instead of the more
+//! complex blending operations. This difference will be resolved for WebGPU in a future release.
+//!
+//! [Shadertoy]: https://www.shadertoy.com/view/slSXRW
+//!
+//! [Unreal Engine Implementation]: https://github.com/sebh/UnrealEngineSkyAtmosphere
+
+mod environment;
+mod node;
+pub mod resources;
+
+use bevy_app::{App, Plugin, Update};
+use bevy_asset::embedded_asset;
+use bevy_camera::Camera3d;
+use crate::render::core_3d::graph::Node3d;
+use bevy_ecs::{
+    component::Component,
+    query::{Changed, QueryItem, With},
+    schedule::IntoScheduleConfigs,
+    system::{lifetimeless::Read, Query},
+};
+use bevy_math::{UVec2, UVec3, Vec3};
+use bevy_reflect::{std_traits::ReflectDefault, Reflect};
+use crate::render::{
+    extract_component::UniformComponentPlugin,
+    render_resource::{DownlevelFlags, ShaderType, SpecializedRenderPipelines},
+    view::Hdr,
+    RenderStartup,
+};
+use crate::render::{
+    extract_component::{ExtractComponent, ExtractComponentPlugin},
+    render_graph::{RenderGraphExt, ViewNodeRunner},
+    render_resource::{TextureFormat, TextureUsages},
+    renderer::RenderAdapter,
+    Render, RenderApp, RenderSystems,
+};
+
+use crate::render::core_3d::graph::Core3d;
+use bevy_shader::load_shader_library;
+use environment::{
+    init_atmosphere_probe_layout, init_atmosphere_probe_pipeline,
+    prepare_atmosphere_probe_bind_groups, prepare_atmosphere_probe_components,
+    prepare_probe_textures, AtmosphereEnvironmentMap, EnvironmentNode,
+};
+use resources::{
+    prepare_atmosphere_transforms, queue_render_sky_pipelines, AtmosphereTransforms,
+    RenderSkyBindGroupLayouts,
+};
+use tracing::warn;
+
+use self::{
+    node::{AtmosphereLutsNode, AtmosphereNode, RenderSkyNode},
+    resources::{
+        prepare_atmosphere_bind_groups, prepare_atmosphere_textures, AtmosphereBindGroupLayouts,
+        AtmosphereLutPipelines, AtmosphereSamplers,
+    },
+};
+
+#[doc(hidden)]
+pub struct AtmospherePlugin;
+
+impl Plugin for AtmospherePlugin {
+    fn build(&self, app: &mut App) {
+        load_shader_library!(app, "types.wgsl");
+        load_shader_library!(app, "functions.wgsl");
+        load_shader_library!(app, "bruneton_functions.wgsl");
+        load_shader_library!(app, "bindings.wgsl");
+
+        embedded_asset!(app, "transmittance_lut.wgsl");
+        embedded_asset!(app, "multiscattering_lut.wgsl");
+        embedded_asset!(app, "sky_view_lut.wgsl");
+        embedded_asset!(app, "aerial_view_lut.wgsl");
+        embedded_asset!(app, "render_sky.wgsl");
+        embedded_asset!(app, "environment.wgsl");
+
+        app.add_plugins((
+            ExtractComponentPlugin::<Atmosphere>::default(),
+            ExtractComponentPlugin::<GpuAtmosphereSettings>::default(),
+            ExtractComponentPlugin::<AtmosphereEnvironmentMap>::default(),
+            UniformComponentPlugin::<Atmosphere>::default(),
+            UniformComponentPlugin::<GpuAtmosphereSettings>::default(),
+        ))
+        .add_systems(Update, prepare_atmosphere_probe_components);
+    }
+
+    fn finish(&self, app: &mut App) {
+        let Some(render_app) = app.get_sub_app_mut(RenderApp) else {
+            return;
+        };
+
+        let render_adapter = render_app.world().resource::<RenderAdapter>();
+
+        if !render_adapter
+            .get_downlevel_capabilities()
+            .flags
+            .contains(DownlevelFlags::COMPUTE_SHADERS)
+        {
+            warn!("AtmospherePlugin not loaded. GPU lacks support for compute shaders.");
+            return;
+        }
+
+        if !render_adapter
+            .get_texture_format_features(TextureFormat::Rgba16Float)
+            .allowed_usages
+            .contains(TextureUsages::STORAGE_BINDING)
+        {
+            warn!("AtmospherePlugin not loaded. GPU lacks support: TextureFormat::Rgba16Float does not support TextureUsages::STORAGE_BINDING.");
+            return;
+        }
+
+        render_app
+            .init_resource::<AtmosphereBindGroupLayouts>()
+            .init_resource::<RenderSkyBindGroupLayouts>()
+            .init_resource::<AtmosphereSamplers>()
+            .init_resource::<AtmosphereLutPipelines>()
+            .init_resource::<AtmosphereTransforms>()
+            .init_resource::<SpecializedRenderPipelines<RenderSkyBindGroupLayouts>>()
+            .add_systems(
+                RenderStartup,
+                (init_atmosphere_probe_layout, init_atmosphere_probe_pipeline).chain(),
+            )
+            .add_systems(
+                Render,
+                (
+                    configure_camera_depth_usages.in_set(RenderSystems::ManageViews),
+                    queue_render_sky_pipelines.in_set(RenderSystems::Queue),
+                    prepare_atmosphere_textures.in_set(RenderSystems::PrepareResources),
+                    prepare_probe_textures
+                        .in_set(RenderSystems::PrepareResources)
+                        .after(prepare_atmosphere_textures),
+                    prepare_atmosphere_probe_bind_groups.in_set(RenderSystems::PrepareBindGroups),
+                    prepare_atmosphere_transforms.in_set(RenderSystems::PrepareResources),
+                    prepare_atmosphere_bind_groups.in_set(RenderSystems::PrepareBindGroups),
+                ),
+            )
+            .add_render_graph_node::<ViewNodeRunner<AtmosphereLutsNode>>(
+                Core3d,
+                AtmosphereNode::RenderLuts,
+            )
+            .add_render_graph_edges(
+                Core3d,
+                (
+                    // END_PRE_PASSES -> RENDER_LUTS -> MAIN_PASS
+                    Node3d::EndPrepasses,
+                    AtmosphereNode::RenderLuts,
+                    Node3d::StartMainPass,
+                ),
+            )
+            .add_render_graph_node::<ViewNodeRunner<RenderSkyNode>>(
+                Core3d,
+                AtmosphereNode::RenderSky,
+            )
+            .add_render_graph_node::<EnvironmentNode>(Core3d, AtmosphereNode::Environment)
+            .add_render_graph_edges(
+                Core3d,
+                (
+                    Node3d::MainOpaquePass,
+                    AtmosphereNode::RenderSky,
+                    Node3d::MainTransparentPass,
+                ),
+            );
+    }
+}
+
+/// This component describes the atmosphere of a planet, and when added to a camera
+/// will enable atmospheric scattering for that camera. This is only compatible with
+/// HDR cameras.
+///
+/// Most atmospheric particles scatter and absorb light in two main ways:
+///
+/// Rayleigh scattering occurs among very small particles, like individual gas
+/// molecules. It's wavelength dependent, and causes colors to separate out as
+/// light travels through the atmosphere. These particles *don't* absorb light.
+///
+/// Mie scattering occurs among slightly larger particles, like dust and sea spray.
+/// These particles *do* absorb light, but Mie scattering and absorption is
+/// *wavelength independent*.
+///
+/// Ozone acts differently from the other two, and is special-cased because
+/// it's very important to the look of Earth's atmosphere. It's wavelength
+/// dependent, but only *absorbs* light. Also, while the density of particles
+/// participating in Rayleigh and Mie scattering falls off roughly exponentially
+/// from the planet's surface, ozone only exists in a band centered at a fairly
+/// high altitude.
+#[derive(Clone, Component, Reflect, ShaderType)]
+#[require(AtmosphereSettings, Hdr)]
+#[reflect(Clone, Default)]
+pub struct Atmosphere {
+    /// Radius of the planet
+    ///
+    /// units: m
+    pub bottom_radius: f32,
+
+    /// Radius at which we consider the atmosphere to 'end' for our
+    /// calculations (from center of planet)
+    ///
+    /// units: m
+    pub top_radius: f32,
+
+    /// An approximation of the average albedo (or color, roughly) of the
+    /// planet's surface. This is used when calculating multiscattering.
+    ///
+    /// units: N/A
+    pub ground_albedo: Vec3,
+
+    /// The rate of falloff of rayleigh particulate with respect to altitude:
+    /// optical density = exp(-rayleigh_density_exp_scale * altitude in meters).
+    ///
+    /// THIS VALUE MUST BE POSITIVE
+    ///
+    /// units: N/A
+    pub rayleigh_density_exp_scale: f32,
+
+    /// The scattering optical density of rayleigh particulate, or how
+    /// much light it scatters per meter
+    ///
+    /// units: m^-1
+    pub rayleigh_scattering: Vec3,
+
+    /// The rate of falloff of mie particulate with respect to altitude:
+    /// optical density = exp(-mie_density_exp_scale * altitude in meters)
+    ///
+    /// THIS VALUE MUST BE POSITIVE
+    ///
+    /// units: N/A
+    pub mie_density_exp_scale: f32,
+
+    /// The scattering optical density of mie particulate, or how much light
+    /// it scatters per meter.
+    ///
+    /// units: m^-1
+    pub mie_scattering: f32,
+
+    /// The absorbing optical density of mie particulate, or how much light
+    /// it absorbs per meter.
+    ///
+    /// units: m^-1
+    pub mie_absorption: f32,
+
+    /// The "asymmetry" of mie scattering, or how much light tends to scatter
+    /// forwards, rather than backwards or to the side.
+    ///
+    /// domain: (-1, 1)
+    /// units: N/A
+    pub mie_asymmetry: f32, //the "asymmetry" value of the phase function, unitless. Domain: (-1, 1)
+
+    /// The altitude at which the ozone layer is centered.
+    ///
+    /// units: m
+    pub ozone_layer_altitude: f32,
+
+    /// The width of the ozone layer
+    ///
+    /// units: m
+    pub ozone_layer_width: f32,
+
+    /// The optical density of ozone, or how much of each wavelength of
+    /// light it absorbs per meter.
+    ///
+    /// units: m^-1
+    pub ozone_absorption: Vec3,
+}
+
+impl Atmosphere {
+    pub const EARTH: Atmosphere = Atmosphere {
+        bottom_radius: 6_360_000.0,
+        top_radius: 6_460_000.0,
+        ground_albedo: Vec3::splat(0.3),
+        rayleigh_density_exp_scale: 1.0 / 8_000.0,
+        rayleigh_scattering: Vec3::new(5.802e-6, 13.558e-6, 33.100e-6),
+        mie_density_exp_scale: 1.0 / 1_200.0,
+        mie_scattering: 3.996e-6,
+        mie_absorption: 0.444e-6,
+        mie_asymmetry: 0.8,
+        ozone_layer_altitude: 25_000.0,
+        ozone_layer_width: 30_000.0,
+        ozone_absorption: Vec3::new(0.650e-6, 1.881e-6, 0.085e-6),
+    };
+
+    pub fn with_density_multiplier(mut self, mult: f32) -> Self {
+        self.rayleigh_scattering *= mult;
+        self.mie_scattering *= mult;
+        self.mie_absorption *= mult;
+        self.ozone_absorption *= mult;
+        self
+    }
+}
+
+impl Default for Atmosphere {
+    fn default() -> Self {
+        Self::EARTH
+    }
+}
+
+impl ExtractComponent for Atmosphere {
+    type QueryData = Read<Atmosphere>;
+
+    type QueryFilter = With<Camera3d>;
+
+    type Out = Atmosphere;
+
+    fn extract_component(item: QueryItem<'_, '_, Self::QueryData>) -> Option<Self::Out> {
+        Some(item.clone())
+    }
+}
+
+/// This component controls the resolution of the atmosphere LUTs, and
+/// how many samples are used when computing them.
+///
+/// The transmittance LUT stores the transmittance from a point in the
+/// atmosphere to the outer edge of the atmosphere in any direction,
+/// parametrized by the point's radius and the cosine of the zenith angle
+/// of the ray.
+///
+/// The multiscattering LUT stores the factor representing luminance scattered
+/// towards the camera with scattering order >2, parametrized by the point's radius
+/// and the cosine of the zenith angle of the sun.
+///
+/// The sky-view lut is essentially the actual skybox, storing the light scattered
+/// towards the camera in every direction with a cubemap.
+///
+/// The aerial-view lut is a 3d LUT fit to the view frustum, which stores the luminance
+/// scattered towards the camera at each point (RGB channels), alongside the average
+/// transmittance to that point (A channel).
+#[derive(Clone, Component, Reflect)]
+#[reflect(Clone, Default)]
+pub struct AtmosphereSettings {
+    /// The size of the transmittance LUT
+    pub transmittance_lut_size: UVec2,
+
+    /// The size of the multiscattering LUT
+    pub multiscattering_lut_size: UVec2,
+
+    /// The size of the sky-view LUT.
+    pub sky_view_lut_size: UVec2,
+
+    /// The size of the aerial-view LUT.
+    pub aerial_view_lut_size: UVec3,
+
+    /// The number of points to sample along each ray when
+    /// computing the transmittance LUT
+    pub transmittance_lut_samples: u32,
+
+    /// The number of rays to sample when computing each
+    /// pixel of the multiscattering LUT
+    pub multiscattering_lut_dirs: u32,
+
+    /// The number of points to sample when integrating along each
+    /// multiscattering ray
+    pub multiscattering_lut_samples: u32,
+
+    /// The number of points to sample along each ray when
+    /// computing the sky-view LUT.
+    pub sky_view_lut_samples: u32,
+
+    /// The number of points to sample for each slice along the z-axis
+    /// of the aerial-view LUT.
+    pub aerial_view_lut_samples: u32,
+
+    /// The maximum distance from the camera to evaluate the
+    /// aerial view LUT. The slices along the z-axis of the
+    /// texture will be distributed linearly from the camera
+    /// to this value.
+    ///
+    /// units: m
+    pub aerial_view_lut_max_distance: f32,
+
+    /// A conversion factor between scene units and meters, used to
+    /// ensure correctness at different length scales.
+    pub scene_units_to_m: f32,
+
+    /// The number of points to sample for each fragment when the using
+    /// ray marching to render the sky
+    pub sky_max_samples: u32,
+
+    /// The rendering method to use for the atmosphere.
+    pub rendering_method: AtmosphereMode,
+}
+
+impl Default for AtmosphereSettings {
+    fn default() -> Self {
+        Self {
+            transmittance_lut_size: UVec2::new(256, 128),
+            transmittance_lut_samples: 40,
+            multiscattering_lut_size: UVec2::new(32, 32),
+            multiscattering_lut_dirs: 64,
+            multiscattering_lut_samples: 20,
+            sky_view_lut_size: UVec2::new(400, 200),
+            sky_view_lut_samples: 16,
+            aerial_view_lut_size: UVec3::new(32, 32, 32),
+            aerial_view_lut_samples: 10,
+            aerial_view_lut_max_distance: 3.2e4,
+            scene_units_to_m: 1.0,
+            sky_max_samples: 16,
+            rendering_method: AtmosphereMode::LookupTexture,
+        }
+    }
+}
+
+#[derive(Clone, Component, Reflect, ShaderType)]
+#[reflect(Default)]
+pub struct GpuAtmosphereSettings {
+    pub transmittance_lut_size: UVec2,
+    pub multiscattering_lut_size: UVec2,
+    pub sky_view_lut_size: UVec2,
+    pub aerial_view_lut_size: UVec3,
+    pub transmittance_lut_samples: u32,
+    pub multiscattering_lut_dirs: u32,
+    pub multiscattering_lut_samples: u32,
+    pub sky_view_lut_samples: u32,
+    pub aerial_view_lut_samples: u32,
+    pub aerial_view_lut_max_distance: f32,
+    pub scene_units_to_m: f32,
+    pub sky_max_samples: u32,
+    pub rendering_method: u32,
+}
+
+impl Default for GpuAtmosphereSettings {
+    fn default() -> Self {
+        AtmosphereSettings::default().into()
+    }
+}
+
+impl From<AtmosphereSettings> for GpuAtmosphereSettings {
+    fn from(s: AtmosphereSettings) -> Self {
+        Self {
+            transmittance_lut_size: s.transmittance_lut_size,
+            multiscattering_lut_size: s.multiscattering_lut_size,
+            sky_view_lut_size: s.sky_view_lut_size,
+            aerial_view_lut_size: s.aerial_view_lut_size,
+            transmittance_lut_samples: s.transmittance_lut_samples,
+            multiscattering_lut_dirs: s.multiscattering_lut_dirs,
+            multiscattering_lut_samples: s.multiscattering_lut_samples,
+            sky_view_lut_samples: s.sky_view_lut_samples,
+            aerial_view_lut_samples: s.aerial_view_lut_samples,
+            aerial_view_lut_max_distance: s.aerial_view_lut_max_distance,
+            scene_units_to_m: s.scene_units_to_m,
+            sky_max_samples: s.sky_max_samples,
+            rendering_method: s.rendering_method as u32,
+        }
+    }
+}
+
+impl ExtractComponent for GpuAtmosphereSettings {
+    type QueryData = Read<AtmosphereSettings>;
+
+    type QueryFilter = (With<Camera3d>, With<Atmosphere>);
+
+    type Out = GpuAtmosphereSettings;
+
+    fn extract_component(item: QueryItem<'_, '_, Self::QueryData>) -> Option<Self::Out> {
+        Some(item.clone().into())
+    }
+}
+
+fn configure_camera_depth_usages(
+    mut cameras: Query<&mut Camera3d, (Changed<Camera3d>, With<Atmosphere>)>,
+) {
+    for mut camera in &mut cameras {
+        camera.depth_texture_usages.0 |= TextureUsages::TEXTURE_BINDING.bits();
+    }
+}
+
+/// Selects how the atmosphere is rendered. Choose based on scene scale and
+/// volumetric shadow quality, and based on performance needs.
+#[repr(u32)]
+#[derive(Clone, Default, Reflect, Copy)]
+pub enum AtmosphereMode {
+    /// High-performance solution tailored to scenes that are mostly inside of the atmosphere.
+    /// Uses a set of lookup textures to approximate scattering integration.
+    /// Slightly less accurate for very long-distance/space views (lighting precision
+    /// tapers as the camera moves far from the scene origin) and for sharp volumetric
+    /// (cloud/fog) shadows.
+    #[default]
+    LookupTexture = 0,
+    /// Slower, more accurate rendering method for any type of scene.
+    /// Integrates the scattering numerically with raymarching and produces sharp volumetric
+    /// (cloud/fog) shadows.
+    /// Best for cinematic shots, planets seen from orbit, and scenes requiring
+    /// accurate long-distance lighting.
+    Raymarched = 1,
+}

crates/libmarathon/src/render/pbr/atmosphere/multiscattering_lut.wgsl

@@ -0,0 +1,139 @@
+#import bevy_pbr::{
+    mesh_view_types::{Lights, DirectionalLight},
+    atmosphere::{
+        types::{Atmosphere, AtmosphereSettings},
+        bindings::{atmosphere, settings},
+        functions::{
+            multiscattering_lut_uv_to_r_mu, sample_transmittance_lut,
+            get_local_r, get_local_up, sample_atmosphere, FRAC_4_PI,
+            max_atmosphere_distance, rayleigh, henyey_greenstein,
+            zenith_azimuth_to_ray_dir,
+        },
+        bruneton_functions::{
+            distance_to_top_atmosphere_boundary, distance_to_bottom_atmosphere_boundary, ray_intersects_ground
+        }
+    }
+}
+
+#import bevy_render::maths::{PI,PI_2}
+
+const PHI_2: vec2<f32> = vec2(1.3247179572447460259609088, 1.7548776662466927600495087);
+
+@group(0) @binding(13) var multiscattering_lut_out: texture_storage_2d<rgba16float, write>;
+
+fn s2_sequence(n: u32) -> vec2<f32> {
+    return fract(0.5 + f32(n) * PHI_2);
+}
+
+// Lambert equal-area projection. 
+fn uv_to_sphere(uv: vec2<f32>) -> vec3<f32> {
+    let phi = PI_2 * uv.y;
+    let sin_lambda = 2 * uv.x - 1;
+    let cos_lambda = sqrt(1 - sin_lambda * sin_lambda);
+
+    return vec3(cos_lambda * cos(phi), cos_lambda * sin(phi), sin_lambda);
+}
+
+// Shared memory arrays for workgroup communication
+var<workgroup> multi_scat_shared_mem: array<vec3<f32>, 64>;
+var<workgroup> l_shared_mem: array<vec3<f32>, 64>;
+
+@compute 
+@workgroup_size(1, 1, 64)
+fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
+    var uv = (vec2<f32>(global_id.xy) + 0.5) / vec2<f32>(settings.multiscattering_lut_size);
+
+    let r_mu = multiscattering_lut_uv_to_r_mu(uv);
+    let light_dir = normalize(vec3(0.0, r_mu.y, -1.0));
+
+    let ray_dir = uv_to_sphere(s2_sequence(global_id.z));
+    let ms_sample = sample_multiscattering_dir(r_mu.x, ray_dir, light_dir);
+    
+    // Calculate the contribution for this sample
+    let sphere_solid_angle = 4.0 * PI;
+    let sample_weight = sphere_solid_angle / 64.0;
+    multi_scat_shared_mem[global_id.z] = ms_sample.f_ms * sample_weight;
+    l_shared_mem[global_id.z] = ms_sample.l_2 * sample_weight;
+
+    workgroupBarrier();
+
+    // Parallel reduction bitshift to the right to divide by 2 each step
+    for (var step = 32u; step > 0u; step >>= 1u) {
+        if global_id.z < step {
+            multi_scat_shared_mem[global_id.z] += multi_scat_shared_mem[global_id.z + step];
+            l_shared_mem[global_id.z] += l_shared_mem[global_id.z + step];
+        }
+        workgroupBarrier();
+    }
+
+    if global_id.z > 0u {
+        return;
+    }
+
+    // Apply isotropic phase function
+    let f_ms = multi_scat_shared_mem[0] * FRAC_4_PI;
+    let l_2 = l_shared_mem[0] * FRAC_4_PI;
+    
+    // Equation 10 from the paper: Geometric series for infinite scattering
+    let psi_ms = l_2 / (1.0 - f_ms);
+    textureStore(multiscattering_lut_out, global_id.xy, vec4<f32>(psi_ms, 1.0));
+}
+
+struct MultiscatteringSample {
+    l_2: vec3<f32>,
+    f_ms: vec3<f32>,
+};
+
+fn sample_multiscattering_dir(r: f32, ray_dir: vec3<f32>, light_dir: vec3<f32>) -> MultiscatteringSample {
+    // get the cosine of the zenith angle of the view direction with respect to the light direction
+    let mu_view = ray_dir.y;
+    let t_max = max_atmosphere_distance(r, mu_view);
+
+    let dt = t_max / f32(settings.multiscattering_lut_samples);
+    var optical_depth = vec3<f32>(0.0);
+
+    var l_2 = vec3(0.0);
+    var f_ms = vec3(0.0);
+    var throughput = vec3(1.0);
+    for (var i: u32 = 0u; i < settings.multiscattering_lut_samples; i++) {
+        let t_i = dt * (f32(i) + 0.5);
+        let local_r = get_local_r(r, mu_view, t_i);
+        let local_up = get_local_up(r, t_i, ray_dir);
+
+        let local_atmosphere = sample_atmosphere(local_r);
+        let sample_optical_depth = local_atmosphere.extinction * dt;
+        let sample_transmittance = exp(-sample_optical_depth);
+        optical_depth += sample_optical_depth;
+
+        let mu_light = dot(light_dir, local_up);
+        let scattering_no_phase = local_atmosphere.rayleigh_scattering + local_atmosphere.mie_scattering;
+
+        let ms = scattering_no_phase;
+        let ms_int = (ms - ms * sample_transmittance) / local_atmosphere.extinction;
+        f_ms += throughput * ms_int;
+
+        let transmittance_to_light = sample_transmittance_lut(local_r, mu_light);
+        let shadow_factor = transmittance_to_light * f32(!ray_intersects_ground(local_r, mu_light));
+
+        let s = scattering_no_phase * shadow_factor * FRAC_4_PI;
+        let s_int = (s - s * sample_transmittance) / local_atmosphere.extinction;
+        l_2 += throughput * s_int;
+
+        throughput *= sample_transmittance;
+        if all(throughput < vec3(0.001)) {
+            break;
+        }
+    }
+
+    //include reflected luminance from planet ground 
+    if ray_intersects_ground(r, mu_view) {
+        let transmittance_to_ground = exp(-optical_depth);
+        let local_up = get_local_up(r, t_max, ray_dir);
+        let mu_light = dot(light_dir, local_up);
+        let transmittance_to_light = sample_transmittance_lut(0.0, mu_light);
+        let ground_luminance = transmittance_to_light * transmittance_to_ground * max(mu_light, 0.0) * atmosphere.ground_albedo;
+        l_2 += ground_luminance;
+    }
+
+    return MultiscatteringSample(l_2, f_ms);
+}

crates/libmarathon/src/render/pbr/atmosphere/node.rs

@@ -0,0 +1,233 @@
+use bevy_ecs::{query::QueryItem, system::lifetimeless::Read, world::World};
+use bevy_math::{UVec2, Vec3Swizzles};
+use crate::render::{
+    diagnostic::RecordDiagnostics,
+    extract_component::DynamicUniformIndex,
+    render_graph::{NodeRunError, RenderGraphContext, RenderLabel, ViewNode},
+    render_resource::{ComputePass, ComputePassDescriptor, PipelineCache, RenderPassDescriptor},
+    renderer::RenderContext,
+    view::{ViewTarget, ViewUniformOffset},
+};
+
+use crate::render::pbr::ViewLightsUniformOffset;
+
+use super::{
+    resources::{
+        AtmosphereBindGroups, AtmosphereLutPipelines, AtmosphereTransformsOffset,
+        RenderSkyPipelineId,
+    },
+    Atmosphere, GpuAtmosphereSettings,
+};
+
+#[derive(PartialEq, Eq, Debug, Copy, Clone, Hash, RenderLabel)]
+pub enum AtmosphereNode {
+    RenderLuts,
+    RenderSky,
+    Environment,
+}
+
+#[derive(Default)]
+pub(super) struct AtmosphereLutsNode {}
+
+impl ViewNode for AtmosphereLutsNode {
+    type ViewQuery = (
+        Read<GpuAtmosphereSettings>,
+        Read<AtmosphereBindGroups>,
+        Read<DynamicUniformIndex<Atmosphere>>,
+        Read<DynamicUniformIndex<GpuAtmosphereSettings>>,
+        Read<AtmosphereTransformsOffset>,
+        Read<ViewUniformOffset>,
+        Read<ViewLightsUniformOffset>,
+    );
+
+    fn run(
+        &self,
+        _graph: &mut RenderGraphContext,
+        render_context: &mut RenderContext,
+        (
+            settings,
+            bind_groups,
+            atmosphere_uniforms_offset,
+            settings_uniforms_offset,
+            atmosphere_transforms_offset,
+            view_uniforms_offset,
+            lights_uniforms_offset,
+        ): QueryItem<Self::ViewQuery>,
+        world: &World,
+    ) -> Result<(), NodeRunError> {
+        let pipelines = world.resource::<AtmosphereLutPipelines>();
+        let pipeline_cache = world.resource::<PipelineCache>();
+        let (
+            Some(transmittance_lut_pipeline),
+            Some(multiscattering_lut_pipeline),
+            Some(sky_view_lut_pipeline),
+            Some(aerial_view_lut_pipeline),
+        ) = (
+            pipeline_cache.get_compute_pipeline(pipelines.transmittance_lut),
+            pipeline_cache.get_compute_pipeline(pipelines.multiscattering_lut),
+            pipeline_cache.get_compute_pipeline(pipelines.sky_view_lut),
+            pipeline_cache.get_compute_pipeline(pipelines.aerial_view_lut),
+        )
+        else {
+            return Ok(());
+        };
+
+        let diagnostics = render_context.diagnostic_recorder();
+
+        let command_encoder = render_context.command_encoder();
+
+        let mut luts_pass = command_encoder.begin_compute_pass(&ComputePassDescriptor {
+            label: Some("atmosphere_luts"),
+            timestamp_writes: None,
+        });
+        let pass_span = diagnostics.pass_span(&mut luts_pass, "atmosphere_luts");
+
+        fn dispatch_2d(compute_pass: &mut ComputePass, size: UVec2) {
+            const WORKGROUP_SIZE: u32 = 16;
+            let workgroups_x = size.x.div_ceil(WORKGROUP_SIZE);
+            let workgroups_y = size.y.div_ceil(WORKGROUP_SIZE);
+            compute_pass.dispatch_workgroups(workgroups_x, workgroups_y, 1);
+        }
+
+        // Transmittance LUT
+
+        luts_pass.set_pipeline(transmittance_lut_pipeline);
+        luts_pass.set_bind_group(
+            0,
+            &bind_groups.transmittance_lut,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+            ],
+        );
+
+        dispatch_2d(&mut luts_pass, settings.transmittance_lut_size);
+
+        // Multiscattering LUT
+
+        luts_pass.set_pipeline(multiscattering_lut_pipeline);
+        luts_pass.set_bind_group(
+            0,
+            &bind_groups.multiscattering_lut,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+            ],
+        );
+
+        luts_pass.dispatch_workgroups(
+            settings.multiscattering_lut_size.x,
+            settings.multiscattering_lut_size.y,
+            1,
+        );
+
+        // Sky View LUT
+
+        luts_pass.set_pipeline(sky_view_lut_pipeline);
+        luts_pass.set_bind_group(
+            0,
+            &bind_groups.sky_view_lut,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+                atmosphere_transforms_offset.index(),
+                view_uniforms_offset.offset,
+                lights_uniforms_offset.offset,
+            ],
+        );
+
+        dispatch_2d(&mut luts_pass, settings.sky_view_lut_size);
+
+        // Aerial View LUT
+
+        luts_pass.set_pipeline(aerial_view_lut_pipeline);
+        luts_pass.set_bind_group(
+            0,
+            &bind_groups.aerial_view_lut,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+                view_uniforms_offset.offset,
+                lights_uniforms_offset.offset,
+            ],
+        );
+
+        dispatch_2d(&mut luts_pass, settings.aerial_view_lut_size.xy());
+
+        pass_span.end(&mut luts_pass);
+
+        Ok(())
+    }
+}
+
+#[derive(Default)]
+pub(super) struct RenderSkyNode;
+
+impl ViewNode for RenderSkyNode {
+    type ViewQuery = (
+        Read<AtmosphereBindGroups>,
+        Read<ViewTarget>,
+        Read<DynamicUniformIndex<Atmosphere>>,
+        Read<DynamicUniformIndex<GpuAtmosphereSettings>>,
+        Read<AtmosphereTransformsOffset>,
+        Read<ViewUniformOffset>,
+        Read<ViewLightsUniformOffset>,
+        Read<RenderSkyPipelineId>,
+    );
+
+    fn run<'w>(
+        &self,
+        _graph: &mut RenderGraphContext,
+        render_context: &mut RenderContext<'w>,
+        (
+            atmosphere_bind_groups,
+            view_target,
+            atmosphere_uniforms_offset,
+            settings_uniforms_offset,
+            atmosphere_transforms_offset,
+            view_uniforms_offset,
+            lights_uniforms_offset,
+            render_sky_pipeline_id,
+        ): QueryItem<'w, '_, Self::ViewQuery>,
+        world: &'w World,
+    ) -> Result<(), NodeRunError> {
+        let pipeline_cache = world.resource::<PipelineCache>();
+        let Some(render_sky_pipeline) =
+            pipeline_cache.get_render_pipeline(render_sky_pipeline_id.0)
+        else {
+            return Ok(());
+        }; //TODO: warning
+
+        let diagnostics = render_context.diagnostic_recorder();
+
+        let mut render_sky_pass =
+            render_context
+                .command_encoder()
+                .begin_render_pass(&RenderPassDescriptor {
+                    label: Some("render_sky"),
+                    color_attachments: &[Some(view_target.get_color_attachment())],
+                    depth_stencil_attachment: None,
+                    timestamp_writes: None,
+                    occlusion_query_set: None,
+                });
+        let pass_span = diagnostics.pass_span(&mut render_sky_pass, "render_sky");
+
+        render_sky_pass.set_pipeline(render_sky_pipeline);
+        render_sky_pass.set_bind_group(
+            0,
+            &atmosphere_bind_groups.render_sky,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+                atmosphere_transforms_offset.index(),
+                view_uniforms_offset.offset,
+                lights_uniforms_offset.offset,
+            ],
+        );
+        render_sky_pass.draw(0..3, 0..1);
+
+        pass_span.end(&mut render_sky_pass);
+
+        Ok(())
+    }
+}

crates/libmarathon/src/render/pbr/atmosphere/render_sky.wgsl

@@ -0,0 +1,82 @@
+enable dual_source_blending;
+
+#import bevy_pbr::atmosphere::{
+    types::{Atmosphere, AtmosphereSettings},
+    bindings::{atmosphere, view, atmosphere_transforms, settings},
+    functions::{
+        sample_transmittance_lut, sample_transmittance_lut_segment,
+        sample_sky_view_lut, direction_world_to_atmosphere,
+        uv_to_ray_direction, uv_to_ndc, sample_aerial_view_lut,
+        sample_sun_radiance, ndc_to_camera_dist, raymarch_atmosphere, 
+        get_view_position, max_atmosphere_distance
+    },
+};
+#import bevy_render::view::View;
+
+#import bevy_core_pipeline::fullscreen_vertex_shader::FullscreenVertexOutput
+
+#ifdef MULTISAMPLED
+@group(0) @binding(13) var depth_texture: texture_depth_multisampled_2d;
+#else
+@group(0) @binding(13) var depth_texture: texture_depth_2d;
+#endif
+
+struct RenderSkyOutput {
+#ifdef DUAL_SOURCE_BLENDING
+    @location(0) @blend_src(0) inscattering: vec4<f32>,
+    @location(0) @blend_src(1) transmittance: vec4<f32>,
+#else
+    @location(0) inscattering: vec4<f32>,
+#endif
+}
+
+@fragment
+fn main(in: FullscreenVertexOutput) -> RenderSkyOutput {
+    let depth = textureLoad(depth_texture, vec2<i32>(in.position.xy), 0);
+
+    let ray_dir_ws = uv_to_ray_direction(in.uv);
+    let world_pos = get_view_position();
+    let r = length(world_pos);
+    let up = normalize(world_pos);
+    let mu = dot(ray_dir_ws, up);
+    let max_samples = settings.sky_max_samples;
+    let should_raymarch = settings.rendering_method == 1u;
+
+    var transmittance: vec3<f32>;
+    var inscattering: vec3<f32>;
+
+    let sun_radiance = sample_sun_radiance(ray_dir_ws);
+
+    if depth == 0.0 {
+        let ray_dir_as = direction_world_to_atmosphere(ray_dir_ws, up);
+        transmittance = sample_transmittance_lut(r, mu);
+        inscattering = sample_sky_view_lut(r, ray_dir_as);
+        if should_raymarch {
+            let t_max = max_atmosphere_distance(r, mu);
+            let result = raymarch_atmosphere(world_pos, ray_dir_ws, t_max, max_samples, in.uv, true);
+            inscattering = result.inscattering;
+            transmittance = result.transmittance;
+        }
+        inscattering += sun_radiance * transmittance;
+    } else {
+        let t = ndc_to_camera_dist(vec3(uv_to_ndc(in.uv), depth));
+        inscattering = sample_aerial_view_lut(in.uv, t);
+        transmittance = sample_transmittance_lut_segment(r, mu, t);
+        if should_raymarch {
+            let result = raymarch_atmosphere(world_pos, ray_dir_ws, t, max_samples, in.uv, false);
+            inscattering = result.inscattering;
+            transmittance = result.transmittance;
+        }
+    }
+
+    // exposure compensation
+    inscattering *= view.exposure;
+    
+#ifdef DUAL_SOURCE_BLENDING
+    return RenderSkyOutput(vec4(inscattering, 0.0), vec4(transmittance, 1.0));
+#else
+    let mean_transmittance = (transmittance.r + transmittance.g + transmittance.b) / 3.0;
+    return RenderSkyOutput(vec4(inscattering, mean_transmittance));
+#endif
+    
+}

crates/libmarathon/src/render/pbr/atmosphere/resources.rs

@@ -0,0 +1,700 @@
+use crate::render::pbr::{GpuLights, LightMeta};
+use bevy_asset::{load_embedded_asset, Handle};
+use bevy_camera::{Camera, Camera3d};
+use crate::render::FullscreenShader;
+use bevy_ecs::{
+    component::Component,
+    entity::Entity,
+    query::With,
+    resource::Resource,
+    system::{Commands, Query, Res, ResMut},
+    world::{FromWorld, World},
+};
+use bevy_image::ToExtents;
+use bevy_math::{Affine3A, Mat4, Vec3A};
+use crate::render::{
+    extract_component::ComponentUniforms,
+    render_resource::{binding_types::*, *},
+    renderer::{RenderDevice, RenderQueue},
+    texture::{CachedTexture, TextureCache},
+    view::{ExtractedView, Msaa, ViewDepthTexture, ViewUniform, ViewUniforms},
+};
+use bevy_shader::Shader;
+use bevy_utils::default;
+
+use super::{Atmosphere, GpuAtmosphereSettings};
+
+#[derive(Resource)]
+pub(crate) struct AtmosphereBindGroupLayouts {
+    pub transmittance_lut: BindGroupLayout,
+    pub multiscattering_lut: BindGroupLayout,
+    pub sky_view_lut: BindGroupLayout,
+    pub aerial_view_lut: BindGroupLayout,
+}
+
+#[derive(Resource)]
+pub(crate) struct RenderSkyBindGroupLayouts {
+    pub render_sky: BindGroupLayout,
+    pub render_sky_msaa: BindGroupLayout,
+    pub fullscreen_shader: FullscreenShader,
+    pub fragment_shader: Handle<Shader>,
+}
+
+impl FromWorld for AtmosphereBindGroupLayouts {
+    fn from_world(world: &mut World) -> Self {
+        let render_device = world.resource::<RenderDevice>();
+        let transmittance_lut = render_device.create_bind_group_layout(
+            "transmittance_lut_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::COMPUTE,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<GpuAtmosphereSettings>(true)),
+                    (
+                        // transmittance lut storage texture
+                        13,
+                        texture_storage_2d(
+                            TextureFormat::Rgba16Float,
+                            StorageTextureAccess::WriteOnly,
+                        ),
+                    ),
+                ),
+            ),
+        );
+
+        let multiscattering_lut = render_device.create_bind_group_layout(
+            "multiscattering_lut_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::COMPUTE,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<GpuAtmosphereSettings>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (
+                        //multiscattering lut storage texture
+                        13,
+                        texture_storage_2d(
+                            TextureFormat::Rgba16Float,
+                            StorageTextureAccess::WriteOnly,
+                        ),
+                    ),
+                ),
+            ),
+        );
+
+        let sky_view_lut = render_device.create_bind_group_layout(
+            "sky_view_lut_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::COMPUTE,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<GpuAtmosphereSettings>(true)),
+                    (2, uniform_buffer::<AtmosphereTransform>(true)),
+                    (3, uniform_buffer::<ViewUniform>(true)),
+                    (4, uniform_buffer::<GpuLights>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (7, texture_2d(TextureSampleType::Float { filterable: true })), //multiscattering lut and sampler
+                    (8, sampler(SamplerBindingType::Filtering)),
+                    (
+                        13,
+                        texture_storage_2d(
+                            TextureFormat::Rgba16Float,
+                            StorageTextureAccess::WriteOnly,
+                        ),
+                    ),
+                ),
+            ),
+        );
+
+        let aerial_view_lut = render_device.create_bind_group_layout(
+            "aerial_view_lut_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::COMPUTE,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<GpuAtmosphereSettings>(true)),
+                    (3, uniform_buffer::<ViewUniform>(true)),
+                    (4, uniform_buffer::<GpuLights>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (7, texture_2d(TextureSampleType::Float { filterable: true })), //multiscattering lut and sampler
+                    (8, sampler(SamplerBindingType::Filtering)),
+                    (
+                        //Aerial view lut storage texture
+                        13,
+                        texture_storage_3d(
+                            TextureFormat::Rgba16Float,
+                            StorageTextureAccess::WriteOnly,
+                        ),
+                    ),
+                ),
+            ),
+        );
+
+        Self {
+            transmittance_lut,
+            multiscattering_lut,
+            sky_view_lut,
+            aerial_view_lut,
+        }
+    }
+}
+
+impl FromWorld for RenderSkyBindGroupLayouts {
+    fn from_world(world: &mut World) -> Self {
+        let render_device = world.resource::<RenderDevice>();
+        let render_sky = render_device.create_bind_group_layout(
+            "render_sky_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::FRAGMENT,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<GpuAtmosphereSettings>(true)),
+                    (2, uniform_buffer::<AtmosphereTransform>(true)),
+                    (3, uniform_buffer::<ViewUniform>(true)),
+                    (4, uniform_buffer::<GpuLights>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (7, texture_2d(TextureSampleType::Float { filterable: true })), //multiscattering lut and sampler
+                    (8, sampler(SamplerBindingType::Filtering)),
+                    (9, texture_2d(TextureSampleType::Float { filterable: true })), //sky view lut and sampler
+                    (10, sampler(SamplerBindingType::Filtering)),
+                    (
+                        // aerial view lut and sampler
+                        11,
+                        texture_3d(TextureSampleType::Float { filterable: true }),
+                    ),
+                    (12, sampler(SamplerBindingType::Filtering)),
+                    (
+                        //view depth texture
+                        13,
+                        texture_2d(TextureSampleType::Depth),
+                    ),
+                ),
+            ),
+        );
+
+        let render_sky_msaa = render_device.create_bind_group_layout(
+            "render_sky_msaa_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::FRAGMENT,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<GpuAtmosphereSettings>(true)),
+                    (2, uniform_buffer::<AtmosphereTransform>(true)),
+                    (3, uniform_buffer::<ViewUniform>(true)),
+                    (4, uniform_buffer::<GpuLights>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (7, texture_2d(TextureSampleType::Float { filterable: true })), //multiscattering lut and sampler
+                    (8, sampler(SamplerBindingType::Filtering)),
+                    (9, texture_2d(TextureSampleType::Float { filterable: true })), //sky view lut and sampler
+                    (10, sampler(SamplerBindingType::Filtering)),
+                    (
+                        // aerial view lut and sampler
+                        11,
+                        texture_3d(TextureSampleType::Float { filterable: true }),
+                    ),
+                    (12, sampler(SamplerBindingType::Filtering)),
+                    (
+                        //view depth texture
+                        13,
+                        texture_2d_multisampled(TextureSampleType::Depth),
+                    ),
+                ),
+            ),
+        );
+
+        Self {
+            render_sky,
+            render_sky_msaa,
+            fullscreen_shader: world.resource::<FullscreenShader>().clone(),
+            fragment_shader: load_embedded_asset!(world, "render_sky.wgsl"),
+        }
+    }
+}
+
+#[derive(Resource)]
+pub struct AtmosphereSamplers {
+    pub transmittance_lut: Sampler,
+    pub multiscattering_lut: Sampler,
+    pub sky_view_lut: Sampler,
+    pub aerial_view_lut: Sampler,
+}
+
+impl FromWorld for AtmosphereSamplers {
+    fn from_world(world: &mut World) -> Self {
+        let render_device = world.resource::<RenderDevice>();
+
+        let base_sampler = SamplerDescriptor {
+            mag_filter: FilterMode::Linear,
+            min_filter: FilterMode::Linear,
+            mipmap_filter: FilterMode::Nearest,
+            ..Default::default()
+        };
+
+        let transmittance_lut = render_device.create_sampler(&SamplerDescriptor {
+            label: Some("transmittance_lut_sampler"),
+            ..base_sampler
+        });
+
+        let multiscattering_lut = render_device.create_sampler(&SamplerDescriptor {
+            label: Some("multiscattering_lut_sampler"),
+            ..base_sampler
+        });
+
+        let sky_view_lut = render_device.create_sampler(&SamplerDescriptor {
+            label: Some("sky_view_lut_sampler"),
+            address_mode_u: AddressMode::Repeat,
+            ..base_sampler
+        });
+
+        let aerial_view_lut = render_device.create_sampler(&SamplerDescriptor {
+            label: Some("aerial_view_lut_sampler"),
+            ..base_sampler
+        });
+
+        Self {
+            transmittance_lut,
+            multiscattering_lut,
+            sky_view_lut,
+            aerial_view_lut,
+        }
+    }
+}
+
+#[derive(Resource)]
+pub(crate) struct AtmosphereLutPipelines {
+    pub transmittance_lut: CachedComputePipelineId,
+    pub multiscattering_lut: CachedComputePipelineId,
+    pub sky_view_lut: CachedComputePipelineId,
+    pub aerial_view_lut: CachedComputePipelineId,
+}
+
+impl FromWorld for AtmosphereLutPipelines {
+    fn from_world(world: &mut World) -> Self {
+        let pipeline_cache = world.resource::<PipelineCache>();
+        let layouts = world.resource::<AtmosphereBindGroupLayouts>();
+
+        let transmittance_lut = pipeline_cache.queue_compute_pipeline(ComputePipelineDescriptor {
+            label: Some("transmittance_lut_pipeline".into()),
+            layout: vec![layouts.transmittance_lut.clone()],
+            shader: load_embedded_asset!(world, "transmittance_lut.wgsl"),
+            ..default()
+        });
+
+        let multiscattering_lut =
+            pipeline_cache.queue_compute_pipeline(ComputePipelineDescriptor {
+                label: Some("multi_scattering_lut_pipeline".into()),
+                layout: vec![layouts.multiscattering_lut.clone()],
+                shader: load_embedded_asset!(world, "multiscattering_lut.wgsl"),
+                ..default()
+            });
+
+        let sky_view_lut = pipeline_cache.queue_compute_pipeline(ComputePipelineDescriptor {
+            label: Some("sky_view_lut_pipeline".into()),
+            layout: vec![layouts.sky_view_lut.clone()],
+            shader: load_embedded_asset!(world, "sky_view_lut.wgsl"),
+            ..default()
+        });
+
+        let aerial_view_lut = pipeline_cache.queue_compute_pipeline(ComputePipelineDescriptor {
+            label: Some("aerial_view_lut_pipeline".into()),
+            layout: vec![layouts.aerial_view_lut.clone()],
+            shader: load_embedded_asset!(world, "aerial_view_lut.wgsl"),
+            ..default()
+        });
+
+        Self {
+            transmittance_lut,
+            multiscattering_lut,
+            sky_view_lut,
+            aerial_view_lut,
+        }
+    }
+}
+
+#[derive(Component)]
+pub(crate) struct RenderSkyPipelineId(pub CachedRenderPipelineId);
+
+#[derive(Copy, Clone, Hash, PartialEq, Eq)]
+pub(crate) struct RenderSkyPipelineKey {
+    pub msaa_samples: u32,
+    pub dual_source_blending: bool,
+}
+
+impl SpecializedRenderPipeline for RenderSkyBindGroupLayouts {
+    type Key = RenderSkyPipelineKey;
+
+    fn specialize(&self, key: Self::Key) -> RenderPipelineDescriptor {
+        let mut shader_defs = Vec::new();
+
+        if key.msaa_samples > 1 {
+            shader_defs.push("MULTISAMPLED".into());
+        }
+        if key.dual_source_blending {
+            shader_defs.push("DUAL_SOURCE_BLENDING".into());
+        }
+
+        let dst_factor = if key.dual_source_blending {
+            BlendFactor::Src1
+        } else {
+            BlendFactor::SrcAlpha
+        };
+
+        RenderPipelineDescriptor {
+            label: Some(format!("render_sky_pipeline_{}", key.msaa_samples).into()),
+            layout: vec![if key.msaa_samples == 1 {
+                self.render_sky.clone()
+            } else {
+                self.render_sky_msaa.clone()
+            }],
+            vertex: self.fullscreen_shader.to_vertex_state(),
+            fragment: Some(FragmentState {
+                shader: self.fragment_shader.clone(),
+                shader_defs,
+                targets: vec![Some(ColorTargetState {
+                    format: TextureFormat::Rgba16Float,
+                    blend: Some(BlendState {
+                        color: BlendComponent {
+                            src_factor: BlendFactor::One,
+                            dst_factor,
+                            operation: BlendOperation::Add,
+                        },
+                        alpha: BlendComponent {
+                            src_factor: BlendFactor::Zero,
+                            dst_factor: BlendFactor::One,
+                            operation: BlendOperation::Add,
+                        },
+                    }),
+                    write_mask: ColorWrites::ALL,
+                })],
+                ..default()
+            }),
+            multisample: MultisampleState {
+                count: key.msaa_samples,
+                ..default()
+            },
+            ..default()
+        }
+    }
+}
+
+pub(super) fn queue_render_sky_pipelines(
+    views: Query<(Entity, &Msaa), (With<Camera>, With<Atmosphere>)>,
+    pipeline_cache: Res<PipelineCache>,
+    layouts: Res<RenderSkyBindGroupLayouts>,
+    mut specializer: ResMut<SpecializedRenderPipelines<RenderSkyBindGroupLayouts>>,
+    render_device: Res<RenderDevice>,
+    mut commands: Commands,
+) {
+    for (entity, msaa) in &views {
+        let id = specializer.specialize(
+            &pipeline_cache,
+            &layouts,
+            RenderSkyPipelineKey {
+                msaa_samples: msaa.samples(),
+                dual_source_blending: render_device
+                    .features()
+                    .contains(WgpuFeatures::DUAL_SOURCE_BLENDING),
+            },
+        );
+        commands.entity(entity).insert(RenderSkyPipelineId(id));
+    }
+}
+
+#[derive(Component)]
+pub struct AtmosphereTextures {
+    pub transmittance_lut: CachedTexture,
+    pub multiscattering_lut: CachedTexture,
+    pub sky_view_lut: CachedTexture,
+    pub aerial_view_lut: CachedTexture,
+}
+
+pub(super) fn prepare_atmosphere_textures(
+    views: Query<(Entity, &GpuAtmosphereSettings), With<Atmosphere>>,
+    render_device: Res<RenderDevice>,
+    mut texture_cache: ResMut<TextureCache>,
+    mut commands: Commands,
+) {
+    for (entity, lut_settings) in &views {
+        let transmittance_lut = texture_cache.get(
+            &render_device,
+            TextureDescriptor {
+                label: Some("transmittance_lut"),
+                size: lut_settings.transmittance_lut_size.to_extents(),
+                mip_level_count: 1,
+                sample_count: 1,
+                dimension: TextureDimension::D2,
+                format: TextureFormat::Rgba16Float,
+                usage: TextureUsages::STORAGE_BINDING | TextureUsages::TEXTURE_BINDING,
+                view_formats: &[],
+            },
+        );
+
+        let multiscattering_lut = texture_cache.get(
+            &render_device,
+            TextureDescriptor {
+                label: Some("multiscattering_lut"),
+                size: lut_settings.multiscattering_lut_size.to_extents(),
+                mip_level_count: 1,
+                sample_count: 1,
+                dimension: TextureDimension::D2,
+                format: TextureFormat::Rgba16Float,
+                usage: TextureUsages::STORAGE_BINDING | TextureUsages::TEXTURE_BINDING,
+                view_formats: &[],
+            },
+        );
+
+        let sky_view_lut = texture_cache.get(
+            &render_device,
+            TextureDescriptor {
+                label: Some("sky_view_lut"),
+                size: lut_settings.sky_view_lut_size.to_extents(),
+                mip_level_count: 1,
+                sample_count: 1,
+                dimension: TextureDimension::D2,
+                format: TextureFormat::Rgba16Float,
+                usage: TextureUsages::STORAGE_BINDING | TextureUsages::TEXTURE_BINDING,
+                view_formats: &[],
+            },
+        );
+
+        let aerial_view_lut = texture_cache.get(
+            &render_device,
+            TextureDescriptor {
+                label: Some("aerial_view_lut"),
+                size: lut_settings.aerial_view_lut_size.to_extents(),
+                mip_level_count: 1,
+                sample_count: 1,
+                dimension: TextureDimension::D3,
+                format: TextureFormat::Rgba16Float,
+                usage: TextureUsages::STORAGE_BINDING | TextureUsages::TEXTURE_BINDING,
+                view_formats: &[],
+            },
+        );
+
+        commands.entity(entity).insert({
+            AtmosphereTextures {
+                transmittance_lut,
+                multiscattering_lut,
+                sky_view_lut,
+                aerial_view_lut,
+            }
+        });
+    }
+}
+
+#[derive(Resource, Default)]
+pub struct AtmosphereTransforms {
+    uniforms: DynamicUniformBuffer<AtmosphereTransform>,
+}
+
+impl AtmosphereTransforms {
+    #[inline]
+    pub fn uniforms(&self) -> &DynamicUniformBuffer<AtmosphereTransform> {
+        &self.uniforms
+    }
+}
+
+#[derive(ShaderType)]
+pub struct AtmosphereTransform {
+    world_from_atmosphere: Mat4,
+}
+
+#[derive(Component)]
+pub struct AtmosphereTransformsOffset {
+    index: u32,
+}
+
+impl AtmosphereTransformsOffset {
+    #[inline]
+    pub fn index(&self) -> u32 {
+        self.index
+    }
+}
+
+pub(super) fn prepare_atmosphere_transforms(
+    views: Query<(Entity, &ExtractedView), (With<Atmosphere>, With<Camera3d>)>,
+    render_device: Res<RenderDevice>,
+    render_queue: Res<RenderQueue>,
+    mut atmo_uniforms: ResMut<AtmosphereTransforms>,
+    mut commands: Commands,
+) {
+    let atmo_count = views.iter().len();
+    let Some(mut writer) =
+        atmo_uniforms
+            .uniforms
+            .get_writer(atmo_count, &render_device, &render_queue)
+    else {
+        return;
+    };
+
+    for (entity, view) in &views {
+        let world_from_view = view.world_from_view.affine();
+        let camera_z = world_from_view.matrix3.z_axis;
+        let camera_y = world_from_view.matrix3.y_axis;
+        let atmo_z = camera_z
+            .with_y(0.0)
+            .try_normalize()
+            .unwrap_or_else(|| camera_y.with_y(0.0).normalize());
+        let atmo_y = Vec3A::Y;
+        let atmo_x = atmo_y.cross(atmo_z).normalize();
+        let world_from_atmosphere =
+            Affine3A::from_cols(atmo_x, atmo_y, atmo_z, world_from_view.translation);
+
+        let world_from_atmosphere = Mat4::from(world_from_atmosphere);
+
+        commands.entity(entity).insert(AtmosphereTransformsOffset {
+            index: writer.write(&AtmosphereTransform {
+                world_from_atmosphere,
+            }),
+        });
+    }
+}
+
+#[derive(Component)]
+pub(crate) struct AtmosphereBindGroups {
+    pub transmittance_lut: BindGroup,
+    pub multiscattering_lut: BindGroup,
+    pub sky_view_lut: BindGroup,
+    pub aerial_view_lut: BindGroup,
+    pub render_sky: BindGroup,
+}
+
+pub(super) fn prepare_atmosphere_bind_groups(
+    views: Query<
+        (Entity, &AtmosphereTextures, &ViewDepthTexture, &Msaa),
+        (With<Camera3d>, With<Atmosphere>),
+    >,
+    render_device: Res<RenderDevice>,
+    layouts: Res<AtmosphereBindGroupLayouts>,
+    render_sky_layouts: Res<RenderSkyBindGroupLayouts>,
+    samplers: Res<AtmosphereSamplers>,
+    view_uniforms: Res<ViewUniforms>,
+    lights_uniforms: Res<LightMeta>,
+    atmosphere_transforms: Res<AtmosphereTransforms>,
+    atmosphere_uniforms: Res<ComponentUniforms<Atmosphere>>,
+    settings_uniforms: Res<ComponentUniforms<GpuAtmosphereSettings>>,
+
+    mut commands: Commands,
+) {
+    if views.iter().len() == 0 {
+        return;
+    }
+
+    let atmosphere_binding = atmosphere_uniforms
+        .binding()
+        .expect("Failed to prepare atmosphere bind groups. Atmosphere uniform buffer missing");
+
+    let transforms_binding = atmosphere_transforms
+        .uniforms()
+        .binding()
+        .expect("Failed to prepare atmosphere bind groups. Atmosphere transforms buffer missing");
+
+    let settings_binding = settings_uniforms.binding().expect(
+        "Failed to prepare atmosphere bind groups. AtmosphereSettings uniform buffer missing",
+    );
+
+    let view_binding = view_uniforms
+        .uniforms
+        .binding()
+        .expect("Failed to prepare atmosphere bind groups. View uniform buffer missing");
+
+    let lights_binding = lights_uniforms
+        .view_gpu_lights
+        .binding()
+        .expect("Failed to prepare atmosphere bind groups. Lights uniform buffer missing");
+
+    for (entity, textures, view_depth_texture, msaa) in &views {
+        let transmittance_lut = render_device.create_bind_group(
+            "transmittance_lut_bind_group",
+            &layouts.transmittance_lut,
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (13, &textures.transmittance_lut.default_view),
+            )),
+        );
+
+        let multiscattering_lut = render_device.create_bind_group(
+            "multiscattering_lut_bind_group",
+            &layouts.multiscattering_lut,
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (5, &textures.transmittance_lut.default_view),
+                (6, &samplers.transmittance_lut),
+                (13, &textures.multiscattering_lut.default_view),
+            )),
+        );
+
+        let sky_view_lut = render_device.create_bind_group(
+            "sky_view_lut_bind_group",
+            &layouts.sky_view_lut,
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (2, transforms_binding.clone()),
+                (3, view_binding.clone()),
+                (4, lights_binding.clone()),
+                (5, &textures.transmittance_lut.default_view),
+                (6, &samplers.transmittance_lut),
+                (7, &textures.multiscattering_lut.default_view),
+                (8, &samplers.multiscattering_lut),
+                (13, &textures.sky_view_lut.default_view),
+            )),
+        );
+
+        let aerial_view_lut = render_device.create_bind_group(
+            "sky_view_lut_bind_group",
+            &layouts.aerial_view_lut,
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (3, view_binding.clone()),
+                (4, lights_binding.clone()),
+                (5, &textures.transmittance_lut.default_view),
+                (6, &samplers.transmittance_lut),
+                (7, &textures.multiscattering_lut.default_view),
+                (8, &samplers.multiscattering_lut),
+                (13, &textures.aerial_view_lut.default_view),
+            )),
+        );
+
+        let render_sky = render_device.create_bind_group(
+            "render_sky_bind_group",
+            if *msaa == Msaa::Off {
+                &render_sky_layouts.render_sky
+            } else {
+                &render_sky_layouts.render_sky_msaa
+            },
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (2, transforms_binding.clone()),
+                (3, view_binding.clone()),
+                (4, lights_binding.clone()),
+                (5, &textures.transmittance_lut.default_view),
+                (6, &samplers.transmittance_lut),
+                (7, &textures.multiscattering_lut.default_view),
+                (8, &samplers.multiscattering_lut),
+                (9, &textures.sky_view_lut.default_view),
+                (10, &samplers.sky_view_lut),
+                (11, &textures.aerial_view_lut.default_view),
+                (12, &samplers.aerial_view_lut),
+                (13, view_depth_texture.view()),
+            )),
+        );
+
+        commands.entity(entity).insert(AtmosphereBindGroups {
+            transmittance_lut,
+            multiscattering_lut,
+            sky_view_lut,
+            aerial_view_lut,
+            render_sky,
+        });
+    }
+}

crates/libmarathon/src/render/pbr/atmosphere/sky_view_lut.wgsl

@@ -0,0 +1,44 @@
+#import bevy_pbr::{
+    mesh_view_types::Lights,
+    atmosphere::{
+        types::{Atmosphere, AtmosphereSettings},
+        bindings::{atmosphere, view, settings},
+        functions::{
+            sample_atmosphere, AtmosphereSample,
+            sample_local_inscattering, get_view_position,
+            max_atmosphere_distance, direction_atmosphere_to_world,
+            sky_view_lut_uv_to_zenith_azimuth, zenith_azimuth_to_ray_dir,
+            MIDPOINT_RATIO, raymarch_atmosphere, EPSILON
+        },
+    }
+}
+
+#import bevy_render::{
+    view::View,
+    maths::HALF_PI,
+}
+#import bevy_core_pipeline::fullscreen_vertex_shader::FullscreenVertexOutput
+
+@group(0) @binding(13) var sky_view_lut_out: texture_storage_2d<rgba16float, write>;
+
+@compute
+@workgroup_size(16, 16, 1)
+fn main(@builtin(global_invocation_id) idx: vec3<u32>) {
+    let uv = vec2<f32>(idx.xy) / vec2<f32>(settings.sky_view_lut_size);
+
+    let cam_pos = get_view_position();
+    let r = length(cam_pos);
+    var zenith_azimuth = sky_view_lut_uv_to_zenith_azimuth(r, uv);
+
+    let ray_dir_as = zenith_azimuth_to_ray_dir(zenith_azimuth.x, zenith_azimuth.y);
+    let ray_dir_ws = direction_atmosphere_to_world(ray_dir_as);
+
+    let world_pos = vec3(0.0, r, 0.0);
+    let up = normalize(world_pos);
+    let mu = dot(ray_dir_ws, up);
+    let t_max = max_atmosphere_distance(r, mu);
+
+    let result = raymarch_atmosphere(world_pos, ray_dir_ws, t_max, settings.sky_view_lut_samples, uv, true);
+
+    textureStore(sky_view_lut_out, idx.xy, vec4(result.inscattering, 1.0));
+}

crates/libmarathon/src/render/pbr/atmosphere/transmittance_lut.wgsl

@@ -0,0 +1,48 @@
+#import bevy_pbr::atmosphere::{
+    types::{Atmosphere, AtmosphereSettings},
+    bindings::{settings, atmosphere},
+    functions::{AtmosphereSample, sample_atmosphere, get_local_r, max_atmosphere_distance, MIDPOINT_RATIO},
+    bruneton_functions::{transmittance_lut_uv_to_r_mu, distance_to_bottom_atmosphere_boundary, distance_to_top_atmosphere_boundary},
+}
+
+
+#import bevy_core_pipeline::fullscreen_vertex_shader::FullscreenVertexOutput
+
+@group(0) @binding(13) var transmittance_lut_out: texture_storage_2d<rgba16float, write>;
+
+@compute 
+@workgroup_size(16, 16, 1)
+fn main(@builtin(global_invocation_id) idx: vec3<u32>) {
+    let uv: vec2<f32> = (vec2<f32>(idx.xy) + 0.5) / vec2<f32>(settings.transmittance_lut_size);
+    // map UV coordinates to view height (r) and zenith cos angle (mu)
+    let r_mu = transmittance_lut_uv_to_r_mu(uv);
+
+    // compute the optical depth from view height r to the top atmosphere boundary
+    let optical_depth = ray_optical_depth(r_mu.x, r_mu.y, settings.transmittance_lut_samples);
+    let transmittance = exp(-optical_depth);
+
+    textureStore(transmittance_lut_out, idx.xy, vec4(transmittance, 1.0));
+}
+
+/// Compute the optical depth of the atmosphere from the ground to the top atmosphere boundary
+/// at a given view height (r) and zenith cos angle (mu)
+fn ray_optical_depth(r: f32, mu: f32, sample_count: u32) -> vec3<f32> {
+    let t_max = max_atmosphere_distance(r, mu);
+    var optical_depth = vec3<f32>(0.0f);
+    var prev_t = 0.0f;
+
+    for (var i = 0u; i < sample_count; i++) {
+        let t_i = t_max * (f32(i) + MIDPOINT_RATIO) / f32(sample_count);
+        let dt = t_i - prev_t;
+        prev_t = t_i;
+
+        let r_i = get_local_r(r, mu, t_i);
+
+        let atmosphere_sample = sample_atmosphere(r_i);
+        let sample_optical_depth = atmosphere_sample.extinction * dt;
+
+        optical_depth += sample_optical_depth;
+    }
+
+    return optical_depth;
+}

crates/libmarathon/src/render/pbr/atmosphere/types.wgsl

@@ -0,0 +1,46 @@
+#define_import_path bevy_pbr::atmosphere::types
+
+struct Atmosphere {
+    // Radius of the planet
+    bottom_radius: f32, // units: m
+
+    // Radius at which we consider the atmosphere to 'end' for out calculations (from center of planet)
+    top_radius: f32, // units: m
+
+    ground_albedo: vec3<f32>,
+
+    rayleigh_density_exp_scale: f32,
+    rayleigh_scattering: vec3<f32>,
+
+    mie_density_exp_scale: f32,
+    mie_scattering: f32, // units: m^-1
+    mie_absorption: f32, // units: m^-1
+    mie_asymmetry: f32, // the "asymmetry" value of the phase function, unitless. Domain: (-1, 1)
+
+    ozone_layer_altitude: f32, // units: m
+    ozone_layer_width: f32, // units: m
+    ozone_absorption: vec3<f32>, // ozone absorption. units: m^-1
+}
+
+struct AtmosphereSettings {
+    transmittance_lut_size: vec2<u32>,
+    multiscattering_lut_size: vec2<u32>,
+    sky_view_lut_size: vec2<u32>,
+    aerial_view_lut_size: vec3<u32>,
+    transmittance_lut_samples: u32,
+    multiscattering_lut_dirs: u32,
+    multiscattering_lut_samples: u32,
+    sky_view_lut_samples: u32,
+    aerial_view_lut_samples: u32,
+    aerial_view_lut_max_distance: f32,
+    scene_units_to_m: f32,
+    sky_max_samples: u32,
+    rendering_method: u32,
+}
+
+
+// "Atmosphere space" is just the view position with y=0 and oriented horizontally,
+// so the horizon stays a horizontal line in our luts
+struct AtmosphereTransforms {
+    world_from_atmosphere: mat4x4<f32>,
+}