Shadows, self and nearby craft

crates/resolved-cli/src/main.rs

@@ -156,6 +156,7 @@ fn render(args: &RenderArgs) -> Result<(), String> {
         appearance: file.appearance(),
         read_radiance: false,
         profile: timings.is_some(),
+        shadow_map_size: None,
     };
     let start = Instant::now();
     let mut renderer = Renderer::new(camera, options).map_err(|e| e.to_string())?;
@@ -217,7 +218,8 @@ fn render(args: &RenderArgs) -> Result<(), String> {
                 v[v.len() / 2]
             };
             println!(
-                "gpu      bodies {:.3} ms, mesh {:.3} ms, sensor {:.3} ms, first to last {:.3} ms (medians)",
+                "gpu      shadow {:.3} ms, bodies {:.3} ms, mesh {:.3} ms, sensor {:.3} ms, first to last {:.3} ms (medians)",
+                median(|t| t.shadow),
                 median(|t| t.bodies),
                 median(|t| t.mesh),
                 median(|t| t.sensor),

crates/resolved/src/fidelity.rs

@@ -19,8 +19,10 @@ pub enum LightingTier {
     #[default]
     L0,
     /// L0, plus Earth's surface from an image (`RendererOptions::earth_texture`, NASA's
-    /// Blue Marble) with a stand-in for the atmosphere's blue. The rest of M2 (shadows,
-    /// Earthshine, the Moon's albedo map) joins this tier as it lands.
+    /// Blue Marble) with a stand-in for the atmosphere's blue, and shadows that targets
+    /// cast on themselves and on each other, with the penumbra the Sun's disk gives
+    /// them. The rest of M2 (Earthshine, the Moon's albedo map) joins this tier as it
+    /// lands.
     L1,
 }
 

crates/resolved/src/gpu/mesh.rs

@@ -4,7 +4,7 @@
 //! per-target uniform block ([`MeshDraw`]) that places it and says how much sunlight
 //! reaches it.
 
-use nalgebra::{Matrix3, Vector3};
+use nalgebra::{Matrix3, Matrix4, Vector3};
 
 use super::bodies::MAX_CHANNELS;
 use super::profile::PassTimestamps;
@@ -18,12 +18,25 @@ use crate::Error;
 
 #[repr(C)]
 #[derive(Clone, Copy, Debug, bytemuck::Pod, bytemuck::Zeroable)]
-struct Vertex {
+pub(crate) struct Vertex {
     position: [f32; 3],
     normal: [f32; 3],
     part: u32,
 }
 
+impl Vertex {
+    /// The vertex buffer layout: position, normal, part.
+    pub(crate) fn layout() -> wgpu::VertexBufferLayout<'static> {
+        const ATTRIBUTES: [wgpu::VertexAttribute; 3] =
+            wgpu::vertex_attr_array![0 => Float32x3, 1 => Float32x3, 2 => Uint32];
+        wgpu::VertexBufferLayout {
+            array_stride: std::mem::size_of::<Vertex>() as u64,
+            step_mode: wgpu::VertexStepMode::Vertex,
+            attributes: &ATTRIBUTES,
+        }
+    }
+}
+
 /// A model on the GPU, ready to draw.
 #[derive(Debug)]
 pub struct GpuMesh {
@@ -31,6 +44,40 @@ pub struct GpuMesh {
     indices: wgpu::Buffer,
     index_count: u32,
     parts: wgpu::BindGroup,
+    /// Bounding sphere center in the body frame, in meters.
+    pub center: Vector3<f64>,
+    /// Bounding sphere radius, in meters.
+    pub radius: f64,
+}
+
+impl GpuMesh {
+    /// Binds the vertex and index buffers and draws every triangle.
+    pub(crate) fn draw(&self, pass: &mut wgpu::RenderPass) {
+        pass.set_vertex_buffer(0, self.vertices.slice(..));
+        pass.set_index_buffer(self.indices.slice(..), wgpu::IndexFormat::Uint32);
+        pass.draw_indexed(0..self.index_count, 0, 0..1);
+    }
+}
+
+/// Where a target's fragments look up their shadow: its own shadow map.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct ShadowDraw {
+    /// Body frame to the map's light coordinates ([`crate::shadow::ShadowPlan`]).
+    pub light_from_body: Matrix4<f64>,
+    /// The map's array layer.
+    pub layer: u32,
+    /// Meters per shadow-map texel.
+    pub texel: f64,
+    /// The map's depth range, in meters.
+    pub depth_range: f64,
+    /// Texels across the map.
+    pub size: u32,
+}
+
+/// The tangent of the Sun's angular diameter at 1 AU: how fast a shadow's edge blurs
+/// with distance from what casts it. 0.93 % of the distance.
+fn sun_spread(sun_distance: f64) -> f64 {
+    (2.0 * (crate::radiometry::SUN_RADIUS / sun_distance).asin()).tan()
 }
 
 /// One target's uniform block; the layout matches `Draw` in `mesh.wgsl`.
@@ -43,6 +90,11 @@ pub struct MeshDraw {
     translation: [f32; 4],
     sun: [f32; 4],
     params: [f32; 4],
+    light_from_body: [[f32; 4]; 4],
+    // Layer, enabled (0 or 1), meters per texel, depth range in meters.
+    shadow: [f32; 4],
+    // Tangent of the Sun's angular diameter, normal offset in meters, map size in texels.
+    shadow_more: [f32; 4],
 }
 
 impl MeshDraw {
@@ -55,6 +107,7 @@ impl MeshDraw {
         near: f64,
         sun: Vector3<f64>,
         sun_visible: f64,
+        shadow: Option<ShadowDraw>,
     ) -> Self {
         let to_camera = camera.attitude.inverse();
         let rotation: Matrix3<f32> = (to_camera * target.attitude)
@@ -80,6 +133,23 @@ impl MeshDraw {
             translation: [translation.x, translation.y, translation.z, 0.0],
             sun: [sun_dir.x, sun_dir.y, sun_dir.z, 0.0],
             params: [(sun_visible / (au * au)) as f32, 0.0, 0.0, 0.0],
+            light_from_body: shadow.map_or([[0.0; 4]; 4], |s| {
+                let m = s.light_from_body.cast::<f32>();
+                std::array::from_fn(|c| [m[(0, c)], m[(1, c)], m[(2, c)], m[(3, c)]])
+            }),
+            shadow: shadow.map_or([0.0; 4], |s| {
+                [s.layer as f32, 1.0, s.texel as f32, s.depth_range as f32]
+            }),
+            // Offsetting the lookup a texel and a half along the normal keeps a surface
+            // from shadowing itself where it is sampled at a slant.
+            shadow_more: shadow.map_or([0.0; 4], |s| {
+                [
+                    sun_spread(to_sun.norm()) as f32,
+                    (1.5 * s.texel) as f32,
+                    s.size as f32,
+                    0.0,
+                ]
+            }),
         }
     }
 }
@@ -90,6 +160,9 @@ pub struct MeshPass {
     pipeline: wgpu::RenderPipeline,
     draw_layout: wgpu::BindGroupLayout,
     parts_layout: wgpu::BindGroupLayout,
+    shadow_layout: wgpu::BindGroupLayout,
+    /// Bound when there are no shadow maps: a 1 x 1 x 1 array the shader does not read.
+    no_shadows: wgpu::BindGroup,
     draws: wgpu::Buffer,
     draws_group: wgpu::BindGroup,
     stride: u64,
@@ -126,11 +199,40 @@ impl MeshPass {
                 count: None,
             }],
         });
+        let shadow_layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
+            label: Some("mesh shadows"),
+            entries: &[
+                wgpu::BindGroupLayoutEntry {
+                    binding: 0,
+                    visibility: wgpu::ShaderStages::FRAGMENT,
+                    ty: wgpu::BindingType::Texture {
+                        sample_type: wgpu::TextureSampleType::Depth,
+                        view_dimension: wgpu::TextureViewDimension::D2Array,
+                        multisampled: false,
+                    },
+                    count: None,
+                },
+                wgpu::BindGroupLayoutEntry {
+                    binding: 1,
+                    visibility: wgpu::ShaderStages::FRAGMENT,
+                    ty: wgpu::BindingType::Sampler(wgpu::SamplerBindingType::Comparison),
+                    count: None,
+                },
+            ],
+        });
         let layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
             label: Some("mesh"),
-            bind_group_layouts: &[Some(&draw_layout), Some(&parts_layout)],
+            bind_group_layouts: &[
+                Some(&draw_layout),
+                Some(&parts_layout),
+                Some(&shadow_layout),
+            ],
             immediate_size: 0,
         });
+        let no_shadows = {
+            let mut placeholder = super::shadow::ShadowPass::new(ctx, 1);
+            placeholder.prepare(ctx, 1, &shadow_layout)
+        };
         let module = device.create_shader_module(wgpu::ShaderModuleDescriptor {
             label: Some("mesh"),
             source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/mesh.wgsl").into()),
@@ -142,11 +244,7 @@ impl MeshPass {
                 module: &module,
                 entry_point: Some("vs_main"),
                 compilation_options: Default::default(),
-                buffers: &[Some(wgpu::VertexBufferLayout {
-                    array_stride: std::mem::size_of::<Vertex>() as u64,
-                    step_mode: wgpu::VertexStepMode::Vertex,
-                    attributes: &wgpu::vertex_attr_array![0 => Float32x3, 1 => Float32x3, 2 => Uint32],
-                })],
+                buffers: &[Some(Vertex::layout())],
             },
             // No culling: thin panels are seen from both sides.
             primitive: wgpu::PrimitiveState::default(),
@@ -174,6 +272,8 @@ impl MeshPass {
             pipeline,
             draw_layout,
             parts_layout,
+            shadow_layout,
+            no_shadows,
             draws,
             draws_group,
             stride,
@@ -234,6 +334,18 @@ impl MeshPass {
             })
             .collect();
         let parts: Vec<[f32; 4]> = mesh.parts.iter().map(|p| lambert[p.material]).collect();
+        let (mut lo, mut hi) = (Vector3::repeat(f64::MAX), Vector3::repeat(f64::MIN));
+        for p in &mesh.positions {
+            let p = Vector3::from(*p).cast::<f64>();
+            lo = lo.inf(&p);
+            hi = hi.sup(&p);
+        }
+        let center = (lo + hi) / 2.0;
+        let radius = mesh
+            .positions
+            .iter()
+            .map(|p| (Vector3::from(*p).cast::<f64>() - center).norm())
+            .fold(0.0, f64::max);
         let vertices: Vec<Vertex> = (0..mesh.positions.len())
             .map(|i| Vertex {
                 position: mesh.positions[i],
@@ -266,6 +378,8 @@ impl MeshPass {
                 wgpu::BufferUsages::INDEX,
             ),
             index_count: mesh.indices.len() as u32,
+            center,
+            radius,
             parts: ctx.device.create_bind_group(&wgpu::BindGroupDescriptor {
                 label: Some("mesh parts"),
                 layout: &self.parts_layout,
@@ -277,14 +391,21 @@ impl MeshPass {
         })
     }
 
+    /// The layout of the shadow-map bind group the shader samples (group 2).
+    pub fn shadow_layout(&self) -> &wgpu::BindGroupLayout {
+        &self.shadow_layout
+    }
+
     /// Uploads each target's uniforms and records the pass into `encoder`, drawing over
-    /// what earlier passes left in `targets`.
+    /// what earlier passes left in `targets`. `shadows` is the shadow maps' bind group
+    /// ([`super::shadow::ShadowPass::prepare`]), or `None` when no draw uses one.
     pub fn encode(
         &mut self,
         ctx: &GpuContext,
         encoder: &mut wgpu::CommandEncoder,
         targets: &FrameTargets,
         draws: &[(&GpuMesh, MeshDraw)],
+        shadows: Option<&wgpu::BindGroup>,
         timestamps: Option<PassTimestamps>,
     ) {
         if draws.len() > self.capacity {
@@ -324,12 +445,11 @@ impl MeshPass {
             ..Default::default()
         });
         pass.set_pipeline(&self.pipeline);
+        pass.set_bind_group(2, shadows.unwrap_or(&self.no_shadows), &[]);
         for (i, (mesh, _)) in draws.iter().enumerate() {
             pass.set_bind_group(0, &self.draws_group, &[(i as u64 * self.stride) as u32]);
             pass.set_bind_group(1, &mesh.parts, &[]);
-            pass.set_vertex_buffer(0, mesh.vertices.slice(..));
-            pass.set_index_buffer(mesh.indices.slice(..), wgpu::IndexFormat::Uint32);
-            pass.draw_indexed(0..mesh.index_count, 0, 0..1);
+            mesh.draw(&mut pass);
         }
     }
 }

crates/resolved/src/gpu/mod.rs

@@ -5,6 +5,7 @@ pub mod context;
 pub mod mesh;
 pub mod profile;
 pub mod sensor;
+pub mod shadow;
 pub mod targets;
 
 pub use context::{adapters, AdapterChoice, GpuContext};

crates/resolved/src/gpu/profile.rs

@@ -8,31 +8,52 @@
 use super::GpuContext;
 
 /// The passes timed, in the order they run.
-pub const PASSES: [&str; 3] = ["bodies", "mesh", "sensor"];
+pub const PASSES: [&str; 4] = ["shadow", "bodies", "mesh", "sensor"];
 
 /// Where one pass writes its begin and end timestamps.
 #[derive(Clone, Copy, Debug)]
 pub struct PassTimestamps<'a> {
     set: &'a wgpu::QuerySet,
     begin: u32,
+    write_begin: bool,
+    write_end: bool,
 }
 
 impl<'a> PassTimestamps<'a> {
+    /// The same pass's timestamps, keeping only the begin, the end, both or neither:
+    /// for one logical pass made of several GPU passes.
+    pub fn split(self, begin: bool, end: bool) -> Self {
+        Self {
+            write_begin: self.write_begin && begin,
+            write_end: self.write_end && end,
+            ..self
+        }
+    }
+
+    fn indices(&self) -> (Option<u32>, Option<u32>) {
+        (
+            self.write_begin.then_some(self.begin),
+            self.write_end.then_some(self.begin + 1),
+        )
+    }
+
     /// For a render pass descriptor.
     pub fn render(self) -> wgpu::RenderPassTimestampWrites<'a> {
+        let (begin, end) = self.indices();
         wgpu::RenderPassTimestampWrites {
             query_set: self.set,
-            beginning_of_pass_write_index: Some(self.begin),
-            end_of_pass_write_index: Some(self.begin + 1),
+            beginning_of_pass_write_index: begin,
+            end_of_pass_write_index: end,
         }
     }
 
     /// For a compute pass descriptor.
     pub fn compute(self) -> wgpu::ComputePassTimestampWrites<'a> {
+        let (begin, end) = self.indices();
         wgpu::ComputePassTimestampWrites {
             query_set: self.set,
-            beginning_of_pass_write_index: Some(self.begin),
-            end_of_pass_write_index: Some(self.begin + 1),
+            beginning_of_pass_write_index: begin,
+            end_of_pass_write_index: end,
         }
     }
 }
@@ -40,6 +61,8 @@ impl<'a> PassTimestamps<'a> {
 /// GPU time spent in each pass of one frame, in milliseconds.
 #[derive(Clone, Copy, Debug, PartialEq)]
 pub struct PassTimings {
+    /// Shadow maps (zero work below lighting L1).
+    pub shadow: f64,
     /// Earth, Moon and Sun.
     pub bodies: f64,
     /// Target meshes.
@@ -105,6 +128,8 @@ impl Profiler {
         PassTimestamps {
             set: &self.set,
             begin: 2 * index as u32,
+            write_begin: true,
+            write_end: true,
         }
     }
 
@@ -118,10 +143,11 @@ impl Profiler {
     pub fn timings(&self, ticks: &[u64]) -> PassTimings {
         let ms = |a: u64, b: u64| b.wrapping_sub(a) as f64 * self.period * 1e-6;
         PassTimings {
-            bodies: ms(ticks[0], ticks[1]),
-            mesh: ms(ticks[2], ticks[3]),
-            sensor: ms(ticks[4], ticks[5]),
-            total: ms(ticks[0], ticks[5]),
+            shadow: ms(ticks[0], ticks[1]),
+            bodies: ms(ticks[2], ticks[3]),
+            mesh: ms(ticks[4], ticks[5]),
+            sensor: ms(ticks[6], ticks[7]),
+            total: ms(ticks[0], ticks[7]),
         }
     }
 }

crates/resolved/src/gpu/shadow.rs

@@ -0,0 +1,267 @@
+//! The pass that renders shadow maps: each target's depth as seen from the Sun.
+//!
+//! One layer of a `Depth32Float` array per target, cleared to 1 (farthest from the Sun)
+//! and drawn with every caster the planner chose for it ([`crate::shadow::plan`]). The
+//! mesh pass then looks each fragment up in its target's layer.
+
+use nalgebra::Matrix4;
+
+use super::mesh::{GpuMesh, Vertex};
+use super::profile::PassTimestamps;
+use super::GpuContext;
+
+/// The shadow maps' depth format.
+pub const SHADOW_FORMAT: wgpu::TextureFormat = wgpu::TextureFormat::Depth32Float;
+
+/// One caster drawn into one layer.
+#[derive(Clone, Copy, Debug)]
+pub struct CasterDraw<'a> {
+    /// Which target's map (array layer).
+    pub layer: u32,
+    /// The caster's model.
+    pub mesh: &'a GpuMesh,
+    /// Body frame to the layer's light coordinates.
+    pub light_from_body: Matrix4<f32>,
+}
+
+/// The shadow pipeline, the map array and the uniforms for each caster draw.
+#[derive(Debug)]
+pub struct ShadowPass {
+    pipeline: wgpu::RenderPipeline,
+    caster_layout: wgpu::BindGroupLayout,
+    casters: wgpu::Buffer,
+    casters_group: wgpu::BindGroup,
+    caster_capacity: usize,
+    stride: u64,
+    size: u32,
+    maps: Option<(wgpu::Texture, u32)>,
+    sampler: wgpu::Sampler,
+}
+
+impl ShadowPass {
+    /// Compiles the shadow pipeline for maps of `size` x `size` texels.
+    pub fn new(ctx: &GpuContext, size: u32) -> Self {
+        let device = &ctx.device;
+        let caster_layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
+            label: Some("shadow caster"),
+            entries: &[wgpu::BindGroupLayoutEntry {
+                binding: 0,
+                visibility: wgpu::ShaderStages::VERTEX,
+                ty: wgpu::BindingType::Buffer {
+                    ty: wgpu::BufferBindingType::Uniform,
+                    has_dynamic_offset: true,
+                    min_binding_size: wgpu::BufferSize::new(64),
+                },
+                count: None,
+            }],
+        });
+        let layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
+            label: Some("shadow"),
+            bind_group_layouts: &[Some(&caster_layout)],
+            immediate_size: 0,
+        });
+        let module = device.create_shader_module(wgpu::ShaderModuleDescriptor {
+            label: Some("shadow"),
+            source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/shadow.wgsl").into()),
+        });
+        let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
+            label: Some("shadow"),
+            layout: Some(&layout),
+            vertex: wgpu::VertexState {
+                module: &module,
+                entry_point: Some("vs_main"),
+                compilation_options: Default::default(),
+                buffers: &[Some(Vertex::layout())],
+            },
+            // No culling: thin panels cast shadows whichever way they face.
+            primitive: wgpu::PrimitiveState::default(),
+            depth_stencil: Some(wgpu::DepthStencilState {
+                format: SHADOW_FORMAT,
+                depth_write_enabled: Some(true),
+                depth_compare: Some(wgpu::CompareFunction::Less),
+                stencil: Default::default(),
+                // Pushes stored depth away from the Sun on slanted surfaces, so a
+                // surface does not shadow itself where its own depth is sampled.
+                bias: wgpu::DepthBiasState {
+                    constant: 2,
+                    slope_scale: 2.0,
+                    clamp: 0.0,
+                },
+            }),
+            multisample: Default::default(),
+            fragment: None,
+            multiview_mask: None,
+            cache: None,
+        });
+        let alignment = device.limits().min_uniform_buffer_offset_alignment as u64;
+        let stride = 64u64.next_multiple_of(alignment);
+        let (casters, casters_group) = Self::caster_slots(ctx, &caster_layout, stride, 1);
+        let sampler = device.create_sampler(&wgpu::SamplerDescriptor {
+            label: Some("shadow compare"),
+            mag_filter: wgpu::FilterMode::Linear,
+            min_filter: wgpu::FilterMode::Linear,
+            compare: Some(wgpu::CompareFunction::LessEqual),
+            ..Default::default()
+        });
+        Self {
+            pipeline,
+            caster_layout,
+            casters,
+            casters_group,
+            caster_capacity: 1,
+            stride,
+            size,
+            maps: None,
+            sampler,
+        }
+    }
+
+    fn caster_slots(
+        ctx: &GpuContext,
+        layout: &wgpu::BindGroupLayout,
+        stride: u64,
+        capacity: usize,
+    ) -> (wgpu::Buffer, wgpu::BindGroup) {
+        let buffer = ctx.device.create_buffer(&wgpu::BufferDescriptor {
+            label: Some("shadow casters"),
+            size: stride * capacity as u64,
+            usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST,
+            mapped_at_creation: false,
+        });
+        let group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor {
+            label: Some("shadow casters"),
+            layout,
+            entries: &[wgpu::BindGroupEntry {
+                binding: 0,
+                resource: wgpu::BindingResource::Buffer(wgpu::BufferBinding {
+                    buffer: &buffer,
+                    offset: 0,
+                    size: wgpu::BufferSize::new(64),
+                }),
+            }],
+        });
+        (buffer, group)
+    }
+
+    /// Texels across each map.
+    pub fn size(&self) -> u32 {
+        self.size
+    }
+
+    /// Makes room for `layers` maps, reallocating the array if it is too small, and
+    /// returns the bind group the mesh pass samples them through (`layout` is its
+    /// shadow group layout). Allocation happens only when the target count grows.
+    pub fn prepare(
+        &mut self,
+        ctx: &GpuContext,
+        layers: u32,
+        layout: &wgpu::BindGroupLayout,
+    ) -> wgpu::BindGroup {
+        let have = self.maps.as_ref().map_or(0, |(_, n)| *n);
+        if have < layers {
+            let count = layers.next_power_of_two();
+            let texture = ctx.device.create_texture(&wgpu::TextureDescriptor {
+                label: Some("shadow maps"),
+                size: wgpu::Extent3d {
+                    width: self.size,
+                    height: self.size,
+                    depth_or_array_layers: count,
+                },
+                mip_level_count: 1,
+                sample_count: 1,
+                dimension: wgpu::TextureDimension::D2,
+                format: SHADOW_FORMAT,
+                usage: wgpu::TextureUsages::RENDER_ATTACHMENT
+                    | wgpu::TextureUsages::TEXTURE_BINDING,
+                view_formats: &[],
+            });
+            self.maps = Some((texture, count));
+        }
+        let (texture, _) = self.maps.as_ref().expect("allocated above");
+        let view = texture.create_view(&wgpu::TextureViewDescriptor {
+            dimension: Some(wgpu::TextureViewDimension::D2Array),
+            ..Default::default()
+        });
+        ctx.device.create_bind_group(&wgpu::BindGroupDescriptor {
+            label: Some("shadow maps"),
+            layout,
+            entries: &[
+                wgpu::BindGroupEntry {
+                    binding: 0,
+                    resource: wgpu::BindingResource::TextureView(&view),
+                },
+                wgpu::BindGroupEntry {
+                    binding: 1,
+                    resource: wgpu::BindingResource::Sampler(&self.sampler),
+                },
+            ],
+        })
+    }
+
+    /// Records one render pass per map layer in `0..layers`, drawing `draws` into their
+    /// layers. `timestamps` brackets all of them; with no layers, an empty pass keeps the
+    /// timestamps written.
+    pub fn encode(
+        &mut self,
+        ctx: &GpuContext,
+        encoder: &mut wgpu::CommandEncoder,
+        layers: u32,
+        draws: &[CasterDraw],
+        timestamps: Option<PassTimestamps>,
+    ) {
+        if draws.len() > self.caster_capacity {
+            self.caster_capacity = draws.len().next_power_of_two();
+            (self.casters, self.casters_group) =
+                Self::caster_slots(ctx, &self.caster_layout, self.stride, self.caster_capacity);
+        }
+        for (i, draw) in draws.iter().enumerate() {
+            let m = draw.light_from_body;
+            let columns: [[f32; 4]; 4] =
+                std::array::from_fn(|c| [m[(0, c)], m[(1, c)], m[(2, c)], m[(3, c)]]);
+            ctx.queue.write_buffer(
+                &self.casters,
+                i as u64 * self.stride,
+                bytemuck::cast_slice(&columns),
+            );
+        }
+        let texture = match &self.maps {
+            Some((texture, _)) if layers > 0 => texture,
+            _ => {
+                // Nothing to draw: an empty pass still writes the timestamps.
+                drop(encoder.begin_compute_pass(&wgpu::ComputePassDescriptor {
+                    label: Some("shadow (none)"),
+                    timestamp_writes: timestamps.map(PassTimestamps::compute),
+                }));
+                return;
+            }
+        };
+        for layer in 0..layers {
+            let view = texture.create_view(&wgpu::TextureViewDescriptor {
+                dimension: Some(wgpu::TextureViewDimension::D2),
+                base_array_layer: layer,
+                array_layer_count: Some(1),
+                ..Default::default()
+            });
+            let stamps = timestamps.map(|t| t.split(layer == 0, layer + 1 == layers));
+            let mut pass = encoder.begin_render_pass(&wgpu::RenderPassDescriptor {
+                label: Some("shadow"),
+                color_attachments: &[],
+                depth_stencil_attachment: Some(wgpu::RenderPassDepthStencilAttachment {
+                    view: &view,
+                    depth_ops: Some(wgpu::Operations {
+                        load: wgpu::LoadOp::Clear(1.0),
+                        store: wgpu::StoreOp::Store,
+                    }),
+                    stencil_ops: None,
+                }),
+                timestamp_writes: stamps.map(PassTimestamps::render),
+                ..Default::default()
+            });
+            pass.set_pipeline(&self.pipeline);
+            for (i, draw) in draws.iter().enumerate().filter(|(_, d)| d.layer == layer) {
+                pass.set_bind_group(0, &self.casters_group, &[(i as u64 * self.stride) as u32]);
+                draw.mesh.draw(&mut pass);
+            }
+        }
+    }
+}

crates/resolved/src/lib.rs

@@ -30,6 +30,7 @@ pub mod model;
 pub mod radiometry;
 pub mod render;
 pub mod scene;
+pub mod shadow;
 mod solar_data;
 pub mod spectrum;
 

crates/resolved/src/render.rs

@@ -18,14 +18,16 @@ use crate::earth::EarthTexture;
 use crate::eclipse::sun_visible_fraction;
 use crate::fidelity::{Fidelity, LightingTier};
 use crate::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry};
-use crate::gpu::mesh::{GpuMesh, MeshDraw, MeshPass};
+use crate::gpu::mesh::{GpuMesh, MeshDraw, MeshPass, ShadowDraw};
 use crate::gpu::profile::{PassTimings, Profiler};
 use crate::gpu::sensor::{SensorPass, SensorUniforms};
+use crate::gpu::shadow::{CasterDraw, ShadowPass};
 use crate::gpu::targets::FrameTargets;
 use crate::gpu::{AdapterChoice, GpuContext};
 use crate::model::{Mesh, ModelId};
 use crate::radiometry::Radiometry;
 use crate::scene::{Pose, Scene};
+use crate::shadow::{self, Bounded};
 use crate::Error;
 
 /// Choices made once, when the renderer is created.
@@ -45,6 +47,10 @@ pub struct RendererOptions {
     /// Time each pass on the GPU and report it in [`Frame::timings`]. Needs an adapter
     /// with timestamp queries; without one, timings are `None`.
     pub profile: bool,
+    /// Texels across each target's shadow map at lighting L1 and up; `None` for 4096.
+    /// Each map is 4 bytes per texel (64 MB at 4096) and covers its target's bounding
+    /// sphere, so a 16 m model gets 4 mm texels.
+    pub shadow_map_size: Option<u32>,
 }
 
 /// One rendered image.
@@ -88,6 +94,7 @@ pub struct Renderer {
     targets: FrameTargets,
     bodies: BodiesPass,
     meshes: MeshPass,
+    shadows: ShadowPass,
     sensor: SensorPass,
     sensor_uniforms: SensorUniforms,
     models: Vec<GpuMesh>,
@@ -140,6 +147,7 @@ impl Renderer {
         let mut renderer = Self {
             bodies: BodiesPass::new(&ctx),
             meshes: MeshPass::new(&ctx),
+            shadows: ShadowPass::new(&ctx, options.shadow_map_size.unwrap_or(4096)),
             sensor: SensorPass::new(&ctx),
             sensor_uniforms: SensorUniforms::new(&camera, &radiometry),
             ctx,
@@ -226,12 +234,58 @@ impl Renderer {
         ];
 
         let bodies = BodiesUniforms::new(scene, &camera, &intrinsics, near, &self.body_radiometry);
-        let draws: Vec<(&GpuMesh, MeshDraw)> = scene
+        // Shadow maps, from lighting L1: one per target, with the targets that can
+        // shadow it drawn in.
+        let bounded: Vec<Bounded> = scene
             .targets
             .iter()
             .map(|t| {
+                let mesh = &self.models[t.model.0];
+                Bounded {
+                    pose: t.pose,
+                    center: mesh.center,
+                    radius: mesh.radius.max(1e-3),
+                }
+            })
+            .collect();
+        let plans = if self.fidelity.lighting >= LightingTier::L1 {
+            shadow::plan(&bounded, scene.sun)
+        } else {
+            Vec::new()
+        };
+        let size = self.shadows.size();
+        let casters: Vec<CasterDraw> = plans
+            .iter()
+            .enumerate()
+            .flat_map(|(layer, plan)| plan.casters.iter().map(move |&j| (layer, plan, j)))
+            .map(|(layer, plan, j)| CasterDraw {
+                layer: layer as u32,
+                mesh: &self.models[scene.targets[j].model.0],
+                light_from_body: plan.light_from_body(&bounded[j]).cast(),
+            })
+            .collect();
+        let draws: Vec<(&GpuMesh, MeshDraw)> = scene
+            .targets
+            .iter()
+            .enumerate()
+            .map(|(i, t)| {
                 let visible = sun_visible_fraction(t.pose.position, scene.sun, &occulters);
-                let draw = MeshDraw::new(&t.pose, &camera, &intrinsics, near, scene.sun, visible);
+                let shadow = plans.get(i).map(|plan| ShadowDraw {
+                    light_from_body: plan.light_from_body(&bounded[i]),
+                    layer: i as u32,
+                    texel: 2.0 * plan.half_width / size as f64,
+                    depth_range: plan.depth_range(),
+                    size,
+                });
+                let draw = MeshDraw::new(
+                    &t.pose,
+                    &camera,
+                    &intrinsics,
+                    near,
+                    scene.sun,
+                    visible,
+                    shadow,
+                );
                 (&self.models[t.model.0], draw)
             })
             .collect();
@@ -239,16 +293,28 @@ impl Renderer {
         let ctx = &self.ctx;
         let mut encoder = ctx.device.create_command_encoder(&Default::default());
         let stamps = |i| self.profiler.as_ref().map(|p| p.pass(i));
+        let shadow_group = (!plans.is_empty()).then(|| {
+            self.shadows
+                .prepare(ctx, plans.len() as u32, self.meshes.shadow_layout())
+        });
+        self.shadows
+            .encode(ctx, &mut encoder, plans.len() as u32, &casters, stamps(0));
         self.bodies
-            .encode(ctx, &mut encoder, &self.targets, &bodies, stamps(0));
-        self.meshes
-            .encode(ctx, &mut encoder, &self.targets, &draws, stamps(1));
+            .encode(ctx, &mut encoder, &self.targets, &bodies, stamps(1));
+        self.meshes.encode(
+            ctx,
+            &mut encoder,
+            &self.targets,
+            &draws,
+            shadow_group.as_ref(),
+            stamps(2),
+        );
         self.sensor.encode(
             ctx,
             &mut encoder,
             &self.targets,
             &self.sensor_uniforms,
-            stamps(2),
+            stamps(3),
         );
         if let Some(profiler) = &self.profiler {
             profiler.resolve(&mut encoder);

crates/resolved/src/shaders/mesh.wgsl

@@ -14,11 +14,19 @@ struct Draw {
     sun: vec4<f32>,
     // x: Sun visibility (eclipse) / (Sun distance in AU)^2.
     params: vec4<f32>,
+    // Body frame to this target's shadow map (shadow.rs).
+    light_from_body: mat4x4<f32>,
+    // Layer, enabled (0 or 1), meters per texel, depth range in meters.
+    shadow: vec4<f32>,
+    // Tangent of the Sun's angular diameter, normal offset in meters, map size in texels.
+    shadow_more: vec4<f32>,
 }
 
 @group(0) @binding(0) var<uniform> draw: Draw;
 // Per part: channel radiance of a surface facing the Sun at 1 AU (albedo S_k / pi).
 @group(1) @binding(0) var<storage, read> parts: array<vec4<f32>>;
+@group(2) @binding(0) var shadow_maps: texture_depth_2d_array;
+@group(2) @binding(1) var shadow_compare: sampler_comparison;
 
 struct VertexIn {
     @location(0) position: vec3<f32>,
@@ -31,19 +39,80 @@ struct VertexOut {
     @location(0) camera_position: vec3<f32>,
     @location(1) normal: vec3<f32>,
     @location(2) @interpolate(flat) part: u32,
+    // Shadow-map coordinates: x and y in [-1, 1], depth in [0, 1].
+    @location(3) light: vec3<f32>,
 }
 
 @vertex
 fn vs_main(v: VertexIn) -> VertexOut {
     let p = draw.camera_from_body * v.position + draw.translation.xyz;
+    let n = draw.camera_from_body * v.normal;
     var out: VertexOut;
     out.clip = draw.clip_from_camera * vec4<f32>(p, 1.0);
     out.camera_position = p;
-    out.normal = draw.camera_from_body * v.normal;
+    out.normal = n;
     out.part = v.part;
+    // Look the shadow up slightly off the surface, on its sunward side, so a surface
+    // does not shadow itself where the map samples it at a slant.
+    let toward_sun = select(-1.0, 1.0, dot(n, draw.sun.xyz) >= 0.0);
+    let offset = v.normal * toward_sun * draw.shadow_more.y;
+    out.light = (draw.light_from_body * vec4<f32>(v.position + offset, 1.0)).xyz;
     return out;
 }
 
+// The fraction of the Sun this point sees past the target's own geometry and its
+// neighbors, from the shadow map, with a penumbra as wide as the Sun's disk makes it.
+//
+// Percentage-closer soft shadows: first the mean depth of whatever blocks the Sun near
+// this point, then a comparison filter whose width is the penumbra that blocker casts,
+// (receiver - blocker distance) x the Sun's angular diameter. Both use a regular grid
+// of samples rather than a random pattern, so the result has no noise.
+fn sunlit(light: vec3<f32>) -> f32 {
+    if draw.shadow.y < 0.5 {
+        return 1.0;
+    }
+    let layer = i32(draw.shadow.x);
+    let texel = draw.shadow.z;
+    let range = draw.shadow.w;
+    let spread = draw.shadow_more.x;
+    let size = draw.shadow_more.z;
+    let uv = vec2<f32>(light.x * 0.5 + 0.5, 0.5 - light.y * 0.5);
+    // Two texels of depth margin against comparing a surface with itself.
+    let depth = light.z - 2.0 * texel / range;
+
+    // Blockers can only be between this point and the map's sunward edge; anything
+    // there casting onto this point lies within that distance times the spread.
+    let search = clamp(light.z * range * spread / texel, 1.0, 32.0);
+    var blockers = 0.0;
+    var blocker_depth = 0.0;
+    for (var i = 0; i < 5; i++) {
+        for (var j = 0; j < 5; j++) {
+            let offset = (vec2<f32>(f32(i), f32(j)) * 0.5 - 1.0) * search;
+            let texel_xy = clamp(vec2<i32>(floor(uv * size + offset)), vec2<i32>(0), vec2<i32>(i32(size) - 1));
+            let d = textureLoad(shadow_maps, texel_xy, layer, 0);
+            if d < depth {
+                blockers += 1.0;
+                blocker_depth += d;
+            }
+        }
+    }
+    if blockers == 0.0 {
+        return 1.0;
+    }
+    let penumbra = (depth - blocker_depth / blockers) * range * spread;
+    // Half the penumbra either side of the edge, and never under half a texel, which
+    // also smooths the map's own stair-steps.
+    let radius = max(0.5 * penumbra / texel, 0.5) / size;
+    var lit = 0.0;
+    for (var i = 0; i < 6; i++) {
+        for (var j = 0; j < 6; j++) {
+            let offset = ((vec2<f32>(f32(i), f32(j)) + 0.5) / 3.0 - 1.0) * radius;
+            lit += textureSampleCompareLevel(shadow_maps, shadow_compare, uv + offset, layer, depth);
+        }
+    }
+    return lit / 36.0;
+}
+
 @fragment
 fn fs_main(in: VertexOut) -> @location(0) vec4<f32> {
     var n = normalize(in.normal);
@@ -52,5 +121,10 @@ fn fs_main(in: VertexOut) -> @location(0) vec4<f32> {
     if dot(n, in.camera_position) > 0.0 {
         n = -n;
     }
-    return parts[in.part] * max(dot(n, draw.sun.xyz), 0.0) * draw.params.x;
+    let lit = max(dot(n, draw.sun.xyz), 0.0);
+    var shade = 0.0;
+    if lit > 0.0 {
+        shade = sunlit(in.light);
+    }
+    return parts[in.part] * lit * shade * draw.params.x;
 }

crates/resolved/src/shaders/shadow.wgsl

@@ -0,0 +1,16 @@
+// Depth from the Sun: draws one caster into one target's shadow map.
+//
+// The matrix takes body-frame vertices straight to the map's light coordinates
+// (shadow.rs): x and y across the map in [-1, 1], depth along the Sun's rays in [0, 1],
+// 0 nearest the Sun. No fragment stage: only depth is written.
+
+struct Caster {
+    light_from_body: mat4x4<f32>,
+}
+
+@group(0) @binding(0) var<uniform> caster: Caster;
+
+@vertex
+fn vs_main(@location(0) position: vec3<f32>) -> @builtin(position) vec4<f32> {
+    return caster.light_from_body * vec4<f32>(position, 1.0);
+}

crates/resolved/src/shadow.rs

@@ -0,0 +1,188 @@
+//! Planning shadow maps: which targets can shadow which, and how each map is framed.
+//!
+//! Each target gets an orthographic shadow map looking along the direction from the
+//! target to the Sun (the Sun is far enough that its rays are parallel across a
+//! spacecraft). The map covers the target's bounding sphere, so every point of the
+//! target lands in it, and holds the depth, along the Sun's direction, of whatever is
+//! nearest the Sun: the target itself or another target sunward of it.
+//!
+//! Light coordinates of a point P for the map of a target with bounding center c and
+//! radius r, along axes u, v (across the Sun's rays) and w (toward the Sun):
+//!
+//! ```text
+//! x = (P - c) . u / r            in [-1, 1] across the map
+//! y = (P - c) . v / r
+//! depth = (near - (P - c) . w) / (near - far)    0 nearest the Sun, 1 farthest
+//! ```
+//!
+//! with `near` the sunward extent of the farthest caster and `far = -r`.
+
+use nalgebra::{Matrix4, Vector3};
+
+use crate::scene::Pose;
+
+/// A target as the planner sees it: where it is and the sphere around its model.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Bounded {
+    /// The target's pose.
+    pub pose: Pose,
+    /// Bounding sphere center in the model's body frame, in meters.
+    pub center: Vector3<f64>,
+    /// Bounding sphere radius, in meters.
+    pub radius: f64,
+}
+
+impl Bounded {
+    /// The bounding sphere's center in the inertial frame.
+    pub fn world_center(&self) -> Vector3<f64> {
+        self.pose.position + self.pose.attitude * self.center
+    }
+}
+
+/// How one target's shadow map is framed and which targets are drawn into it.
+#[derive(Clone, Debug, PartialEq)]
+pub struct ShadowPlan {
+    /// The receiver's bounding center, inertial.
+    pub center: Vector3<f64>,
+    /// Half the map's width, in meters: the receiver's bounding radius.
+    pub half_width: f64,
+    /// Unit vectors across the map (u, v) and toward the Sun (w), inertial.
+    pub axes: [Vector3<f64>; 3],
+    /// Distance from `center` toward the Sun of the map's near plane, in meters.
+    pub near: f64,
+    /// The same for the far plane (negative: away from the Sun).
+    pub far: f64,
+    /// Indices of the targets drawn into this map, the receiver included.
+    pub casters: Vec<usize>,
+}
+
+impl ShadowPlan {
+    /// The map's depth range in meters.
+    pub fn depth_range(&self) -> f64 {
+        self.near - self.far
+    }
+
+    /// The matrix taking a point in `target`'s body frame to this map's light
+    /// coordinates (x, y in [-1, 1], depth in [0, 1], w = 1).
+    pub fn light_from_body(&self, target: &Bounded) -> Matrix4<f64> {
+        let [u, v, w] = self.axes;
+        let range = self.depth_range();
+        // Rows act on (P - c); P = R p + t for body point p.
+        let rows = [u / self.half_width, v / self.half_width, -w / range];
+        let rotation = target.pose.attitude.to_rotation_matrix();
+        let offset = target.pose.position - self.center;
+        let mut m = Matrix4::identity();
+        for (i, row) in rows.iter().enumerate() {
+            let r = rotation.matrix().transpose() * row;
+            m[(i, 0)] = r.x;
+            m[(i, 1)] = r.y;
+            m[(i, 2)] = r.z;
+            m[(i, 3)] = row.dot(&offset);
+        }
+        // depth = (near - along) / range: the -along / range is in the row; add near / range.
+        m[(2, 3)] += self.near / range;
+        m
+    }
+}
+
+/// One shadow map per target, lit from the Sun at `sun` (inertial, meters).
+///
+/// A target shadows another when its bounding sphere reaches the cylinder swept by the
+/// other's sphere toward the Sun.
+pub fn plan(targets: &[Bounded], sun: Vector3<f64>) -> Vec<ShadowPlan> {
+    targets
+        .iter()
+        .enumerate()
+        .map(|(i, receiver)| {
+            let center = receiver.world_center();
+            let w = (sun - center).normalize();
+            let u = w
+                .cross(&Vector3::z())
+                .try_normalize(1e-6)
+                .unwrap_or_else(|| w.cross(&Vector3::x()).normalize());
+            let v = w.cross(&u);
+            let r = receiver.radius;
+            let mut near = r;
+            let mut casters = vec![i];
+            for (j, other) in targets.iter().enumerate().filter(|(j, _)| *j != i) {
+                let d = other.world_center() - center;
+                let along = d.dot(&w);
+                let across = (d - along * w).norm();
+                let reach = r + other.radius;
+                if across < reach && along > -reach {
+                    casters.push(j);
+                    near = near.max(along + other.radius);
+                }
+            }
+            ShadowPlan {
+                center,
+                half_width: r,
+                axes: [u, v, w],
+                near,
+                far: -r,
+                casters,
+            }
+        })
+        .collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+    use nalgebra::{UnitQuaternion, Vector4};
+
+    fn at(x: f64, y: f64, z: f64, radius: f64) -> Bounded {
+        Bounded {
+            pose: Pose::new(Vector3::new(x, y, z), UnitQuaternion::identity()),
+            center: Vector3::zeros(),
+            radius,
+        }
+    }
+
+    const SUN: Vector3<f64> = Vector3::new(1.496e11, 0.0, 0.0);
+
+    #[test]
+    fn sunward_neighbor_casts_and_side_neighbor_does_not() {
+        // Target 0 at the origin; 1 is 20 m sunward (+X); 2 is 20 m off to the side.
+        let targets = [
+            at(0.0, 0.0, 0.0, 5.0),
+            at(20.0, 2.0, 0.0, 3.0),
+            at(0.0, 20.0, 0.0, 3.0),
+        ];
+        let plans = plan(&targets, SUN);
+        assert_eq!(plans[0].casters, [0, 1]);
+        assert_relative_eq!(plans[0].near, 23.0, epsilon = 1e-6);
+        assert_relative_eq!(plans[0].far, -5.0);
+        // Nothing is sunward of target 1.
+        assert_eq!(plans[1].casters, [1]);
+        // Target 2 is beside 0, not sunward of it: it does not shadow 0, nor 0 it.
+        assert_eq!(plans[2].casters, [2]);
+    }
+
+    #[test]
+    fn light_coordinates_span_the_bounding_sphere() {
+        let target = Bounded {
+            pose: Pose::new(
+                Vector3::new(7.0e6, 1.0, 2.0),
+                UnitQuaternion::from_euler_angles(0.3, -0.5, 1.1),
+            ),
+            center: Vector3::new(0.5, -1.0, 2.0),
+            radius: 4.0,
+        };
+        let plan = &plan(&[target], SUN)[0];
+        let m = plan.light_from_body(&target);
+        let light = |body: Vector3<f64>| m * Vector4::new(body.x, body.y, body.z, 1.0);
+        // The bounding center maps to the middle of the map, at the depth of `center`.
+        let c = light(target.center);
+        assert_relative_eq!(c.xy(), Vector3::zeros().xy(), epsilon = 1e-9);
+        assert_relative_eq!(c.z, plan.near / plan.depth_range(), epsilon = 1e-9);
+        // A body point 4 m toward the Sun is at depth 0 (the near plane); 4 m away, 1.
+        let toward = target.pose.attitude.inverse() * plan.axes[2];
+        assert_relative_eq!(light(target.center + 4.0 * toward).z, 0.0, epsilon = 1e-9);
+        assert_relative_eq!(light(target.center - 4.0 * toward).z, 1.0, epsilon = 1e-9);
+        // 4 m across along u is the map's edge.
+        let across = target.pose.attitude.inverse() * plan.axes[0];
+        assert_relative_eq!(light(target.center + 4.0 * across).x, 1.0, epsilon = 1e-9);
+    }
+}

crates/resolved/tests/mesh.rs

@@ -112,13 +112,21 @@ fn render(ctx: &GpuContext, name: &str, scene: &Scene, camera: &Pose, mesh: &Mes
             );
             (
                 &gpu_mesh,
-                MeshDraw::new(&t.pose, camera, &k, config.near_clip, scene.sun, visible),
+                MeshDraw::new(
+                    &t.pose,
+                    camera,
+                    &k,
+                    config.near_clip,
+                    scene.sun,
+                    visible,
+                    None,
+                ),
             )
         })
         .collect();
     let mut encoder = ctx.device.create_command_encoder(&Default::default());
     bodies.encode(ctx, &mut encoder, &targets, &uniforms, None);
-    meshes.encode(ctx, &mut encoder, &targets, &draws, None);
+    meshes.encode(ctx, &mut encoder, &targets, &draws, None, None);
     ctx.queue.submit([encoder.finish()]);
     let radiance = targets.read_radiance(ctx);
     common::save(name, &ctx.info, &radiance);

crates/resolved/tests/shadows.rs

@@ -0,0 +1,384 @@
+//! Shadows on targets at lighting L1, on the GPU and on the CPU (lavapipe): where they
+//! fall, how wide their penumbra is, and that surfaces do not shadow themselves.
+
+mod common;
+
+use std::io::Cursor;
+use std::sync::Arc;
+
+use resolved::gpu::AdapterChoice;
+use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3};
+use resolved::radiometry::{Radiometry, AU};
+use resolved::{
+    CameraConfig, EarthTexture, Fidelity, LightingTier, LoadOptions, Material, Mesh, Pose,
+    Renderer, RendererOptions, Scene, Target,
+};
+
+const ALBEDO: f64 = 0.5;
+const TARGET: Vector3<f64> = Vector3::new(7.0e6, 0.0, 0.0);
+/// Straight up from the target: body +X faces the Sun with the identity attitude.
+const SUN: Vector3<f64> = Vector3::new(AU, 0.0, 0.0);
+/// Half the width of the plate the shadows fall on, in meters.
+const PLATE: f32 = 4.0;
+/// The Sun's angular diameter seen from 1 AU, as a tangent: 2 x 696,000 km / 1 AU.
+const SUN_SPREAD: f64 = 0.0093;
+
+/// Four-cornered faces, each split into two triangles, as one gray mesh.
+fn quads(name: &str, faces: &[[[f32; 3]; 4]]) -> Mesh {
+    let mut s = format!("solid {name}\n");
+    for f in faces {
+        for tri in [[f[0], f[1], f[2]], [f[0], f[2], f[3]]] {
+            s += "facet normal 0 0 0\nouter loop\n";
+            for [x, y, z] in tri {
+                s += &format!("vertex {x} {y} {z}\n");
+            }
+            s += "endloop\nendfacet\n";
+        }
+    }
+    s += &format!("endsolid {name}\n");
+    let options = LoadOptions {
+        default_material: Material::gray("gray", ALBEDO),
+        ..Default::default()
+    };
+    Mesh::from_stl(&mut Cursor::new(s), name, &options).unwrap()
+}
+
+/// A square in the body YZ plane at height `x`, facing +X.
+fn square(x: f32, half: f32) -> [[f32; 3]; 4] {
+    [
+        [x, -half, -half],
+        [x, half, -half],
+        [x, half, half],
+        [x, -half, half],
+    ]
+}
+
+/// The plate alone.
+fn plate() -> Mesh {
+    quads("plate", &[square(0.0, PLATE)])
+}
+
+/// The plate with a 2 m square occluder `height` meters above its center.
+fn plate_under(height: f32) -> Mesh {
+    quads("plate", &[square(0.0, PLATE), square(height, 1.0)])
+}
+
+/// A UV sphere of radius 3 m, flat-faceted: 48 segments around, 24 from pole to pole.
+fn sphere() -> Mesh {
+    let (around, down) = (48, 24);
+    let point = |i: usize, j: usize| {
+        let lon = i as f32 / around as f32 * std::f32::consts::TAU;
+        let lat = (j as f32 / down as f32 - 0.5) * std::f32::consts::PI;
+        [lat.cos() * lon.cos(), lat.cos() * lon.sin(), lat.sin()].map(|c| 3.0 * c)
+    };
+    let faces: Vec<_> = (0..around)
+        .flat_map(|i| (0..down).map(move |j| (i, j)))
+        .map(|(i, j)| {
+            [
+                point(i, j),
+                point(i + 1, j),
+                point(i + 1, j + 1),
+                point(i, j + 1),
+            ]
+        })
+        .collect();
+    quads("sphere", &faces)
+}
+
+fn renderer(adapter: &AdapterChoice, lighting: LightingTier) -> Renderer {
+    let gray = EarthTexture::from_rgba8(1, 1, vec![128, 128, 128, 255]).unwrap();
+    let options = RendererOptions {
+        adapter: adapter.clone(),
+        fidelity: Fidelity {
+            lighting,
+            ..Default::default()
+        },
+        earth_texture: Some(Arc::new(gray)),
+        read_radiance: true,
+        ..Default::default()
+    };
+    Renderer::new(CameraConfig::example(), options).unwrap()
+}
+
+/// 30 m from the target along `direction` (target body frame), looking at it.
+fn observer(direction: Vector3<f64>, up: Vector3<f64>) -> Pose {
+    Pose::looking_at(TARGET + 30.0 * direction.normalize(), TARGET, up).unwrap()
+}
+
+/// Looking down at the plate 50 degrees from overhead, so an occluder above the plate
+/// appears well clear of the shadow it casts.
+fn slanted() -> Pose {
+    let a = 50f64.to_radians();
+    observer(Vector3::new(a.cos(), 0.0, a.sin()), Vector3::z())
+}
+
+/// 10 m from the occluder's shadow edge (body z = 1 m), 40 degrees from overhead,
+/// where a pixel spans 1.4 mm of plate across the edge, a seventh of a 1 m occluder's
+/// penumbra.
+fn close_to_edge() -> Pose {
+    let a = 40f64.to_radians();
+    let edge = TARGET + Vector3::z();
+    let eye = edge + 10.0 * Vector3::new(a.cos(), 0.0, a.sin());
+    Pose::looking_at(eye, edge, Vector3::z()).unwrap()
+}
+
+fn scene(targets: Vec<Target>) -> Scene {
+    Scene {
+        earth_attitude: UnitQuaternion::identity(),
+        sun: SUN,
+        moon: Vector3::new(0.0, -3.844e8, 0.0),
+        targets,
+    }
+}
+
+fn at(model: resolved::ModelId, offset: Vector3<f64>) -> Target {
+    Target {
+        model,
+        pose: Pose::new(TARGET + offset, UnitQuaternion::identity()),
+    }
+}
+
+/// Green-channel radiance of a Lambertian surface at `cos_i` to the Sun, at 1 AU.
+fn lit(cos_i: f64) -> f64 {
+    Radiometry::new(&CameraConfig::example()).lambert_radiance(
+        1,
+        ALBEDO,
+        cos_i,
+        (SUN - TARGET).norm(),
+    )
+}
+
+/// The plate's (y, z), body frame, and the green radiance there, for each pixel of
+/// image column 960, the column through the plate's center line. The plate sits at the
+/// identity attitude on `TARGET`.
+fn column(observer: &Pose, radiance: &[[f32; 4]]) -> Vec<(f64, f64)> {
+    let intrinsics = CameraConfig::example().intrinsics();
+    let mut samples: Vec<(f64, f64)> = (0..1080)
+        .filter_map(|row| {
+            let ray = observer.attitude * intrinsics.ray(&Vector2::new(960.5, row as f64 + 0.5));
+            // Where the ray meets the plate's plane, body x = 0.
+            let t = (TARGET.x - observer.position.x) / ray.x;
+            let p = observer.position + t * ray - TARGET;
+            (p.y.abs() < PLATE as f64 && p.z.abs() < PLATE as f64)
+                .then(|| (p.z, radiance[row * 1920 + 960][1] as f64))
+        })
+        .collect();
+    samples.sort_by(|a, b| a.0.total_cmp(&b.0));
+    samples
+}
+
+/// Where, going up through `samples` from `from`, the radiance first reaches `level`,
+/// interpolated between pixels.
+fn crossing(samples: &[(f64, f64)], from: f64, level: f64) -> f64 {
+    let i = samples
+        .iter()
+        .position(|&(z, r)| z > from && r >= level)
+        .expect("never reaches the level");
+    let ((z0, r0), (z1, r1)) = (samples[i - 1], samples[i]);
+    z0 + (level - r0) / (r1 - r0) * (z1 - z0)
+}
+
+#[test]
+fn occluder_shadows_the_plate_below_it() {
+    let _devices = common::device_lock();
+    let full = lit(1.0);
+    for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
+        let mut renderer = renderer(&adapter, LightingTier::L1);
+        let name = renderer.adapter_info().name.clone();
+        let model = renderer.load_model(&plate_under(2.0)).unwrap();
+        let frame = renderer
+            .submit(&scene(vec![at(model, Vector3::zeros())]), &slanted())
+            .unwrap()
+            .wait()
+            .unwrap();
+        let radiance = frame.radiance.unwrap();
+        common::save("shadow-occluder", renderer.adapter_info(), &radiance);
+        let samples = column(&slanted(), &radiance);
+        for &(z, r) in &samples {
+            // Under the occluder, 2 cm clear of the penumbra: no direct sunlight at all.
+            if z.abs() < 0.97 {
+                assert_eq!(r, 0.0, "{name}: z = {z}");
+            }
+            // Beyond it, and short of where the occluder hides the plate from the
+            // camera: plain Lambert.
+            if z > 1.03 {
+                assert!(
+                    (r - full).abs() < 1e-4 * full,
+                    "{name}: z = {z}: {r} vs {full}"
+                );
+            }
+        }
+        // The edge is the occluder's edge, z = 1 m, to within a texel of the 4096 map
+        // across the plate's 11.3 m bounding sphere (2.8 mm).
+        let edge = crossing(&samples, 0.0, 0.5 * full);
+        assert!((edge - 1.0).abs() < 2.8e-3, "{name}: edge at {edge}");
+    }
+}
+
+#[test]
+fn penumbra_widens_with_occluder_height() {
+    let _devices = common::device_lock();
+    let full = lit(1.0);
+    for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
+        let mut renderer = renderer(&adapter, LightingTier::L1);
+        let name = renderer.adapter_info().name.clone();
+        for height in [1.0, 2.0, 3.0] {
+            let model = renderer.load_model(&plate_under(height)).unwrap();
+            let frame = renderer
+                .submit(&scene(vec![at(model, Vector3::zeros())]), &close_to_edge())
+                .unwrap()
+                .wait()
+                .unwrap();
+            let samples = column(&close_to_edge(), &frame.radiance.unwrap());
+            // The Sun's disk, seen past a straight edge, uncovers linearly in the
+            // filter this models, so 10 % to 90 % is 0.8 of the full width.
+            let width =
+                (crossing(&samples, 0.0, 0.9 * full) - crossing(&samples, 0.0, 0.1 * full)) / 0.8;
+            let want = height as f64 * SUN_SPREAD;
+            assert!(
+                (width - want).abs() < 0.05 * want,
+                "{name}: occluder {height} m up: penumbra {width} m, want {want} m"
+            );
+        }
+    }
+}
+
+/// Pixels lit at L0 but darker at L1, among those at least 10 % of `full`: a surface
+/// that nothing shadows must not darken (shadow acne).
+fn darkened(l0: &[[f32; 4]], l1: &[[f32; 4]], full: f64) -> usize {
+    l0.iter()
+        .zip(l1)
+        .filter(|(a, b)| a[1] as f64 > 0.1 * full && b[1] < 0.999 * a[1])
+        .count()
+}
+
+#[test]
+fn unshadowed_surfaces_do_not_shadow_themselves() {
+    let _devices = common::device_lock();
+    for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
+        let mut renderers =
+            [LightingTier::L0, LightingTier::L1].map(|lighting| renderer(&adapter, lighting));
+        let name = renderers[1].adapter_info().name.clone();
+        // The plate tipped 80 and 87 degrees from the Sun, viewed face on, and a sphere
+        // of flat facets, whose terminator crosses facets at every angle to the Sun,
+        // seen against space so only the target can differ.
+        let a = |deg: f64| deg.to_radians();
+        let tipped = |deg: f64| {
+            let attitude = UnitQuaternion::from_axis_angle(&Vector3::y_axis(), -a(deg));
+            let normal = attitude * Vector3::x();
+            (
+                plate(),
+                attitude,
+                observer(normal, Vector3::y()),
+                format!("plate at {deg}"),
+            )
+        };
+        let cases = [
+            tipped(80.0),
+            tipped(87.0),
+            (
+                sphere(),
+                UnitQuaternion::identity(),
+                observer(Vector3::new(-0.3, 1.0, 0.3), Vector3::z()),
+                "sphere".to_string(),
+            ),
+        ];
+        for (mesh, attitude, observer, case) in cases {
+            let frames = renderers.each_mut().map(|renderer| {
+                let model = renderer.load_model(&mesh).unwrap();
+                let target = Target {
+                    model,
+                    pose: Pose::new(TARGET, attitude),
+                };
+                renderer
+                    .submit(&scene(vec![target]), &observer)
+                    .unwrap()
+                    .wait()
+                    .unwrap()
+                    .radiance
+                    .unwrap()
+            });
+            let file = format!("shadow-acne-{}", case.replace(' ', "-"));
+            common::save(&file, renderers[1].adapter_info(), &frames[1]);
+            let count = darkened(&frames[0], &frames[1], lit(1.0));
+            assert_eq!(count, 0, "{name}: {case}: {count} pixels darkened");
+        }
+    }
+}
+
+#[test]
+fn sunward_target_shadows_and_side_target_does_not() {
+    let _devices = common::device_lock();
+    let full = lit(1.0);
+    for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
+        let mut renderer = renderer(&adapter, LightingTier::L1);
+        let name = renderer.adapter_info().name.clone();
+        let plate = renderer.load_model(&plate()).unwrap();
+        let panel = renderer
+            .load_model(&quads("panel", &[square(0.0, 1.0)]))
+            .unwrap();
+        // A 2 m panel tipped 45 degrees, so its +Z half is the nearer the Sun; that
+        // half's shadow covers plate z from 0 to 0.71 m.
+        let tipped = |offset: Vector3<f64>| Target {
+            model: panel,
+            pose: Pose::new(
+                TARGET + offset,
+                UnitQuaternion::from_axis_angle(&Vector3::y_axis(), 45f64.to_radians()),
+            ),
+        };
+        let mut at_z = |targets, want: f64| {
+            let frame = renderer
+                .submit(&scene(targets), &slanted())
+                .unwrap()
+                .wait()
+                .unwrap();
+            let samples = column(&slanted(), &frame.radiance.unwrap());
+            let (_, r) = *samples
+                .iter()
+                .min_by(|a, b| (a.0 - want).abs().total_cmp(&(b.0 - want).abs()))
+                .unwrap();
+            r
+        };
+        // The panel 6 m sunward of the plate's center, listed before or after the
+        // plate, shadows it on both sides of center.
+        let sunward = Vector3::new(6.0, 0.0, 0.0);
+        for z in [-0.4, 0.4] {
+            let before = at_z(vec![tipped(sunward), at(plate, Vector3::zeros())], z);
+            let after = at_z(vec![at(plate, Vector3::zeros()), tipped(sunward)], z);
+            assert_eq!((before, after), (0.0, 0.0), "{name}: z = {z}");
+        }
+        // The same panel 12 m to the side does not.
+        let side = at_z(
+            vec![
+                at(plate, Vector3::zeros()),
+                tipped(Vector3::new(0.0, -12.0, 0.0)),
+            ],
+            0.0,
+        );
+        assert!(
+            (side - full).abs() < 1e-4 * full,
+            "{name}: {side} vs {full}"
+        );
+    }
+}
+
+#[test]
+fn lighting_l0_draws_no_shadows() {
+    let _devices = common::device_lock();
+    let full = lit(1.0);
+    for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
+        let mut renderer = renderer(&adapter, LightingTier::L0);
+        let name = renderer.adapter_info().name.clone();
+        let model = renderer.load_model(&plate_under(2.0)).unwrap();
+        let frame = renderer
+            .submit(&scene(vec![at(model, Vector3::zeros())]), &slanted())
+            .unwrap()
+            .wait()
+            .unwrap();
+        for (z, r) in column(&slanted(), &frame.radiance.unwrap()) {
+            if z > -1.3 {
+                assert!((r - full).abs() < 1e-4 * full, "{name}: z = {z}: {r}");
+            }
+        }
+    }
+}