Limb darkening for the sun

crates/resolved-cli/src/main.rs

@@ -46,7 +46,8 @@ struct RenderArgs {
     #[arg(long)]
     preview: Option<PathBuf>,
     /// Lighting fidelity tier: L0, or L1 (adds Earth's surface, shadows on targets,
-    /// light from Earth and the Moon, and eclipses on them; needs Blue Marble data).
+    /// light from Earth and the Moon, eclipses on them, and the Sun's limb darkening;
+    /// needs Blue Marble data).
     #[arg(long, default_value = "L0")]
     lighting: LightingTier,
     /// Camera fidelity tier: C0.

crates/resolved/src/eclipse.rs

@@ -9,6 +9,7 @@
 use nalgebra::Vector3;
 
 use crate::radiometry::SUN_RADIUS;
+use crate::sun::LimbDarkening;
 
 /// The fraction of the solar disk visible from `point`, between 0 and 1, with each
 /// occulter given as (center, radius). Positions are in meters in any common frame.
@@ -27,6 +28,44 @@ pub fn sun_visible_fraction(
     (1.0 - covered).clamp(0.0, 1.0)
 }
 
+/// Per channel, the fraction of the Sun's light that reaches `point` past the
+/// occulters, with the solar disk's brightness spread across it as `limb` gives: Earth
+/// or the Moon over the bright center of the disk hides more than its share of the
+/// disk's area. With [`LimbDarkening::uniform`] it is [`sun_visible_fraction`].
+pub fn sun_visible_light(
+    point: Vector3<f64>,
+    sun: Vector3<f64>,
+    occulters: &[(Vector3<f64>, f64)],
+    limb: &LimbDarkening,
+) -> [f64; 4] {
+    let mut covered = [0.0; 4];
+    for (center, radius) in occulters {
+        let (a, b, c) = angles(point, sun, *center, *radius);
+        let part = limb.covered(c / a, b / a);
+        for (sum, x) in covered.iter_mut().zip(part) {
+            *sum += x;
+        }
+    }
+    covered.map(|x| (1.0 - x).clamp(0.0, 1.0))
+}
+
+/// Seen from `point`: the Sun's angular radius, the occulter's, and the angle between
+/// their centers.
+fn angles(
+    point: Vector3<f64>,
+    sun: Vector3<f64>,
+    center: Vector3<f64>,
+    radius: f64,
+) -> (f64, f64, f64) {
+    let to_sun = sun - point;
+    let to_body = center - point;
+    let a = (SUN_RADIUS / to_sun.norm()).min(1.0).asin();
+    let b = (radius / to_body.norm()).min(1.0).asin();
+    // atan2 keeps precision near zero separation, where acos of a dot product loses it.
+    let c = to_sun.cross(&to_body).norm().atan2(to_sun.dot(&to_body));
+    (a, b, c)
+}
+
 /// The fraction of the solar disk one occulter covers, seen from `point`.
 fn covered_fraction(
     point: Vector3<f64>,
@@ -34,12 +73,7 @@ fn covered_fraction(
     center: Vector3<f64>,
     radius: f64,
 ) -> f64 {
-    let to_sun = sun - point;
-    let to_body = center - point;
-    let a = (SUN_RADIUS / to_sun.norm()).min(1.0).asin();
-    let b = (radius / to_body.norm()).min(1.0).asin();
-    // atan2 keeps precision near zero separation, where acos of a dot product loses it.
-    let c = to_sun.cross(&to_body).norm().atan2(to_sun.dot(&to_body));
+    let (a, b, c) = angles(point, sun, center, radius);
     if c >= a + b {
         0.0
     } else if c <= b - a {
@@ -73,11 +107,21 @@ mod tests {
     /// Independent check: the fraction of points spread evenly over the solar disk
     /// (perpendicular to the line of sight) whose line to `point` misses Earth.
     fn sampled(point: Vector3<f64>, sun: Vector3<f64>) -> f64 {
+        sampled_light(point, sun, 150, |_| [1.0; 4])[0]
+    }
+
+    /// The same with `n` x `n` points, each weighted by the disk's brightness there,
+    /// given `mu`.
+    fn sampled_light(
+        point: Vector3<f64>,
+        sun: Vector3<f64>,
+        n: usize,
+        brightness: impl Fn(f64) -> [f64; 4],
+    ) -> [f64; 4] {
         let axis = (sun - point).normalize();
         let u = axis.cross(&Vector3::z()).normalize();
         let v = axis.cross(&u);
-        let n = 150;
-        let (mut inside, mut visible) = (0, 0);
+        let (mut inside, mut visible) = ([0.0; 4], [0.0; 4]);
         for i in 0..n {
             for j in 0..n {
                 let (s, t) = (
@@ -87,17 +131,24 @@ mod tests {
                 if s * s + t * t > 1.0 {
                     continue;
                 }
-                inside += 1;
+                let w = brightness((1.0 - s * s - t * t).sqrt());
+                for k in 0..4 {
+                    inside[k] += w[k];
+                }
                 let target = sun + SUN_RADIUS * (s * u + t * v);
                 let d = (target - point).normalize();
                 // Does the ray from `point` toward `target` hit Earth?
                 let b = point.dot(&d);
                 let q = point - b * d;
                 let blocked = b < 0.0 && q.norm() < RE;
-                visible += usize::from(!blocked);
+                if !blocked {
+                    for k in 0..4 {
+                        visible[k] += w[k];
+                    }
+                }
             }
         }
-        visible as f64 / inside as f64
+        std::array::from_fn(|k| visible[k] / inside[k])
     }
 
     #[test]
@@ -135,6 +186,45 @@ mod tests {
         );
     }
 
+    #[test]
+    fn limb_darkened_penumbra_matches_sampling_the_disk() {
+        // Through Earth's penumbra again, now weighting the disk by its brightness per
+        // channel (example camera). Near the umbra only a sliver of the dim limb is
+        // left, so less light gets through than the disk's area says (4.5 % against
+        // 6.1 % in green 25 km inside); near the outside it is the dim limb that is
+        // covered, so more gets through (0.882 against 0.860). Halfway they agree. A 150 x 150 grid of samples is
+        // good to about 5e-3 here (the uniform test above allows that much); 400 x 400
+        // agrees to 7e-4, the same in every channel, so the grid's geometry rather than
+        // the brightness.
+        let r = crate::radiometry::Radiometry::new(&crate::camera::CameraConfig::example());
+        let limb = LimbDarkening::new(&r);
+        let sun = Vector3::new(AU, 0.0, 0.0);
+        for dy in [-20e3, -5e3, 0.0, 5e3, 20e3] {
+            let point = Vector3::new(-7.0e6, RE + dy, 0.0);
+            let got = sun_visible_light(point, sun, &earth(), &limb);
+            let area = sun_visible_fraction(point, sun, &earth());
+            let want = sampled_light(point, sun, 400, |mu| limb.profile(mu));
+            for k in 0..3 {
+                assert_relative_eq!(got[k], want[k], epsilon = 1e-3);
+            }
+            // Darker than the area says toward the umbra, brighter toward the outside.
+            if dy < -10e3 {
+                assert!(got[1] < area, "{dy}: {got:?} {area}");
+            } else if dy > 10e3 {
+                assert!(got[1] > area, "{dy}: {got:?} {area}");
+            }
+        }
+        // A uniform disk gives the area overlap (measured: to 6e-7, with Earth 250 solar
+        // radii across here).
+        let point = Vector3::new(-7.0e6, RE, 0.0);
+        let uniform = sun_visible_light(point, sun, &earth(), &LimbDarkening::uniform());
+        assert_relative_eq!(
+            uniform[0],
+            sun_visible_fraction(point, sun, &earth()),
+            epsilon = 2e-6
+        );
+    }
+
     #[test]
     fn annular_takes_a_bite() {
         // A body small enough to sit inside the solar disk: covered fraction is the

crates/resolved/src/fidelity.rs

@@ -22,9 +22,10 @@ pub enum LightingTier {
     /// Blue Marble) with a stand-in for the atmosphere's blue; shadows that targets
     /// cast on themselves and on each other, with the penumbra the Sun's disk gives
     /// them; sunlight reflected onto targets by Earth and the Moon (not shadowed by the
-    /// targets); and Earth and the Moon shadowing each other, lunar and solar eclipses.
-    /// The rest of M2 (the Moon's albedo map, limb darkening) joins this tier as it
-    /// lands.
+    /// targets); Earth and the Moon shadowing each other, lunar and solar eclipses; and
+    /// a Sun brighter at the center of its disk than at the edge (limb darkening), in
+    /// its image and in every eclipse. The rest of M2 (the Moon's albedo map) joins
+    /// this tier as it lands.
     L1,
 }
 

crates/resolved/src/gpu/bodies.rs

@@ -15,11 +15,54 @@ use crate::earth::EarthTexture;
 use crate::model::rgb_spectrum;
 use crate::radiometry::{Radiometry, AU, SUN_RADIUS};
 use crate::scene::{Pose, Scene};
+use crate::sun::{ring_rule, LimbDarkening, RINGS, SAMPLES};
 use crate::Error;
 
 /// Channels one radiance target holds.
 pub const MAX_CHANNELS: usize = 4;
 
+/// The solar disk's brightness spread, as `SunDisk` in `common.wgsl` takes it: the
+/// tables of [`LimbDarkening`] and its rule for partly covered circles.
+#[repr(C)]
+#[derive(Clone, Copy, Debug, PartialEq, bytemuck::Pod, bytemuck::Zeroable)]
+pub struct SunDisk {
+    profile: [[f32; 4]; SAMPLES],
+    inside: [[f32; 4]; SAMPLES],
+    rule: [[f32; 4]; RINGS],
+}
+
+impl SunDisk {
+    /// The uniform for `limb`.
+    pub fn new(limb: &LimbDarkening) -> Self {
+        let [profile, inside] = limb.tables();
+        // The cosine, and the sine times the weight over pi: all the shader needs of each.
+        let rule = ring_rule().map(|(angle, weight)| {
+            [
+                angle.cos() as f32,
+                (angle.sin() * weight / std::f64::consts::PI) as f32,
+                0.0,
+                0.0,
+            ]
+        });
+        Self {
+            profile,
+            inside,
+            rule,
+        }
+    }
+
+    /// A buffer holding it, for a pass to bind.
+    pub fn buffer(ctx: &GpuContext, limb: &LimbDarkening) -> wgpu::Buffer {
+        use wgpu::util::DeviceExt;
+        ctx.device
+            .create_buffer_init(&wgpu::util::BufferInitDescriptor {
+                label: Some("sun disk"),
+                contents: bytemuck::bytes_of(&Self::new(limb)),
+                usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST,
+            })
+    }
+}
+
 /// Per-channel constants for the bodies that do not change from frame to frame.
 #[derive(Clone, Copy, Debug, PartialEq)]
 pub struct BodyRadiometry {
@@ -177,6 +220,7 @@ impl BodiesUniforms {
 pub struct BodiesPass {
     pipeline: wgpu::RenderPipeline,
     uniforms: wgpu::Buffer,
+    sun_disk: wgpu::Buffer,
     sampler: wgpu::Sampler,
     earth: wgpu::TextureView,
     earth_height: u32,
@@ -244,10 +288,13 @@ impl BodiesPass {
         });
         // A 1 x 1 placeholder until a texture is set; the shader ignores it then.
         let placeholder = EarthTexture::from_rgba8(1, 1, vec![0, 0, 0, 255]).expect("1 x 1");
-        let (earth, bind_group) = Self::bind(ctx, &pipeline, &uniforms, &sampler, &placeholder);
+        let sun_disk = SunDisk::buffer(ctx, &LimbDarkening::uniform());
+        let (earth, bind_group) =
+            Self::bind(ctx, &pipeline, &uniforms, &sun_disk, &sampler, &placeholder);
         Self {
             pipeline,
             uniforms,
+            sun_disk,
             sampler,
             earth,
             earth_height: placeholder.height(),
@@ -257,11 +304,24 @@ impl BodiesPass {
 
     /// Uploads `texture` with its mip chain and binds it for later frames.
     pub fn set_earth_texture(&mut self, ctx: &GpuContext, texture: &EarthTexture) {
-        (self.earth, self.bind_group) =
-            Self::bind(ctx, &self.pipeline, &self.uniforms, &self.sampler, texture);
+        (self.earth, self.bind_group) = Self::bind(
+            ctx,
+            &self.pipeline,
+            &self.uniforms,
+            &self.sun_disk,
+            &self.sampler,
+            texture,
+        );
         self.earth_height = texture.height();
     }
 
+    /// Draws the Sun with `limb`'s spread of brightness, and eclipses on Earth and the
+    /// Moon with it, from now on. Until called, the disk is uniform.
+    pub fn set_limb_darkening(&self, ctx: &GpuContext, limb: &LimbDarkening) {
+        ctx.queue
+            .write_buffer(&self.sun_disk, 0, bytemuck::bytes_of(&SunDisk::new(limb)));
+    }
+
     /// Earth's texture (a 1 x 1 placeholder until one is set), the sampler that reads
     /// it and its height in texels, for the Earthshine pass.
     pub fn earth_texture(&self) -> (&wgpu::TextureView, &wgpu::Sampler, u32) {
@@ -272,6 +332,7 @@ impl BodiesPass {
         ctx: &GpuContext,
         pipeline: &wgpu::RenderPipeline,
         uniforms: &wgpu::Buffer,
+        sun_disk: &wgpu::Buffer,
         sampler: &wgpu::Sampler,
         texture: &EarthTexture,
     ) -> (wgpu::TextureView, wgpu::BindGroup) {
@@ -329,6 +390,10 @@ impl BodiesPass {
                     binding: 2,
                     resource: wgpu::BindingResource::Sampler(sampler),
                 },
+                wgpu::BindGroupEntry {
+                    binding: 3,
+                    resource: sun_disk.as_entire_binding(),
+                },
             ],
         });
         (view, group)

crates/resolved/src/gpu/earthshine.rs

@@ -7,12 +7,13 @@
 
 use nalgebra::{UnitQuaternion, Vector3};
 
-use super::bodies::BodyRadiometry;
+use super::bodies::{BodyRadiometry, SunDisk};
 use super::profile::PassTimestamps;
 use super::GpuContext;
 use crate::bodies::{EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS};
 use crate::radiometry::SUN_RADIUS;
 use crate::scene::{Pose, Scene};
+use crate::sun::LimbDarkening;
 
 /// Rings of cells over each body's cone, from its center out to its limb.
 pub const RINGS: u32 = 32;
@@ -112,6 +113,7 @@ pub struct EarthshinePass {
     irradiance: wgpu::ComputePipeline,
     layout: wgpu::BindGroupLayout,
     globals: wgpu::Buffer,
+    sun_disk: wgpu::Buffer,
     lights: wgpu::Buffer,
     samples: wgpu::Buffer,
     cube: wgpu::Texture,
@@ -165,6 +167,7 @@ impl EarthshinePass {
                         view_dimension: wgpu::TextureViewDimension::D2Array,
                     },
                 ),
+                entry(6, buffer(wgpu::BufferBindingType::Uniform)),
             ],
         });
         let pipeline_layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
@@ -204,6 +207,7 @@ impl EarthshinePass {
             irradiance: pipeline("irradiance"),
             layout,
             globals,
+            sun_disk: SunDisk::buffer(ctx, &LimbDarkening::uniform()),
             lights,
             samples,
             cube,
@@ -249,6 +253,13 @@ impl EarthshinePass {
         (lights, samples, cube)
     }
 
+    /// Shades eclipses on Earth and the Moon with `limb`'s spread of the Sun's
+    /// brightness from now on. Until called, the disk is uniform.
+    pub fn set_limb_darkening(&self, ctx: &GpuContext, limb: &LimbDarkening) {
+        ctx.queue
+            .write_buffer(&self.sun_disk, 0, bytemuck::bytes_of(&SunDisk::new(limb)));
+    }
+
     /// Makes room for `targets` targets, reallocating only when the count grows, and
     /// returns the cubes as a cube array for the mesh pass to sample (layer i is
     /// target i).
@@ -336,6 +347,10 @@ impl EarthshinePass {
                     binding: 5,
                     resource: wgpu::BindingResource::TextureView(&cube),
                 },
+                wgpu::BindGroupEntry {
+                    binding: 6,
+                    resource: self.sun_disk.as_entire_binding(),
+                },
             ],
         });
         let n = lights.len() as u32;

crates/resolved/src/gpu/mesh.rs

@@ -95,11 +95,13 @@ pub struct MeshDraw {
     shadow: [f32; 4],
     // Tangent of the Sun's angular diameter, normal offset in meters, map size in texels.
     shadow_more: [f32; 4],
+    sunlight: [f32; 4],
 }
 
 impl MeshDraw {
     /// The uniforms to draw a model at `target` seen from `camera`, lit by the Sun at
-    /// `sun` (inertial, meters) with fraction `sun_visible` of its disk unobstructed,
+    /// `sun` (inertial, meters) with fraction `sun_visible` of its light, per channel,
+    /// past Earth and the Moon,
     /// shadowed through `shadow` and lit by Earth and the Moon through layer
     /// `earthshine` of the irradiance cubes, if given.
     #[allow(clippy::too_many_arguments)]
@@ -109,7 +111,7 @@ impl MeshDraw {
         intrinsics: &Intrinsics,
         near: f64,
         sun: Vector3<f64>,
-        sun_visible: f64,
+        sun_visible: [f64; 4],
         shadow: Option<ShadowDraw>,
         earthshine: Option<u32>,
     ) -> Self {
@@ -137,7 +139,7 @@ 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,
                 (1.0 / (au * au)) as f32,
                 f32::from(u8::from(earthshine.is_some())),
                 earthshine.unwrap_or(0) as f32,
@@ -159,6 +161,7 @@ impl MeshDraw {
                     0.0,
                 ]
             }),
+            sunlight: sun_visible.map(|v| (v / (au * au)) as f32),
         }
     }
 }

crates/resolved/src/lib.rs

@@ -33,6 +33,7 @@ pub mod scene;
 pub mod shadow;
 mod solar_data;
 pub mod spectrum;
+pub mod sun;
 
 pub use camera::{CameraConfig, Channel, Detector, Intrinsics};
 pub use earth::EarthTexture;

crates/resolved/src/radiometry.rs

@@ -155,8 +155,20 @@ impl Radiometry {
             / (std::f64::consts::PI * au * au)
     }
 
-    /// Channel radiance of the solar disk, treated as uniformly bright. It does not depend
-    /// on distance.
+    /// Sunlight through channel `channel`, weighted by `weight` at each wavelength:
+    /// `∫ E_sun τ QE_k λ / (h c) weight(λ) dλ`. With a weight of 1 it is `solar_flux[k]`.
+    pub(crate) fn channel_integral(&self, channel: usize, weight: impl Fn(f64) -> f64) -> f64 {
+        let spectra = [
+            solar_spectrum(),
+            &self.transmission,
+            &self.response[channel],
+        ];
+        integrate_weighted(&spectra, |l| weight(l) * l / (PLANCK * SPEED_OF_LIGHT))
+            .expect("solar spectrum is bounded")
+    }
+
+    /// Channel radiance of the solar disk, its mean over the disk. It does not depend on
+    /// distance. [`crate::sun::LimbDarkening`] gives how it varies across the disk.
     pub fn sun_radiance(&self, channel: usize) -> f64 {
         self.solar_flux[channel] / sun_solid_angle_at_1au()
     }

crates/resolved/src/render.rs

@@ -15,7 +15,7 @@ use crate::camera::CameraConfig;
 use std::sync::Arc;
 
 use crate::earth::EarthTexture;
-use crate::eclipse::sun_visible_fraction;
+use crate::eclipse::{sun_visible_fraction, sun_visible_light};
 use crate::fidelity::{Fidelity, LightingTier};
 use crate::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry};
 use crate::gpu::earthshine::{EarthshinePass, TargetLight};
@@ -29,6 +29,7 @@ use crate::model::{Mesh, ModelId};
 use crate::radiometry::Radiometry;
 use crate::scene::{Pose, Scene};
 use crate::shadow::{self, Bounded};
+use crate::sun::LimbDarkening;
 use crate::Error;
 
 /// Choices made once, when the renderer is created.
@@ -103,6 +104,8 @@ pub struct Renderer {
     dn_staging: wgpu::Buffer,
     radiance_staging: Option<(wgpu::Buffer, u32)>,
     profiler: Option<Profiler>,
+    /// How the Sun's brightness is spread across its disk: uniform below lighting L1.
+    limb: LimbDarkening,
 }
 
 impl Renderer {
@@ -164,7 +167,17 @@ impl Renderer {
             dn_staging,
             radiance_staging,
             profiler,
+            limb: LimbDarkening::uniform(),
         };
+        if options.fidelity.lighting >= LightingTier::L1 {
+            renderer.limb = LimbDarkening::new(&renderer.radiometry);
+            renderer
+                .bodies
+                .set_limb_darkening(&renderer.ctx, &renderer.limb);
+            renderer
+                .earthshine
+                .set_limb_darkening(&renderer.ctx, &renderer.limb);
+        }
         if let Some(texture) = &options.earth_texture {
             renderer.set_earth_texture(texture);
         }
@@ -216,7 +229,8 @@ impl Renderer {
     /// returns without waiting for it.
     ///
     /// The camera's pose is `observer.then(&camera.mount)`. Each target is lit by the
-    /// fraction of the solar disk Earth and the Moon leave visible from its center, and
+    /// fraction of the Sun's light Earth and the Moon leave from its center (from
+    /// lighting L1 weighted by the disk's limb darkening, per channel), and
     /// from lighting L1 also by the light Earth and the Moon reflect onto it, computed
     /// at its bounding center.
     pub fn submit(&mut self, scene: &Scene, observer: &Pose) -> Result<Pending<'_>, Error> {
@@ -285,7 +299,12 @@ impl Renderer {
             .iter()
             .enumerate()
             .map(|(i, t)| {
-                let visible = sun_visible_fraction(t.pose.position, scene.sun, &occulters);
+                // From lighting L1, through the limb-darkened disk, per channel.
+                let visible = if l1 {
+                    sun_visible_light(t.pose.position, scene.sun, &occulters, &self.limb)
+                } else {
+                    [sun_visible_fraction(t.pose.position, scene.sun, &occulters); 4]
+                };
                 let shadow = plans.get(i).map(|plan| ShadowDraw {
                     light_from_body: plan.light_from_body(&bounded[i]),
                     layer: i as u32,

crates/resolved/src/shaders/bodies.wgsl

@@ -34,7 +34,7 @@ struct Uniforms {
     moon_lambert: vec4<f32>,
     // xyz: Sun center. w: radius.
     sun_center: vec4<f32>,
-    // Per channel: channel radiance of the solar disk.
+    // Per channel: channel radiance of the solar disk, its mean over the disk.
     sun_radiance: vec4<f32>,
 }
 
@@ -42,6 +42,7 @@ struct Uniforms {
 // Equirectangular reflectance: longitude -180..180 left to right, latitude +90 at the top.
 @group(0) @binding(1) var earth_texture: texture_2d<f32>;
 @group(0) @binding(2) var earth_sampler: sampler;
+@group(0) @binding(3) var<uniform> sun_disk: SunDisk;
 
 struct Out {
     @location(0) radiance: vec4<f32>,
@@ -124,7 +125,13 @@ fn fs_main(@builtin(position) pos: vec4<f32>) -> Out {
     let ts = sphere_hit(-u.sun_center.xyz, d, u.sun_center.w);
     if ts > 0.0 && ts < best {
         best = ts;
-        radiance = u.sun_radiance;
+        // The Sun's surface seen at `mu`, from the ray's miss distance from its center:
+        // mu = sqrt(r^2 - miss^2) / r, which keeps its precision 1 AU away, where the hit
+        // point's normal would not.
+        let o = -u.sun_center.xyz;
+        let miss = length(o - dot(o, d) * d);
+        let r = u.sun_center.w;
+        radiance = u.sun_radiance * disk_profile(sqrt(max((r - miss) * (r + miss), 0.0)) / r);
     }
 
     if best > 3.0e38 {

crates/resolved/src/shaders/common.wgsl

@@ -58,32 +58,80 @@ fn textured_albedo(rgb: vec3<f32>, rgb_to_channels: mat3x4<f32>, haze: vec4<f32>
     return rgb_to_channels * rgb + haze;
 }
 
-// The fraction of the solar disk visible from `p` past one spherical occulter, with
-// the Sun and the occulter given as (center, radius), all relative to the same origin.
-//
-// The overlap of two circles on the sky, as in eclipse.rs (Montenbruck and Gill). The
-// lens area uses atan2 where eclipse.rs uses acos: acos loses its precision in f32
-// near 1, which is where a sliver of the Sun is covered.
-fn sun_visible(p: vec3<f32>, sun: vec4<f32>, occulter: vec4<f32>) -> f32 {
+// How the solar disk's brightness is spread, per channel (sun.rs, `LimbDarkening`). A
+// shader that includes this file declares a uniform `sun_disk` of this type.
+struct SunDisk {
+    // Radiance over the disk's mean, at mu = (i / 63)^4.
+    profile: array<vec4<f32>, 64>,
+    // The fraction of the disk's light inside radius sqrt(i / 63), in solar radii.
+    inside: array<vec4<f32>, 64>,
+    // A Gauss-Legendre rule on [0, pi] for the partly covered circles: x the cosine of
+    // each point's angle, y its sine times its weight over pi.
+    rule: array<vec4<f32>, 8>,
+}
+
+// Per channel, the radiance where the line of sight meets the Sun at `mu`, over the
+// disk's mean.
+fn disk_profile(mu: f32) -> vec4<f32> {
+    let x = sqrt(sqrt(clamp(mu, 0.0, 1.0))) * 63.0;
+    let i = min(u32(x), 62u);
+    return mix(sun_disk.profile[i], sun_disk.profile[i + 1u], x - f32(i));
+}
+
+// Per channel, the fraction of the disk's light from inside `r` solar radii.
+fn disk_inside(r: f32) -> vec4<f32> {
+    let x = clamp(r * r, 0.0, 1.0) * 63.0;
+    let i = min(u32(x), 62u);
+    return mix(sun_disk.inside[i], sun_disk.inside[i + 1u], x - f32(i));
+}
+
+// Half the angle, from the Sun's center, of the arc of the circle of radius `r` around
+// it that lies inside a disk of radius `b` centered `c` away. The sine is a product of
+// four factors, so neither it nor the cosine is a small difference of large numbers.
+fn covered_angle(r: f32, c: f32, b: f32) -> f32 {
+    let sin2 = (b - r + c) * (b + r - c) * (r + c - b) * (r + c + b);
+    return atan2(sqrt(max(sin2, 0.0)), r * r + c * c - b * b);
+}
+
+// Per channel, the fraction of the disk's light hidden by a disk of radius `b` centered
+// `c` from the Sun's, in solar radii: `LimbDarkening::covered` in sun.rs, which explains
+// the method.
+fn disk_covered(c: f32, b: f32) -> vec4<f32> {
+    if c >= 1.0 + b {
+        return vec4<f32>(0.0);
+    }
+    if b >= c + 1.0 {
+        return vec4<f32>(1.0);
+    }
+    var covered = vec4<f32>(0.0);
+    if b > c {
+        covered = disk_inside(b - c);
+    }
+    let lo = abs(c - b);
+    let hi = min(c + b, 1.0);
+    if hi > lo {
+        let mid = (hi + lo) * 0.5;
+        let half = (hi - lo) * 0.5;
+        for (var j = 0u; j < 8u; j++) {
+            let node = sun_disk.rule[j];
+            let r = mid - half * node.x;
+            let arc = 2.0 * covered_angle(r, c, b) * r;
+            let mu = sqrt(max(1.0 - r * r, 0.0));
+            covered += disk_profile(mu) * (arc * half * node.y);
+        }
+    }
+    return clamp(covered, vec4<f32>(0.0), vec4<f32>(1.0));
+}
+
+// Per channel, the fraction of the Sun's light that reaches `p` past one spherical
+// occulter, with the Sun and the occulter given as (center, radius), all relative to the
+// same origin: eclipse.rs's `sun_visible_light`.
+fn sun_visible(p: vec3<f32>, sun: vec4<f32>, occulter: vec4<f32>) -> vec4<f32> {
     let to_sun = sun.xyz - p;
     let to_body = occulter.xyz - p;
     let a = asin(min(sun.w / length(to_sun), 1.0));
     let b = asin(min(occulter.w / length(to_body), 1.0));
     // atan2 keeps precision near zero separation, where acos of a dot product loses it.
     let c = atan2(length(cross(to_sun, to_body)), dot(to_sun, to_body));
-    if c >= a + b {
-        return 1.0;
-    }
-    if c <= b - a {
-        return 0.0;
-    }
-    if c <= a - b {
-        return 1.0 - (b / a) * (b / a);
-    }
-    // A lens-shaped overlap: the circles cross at distance x along the line between
-    // their centers from the Sun's, y either side of it.
-    let x = (c * c + (a - b) * (a + b)) / (2.0 * c);
-    let y = sqrt(max((a - x) * (a + x), 0.0));
-    let area = a * a * atan2(y, x) + b * b * atan2(y, c - x) - c * y;
-    return clamp(1.0 - area / (PI * a * a), 0.0, 1.0);
+    return vec4<f32>(1.0) - disk_covered(c / a, b / a);
 }

crates/resolved/src/shaders/earthshine.wgsl

@@ -70,6 +70,7 @@ struct Sample {
 @group(0) @binding(3) var earth_texture: texture_2d<f32>;
 @group(0) @binding(4) var earth_sampler: sampler;
 @group(0) @binding(5) var cube: texture_storage_2d_array<rgba16float, write>;
+@group(0) @binding(6) var<uniform> sun_disk: SunDisk;
 
 // Two unit vectors perpendicular to unit `a` and to each other.
 fn basis(a: vec3<f32>) -> mat2x3<f32> {
@@ -133,7 +134,7 @@ fn sources(@builtin(global_invocation_id) id: vec3<u32>) {
     var albedo = vec4<f32>(0.0);
     var normal = vec3<f32>(0.0);
     // Each body in the other's shadow (eclipses) as the bodies pass draws it.
-    var visible = 1.0;
+    var visible = vec4<f32>(1.0);
     if !is_moon && te > 0.0 && !(tm > 0.0 && tm < te) {
         let surface = earth_surface(t, d, te);
         albedo = surface[0];
@@ -146,7 +147,7 @@ fn sources(@builtin(global_invocation_id) id: vec3<u32>) {
     }
     let lit = max(dot(normal, source.sun.xyz), 0.0) * visible;
     var out: Sample;
-    out.light = albedo * (lit * source.sun.w * source.cone.y);
+    out.light = albedo * lit * (source.sun.w * source.cone.y);
     out.direction = vec4<f32>(d, 0.0);
     samples[target_index * 2u * cells + id.x] = out;
 }

crates/resolved/src/shaders/mesh.wgsl

@@ -12,8 +12,8 @@ struct Draw {
     translation: vec4<f32>,
     // xyz: unit vector from the target toward the Sun, camera frame.
     sun: vec4<f32>,
-    // x: Sun visibility (eclipse) / (Sun distance in AU)^2. y: 1 / (Sun distance in
-    // AU)^2. z: 1 if lit by Earth and the Moon. w: this target's irradiance cube.
+    // y: 1 / (Sun distance in AU)^2. z: 1 if lit by Earth and the Moon. w: this
+    // target's irradiance cube.
     params: vec4<f32>,
     // Body frame to this target's shadow map (shadow.rs).
     light_from_body: mat4x4<f32>,
@@ -21,6 +21,9 @@ struct Draw {
     shadow: vec4<f32>,
     // Tangent of the Sun's angular diameter, normal offset in meters, map size in texels.
     shadow_more: vec4<f32>,
+    // Per channel: the fraction of sunlight Earth and the Moon leave (eclipses) /
+    // (Sun distance in AU)^2.
+    sunlight: vec4<f32>,
 }
 
 @group(0) @binding(0) var<uniform> draw: Draw;
@@ -140,5 +143,5 @@ fn fs_main(in: VertexOut) -> @location(0) vec4<f32> {
         let stored = textureSampleLevel(earthshine, earthshine_sampler, n, i32(draw.params.w), 0.0);
         bodies = stored / (EARTHSHINE_SCALE * PI);
     }
-    return parts[in.part] * (lit * shade * draw.params.x + bodies * draw.params.y);
+    return parts[in.part] * (lit * shade * draw.sunlight + bodies * draw.params.y);
 }

crates/resolved/src/sun.rs

@@ -0,0 +1,410 @@
+//! The Sun's disk is brighter at its center than at its edge: limb darkening.
+//!
+//! Toward the edge of the disk, the line of sight meets the Sun's surface at a slant
+//! and sees higher, cooler gas, which shines less. At a point of the disk where the
+//! line of sight makes an angle with the Sun's surface normal whose cosine is `mu`
+//! (1 at the center, 0 at the edge; `mu = sqrt(1 - r^2)` at `r` solar radii from the
+//! center), the radiance relative to the center's is close to `mu^alpha`, with `alpha`
+//! depending on wavelength: about 0.94 at 303 nm, 0.48 at 580 nm and 0.26 at 1099 nm.
+//! Hestroffer and Magnan (1998, A&A 333, 338) fit that law to the limb darkening
+//! Neckel and Labs (1994, Solar Physics 153, 91) measured at 30 wavelengths from 303 to
+//! 1099 nm, and it reproduces the measured profiles to about 1 % (their fit excludes
+//! the outermost 10 % of the radius, where measurements are least certain).
+//!
+//! Over the disk, `mu^alpha` averages `2 / (alpha + 2)`, so a wavelength's radiance at
+//! `mu` is its disk-averaged radiance times `(alpha + 2) / 2 * mu^alpha`, and the share
+//! of its light from inside radius `r` is `1 - mu^(alpha + 2)`. Both are folded into
+//! each camera channel exactly, weighting wavelengths as the channel's solar flux does,
+//! and tabulated. The disk-averaged radiance stays the measured solar spectrum's.
+
+use crate::radiometry::Radiometry;
+
+/// Neckel and Labs' limb darkening as exponents of `mu`, by wavelength in nm: Table 2 of
+/// Hestroffer and Magnan (1998), the columns fitted to Neckel and Labs' data.
+const NECKEL_LABS: [(f64, f64); 30] = [
+    (303.327, 0.939),
+    (310.843, 0.914),
+    (320.468, 0.871),
+    (329.897, 0.817),
+    (349.947, 0.763),
+    (365.875, 0.695),
+    (374.086, 0.750),
+    (390.915, 0.791),
+    (401.970, 0.751),
+    (416.319, 0.724),
+    (427.930, 0.682),
+    (443.885, 0.649),
+    (445.125, 0.646),
+    (457.345, 0.628),
+    (477.427, 0.594),
+    (492.905, 0.570),
+    (519.930, 0.538),
+    (541.760, 0.514),
+    (559.950, 0.496),
+    (579.880, 0.477),
+    (610.975, 0.447),
+    (640.970, 0.428),
+    (669.400, 0.407),
+    (700.875, 0.386),
+    (748.710, 0.361),
+    (811.760, 0.333),
+    (869.600, 0.311),
+    (948.850, 0.292),
+    (1046.600, 0.274),
+    (1098.950, 0.260),
+];
+
+/// The exponent `alpha` of the limb darkening `mu^alpha` at `wavelength` (meters),
+/// interpolated linearly between Neckel and Labs' wavelengths and held at the end values
+/// outside 303 to 1099 nm.
+pub fn limb_exponent(wavelength: f64) -> f64 {
+    let nm = wavelength * 1e9;
+    let (first, last) = (NECKEL_LABS[0], NECKEL_LABS[NECKEL_LABS.len() - 1]);
+    if nm <= first.0 {
+        return first.1;
+    }
+    if nm >= last.0 {
+        return last.1;
+    }
+    let i = NECKEL_LABS.partition_point(|&(l, _)| l <= nm);
+    let ((l0, a0), (l1, a1)) = (NECKEL_LABS[i - 1], NECKEL_LABS[i]);
+    a0 + (a1 - a0) * (nm - l0) / (l1 - l0)
+}
+
+/// Entries in each table, spaced evenly in a variable chosen so the table is smooth
+/// everywhere across the disk. The profile's sample `i` is at `mu = s^4`, `s = i /
+/// (SAMPLES - 1)`: in `s`, `mu^alpha` is `s^(4 alpha)` with `4 alpha` at least 1, smooth
+/// right to the limb, where it is steepest in `mu`. The light inside a circle is
+/// sampled evenly in its radius squared, `r^2 = i / (SAMPLES - 1)`, in which a uniform
+/// disk's is exactly linear.
+pub const SAMPLES: usize = 64;
+
+/// Points across the partly covered rings of the disk in [`LimbDarkening::covered`].
+/// Measured against 32 points over occulters from 0.2 to 250 solar radii: 2.4e-5 at
+/// most (16 points: 1e-5 or so, at twice the cost on the GPU). A midpoint rule, tried
+/// first, was 2e-4 off at 16.
+pub const RINGS: usize = 8;
+
+/// Gauss-Legendre nodes and weights for [`RINGS`] points on [0, pi], as (angle, weight).
+/// Found once by Newton's method on the Legendre polynomial's roots.
+pub fn ring_rule() -> &'static [(f64, f64); RINGS] {
+    static RULE: std::sync::OnceLock<[(f64, f64); RINGS]> = std::sync::OnceLock::new();
+    RULE.get_or_init(|| {
+        let n = RINGS;
+        std::array::from_fn(|i| {
+            // Root i of P_n, from the usual starting guess.
+            let mut x = (std::f64::consts::PI * (i as f64 + 0.75) / (n as f64 + 0.5)).cos();
+            let mut derivative = 0.0;
+            for _ in 0..100 {
+                let (mut p0, mut p1) = (1.0, x);
+                for k in 2..=n {
+                    let p2 = ((2 * k - 1) as f64 * x * p1 - (k - 1) as f64 * p0) / k as f64;
+                    p0 = p1;
+                    p1 = p2;
+                }
+                derivative = n as f64 * (x * p1 - p0) / (x * x - 1.0);
+                let step = p1 / derivative;
+                x -= step;
+                if step.abs() < 1e-15 {
+                    break;
+                }
+            }
+            let weight = 2.0 / ((1.0 - x * x) * derivative * derivative);
+            // From [-1, 1] to [0, pi].
+            let half = std::f64::consts::PI / 2.0;
+            (half * (x + 1.0), half * weight)
+        })
+    })
+}
+
+/// How each camera channel's sunlight is spread over the solar disk.
+#[derive(Clone, Debug, PartialEq)]
+pub struct LimbDarkening {
+    /// Per sample, per channel: radiance over the disk's mean.
+    profile: [[f64; 4]; SAMPLES],
+    /// Per sample, per channel: the fraction of the disk's light from inside the
+    /// sample's circle.
+    inside: [[f64; 4]; SAMPLES],
+}
+
+/// Where the profile's sample `i` is.
+fn mu_at(i: usize) -> f64 {
+    (i as f64 / (SAMPLES - 1) as f64).powi(4)
+}
+
+/// Where the light-inside table's sample `i` is: the `mu` at the edge of its circle.
+fn inside_mu_at(i: usize) -> f64 {
+    (1.0 - i as f64 / (SAMPLES - 1) as f64).sqrt()
+}
+
+impl LimbDarkening {
+    /// A uniformly bright disk, in every channel.
+    pub fn uniform() -> Self {
+        Self {
+            profile: [[1.0; 4]; SAMPLES],
+            inside: std::array::from_fn(|i| [i as f64 / (SAMPLES - 1) as f64; 4]),
+        }
+    }
+
+    /// Neckel and Labs' limb darkening, folded through `radiometry`'s channels. Integrates
+    /// over wavelength, so it belongs in setup.
+    pub fn new(radiometry: &Radiometry) -> Self {
+        let mut out = Self::uniform();
+        for (k, flux) in radiometry.solar_flux.iter().enumerate() {
+            for i in 0..SAMPLES {
+                let mu = mu_at(i);
+                out.profile[i][k] = radiometry.channel_integral(k, |l| {
+                    let a = limb_exponent(l);
+                    (a + 2.0) / 2.0 * mu.powf(a)
+                }) / flux;
+                let edge = inside_mu_at(i);
+                out.inside[i][k] = radiometry
+                    .channel_integral(k, |l| 1.0 - edge.powf(limb_exponent(l) + 2.0))
+                    / flux;
+            }
+        }
+        out
+    }
+
+    /// Linear interpolation of `table` at `x` in [0, 1] of its spacing variable.
+    fn lookup(table: &[[f64; 4]; SAMPLES], x: f64) -> [f64; 4] {
+        let x = x.clamp(0.0, 1.0) * (SAMPLES - 1) as f64;
+        let i = (x.floor() as usize).min(SAMPLES - 2);
+        let t = x - i as f64;
+        std::array::from_fn(|k| table[i][k] + t * (table[i + 1][k] - table[i][k]))
+    }
+
+    /// Per channel, the radiance where the line of sight meets the Sun's surface at
+    /// `mu`, over the disk's mean.
+    pub fn profile(&self, mu: f64) -> [f64; 4] {
+        Self::lookup(&self.profile, mu.max(0.0).sqrt().sqrt())
+    }
+
+    /// Per channel, the fraction of the disk's light from inside `r` solar radii of its
+    /// center.
+    pub fn inside(&self, r: f64) -> [f64; 4] {
+        Self::lookup(&self.inside, r * r)
+    }
+
+    /// Per channel, the fraction of the disk's light hidden by a disk of radius `b` whose
+    /// center is `c` from the Sun's, both in solar radii (the Sun's disk being small
+    /// enough on the sky to treat as flat).
+    ///
+    /// Whatever circles around the Sun's center the occulter covers whole (those inside
+    /// radius `b - c`) count through [`Self::inside`]. The circles it covers in part, from
+    /// `|c - b|` out to `c + b` or the limb, count by a sum over [`RINGS`] radii, each
+    /// circle's brightness times the arc of it inside the occulter. The sum is a
+    /// Gauss-Legendre rule in an angle `phi` with the radius at `cos(phi)` across the
+    /// range, which removes the square-root behavior the arc has at both ends of it.
+    pub fn covered(&self, c: f64, b: f64) -> [f64; 4] {
+        if c >= 1.0 + b {
+            return [0.0; 4];
+        }
+        if b >= c + 1.0 {
+            return [1.0; 4];
+        }
+        let mut covered = if b > c { self.inside(b - c) } else { [0.0; 4] };
+        let (lo, hi) = ((c - b).abs(), (c + b).min(1.0));
+        if hi > lo {
+            let (mid, half) = ((hi + lo) / 2.0, (hi - lo) / 2.0);
+            for &(phi, step) in ring_rule() {
+                let r = mid - half * phi.cos();
+                let arc = 2.0 * covered_angle(r, c, b) * r;
+                let mu = (1.0 - r * r).max(0.0).sqrt();
+                let p = self.profile(mu);
+                let weight = arc * half * phi.sin() * step / std::f64::consts::PI;
+                for k in 0..4 {
+                    covered[k] += p[k] * weight;
+                }
+            }
+        }
+        covered.map(|x| x.clamp(0.0, 1.0))
+    }
+
+    /// The tables as the GPU takes them (`SunDisk` in `common.wgsl`): the profile, then
+    /// the fraction inside, each `SAMPLES` entries of four channels.
+    pub fn tables(&self) -> [[[f32; 4]; SAMPLES]; 2] {
+        let cast = |t: &[[f64; 4]; SAMPLES]| t.map(|row| row.map(|x| x as f32));
+        [cast(&self.profile), cast(&self.inside)]
+    }
+}
+
+/// Half the angle, seen from the Sun's center, of the arc of the circle of radius `r`
+/// around it that lies inside a disk of radius `b` centered `c` away.
+///
+/// From `cos = (r^2 + c^2 - b^2) / (2 r c)`, with the sine as a product of four factors
+/// so that neither is the small difference of large numbers near the ends of the arc.
+fn covered_angle(r: f64, c: f64, b: f64) -> f64 {
+    let sin2 = (b - r + c) * (b + r - c) * (r + c - b) * (r + c + b);
+    sin2.max(0.0).sqrt().atan2(r * r + c * c - b * b)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::camera::CameraConfig;
+    use std::f64::consts::PI;
+
+    /// The example camera's radiometry and limb darkening, built once: folding takes
+    /// 0.15 s in a debug build.
+    fn example() -> &'static (Radiometry, LimbDarkening) {
+        static EXAMPLE: std::sync::OnceLock<(Radiometry, LimbDarkening)> =
+            std::sync::OnceLock::new();
+        EXAMPLE.get_or_init(|| {
+            let r = Radiometry::new(&CameraConfig::example());
+            let limb = LimbDarkening::new(&r);
+            (r, limb)
+        })
+    }
+
+    #[test]
+    fn exponents_follow_the_table() {
+        assert_eq!(limb_exponent(579.88e-9), 0.477);
+        assert_eq!(limb_exponent(200e-9), 0.939);
+        assert_eq!(limb_exponent(2000e-9), 0.260);
+        assert!((limb_exponent(590e-9) - 0.4672).abs() < 1e-4);
+    }
+
+    #[test]
+    fn every_channel_keeps_its_light() {
+        // The whole disk holds all of it, and the profile averages 1 over the disk
+        // (measured: to 1e-4, the profile table's interpolation).
+        let (_, limb) = example();
+        let n = 4000;
+        for k in 0..3 {
+            assert!((limb.inside(1.0)[k] - 1.0).abs() < 1e-12);
+            assert!((limb.inside(0.0)[k]).abs() < 1e-12);
+            let mean: f64 = (0..n)
+                .map(|i| {
+                    let r = (i as f64 + 0.5) / n as f64;
+                    limb.profile((1.0 - r * r).sqrt())[k] * 2.0 * r / n as f64
+                })
+                .sum();
+            assert!((mean - 1.0).abs() < 2e-4, "channel {k}: {mean}");
+        }
+    }
+
+    #[test]
+    fn tables_match_folding_each_mu_directly() {
+        // Between the table's samples, against the integral over wavelength done at
+        // that mu. Measured: 1.1e-4 of the disk's mean radiance at most, in blue.
+        let (radiometry, limb) = example();
+        for mu in [0.003, 0.02, 0.1, 0.33, 0.5, 0.77, 0.95] {
+            let got = limb.profile(mu);
+            for (k, flux) in radiometry.solar_flux.iter().enumerate() {
+                let want = radiometry.channel_integral(k, |l| {
+                    let a = limb_exponent(l);
+                    (a + 2.0) / 2.0 * mu.powf(a)
+                }) / flux;
+                assert!(
+                    (got[k] - want).abs() < 2e-4,
+                    "mu {mu}, channel {k}: {got:?} {want}"
+                );
+            }
+        }
+    }
+
+    #[test]
+    fn darkening_matches_the_measured_size() {
+        // Example camera, green channel (540 nm): alpha about 0.52, so the center is
+        // (alpha + 2) / 2 = 1.26 times the mean, and at mu = 0.5 the radiance is
+        // 0.5^0.52 = 0.70 of the center's.
+        let (_, limb) = example();
+        let center = limb.profile(1.0)[1];
+        assert!((center - 1.26).abs() < 0.01, "{center}");
+        assert!((limb.profile(0.5)[1] / center - 0.70).abs() < 0.01);
+        // Blue darkens more than red.
+        let [r, _, b] = [0, 1, 2].map(|k| limb.profile(0.3)[k] / limb.profile(1.0)[k]);
+        assert!(b < r, "{r} {b}");
+    }
+
+    /// The fraction of `limb`'s channel-`k` light that a disk of radius `b` centered `c`
+    /// away hides, by summing brightness over a fine grid on the Sun's disk.
+    fn sampled_covered(limb: &LimbDarkening, k: usize, c: f64, b: f64) -> f64 {
+        let n = 600;
+        let (mut hidden, mut total) = (0.0, 0.0);
+        for i in 0..n {
+            for j in 0..n {
+                let x = 2.0 * (i as f64 + 0.5) / n as f64 - 1.0;
+                let y = 2.0 * (j as f64 + 0.5) / n as f64 - 1.0;
+                let r2 = x * x + y * y;
+                if r2 > 1.0 {
+                    continue;
+                }
+                let p = limb.profile((1.0 - r2).sqrt())[k];
+                total += p;
+                if (x - c).powi(2) + y * y < b * b {
+                    hidden += p;
+                }
+            }
+        }
+        hidden / total
+    }
+
+    #[test]
+    fn covered_matches_summing_over_the_disk() {
+        // Partial, central, annular, and the Moon's and Earth's sizes as seen in eclipses.
+        let (_, limb) = example();
+        for (c, b) in [
+            (0.0, 0.3),
+            (0.2, 0.5),
+            (0.7, 0.5),
+            (1.2, 0.5),
+            (0.5, 0.98),
+            (0.05, 0.97),
+            (2.0, 3.6),
+            (4.0, 3.6),
+        ] {
+            let got = limb.covered(c, b)[1];
+            let want = sampled_covered(limb, 1, c, b);
+            // Measured: 7e-5 at most with this 600 x 600 grid, 2.6e-5 with 1500 x 1500:
+            // the grid's own precision.
+            assert!(
+                (got - want).abs() < 1e-4,
+                "c {c}, b {b}: {got} against {want}"
+            );
+        }
+    }
+
+    #[test]
+    fn a_uniform_disk_gives_the_area_overlap() {
+        // Measured: 7e-10 at most at these sizes.
+        let limb = LimbDarkening::uniform();
+        for (c, b) in [
+            (0.0, 0.3),
+            (0.2, 0.5),
+            (0.7, 0.5),
+            (1.2, 0.5),
+            (0.5, 0.98),
+            (2.0, 3.6),
+        ] {
+            let got = limb.covered(c, b)[0];
+            // Area of the lens, over the Sun's.
+            let want = if c >= 1.0 + b {
+                0.0
+            } else if c <= b - 1.0 {
+                1.0
+            } else if c <= 1.0 - b {
+                b * b
+            } else {
+                let x = (c * c + 1.0 - b * b) / (2.0 * c);
+                let y = (1.0 - x * x).sqrt();
+                ((x).acos() + b * b * ((c - x) / b).acos() - c * y) / PI
+            };
+            assert!(
+                (got - want).abs() < 2e-6,
+                "c {c}, b {b}: {got} against {want}"
+            );
+        }
+    }
+
+    #[test]
+    fn the_bright_center_counts_for_more() {
+        // A small occulter over the center hides more than its share of the area, and
+        // one at the limb less.
+        let (_, limb) = example();
+        let b = 0.1;
+        assert!(limb.covered(0.0, b)[1] > 1.2 * b * b);
+        assert!(limb.covered(0.95, b)[1] < 0.6 * limb.covered(0.0, b)[1]);
+    }
+}

crates/resolved/tests/bodies.rs

@@ -10,6 +10,7 @@ use resolved::gpu::targets::FrameTargets;
 use resolved::gpu::GpuContext;
 use resolved::nalgebra::{UnitQuaternion, Vector3};
 use resolved::radiometry::{Radiometry, AU, SUN_RADIUS};
+use resolved::sun::LimbDarkening;
 use resolved::{CameraConfig, Pose, Scene};
 
 /// Camera distance from Earth's center: far enough that Earth fits in the example
@@ -44,6 +45,18 @@ fn render_with(
     scene: &Scene,
     camera: &Pose,
     config: &CameraConfig,
+) -> Frame {
+    render_limb(ctx, name, scene, camera, config, &LimbDarkening::uniform())
+}
+
+/// [`render_with`], with the Sun's brightness spread across its disk as `limb` says.
+fn render_limb(
+    ctx: &GpuContext,
+    name: &str,
+    scene: &Scene,
+    camera: &Pose,
+    config: &CameraConfig,
+    limb: &LimbDarkening,
 ) -> Frame {
     let targets = FrameTargets::new(
         ctx,
@@ -52,6 +65,7 @@ fn render_with(
         config.channels.len() as u32,
     );
     let pass = BodiesPass::new(ctx);
+    pass.set_limb_darkening(ctx, limb);
     let uniforms = BodiesUniforms::new(
         scene,
         camera,
@@ -175,6 +189,64 @@ fn solar_disk_size_and_radiance() {
     }
 }
 
+#[test]
+fn limb_darkened_sun_keeps_its_light() {
+    // Through a 250 mm lens the Sun is 420 pixels across. Each pixel of the disk must
+    // read the mean radiance times the profile at the mu its line of sight meets, and
+    // the disk's light summed over pixels must still be the channel's solar flux.
+    let sun = Vector3::new(0.0, AU, 0.0);
+    let scene = scene(sun, moon_away());
+    let position = Vector3::new(RANGE, 0.0, 0.0);
+    let camera = Pose::looking_at(position, sun, Vector3::z()).unwrap();
+    let config = CameraConfig {
+        focal_length: 0.25,
+        ..CameraConfig::example()
+    };
+    let radiometry = Radiometry::new(&config);
+    let limb = LimbDarkening::new(&radiometry);
+    let k = config.intrinsics();
+    let mean = body_radiometry().sun_radiance;
+    for ctx in common::contexts() {
+        let frame = render_limb(ctx, "sun-limb", &scene, &camera, &config, &limb);
+        let mut total = [0.0f64; 3];
+        let mut worst: f64 = 0.0;
+        for v in 0..config.height as usize {
+            for u in 0..config.width as usize {
+                let ray = k.ray(&resolved::nalgebra::Vector2::new(
+                    u as f64 + 0.5,
+                    v as f64 + 0.5,
+                ));
+                let d = camera.attitude * ray;
+                let to_sun = sun - position;
+                let miss = (to_sun - to_sun.dot(&d) * d).norm() / SUN_RADIUS;
+                let got = frame.at(u, v).0;
+                // Pixel solid angle: (1 / f)^2 in pixels, times cos^3 off axis.
+                let solid = ray.z.powi(3) / (k.fx * k.fy);
+                for c in 0..3 {
+                    total[c] += got[c] as f64 * solid;
+                }
+                // Clear of the limb, where a pixel's center and the shader's ray may
+                // fall either side of it.
+                if miss < 0.999 {
+                    let want = limb.profile((1.0 - miss * miss).sqrt());
+                    for c in 0..3 {
+                        let expect = mean[c] as f64 * want[c];
+                        worst = worst.max((got[c] as f64 / expect - 1.0).abs());
+                    }
+                }
+            }
+        }
+        // Measured: 5.7e-5 on the RTX 3060, 4.1e-5 on lavapipe, where f32 places a ray
+        // near the limb.
+        assert!(worst < 2e-4, "{}: {worst}", ctx.info.name);
+        for (total, flux) in total.iter().zip(&radiometry.solar_flux) {
+            let flux = flux * (AU / (sun - position).norm()).powi(2);
+            // Pixels cut the disk's edge into steps, but the edge is dim: measured 7e-5.
+            assert_relative_eq!(*total, flux, max_relative = 5e-4);
+        }
+    }
+}
+
 #[test]
 fn moon_in_front_of_earth_wins() {
     let moon = Vector3::new(RANGE / 2.0, 0.0, 0.0);

crates/resolved/tests/earthshine.rs

@@ -8,17 +8,24 @@ use std::io::Cursor;
 use std::sync::Arc;
 
 use resolved::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS};
-use resolved::eclipse::sun_visible_fraction;
+use resolved::eclipse::sun_visible_light;
 use resolved::gpu::bodies::BodyRadiometry;
 use resolved::gpu::AdapterChoice;
 use resolved::nalgebra::{Matrix3, UnitQuaternion, Vector3};
 use resolved::radiometry::{Radiometry, AU};
+use resolved::sun::LimbDarkening;
 use resolved::{
     CameraConfig, EarthTexture, Fidelity, LightingTier, LoadOptions, Material, Mesh, Pose,
     Renderer, RendererOptions, Scene, Spectrum, Target,
 };
 
 const ALBEDO: f64 = 0.5;
+
+/// The example camera's limb darkening, as the renderer folds it at L1.
+fn limb() -> &'static LimbDarkening {
+    static LIMB: std::sync::OnceLock<LimbDarkening> = std::sync::OnceLock::new();
+    LIMB.get_or_init(|| LimbDarkening::new(&Radiometry::new(&CameraConfig::example())))
+}
 const MOON_DISTANCE: f64 = 3.844e8;
 
 /// An 8 m square plate facing body +X, gray.
@@ -173,20 +180,23 @@ impl Case {
         Pose::new(self.target - 20.0 * (attitude * Vector3::z()), attitude)
     }
 
-    /// Direct sunlight on the plate, as a fraction of the Sun's irradiance there.
-    fn direct(&self) -> f64 {
-        let visible = sun_visible_fraction(
+    /// Direct sunlight on the plate at lighting L1, per channel, as a fraction of the
+    /// Sun's irradiance there.
+    fn direct(&self) -> [f64; 4] {
+        let visible = sun_visible_light(
             self.target,
             self.sun,
             &[
                 (Vector3::zeros(), EARTH_EQUATORIAL_RADIUS),
                 (self.moon, MOON_RADIUS),
             ],
+            limb(),
         );
-        self.normal
+        let cos = self
+            .normal
             .dot(&(self.sun - self.target).normalize())
-            .max(0.0)
-            * visible
+            .max(0.0);
+        visible.map(|v| cos * v)
     }
 }
 
@@ -293,10 +303,11 @@ fn quadrature(case: &Case, surfaces: &Surfaces, rings: usize, sectors: usize) ->
                 }
                 // The other body's shadow on this one: a lunar or a solar eclipse.
                 let point = p + t * d;
-                let shadow = sun_visible_fraction(point, case.sun, &[(other.center, other.radius)]);
-                let lit = body.normal(point).dot(&to_sun).max(0.0) * shadow;
-                for (s, albedo) in sum.iter_mut().zip(body.albedo) {
-                    *s += albedo * lit / PI * ratio * cell * facing;
+                let shadow =
+                    sun_visible_light(point, case.sun, &[(other.center, other.radius)], limb());
+                let lit = body.normal(point).dot(&to_sun).max(0.0);
+                for k in 0..3 {
+                    sum[k] += body.albedo[k] * lit * shadow[k] / PI * ratio * cell * facing;
                 }
             }
             sum
@@ -346,7 +357,7 @@ fn rendered(renderer: &mut Renderer, cases: &[Case]) -> Vec<[f64; 3]> {
             let distance = (case.sun - case.target).norm();
             [0, 1, 2].map(|k| {
                 pixel[k] as f64 / radiometry.lambert_radiance(k, ALBEDO, 1.0, distance)
-                    - case.direct()
+                    - case.direct()[k]
             })
         })
         .collect()
@@ -413,6 +424,14 @@ fn low_orbit() {
             Case::new("leo-side-sun100", LEO, Vector3::new(-1.0, 1.0, 0.0)).sun_at(100.0),
             // In Earth's shadow, 500 km up only the night side is in view.
             Case::new("leo-down-umbra", LEO, down).sun_at(150.0),
+            // Halfway through Earth's penumbra, 7000 km behind it, facing the Sun: the
+            // direct light is the limb-darkened disk's, per channel, and Earth's night
+            // side and lit edge are in view.
+            Case::new(
+                "leo-penumbra",
+                Vector3::new(-7.0e6, EARTH_EQUATORIAL_RADIUS, 0.0),
+                Vector3::x(),
+            ),
             // Earth's axes turned: the flattened ellipsoid moves with them.
             Case {
                 earth_attitude: tilted,

crates/resolved/tests/eclipses.rs

@@ -1,21 +1,30 @@
 //! Eclipses seen on the bodies at lighting L1, on the GPU and on the CPU (lavapipe):
 //! Earth's shadow on the Moon and the Moon's on Earth, pixel by pixel against
-//! `eclipse::sun_visible_fraction` at each pixel's surface point.
+//! `eclipse::sun_visible_light` (the limb-darkened disk, per channel) at each pixel's
+//! surface point.
 
 mod common;
 
 use std::sync::Arc;
 
 use resolved::bodies::{EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS};
-use resolved::eclipse::sun_visible_fraction;
+use resolved::eclipse::sun_visible_light;
 use resolved::gpu::AdapterChoice;
 use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3};
+use resolved::radiometry::Radiometry;
 use resolved::radiometry::AU;
+use resolved::sun::LimbDarkening;
 use resolved::{
     CameraConfig, EarthTexture, Fidelity, LightingTier, Pose, Renderer, RendererOptions, Scene,
 };
 
 const SUN: Vector3<f64> = Vector3::new(AU, 0.0, 0.0);
+
+/// The example camera's limb darkening, as the renderer folds it at L1.
+fn limb() -> &'static LimbDarkening {
+    static LIMB: std::sync::OnceLock<LimbDarkening> = std::sync::OnceLock::new();
+    LIMB.get_or_init(|| LimbDarkening::new(&Radiometry::new(&CameraConfig::example())))
+}
 const MOON_DISTANCE: f64 = 3.844e8;
 
 fn renderer(adapter: &AdapterChoice, lighting: LightingTier) -> Renderer {
@@ -97,8 +106,8 @@ fn sampled() -> impl Iterator<Item = (usize, usize)> {
 /// both devices, on every core.
 fn expected(
     point: impl Fn(usize, usize) -> Option<Vector3<f64>> + Sync,
-    visible: impl Fn(Vector3<f64>) -> f64 + Sync,
-) -> Vec<Option<f64>> {
+    visible: impl Fn(Vector3<f64>) -> [f64; 4] + Sync,
+) -> Vec<Option<[f64; 4]>> {
     let pixels: Vec<_> = sampled().collect();
     let threads = std::thread::available_parallelism().map_or(1, |n| n.get());
     let chunk = pixels.len().div_ceil(threads);
@@ -122,10 +131,10 @@ fn expected(
 }
 
 /// Compares `eclipsed` with `clear` (the same pixels without the shadow) against
-/// `expected` at every [`sampled`] pixel that has a value. Pixels at the body's edge,
+/// `expected` at every [`sampled`] pixel that has a value, channel by channel. Pixels at the body's edge,
 /// where `clear` is under 1 % of its largest value, are skipped: a rendered edge and
 /// an f64 ray need not agree on which pixels they cover.
-fn tally(eclipsed: &[[f32; 4]], clear: &[[f32; 4]], expected: &[Option<f64>]) -> Tally {
+fn tally(eclipsed: &[[f32; 4]], clear: &[[f32; 4]], expected: &[Option<[f64; 4]>]) -> Tally {
     let peak = clear.iter().map(|p| p[1]).fold(0.0, f32::max);
     let mut t = Tally {
         dark: 0,
@@ -142,20 +151,28 @@ fn tally(eclipsed: &[[f32; 4]], clear: &[[f32; 4]], expected: &[Option<f64>]) ->
         }
         for k in 0..3 {
             let ratio = got[k] as f64 / base[k] as f64;
+            let want = want[k];
+            // Fully in or out of the shadow: to within f32 rounding of the sum over
+            // the disk, which near an edge of the shadow can land on the other side of
+            // it (measured: 1.4e-7 of the unshadowed value).
             if want == 0.0 {
-                assert_eq!(got[k], 0.0, "pixel ({u}, {v}) is in the umbra");
+                assert!(ratio < 1e-6, "pixel ({u}, {v}) is in the umbra: {ratio}");
             } else if want == 1.0 {
-                assert_eq!(got[k], base[k], "pixel ({u}, {v}) is outside the shadow");
+                assert!(
+                    (ratio - 1.0).abs() < 1e-6,
+                    "pixel ({u}, {v}) is in sunlight: {ratio}"
+                );
             } else {
                 t.worst = t.worst.max((ratio - want).abs());
             }
         }
-        match want {
+        // Counted by the green channel.
+        match want[1] {
             0.0 => t.dark += 1,
             1.0 => t.lit += 1,
-            _ => {
+            green => {
                 t.partial += 1;
-                t.least = t.least.min(want);
+                t.least = t.least.min(green);
             }
         }
     }
@@ -179,7 +196,7 @@ fn earth_shadows_the_moon() {
         let camera = moon_camera(moon);
         let want = expected(
             |u, v| hit(&camera, u, v, moon, MOON_RADIUS, 1.0),
-            |p| sun_visible_fraction(p, SUN, &earth),
+            |p| sun_visible_light(p, SUN, &earth, limb()),
         );
         (label, moon, camera, want)
     });
@@ -240,7 +257,7 @@ fn moon_shadows_earth() {
                     flattening,
                 )
             },
-            |p| sun_visible_fraction(p, SUN, &[(moon, MOON_RADIUS)]),
+            |p| sun_visible_light(p, SUN, &[(moon, MOON_RADIUS)], limb()),
         );
         (label, moon, want)
     });
@@ -259,8 +276,9 @@ fn moon_shadows_earth() {
             let t = tally(&eclipsed, &clear, want);
             match *label {
                 "total" => assert!(t.dark > 10 && t.partial > 10_000, "{name}: {t:?}"),
-                // In the antumbra the Moon sits inside the Sun's disk and 1 - (b/a)^2,
-                // about 0.024, of it shows; a partial overlap never gets that low.
+                // In the antumbra the Moon sits inside the Sun's disk and covers its
+                // bright center: 1 % of the light gets past (2.4 % of the disk's area),
+                // lower than a partial overlap reaches.
                 _ => assert!(
                     t.dark == 0 && t.partial > 10_000 && t.least < 0.03,
                     "{name}: {t:?}"

crates/resolved/tests/mesh.rs

@@ -118,7 +118,7 @@ fn render(ctx: &GpuContext, name: &str, scene: &Scene, camera: &Pose, mesh: &Mes
                     &k,
                     config.near_clip,
                     scene.sun,
-                    visible,
+                    [visible; 4],
                     None,
                     None,
                 ),

crates/resolved/tests/renderer.rs

@@ -3,12 +3,12 @@
 mod common;
 
 use std::io::Cursor;
+use std::sync::Arc;
 
+use approx::assert_relative_eq;
 use resolved::gpu::AdapterChoice;
 use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3};
 use resolved::radiometry::{Radiometry, AU};
-use std::sync::Arc;
-
 use resolved::{
     CameraConfig, EarthTexture, Fidelity, LightingTier, LoadOptions, Material, Mesh, ModelId, Pose,
     Renderer, RendererOptions, Scene, Target,
@@ -274,3 +274,47 @@ fn profiling_times_every_pass_at_l1_with_several_targets() {
         }
     }
 }
+
+#[test]
+fn the_sun_is_uniform_at_l0_and_limb_darkened_at_l1() {
+    // Looking straight at the Sun's center: L0 draws the disk's mean radiance, L1 the
+    // brighter center of a limb-darkened disk.
+    let _devices = common::device_lock();
+    let radiometry = Radiometry::new(&CameraConfig::example());
+    let limb = resolved::sun::LimbDarkening::new(&radiometry);
+    let camera = Pose::looking_at(TARGET, SUN, Vector3::z()).unwrap();
+    for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
+        for lighting in [LightingTier::L0, LightingTier::L1] {
+            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()
+            };
+            let mut renderer = Renderer::new(CameraConfig::example(), options).unwrap();
+            let frame = renderer
+                .submit(&scene(vec![]), &camera)
+                .unwrap()
+                .wait()
+                .unwrap();
+            let center = frame.radiance.unwrap()[540 * 1920 + 960];
+            let profile = limb.profile(1.0);
+            for k in 0..3 {
+                let mean = radiometry.sun_radiance(k) as f32;
+                match lighting {
+                    LightingTier::L0 => assert_eq!(center[k], mean),
+                    _ => assert_relative_eq!(
+                        center[k] as f64,
+                        mean as f64 * profile[k],
+                        max_relative = 1e-4
+                    ),
+                }
+            }
+        }
+    }
+}