sunlight through atmosphere
This commit is contained in:
connor committed
1 parent
8c0a1439d7
commit
8e57375059
18 files changed
+1921
-134
No files matched your search
@@ -189,6 +189,9 @@ fn render(args: &RenderArgs) -> Result<(), String> {
|
||||
let file = SceneFile::read(path)?;
|
||||
let dir = path.parent().unwrap_or(Path::new("."));
|
||||
let mut camera = file.camera(dir)?;
|
||||
if let (LightingTier::L0, Some(exposure)) = (args.lighting, file.camera.exposure_l0) {
|
||||
camera.exposure = exposure;
|
||||
}
|
||||
if args.noiseless {
|
||||
camera.detector = camera.detector.noiseless();
|
||||
}
|
||||
|
||||
@@ -63,6 +63,10 @@ pub struct CameraFile {
|
||||
pub principal_point: Option<[f64; 2]>,
|
||||
pub f_number: Option<f64>,
|
||||
pub exposure: Option<f64>,
|
||||
/// The exposure at lighting L0 in place of `exposure`, for a scene whose bodies L0
|
||||
/// draws far brighter than L1 does (a lunar eclipse, which L0 draws without Earth's
|
||||
/// shadow).
|
||||
pub exposure_l0: Option<f64>,
|
||||
pub transmission: Option<f64>,
|
||||
pub gain: Option<f64>,
|
||||
pub bias: Option<f64>,
|
||||
@@ -540,7 +544,7 @@ mod tests {
|
||||
#[test]
|
||||
fn camera_overrides_apply() {
|
||||
let file: SceneFile = toml::from_str(&format!(
|
||||
"{MINIMAL}\n[camera]\nwidth = 640\nexposure = 2e-3\nadc_bits = 10\nread_noise = 2.5\n\
|
||||
"{MINIMAL}\n[camera]\nwidth = 640\nexposure = 2e-3\nexposure_l0 = 5e-4\nadc_bits = 10\nread_noise = 2.5\n\
|
||||
dark_current = 0\nshot_noise = false\nprnu = 0.02\nhot_pixels = [1e-3, 20]\n\
|
||||
pattern_seed = 9\nblooming = {{ column = 0.25 }}\n"
|
||||
))
|
||||
@@ -548,6 +552,7 @@ mod tests {
|
||||
let camera = file.camera(Path::new(".")).unwrap();
|
||||
assert_eq!((camera.width, camera.height), (640, 1080));
|
||||
assert_eq!(camera.exposure, 2e-3);
|
||||
assert_eq!(file.camera.exposure_l0, Some(5e-4));
|
||||
assert_eq!(camera.detector.adc_bits, 10);
|
||||
let d = camera.detector;
|
||||
assert_eq!(
|
||||
|
||||
@@ -185,21 +185,24 @@ fn pretty_previews_lift_the_darks_and_leave_space_black() {
|
||||
|
||||
#[cfg(feature = "astrokit")]
|
||||
#[test]
|
||||
fn lunar_eclipse_darkens_the_moon() {
|
||||
// Mid-totality: at L1 the whole Moon is in Earth's umbra and reads the bias. At L0
|
||||
// there are no eclipses on the bodies and it is a full Moon.
|
||||
fn lunar_eclipse_darkens_and_reddens_the_moon() {
|
||||
// Mid-totality: at L1 the whole Moon is in Earth's umbra, lit only through Earth's
|
||||
// atmosphere, red, and brighter at its northern edge, farther from the shadow's
|
||||
// axis, than at its southern. At L0 there are no eclipses on the bodies, and the
|
||||
// full Moon at the scene's L0 exposure, 300 times shorter, is bright but short of
|
||||
// saturation.
|
||||
let scene = "examples/scenes/lunar-eclipse.toml";
|
||||
let l0 = render_at(scene, "gpu", "L0");
|
||||
let l1 = render_at(scene, "gpu", "L1");
|
||||
assert!(
|
||||
l0.get_pixel(960, 540).0.iter().all(|&v| v > 150),
|
||||
"{:?}",
|
||||
l0.get_pixel(960, 540)
|
||||
);
|
||||
assert!(
|
||||
l1.pixels().all(|p| p.0 == [64; 3]),
|
||||
"every pixel at the bias"
|
||||
);
|
||||
let full = l0.get_pixel(960, 540).0;
|
||||
assert!(full.iter().all(|&v| v > 1000 && v < 4000), "{full:?}");
|
||||
let [north, center, south] = [260, 540, 820].map(|y| l1.get_pixel(960, y).0);
|
||||
// Measured: 695, 287, 145 DN at the north, 592, 164, 77 at the center, 208, 82, 65
|
||||
// at the south, over a bias of 64.
|
||||
for p in [north, center, south] {
|
||||
assert!(p[0] > p[1] && p[1] >= p[2] && p[2] >= 64, "{p:?}");
|
||||
}
|
||||
assert!(north[0] > center[0] && center[0] > south[0] && south[0] > 100);
|
||||
}
|
||||
|
||||
#[cfg(feature = "astrokit")]
|
||||
|
||||
@@ -29,6 +29,11 @@ use std::sync::{Arc, Mutex};
|
||||
use crate::bodies::EARTH_EQUATORIAL_RADIUS;
|
||||
use crate::radiometry::{Radiometry, AU, SUN_RADIUS};
|
||||
|
||||
mod grazing;
|
||||
pub use grazing::{
|
||||
EarthShadow, Grazing, GRAZING_BOTTOM, GRAZING_SIZE, MAX_BENDING, REFRACTIVITY, SHADOW_SIZE,
|
||||
};
|
||||
|
||||
/// The ground's radius, in meters.
|
||||
pub const BOTTOM_RADIUS: f64 = EARTH_EQUATORIAL_RADIUS;
|
||||
/// Height of the top of the atmosphere above the ground, in meters. Above it the air is
|
||||
@@ -123,6 +128,9 @@ pub struct AtmosphereTables {
|
||||
/// Irradiance from the sky on a horizontal surface, as a fraction of the Sun's,
|
||||
/// [`IRRADIANCE_SIZE`].
|
||||
pub irradiance: Vec<[f32; 4]>,
|
||||
/// Rays that pass through the atmosphere from space to space: their bending and
|
||||
/// transmittance.
|
||||
pub grazing: Grazing,
|
||||
}
|
||||
|
||||
impl AtmosphereTables {
|
||||
@@ -161,6 +169,7 @@ impl AtmosphereTables {
|
||||
rayleigh,
|
||||
mie,
|
||||
irradiance: Vec::new(),
|
||||
grazing: Grazing::new(spectrum),
|
||||
};
|
||||
tables.irradiance = sky_irradiance(&tables);
|
||||
tables
|
||||
@@ -196,9 +205,15 @@ impl AtmosphereTables {
|
||||
/// as the shaders see it: in Earth-fixed axes stretched along the polar axis, where
|
||||
/// Earth's ellipsoid is a sphere of the equatorial radius.
|
||||
pub fn altitude(scene: &crate::Scene, p: nalgebra::Vector3<f64>) -> f64 {
|
||||
stretched(scene, p).norm() - BOTTOM_RADIUS
|
||||
}
|
||||
|
||||
/// Inertial point `p` in Earth-fixed axes stretched along the polar axis, the frame the
|
||||
/// atmosphere is spherical in, in meters.
|
||||
pub fn stretched(scene: &crate::Scene, p: nalgebra::Vector3<f64>) -> nalgebra::Vector3<f64> {
|
||||
let fixed = scene.earth_attitude.inverse() * p;
|
||||
let stretch = BOTTOM_RADIUS / crate::bodies::EARTH_POLAR_RADIUS;
|
||||
nalgebra::Vector3::new(fixed.x, fixed.y, fixed.z * stretch).norm() - BOTTOM_RADIUS
|
||||
nalgebra::Vector3::new(fixed.x, fixed.y, fixed.z * stretch)
|
||||
}
|
||||
|
||||
/// The constants `shaders/atmosphere.wgsl` expects: lengths in kilometers, and the
|
||||
@@ -215,12 +230,16 @@ pub(crate) fn wgsl_constants() -> String {
|
||||
const ATMOSPHERE_SUN_ANGULAR_RADIUS: f32 = {:?};\n\
|
||||
const ATMOSPHERE_TRANSMITTANCE_SIZE = vec2<f32>({tw}.0, {th}.0);\n\
|
||||
const ATMOSPHERE_SCATTERING_SIZE = vec4<f32>({nu}.0, {mu_s}.0, {mu}.0, {r}.0);\n\
|
||||
const ATMOSPHERE_IRRADIANCE_SIZE = vec2<f32>({iw}.0, {ih}.0);\n",
|
||||
const ATMOSPHERE_IRRADIANCE_SIZE = vec2<f32>({iw}.0, {ih}.0);\n\
|
||||
const ATMOSPHERE_GRAZING_BOTTOM: f32 = {:?};\n\
|
||||
const ATMOSPHERE_GRAZING_SIZE: u32 = {GRAZING_SIZE}u;\n\
|
||||
const ATMOSPHERE_SHADOW_SIZE: u32 = {SHADOW_SIZE}u;\n",
|
||||
BOTTOM_RADIUS / 1000.0,
|
||||
TOP_RADIUS / 1000.0,
|
||||
MU_S_MIN,
|
||||
MIE_G,
|
||||
SUN_ANGULAR_RADIUS,
|
||||
GRAZING_BOTTOM / 1000.0,
|
||||
)
|
||||
}
|
||||
|
||||
@@ -307,6 +326,10 @@ impl Spectral {
|
||||
/// Per channel, the sum over wavelengths of `value(wavelength index)` with the
|
||||
/// channels' weights.
|
||||
fn fold(&self, value: impl Fn(usize) -> f64) -> [f32; 4] {
|
||||
self.fold_f64(value).map(|x| x as f32)
|
||||
}
|
||||
|
||||
fn fold_f64(&self, value: impl Fn(usize) -> f64) -> [f64; 4] {
|
||||
let mut out = [0.0f64; 4];
|
||||
for (l, w) in self.weights.iter().enumerate() {
|
||||
let v = value(l);
|
||||
@@ -314,7 +337,7 @@ impl Spectral {
|
||||
out[k] += w[k] * v;
|
||||
}
|
||||
}
|
||||
out.map(|x| x as f32)
|
||||
out
|
||||
}
|
||||
|
||||
/// Optical depth at wavelength index `l` of the column densities `c` (Rayleigh, Mie,
|
||||
|
||||
@@ -0,0 +1,774 @@
|
||||
//! Rays that pass through Earth's atmosphere from space to space: the sunlight that
|
||||
//! reaches a point in space through the limb ([`Grazing::sunlight`]), and how the Sun and
|
||||
//! the Moon look through it (`shaders/bodies.wgsl`, `through_air`).
|
||||
//!
|
||||
//! Air bends light toward the denser air below it, so a ray that passes the limb comes
|
||||
//! out turned toward Earth, by more the lower it passes: 1.25 degrees for a ray that
|
||||
//! skims the ground, against the Sun's diameter of 0.53. Seen through the limb, the Sun
|
||||
//! is raised, and flattened, since its lower edge is raised more than its upper. Rays
|
||||
//! that pass at different heights come out spread over a wider angle than they went in,
|
||||
//! which dims what lies beyond (refractive dilution), and the air along a ray scatters
|
||||
//! and absorbs its light, blue most.
|
||||
//!
|
||||
//! The atmosphere is spherical and its refractivity is in proportion to the air's
|
||||
//! density, so a ray's path depends only on its impact parameter b, the distance of its
|
||||
//! straight line outside the atmosphere from Earth's center. The tables are functions of
|
||||
//! b: the ray's total bending, and its transmittance along the whole curved path. Air's
|
||||
//! refractivity changes by about 2 % across the visible spectrum; it is taken at 550 nm
|
||||
//! for every wavelength, so the Sun's blue and red images at the limb do not separate
|
||||
//! (the green flash).
|
||||
|
||||
use std::f64::consts::PI;
|
||||
|
||||
use nalgebra::Vector3;
|
||||
|
||||
use super::{densities, parallel, Spectral, BOTTOM_RADIUS, RAYLEIGH_SCALE_HEIGHT, TOP_RADIUS};
|
||||
use crate::radiometry::SUN_RADIUS;
|
||||
use crate::sun::{legendre_rule, ring_rule, LimbDarkening};
|
||||
|
||||
/// Air's refractivity, n - 1, at the ground, for 550 nm: dry air at 15 C and
|
||||
/// 101,325 Pa (Edlen 1966). It falls with the air's density, at Rayleigh's scale height.
|
||||
pub const REFRACTIVITY: f64 = 2.778e-4;
|
||||
|
||||
/// The impact parameter of the lowest ray that clears the ground, in meters: its lowest
|
||||
/// point is at the ground, where n r = b.
|
||||
pub const GRAZING_BOTTOM: f64 = (1.0 + REFRACTIVITY) * BOTTOM_RADIUS;
|
||||
|
||||
/// More than any ray's bending, in radians: the lowest ray's is 0.0218.
|
||||
pub const MAX_BENDING: f64 = 0.025;
|
||||
|
||||
/// Entries in the tables, evenly spaced in impact parameter from [`GRAZING_BOTTOM`] to
|
||||
/// the top of the atmosphere: 96 m apart.
|
||||
pub const GRAZING_SIZE: usize = 1024;
|
||||
|
||||
/// Steps along each half of a ray, evenly spaced in the square root of the height above
|
||||
/// its lowest point.
|
||||
const STEPS: usize = 2000;
|
||||
|
||||
/// Rows across the Sun's disk in [`Grazing::sunlight`]. Measured against 512 rows: 2e-5
|
||||
/// at most through a sunset seen from 700 km.
|
||||
const ROWS: usize = 64;
|
||||
|
||||
/// The bending and transmittance of rays through the atmosphere, by impact parameter
|
||||
/// ([`GRAZING_SIZE`] entries from [`GRAZING_BOTTOM`] to the top of the atmosphere).
|
||||
#[derive(Debug)]
|
||||
pub struct Grazing {
|
||||
/// Per entry, the ray's total bending toward Earth, in radians.
|
||||
pub bending: Vec<f64>,
|
||||
/// Per entry, per channel: the transmittance along the whole ray, weighted by
|
||||
/// sunlight over each channel's band, as the transmittance table is.
|
||||
pub transmittance: Vec<[f64; 4]>,
|
||||
}
|
||||
|
||||
/// The impact parameter of entry `i`.
|
||||
fn impact(i: usize) -> f64 {
|
||||
GRAZING_BOTTOM + (TOP_RADIUS - GRAZING_BOTTOM) * i as f64 / (GRAZING_SIZE - 1) as f64
|
||||
}
|
||||
|
||||
/// The path of the ray of impact parameter `b` through an atmosphere whose refractivity
|
||||
/// at the ground is `refractivity`: its total bending, and the column densities along
|
||||
/// it (meters of air, aerosols and ozone at ground or peak density, as `Columns`).
|
||||
///
|
||||
/// Along a ray in a spherical atmosphere n r sin(z) = b, where z is the angle from the
|
||||
/// zenith, so dr / ds = sqrt(n^2 r^2 - b^2) / (n r). The bending is
|
||||
/// -2 b times the integral of (dn/dr) / (n sqrt(n^2 r^2 - b^2)) from the lowest point,
|
||||
/// where n r = b, to the top. Both integrands go as 1 / sqrt(r - r0) at the lowest point
|
||||
/// r0; with r = r0 + x^2 they are smooth in x.
|
||||
pub(crate) fn path(b: f64, refractivity: f64) -> (f64, [f64; 3]) {
|
||||
let n = |r: f64| 1.0 + refractivity * densities(r - BOTTOM_RADIUS)[0];
|
||||
let dn = |r: f64| -refractivity * densities(r - BOTTOM_RADIUS)[0] / RAYLEIGH_SCALE_HEIGHT;
|
||||
// n r grows with r: its derivative, n + r dn/dr, is 0.78 at the ground and more
|
||||
// above, so the lowest point is unique (no ducting).
|
||||
let mut r0 = b;
|
||||
for _ in 0..50 {
|
||||
let step = (n(r0) * r0 - b) / (n(r0) + r0 * dn(r0));
|
||||
r0 -= step;
|
||||
if step.abs() < 1e-9 {
|
||||
break;
|
||||
}
|
||||
}
|
||||
let dx = (TOP_RADIUS - r0).max(0.0).sqrt() / STEPS as f64;
|
||||
let mut bending = 0.0;
|
||||
let mut column = [0.0; 3];
|
||||
for step in 0..=STEPS {
|
||||
let x = step as f64 * dx;
|
||||
let r = r0 + x * x;
|
||||
let nr = n(r) * r;
|
||||
// dr / sqrt(n^2 r^2 - b^2) per dx, and its limit at the lowest point, where
|
||||
// n^2 r^2 - b^2 = 2 b (n + r dn/dr) x^2 to first order.
|
||||
let slope = if step == 0 {
|
||||
2.0 / (2.0 * b * (n(r0) + r0 * dn(r0))).sqrt()
|
||||
} else {
|
||||
2.0 * x / ((nr - b) * (nr + b)).sqrt()
|
||||
};
|
||||
// The trapezoid rule, over both halves of the ray.
|
||||
let weight = if step == 0 || step == STEPS { 0.5 } else { 1.0 } * dx * 2.0;
|
||||
bending -= b * dn(r) / n(r) * slope * weight;
|
||||
let rho = densities(r - BOTTOM_RADIUS);
|
||||
for s in 0..3 {
|
||||
column[s] += rho[s] * nr * slope * weight;
|
||||
}
|
||||
}
|
||||
(bending, column)
|
||||
}
|
||||
|
||||
impl Grazing {
|
||||
pub(super) fn new(spectrum: &Spectral) -> Self {
|
||||
let paths = parallel(GRAZING_SIZE, |i| path(impact(i), REFRACTIVITY));
|
||||
Self {
|
||||
bending: paths.iter().map(|(bending, _)| *bending).collect(),
|
||||
transmittance: paths
|
||||
.iter()
|
||||
.map(|(_, c)| spectrum.fold_f64(|l| spectrum.transmittance(l, c)))
|
||||
.collect(),
|
||||
}
|
||||
}
|
||||
|
||||
/// The ray of impact parameter `b`, in meters: its bending and its transmittance per
|
||||
/// channel, linear between entries; `None` if it meets the ground. A ray above the
|
||||
/// atmosphere passes straight and whole.
|
||||
pub fn ray(&self, b: f64) -> Option<(f64, [f64; 4])> {
|
||||
if b < GRAZING_BOTTOM {
|
||||
return None;
|
||||
}
|
||||
if b >= TOP_RADIUS {
|
||||
return Some((0.0, [1.0; 4]));
|
||||
}
|
||||
let x = (b - GRAZING_BOTTOM) / (TOP_RADIUS - GRAZING_BOTTOM) * (GRAZING_SIZE - 1) as f64;
|
||||
let i = (x as usize).min(GRAZING_SIZE - 2);
|
||||
let f = x - i as f64;
|
||||
let (t0, t1) = (self.transmittance[i], self.transmittance[i + 1]);
|
||||
Some((
|
||||
self.bending[i] + f * (self.bending[i + 1] - self.bending[i]),
|
||||
std::array::from_fn(|k| t0[k] + f * (t1[k] - t0[k])),
|
||||
))
|
||||
}
|
||||
|
||||
/// Per channel, the fraction of the Sun's light that reaches `point` past Earth and
|
||||
/// through its atmosphere, with the brightness spread over the Sun's disk as `limb`
|
||||
/// gives. `point` and `sun` are relative to Earth's center, in Earth-fixed axes
|
||||
/// stretched along the polar axis ([`super::stretched`]), in meters; `point` must be
|
||||
/// above the atmosphere.
|
||||
///
|
||||
/// A ray seen from `point` at angle `theta` from Earth's center passes it at impact
|
||||
/// parameter `|point| sin(theta)` and arrives from `theta` less its bending, keeping
|
||||
/// its azimuth around the line to Earth's center. Radiance is unchanged along a ray
|
||||
/// but for the transmittance, so the light is the integral of the Sun's radiance
|
||||
/// where each ray comes from, times its transmittance, over the solid angle the
|
||||
/// rays arrive in. That is taken in rows across the Sun's image at constant `theta`,
|
||||
/// each an arc of the circle around Earth's center.
|
||||
pub fn sunlight(
|
||||
&self,
|
||||
point: Vector3<f64>,
|
||||
sun: Vector3<f64>,
|
||||
limb: &LimbDarkening,
|
||||
) -> [f64; 4] {
|
||||
static RULE: std::sync::OnceLock<Vec<(f64, f64)>> = std::sync::OnceLock::new();
|
||||
self.sunlight_by(point, sun, limb, RULE.get_or_init(|| legendre_rule(ROWS)))
|
||||
}
|
||||
|
||||
/// [`Self::sunlight`] with `rule`'s rows across the Sun's image.
|
||||
fn sunlight_by(
|
||||
&self,
|
||||
point: Vector3<f64>,
|
||||
sun: Vector3<f64>,
|
||||
limb: &LimbDarkening,
|
||||
rule: &[(f64, f64)],
|
||||
) -> [f64; 4] {
|
||||
let distance = point.norm();
|
||||
let to_sun = sun - point;
|
||||
let a = (SUN_RADIUS / to_sun.norm()).min(1.0).asin();
|
||||
let down = -point / distance;
|
||||
// The angle from Earth's center to the Sun's.
|
||||
let c = to_sun.cross(&down).norm().atan2(to_sun.dot(&down));
|
||||
let top = (TOP_RADIUS / distance).asin();
|
||||
if c - a >= top {
|
||||
return [1.0; 4];
|
||||
}
|
||||
let ground = (GRAZING_BOTTOM / distance).asin();
|
||||
// Where a ray seen at `theta`, at or above `ground`, comes from: `theta` less its
|
||||
// bending. It grows with `theta`. (At `ground` itself, the impact parameter may
|
||||
// round to below the lowest ray's.)
|
||||
let source = |theta: f64| {
|
||||
theta
|
||||
- self
|
||||
.ray((distance * theta.sin()).max(GRAZING_BOTTOM))
|
||||
.map_or(0.0, |(bend, _)| bend)
|
||||
};
|
||||
// Far enough from Earth, rays that pass low are bent past its center, and come
|
||||
// from the other side of it: a ray seen at azimuth phi from `source` below zero
|
||||
// comes from -`source` at phi + pi. That is how sunlight reaches the axis of
|
||||
// Earth's shadow beyond 290,000 km, where the lowest rays cross it.
|
||||
let lowest = source(ground);
|
||||
// Where rays may come from, as `source`: the circles around Earth's center that
|
||||
// cross the Sun's disk, on this side and past the center.
|
||||
let near = (c - a).max(0.0);
|
||||
let sides = [
|
||||
(near.max(lowest), c + a),
|
||||
(lowest.max(-c - a), (-near).min(0.0)),
|
||||
];
|
||||
if sides.iter().all(|(lo, hi)| lo >= hi) {
|
||||
return [0.0; 4];
|
||||
}
|
||||
// The angle a ray must be seen at to come from `theta`, by bisection.
|
||||
let seen = |theta: f64| {
|
||||
if theta >= top {
|
||||
return theta;
|
||||
}
|
||||
let (mut lo, mut hi) = (ground, top);
|
||||
for _ in 0..60 {
|
||||
let mid = 0.5 * (lo + hi);
|
||||
if source(mid) < theta {
|
||||
lo = mid;
|
||||
} else {
|
||||
hi = mid;
|
||||
}
|
||||
}
|
||||
0.5 * (lo + hi)
|
||||
};
|
||||
let mut light = [0.0; 4];
|
||||
for (from, to) in sides {
|
||||
if from >= to {
|
||||
continue;
|
||||
}
|
||||
let lo = if from <= lowest { ground } else { seen(from) };
|
||||
let hi = seen(to);
|
||||
// A Gauss-Legendre rule in an angle with the rows at its cosine across the
|
||||
// image, which takes out the square root of the disk's width at its ends.
|
||||
let (mid, half) = (0.5 * (hi + lo), 0.5 * (hi - lo));
|
||||
for &(phi, weight) in rule {
|
||||
let theta = mid - half * phi.cos();
|
||||
let Some((bend, transmittance)) = self.ray(distance * theta.sin()) else {
|
||||
continue;
|
||||
};
|
||||
// A circle past the center is the circle as far on this side.
|
||||
let arc = arc(limb, (theta - bend).abs(), c, a);
|
||||
// The row's solid angle is sin(theta) d(theta) d(azimuth); the whole
|
||||
// disk's light is pi a^2 times its mean radiance.
|
||||
let span = theta.sin() * half * phi.sin() * weight / (PI * a * a);
|
||||
for k in 0..4 {
|
||||
light[k] += arc[k] * transmittance[k] * span;
|
||||
}
|
||||
}
|
||||
}
|
||||
light.map(|x| x.clamp(0.0, 1.0))
|
||||
}
|
||||
}
|
||||
|
||||
/// Entries in each row of an [`EarthShadow`].
|
||||
pub const SHADOW_SIZE: usize = 512;
|
||||
|
||||
/// The sunlight that reaches points around the Moon past Earth and through its
|
||||
/// atmosphere ([`Grazing::sunlight`]), as a table: the GPU lights the eclipsed Moon from
|
||||
/// it.
|
||||
///
|
||||
/// Earth's shadow is symmetric about the line from the Sun through Earth's center (in
|
||||
/// the stretched frame, where Earth is a sphere), so the light at a point depends only
|
||||
/// on its distance from Earth's center and the angle c between the Sun's center and
|
||||
/// Earth's seen from it. The table has two rows, at distances `near` and `far` either
|
||||
/// side of the Moon, of [`SHADOW_SIZE`] entries evenly spaced in c less the angular
|
||||
/// radius of the atmosphere's top seen from that distance, which lines the shadow's
|
||||
/// edges up in both rows; between the rows it is linear in distance. It holds the
|
||||
/// light's natural log, which follows the umbra's falloff, a factor of 14 over 400 km
|
||||
/// at the Moon, more closely than the light itself does. Measured against
|
||||
/// [`Grazing::sunlight`] over the Moon's distances: within 0.5 %.
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct EarthShadow {
|
||||
/// The near row's distance from Earth's center, in meters.
|
||||
pub near: f64,
|
||||
/// The far row's distance from Earth's center, in meters.
|
||||
pub far: f64,
|
||||
/// The first entry's c less the atmosphere's angular radius, in radians.
|
||||
pub start: f64,
|
||||
/// The step from one entry to the next, in radians.
|
||||
pub step: f64,
|
||||
/// The Sun's angular radius the table was made for, in radians.
|
||||
pub sun_radius: f64,
|
||||
/// The near row's entries, then the far row's: per channel, the natural log of the
|
||||
/// fraction of the Sun's light.
|
||||
pub light: Vec<[f64; 4]>,
|
||||
}
|
||||
|
||||
impl EarthShadow {
|
||||
/// The table between `near` and `far` meters from Earth's center, where the Sun's
|
||||
/// angular radius is `sun_radius`, with the brightness spread over the Sun's disk as
|
||||
/// `limb` gives.
|
||||
pub fn new(
|
||||
grazing: &Grazing,
|
||||
limb: &LimbDarkening,
|
||||
near: f64,
|
||||
far: f64,
|
||||
sun_radius: f64,
|
||||
) -> Self {
|
||||
// The light's slope jumps where the Sun's far edge reaches the atmosphere's top,
|
||||
// at -sun_radius, and it is whole from where the near edge clears it, at
|
||||
// +sun_radius; entries fall on both, and reach back to the shadow's axis.
|
||||
let last = (SHADOW_SIZE - 1) as f64;
|
||||
let reach = (TOP_RADIUS / near).asin() + sun_radius;
|
||||
let step = 2.0 * sun_radius / (2.0 * sun_radius * last / reach).floor();
|
||||
let start = sun_radius - step * last;
|
||||
let to_sun = SUN_RADIUS / sun_radius.sin();
|
||||
let light = super::parallel(2 * SHADOW_SIZE, |i| {
|
||||
let distance = if i < SHADOW_SIZE { near } else { far };
|
||||
let c = start + step * (i % SHADOW_SIZE) as f64 + (TOP_RADIUS / distance).asin();
|
||||
// Below zero, c is past the axis: the light there is the light as far on
|
||||
// this side.
|
||||
let c = c.abs();
|
||||
let point = Vector3::new(-distance, 0.0, 0.0);
|
||||
let sun = point + to_sun * Vector3::new(c.cos(), c.sin(), 0.0);
|
||||
grazing
|
||||
.sunlight(point, sun, limb)
|
||||
.map(|x| x.max(1e-30).ln())
|
||||
});
|
||||
Self {
|
||||
near,
|
||||
far,
|
||||
start,
|
||||
step,
|
||||
sun_radius,
|
||||
light,
|
||||
}
|
||||
}
|
||||
|
||||
/// Whether the table holds points from `near` to `far` meters from Earth's center,
|
||||
/// where the Sun's angular radius is `sun_radius`, to 1e-4 of it.
|
||||
pub fn covers(&self, near: f64, far: f64, sun_radius: f64) -> bool {
|
||||
self.near <= near && far <= self.far && (sun_radius / self.sun_radius - 1.0).abs() < 1e-4
|
||||
}
|
||||
|
||||
/// Per channel, the fraction of the Sun's light at `distance` meters from Earth's
|
||||
/// center, where the Sun's center is `c` radians from Earth's, from the table, as
|
||||
/// the GPU reads it (`bending.wgsl`, `earth_shadow`).
|
||||
pub fn light(&self, distance: f64, c: f64) -> [f64; 4] {
|
||||
let last = (SHADOW_SIZE - 1) as f64;
|
||||
let x = (c - (TOP_RADIUS / distance).asin() - self.start) / self.step;
|
||||
if x >= last {
|
||||
return [1.0; 4];
|
||||
}
|
||||
let x = x.max(0.0);
|
||||
let i = (x as usize).min(SHADOW_SIZE - 2);
|
||||
let f = x - i as f64;
|
||||
let w = ((distance - self.near) / (self.far - self.near)).clamp(0.0, 1.0);
|
||||
let row = |r: usize, k: usize| {
|
||||
let (a, b) = (self.light[r + i][k], self.light[r + i + 1][k]);
|
||||
a + f * (b - a)
|
||||
};
|
||||
std::array::from_fn(|k| {
|
||||
let (a, b) = (row(0, k), row(SHADOW_SIZE, k));
|
||||
(a + w * (b - a)).exp()
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
/// Per channel, the integral over azimuth around Earth's center of the Sun's radiance
|
||||
/// (over its mean) along the circle `theta` from Earth's center, where the Sun's
|
||||
/// center is `c` from it and its radius is `a`.
|
||||
///
|
||||
/// Angles on the sky from haversines, which keep their precision for a disk this small:
|
||||
/// a point at azimuth `phi` from the Sun's on the circle is `rho` from its center, with
|
||||
/// hav(rho) = hav(theta - c) + sin(theta) sin(c) hav(phi).
|
||||
fn arc(limb: &LimbDarkening, theta: f64, c: f64, a: f64) -> [f64; 4] {
|
||||
let hav = |x: f64| (0.5 * x).sin().powi(2);
|
||||
let across = (theta.sin() * c.sin()).max(1e-300);
|
||||
let edge = (hav(a) - hav(theta - c)) / across;
|
||||
if edge <= 0.0 {
|
||||
return [0.0; 4];
|
||||
}
|
||||
// The arc's half-width in azimuth (all the way around if the disk covers Earth's
|
||||
// center).
|
||||
let width = if edge >= 1.0 {
|
||||
PI
|
||||
} else {
|
||||
2.0 * edge.sqrt().asin()
|
||||
};
|
||||
let mut sum = [0.0; 4];
|
||||
// phi = width cos(psi), which takes out the square root of the profile at the
|
||||
// limb.
|
||||
for &(psi, weight) in ring_rule() {
|
||||
let phi = width * psi.cos();
|
||||
let rho = 2.0 * (hav(theta - c) + across * hav(phi)).sqrt().min(1.0).asin();
|
||||
let mu = (1.0 - (rho / a).powi(2)).max(0.0).sqrt();
|
||||
let p = limb.profile(mu);
|
||||
for k in 0..4 {
|
||||
sum[k] += p[k] * width * psi.sin() * weight;
|
||||
}
|
||||
}
|
||||
sum
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::camera::CameraConfig;
|
||||
use crate::radiometry::{Radiometry, AU};
|
||||
use approx::assert_relative_eq;
|
||||
|
||||
fn radiometry() -> Radiometry {
|
||||
Radiometry::new(&CameraConfig::example())
|
||||
}
|
||||
|
||||
fn tables() -> Grazing {
|
||||
Grazing::new(&Spectral::new(&radiometry()))
|
||||
}
|
||||
|
||||
/// Column densities along the straight line `b` from Earth's center, by a fine
|
||||
/// direct integral.
|
||||
fn straight(b: f64) -> [f64; 3] {
|
||||
let steps = 200_000;
|
||||
let half = (TOP_RADIUS * TOP_RADIUS - b * b).sqrt();
|
||||
let dx = half / steps as f64;
|
||||
let mut sum = [0.0; 3];
|
||||
for i in 0..steps {
|
||||
let x = (i as f64 + 0.5) * dx;
|
||||
let rho = densities((x * x + b * b).sqrt() - BOTTOM_RADIUS);
|
||||
for s in 0..3 {
|
||||
sum[s] += 2.0 * rho[s] * dx;
|
||||
}
|
||||
}
|
||||
sum
|
||||
}
|
||||
|
||||
/// A ray of impact parameter `b` traced through the atmosphere by its equation of
|
||||
/// motion, d(n t)/ds = grad n with t its unit direction, in steps of `ds` meters:
|
||||
/// its bending and the column of air along it.
|
||||
fn traced(b: f64, ds: f64) -> (f64, f64) {
|
||||
let n = |p: [f64; 2]| 1.0 + REFRACTIVITY * densities((p[0].hypot(p[1])) - BOTTOM_RADIUS)[0];
|
||||
// d/ds of (position, n t).
|
||||
let derivative = |s: [f64; 4]| {
|
||||
let p = [s[0], s[1]];
|
||||
let r = p[0].hypot(p[1]);
|
||||
let index = n(p);
|
||||
let grad = -REFRACTIVITY * densities(r - BOTTOM_RADIUS)[0] / RAYLEIGH_SCALE_HEIGHT;
|
||||
[s[2] / index, s[3] / index, grad * p[0] / r, grad * p[1] / r]
|
||||
};
|
||||
let x0 = -(TOP_RADIUS * TOP_RADIUS - b * b).sqrt();
|
||||
let mut s = [x0, b, n([x0, b]), 0.0];
|
||||
let mut air = 0.0;
|
||||
loop {
|
||||
let k1 = derivative(s);
|
||||
let at = |k: [f64; 4], h: f64| std::array::from_fn(|i| s[i] + h * k[i]);
|
||||
let k2 = derivative(at(k1, ds / 2.0));
|
||||
let k3 = derivative(at(k2, ds / 2.0));
|
||||
let k4 = derivative(at(k3, ds));
|
||||
let r_before = s[0].hypot(s[1]);
|
||||
s = std::array::from_fn(|i| {
|
||||
s[i] + ds / 6.0 * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i])
|
||||
});
|
||||
let r = s[0].hypot(s[1]);
|
||||
air += ds
|
||||
* 0.5
|
||||
* (densities(r_before - BOTTOM_RADIUS)[0] + densities(r - BOTTOM_RADIUS)[0]);
|
||||
if r > TOP_RADIUS && s[0] > 0.0 {
|
||||
break;
|
||||
}
|
||||
}
|
||||
((-s[3]).atan2(s[2]), air)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn straight_rays_match_a_direct_integral() {
|
||||
// Without refraction the path is the straight line.
|
||||
for h in [2e3, 10e3, 30e3, 60e3] {
|
||||
let b = BOTTOM_RADIUS + h;
|
||||
let (bending, column) = path(b, 0.0);
|
||||
assert_eq!(bending, 0.0);
|
||||
let want = straight(b);
|
||||
for s in 0..3 {
|
||||
assert_relative_eq!(column[s], want[s], max_relative = 1e-5);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn bending_matches_a_traced_ray() {
|
||||
// Measured: within 2e-5 in the bending and 3e-5 in the air.
|
||||
for h in [3e3, 15e3, 40e3] {
|
||||
let b = GRAZING_BOTTOM + h;
|
||||
let (bending, column) = path(b, REFRACTIVITY);
|
||||
let (traced_bending, traced_air) = traced(b, 50.0);
|
||||
assert_relative_eq!(bending, traced_bending, max_relative = 1e-4);
|
||||
assert_relative_eq!(column[0], traced_air, max_relative = 1e-4);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_ray_skimming_the_ground_bends_more_than_two_suns() {
|
||||
// Astronomical refraction at the horizon is about 35 arcminutes for an observer
|
||||
// at the ground: half of the bending of a ray from space to space that skims
|
||||
// it, so about 70. The model's 8 km scale height gives 74.8: the real air's
|
||||
// density falls more slowly near the ground, over about 10 km, and bends less.
|
||||
let bending = tables().bending[0];
|
||||
let arcmin = bending.to_degrees() * 60.0;
|
||||
assert!((70.0..80.0).contains(&arcmin), "{arcmin}");
|
||||
assert!(bending < MAX_BENDING);
|
||||
}
|
||||
|
||||
/// A point `altitude` above the ground in the plane of the Sun (along +X), `past`
|
||||
/// radians on in its orbit from where the straight line toward the Sun's center
|
||||
/// grazes the ground (Earth's geometric shadow, taken as a cylinder).
|
||||
fn orbit(altitude: f64, past: f64) -> Vector3<f64> {
|
||||
let r = BOTTOM_RADIUS + altitude;
|
||||
let alpha = PI - (BOTTOM_RADIUS / r).asin() + past;
|
||||
r * Vector3::new(alpha.cos(), alpha.sin(), 0.0)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn whole_above_the_air_and_dark_in_the_umbra() {
|
||||
let grazing = tables();
|
||||
let limb = LimbDarkening::new(&radiometry());
|
||||
let sun = Vector3::new(AU, 0.0, 0.0);
|
||||
let r = BOTTOM_RADIUS + 700e3;
|
||||
// Over the day side, and deep in the night.
|
||||
assert_eq!(
|
||||
grazing.sunlight(Vector3::new(0.0, r, 0.0), sun, &limb),
|
||||
[1.0; 4]
|
||||
);
|
||||
assert_eq!(
|
||||
grazing.sunlight(Vector3::new(-r, 0.0, 0.0), sun, &limb),
|
||||
[0.0; 4]
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn airless_limb_matches_the_eclipse_model() {
|
||||
// No bending and no air: Earth's hard edge, at the lowest ray's height. The
|
||||
// eclipse model draws that edge as a circle on a flat sky, and measured 1.7e-4
|
||||
// brighter than these rows; a 4000 x 4000 sampling of the disk agreed with the
|
||||
// rows to 3.4e-5.
|
||||
let airless = Grazing {
|
||||
bending: vec![0.0; GRAZING_SIZE],
|
||||
transmittance: vec![[1.0; 4]; GRAZING_SIZE],
|
||||
};
|
||||
let limb = LimbDarkening::new(&radiometry());
|
||||
let sun = Vector3::new(AU, 0.0, 0.0);
|
||||
for past in [-2e-3, -1e-3, 0.0, 1e-3, 2e-3] {
|
||||
let p = orbit(700e3, past);
|
||||
let got = airless.sunlight(p, sun, &limb);
|
||||
let want = crate::eclipse::sun_visible_light(
|
||||
p,
|
||||
sun,
|
||||
&[(Vector3::zeros(), GRAZING_BOTTOM)],
|
||||
&limb,
|
||||
);
|
||||
for k in 0..3 {
|
||||
assert_relative_eq!(got[k], want[k], epsilon = 3e-4);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn more_rows_change_nothing() {
|
||||
let grazing = tables();
|
||||
let limb = LimbDarkening::new(&radiometry());
|
||||
let sun = Vector3::new(AU, 0.0, 0.0);
|
||||
let fine = legendre_rule(512);
|
||||
let mut worst = 0.0f64;
|
||||
for i in 0..80 {
|
||||
let p = orbit(700e3, -0.04 + 1e-3 * i as f64);
|
||||
let a = grazing.sunlight(p, sun, &limb);
|
||||
let b = grazing.sunlight_by(p, sun, &limb, &fine);
|
||||
for k in 0..3 {
|
||||
worst = worst.max((a[k] - b[k]).abs());
|
||||
}
|
||||
}
|
||||
// Measured: 2e-5, in the middle of the sunset.
|
||||
assert!(worst < 3e-5, "{worst}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_point_like_sun_is_dimmed_by_the_air_and_by_dilution() {
|
||||
// A Sun 10^4 times farther, 9e-7 rad across, seen through the limb at impact
|
||||
// parameter b: the transmittance there times the refractive dilution, the solid
|
||||
// angle the Sun covers over the solid angle its image covers.
|
||||
let grazing = tables();
|
||||
let limb = LimbDarkening::uniform();
|
||||
let p = Vector3::new(0.0, BOTTOM_RADIUS + 700e3, 0.0);
|
||||
let distance = p.norm();
|
||||
for h in [8e3, 15e3, 25e3] {
|
||||
let b = GRAZING_BOTTOM + h;
|
||||
let seen = (b / distance).asin();
|
||||
// Seen at `seen` below the horizontal (Earth's center is along -Y).
|
||||
let bend = |b: f64| path(b, REFRACTIVITY).0;
|
||||
let from = seen - bend(b);
|
||||
// d(from)/d(seen), from the bending's derivative in b, directly.
|
||||
let db = 1.0;
|
||||
let slope = 1.0 - (bend(b + db) - bend(b - db)) / (2.0 * db) * distance * seen.cos();
|
||||
let dilution = seen.sin() / from.sin() / slope;
|
||||
let direction = Vector3::new(from.sin(), -from.cos(), 0.0);
|
||||
let sun = p + 1e4 * AU * direction;
|
||||
let got = grazing.sunlight(p, sun, &limb);
|
||||
let (_, transmittance) = grazing.ray(b).unwrap();
|
||||
for k in 0..3 {
|
||||
let want = transmittance[k] * dilution;
|
||||
assert_relative_eq!(got[k], want, max_relative = 1e-2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Per channel, the sunlight at `point` by sampling the sky around Earth, each
|
||||
/// sample's ray followed back to where it came from: on an `n` x `n` grid in the
|
||||
/// angle from Earth's center and the azimuth around it, over the rays through the
|
||||
/// atmosphere and again over the straight rays above it, out to the Sun's far edge.
|
||||
fn sampled(
|
||||
grazing: &Grazing,
|
||||
point: Vector3<f64>,
|
||||
sun: Vector3<f64>,
|
||||
limb: &LimbDarkening,
|
||||
n: usize,
|
||||
) -> [f64; 4] {
|
||||
let distance = point.norm();
|
||||
let down = -point / distance;
|
||||
let e1 = down.cross(&Vector3::z()).normalize();
|
||||
let e2 = down.cross(&e1);
|
||||
let to_sun = (sun - point).normalize();
|
||||
let a = (SUN_RADIUS / (sun - point).norm()).asin();
|
||||
let c = to_sun.cross(&down).norm().atan2(to_sun.dot(&down));
|
||||
let (ground, top) = (
|
||||
(GRAZING_BOTTOM / distance).asin(),
|
||||
(TOP_RADIUS / distance).asin(),
|
||||
);
|
||||
let mut light = [0.0; 4];
|
||||
for (lo, hi) in [(ground, top), (top, top.max(c + a))] {
|
||||
let (dt, dp) = ((hi - lo) / n as f64, 2.0 * PI / n as f64);
|
||||
for i in 0..n {
|
||||
let theta = lo + (i as f64 + 0.5) * dt;
|
||||
let (bend, transmittance) = grazing.ray(distance * theta.sin()).unwrap();
|
||||
let from = theta - bend;
|
||||
for j in 0..n {
|
||||
let phi = (j as f64 + 0.5) * dp;
|
||||
let d = from.cos() * down + from.sin() * (phi.cos() * e1 + phi.sin() * e2);
|
||||
let rho = d.cross(&to_sun).norm().atan2(d.dot(&to_sun));
|
||||
if rho < a {
|
||||
let p = limb.profile((1.0 - (rho / a).powi(2)).sqrt());
|
||||
for k in 0..4 {
|
||||
light[k] +=
|
||||
p[k] * transmittance[k] * theta.sin() * dt * dp / (PI * a * a);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
light
|
||||
}
|
||||
|
||||
/// A point `distance` from Earth's center, `off` meters from the axis of its shadow
|
||||
/// (with the Sun along +X).
|
||||
fn behind(distance: f64, off: f64) -> Vector3<f64> {
|
||||
Vector3::new(-(distance * distance - off * off).sqrt(), off, 0.0)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_shadow_at_the_moon_matches_sampling_the_sky() {
|
||||
// At the Moon's distance, across the umbra, where the lowest rays are bent past
|
||||
// Earth's center, and the penumbra. Measured: within 5e-4 of 1000 x 1000
|
||||
// samples; 3000 x 3000 agree to 2e-4.
|
||||
let grazing = tables();
|
||||
let limb = LimbDarkening::new(&radiometry());
|
||||
let sun = Vector3::new(AU, 0.0, 0.0);
|
||||
for off in [1000e3, 3000e3, 4600e3, 5000e3, 6000e3] {
|
||||
let p = behind(3.844e8, off);
|
||||
let got = grazing.sunlight(p, sun, &limb);
|
||||
let want = sampled(&grazing, p, sun, &limb, 1000);
|
||||
for k in 0..3 {
|
||||
assert_relative_eq!(got[k], want[k], max_relative = 1e-3);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_eclipsed_moon_is_dim_and_red() {
|
||||
// Across the umbra at the Moon's distance. Observed in June 2011 (Ugolnikov et
|
||||
// al., arXiv:1106.6178, with a model of clear air like this one), from the
|
||||
// center out to 0.6 degrees (4000 km): 1e-4 to 1e-3 of the full Moon at 677 nm,
|
||||
// falling to 1e-6 at 503 nm. Measured here, red: 5.9e-5 at the center, 3.1e-4
|
||||
// at 3000 km, 6.5e-4 at 4000 km; green 10 to 2 times less, blue 170 to 3 times
|
||||
// less. The center is darker than the edge: the Moon is too far beyond the
|
||||
// lowest rays' focus, 290,000 km behind Earth, to see them brighten its axis.
|
||||
let grazing = tables();
|
||||
let limb = LimbDarkening::new(&radiometry());
|
||||
let sun = Vector3::new(AU, 0.0, 0.0);
|
||||
let mut last = 0.0;
|
||||
for off in [0.0, 1000e3, 2000e3, 3000e3, 4000e3] {
|
||||
let light = grazing.sunlight(behind(3.844e8, off), sun, &limb);
|
||||
assert!((3e-5..1e-3).contains(&light[0]), "{off}: {light:?}");
|
||||
assert!(
|
||||
light[0] > 2.0 * light[1] && light[1] > light[2],
|
||||
"{light:?}"
|
||||
);
|
||||
assert!(light[0] > last, "{off}: {light:?}");
|
||||
last = light[0];
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_shadow_table_follows_the_sunlight() {
|
||||
// Over the Moon's depth and across the umbra and the penumbra at 380,000 km.
|
||||
// Measured: within 4.9e-3, worst where the Sun's far edge nears the atmosphere's
|
||||
// top and the light rises steeply; 256 entries a row were within 1.7e-2.
|
||||
let grazing = tables();
|
||||
let limb = LimbDarkening::new(&radiometry());
|
||||
let (center, moon) = (3.8e8, 1.7374e6);
|
||||
let sun_radius = (SUN_RADIUS / AU).asin();
|
||||
let table = EarthShadow::new(
|
||||
&grazing,
|
||||
&limb,
|
||||
center - 1.5 * moon,
|
||||
center + 1.5 * moon,
|
||||
sun_radius,
|
||||
);
|
||||
assert!(table.covers(center - moon, center + moon, sun_radius * (1.0 + 5e-5)));
|
||||
assert!(!table.covers(center - 2.0 * moon, center + moon, sun_radius));
|
||||
let mut worst = 0.0f64;
|
||||
for i in 0..=20 {
|
||||
let distance = center - moon + 0.1 * moon * i as f64;
|
||||
let top = (TOP_RADIUS / distance).asin();
|
||||
for j in 0..=400 {
|
||||
let c = (top + 1.1 * sun_radius) * j as f64 / 400.0;
|
||||
let point = Vector3::new(-distance, 0.0, 0.0);
|
||||
let sun = point + AU * Vector3::new(c.cos(), c.sin(), 0.0);
|
||||
let want = grazing.sunlight(point, sun, &limb);
|
||||
let got = table.light(distance, c);
|
||||
for k in 0..3 {
|
||||
worst = worst.max((got[k] / want[k] - 1.0).abs());
|
||||
}
|
||||
}
|
||||
}
|
||||
assert!(worst < 6e-3, "{worst}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sunset_dims_blue_first_and_red_last() {
|
||||
// Through a sunset seen from 700 km, 0.02 rad (19 s) of orbit to the geometric
|
||||
// shadow's edge. Each channel dims steadily. Channels 2, 1, 0 are blue, green,
|
||||
// red. Blue dims first, where air scatters it along the ray. Through the ozone
|
||||
// layer, 20 to 30 km up, ozone's Chappuis band absorbs green and red more than
|
||||
// air scatters blue, and blue leads for a while (the turquoise edge of Earth's
|
||||
// shadow in lunar eclipses). Below it red leads again. Refraction still brings
|
||||
// in a twentieth of the red light at the geometric edge.
|
||||
let grazing = tables();
|
||||
let limb = LimbDarkening::new(&radiometry());
|
||||
let sun = Vector3::new(AU, 0.0, 0.0);
|
||||
let mut last = [1.0; 4];
|
||||
let mut blue_led = false;
|
||||
for i in 0..=200 {
|
||||
let past = -0.02 + 1e-4 * i as f64;
|
||||
let light = grazing.sunlight(orbit(700e3, past), sun, &limb);
|
||||
for k in 0..3 {
|
||||
assert!(
|
||||
light[k] <= last[k] + 1e-12,
|
||||
"{past}: {light:?} after {last:?}"
|
||||
);
|
||||
}
|
||||
if light[0] > 0.9 {
|
||||
assert!(
|
||||
light[2] < light[1] && light[1] < light[0],
|
||||
"{past}: {light:?}"
|
||||
);
|
||||
}
|
||||
blue_led |= light[2] > light[0] + 0.03;
|
||||
last = light;
|
||||
}
|
||||
assert!(blue_led);
|
||||
assert!(last[0] > 0.04 && last[2] < 0.5 * last[0], "{last:?}");
|
||||
}
|
||||
}
|
||||
@@ -8,8 +8,11 @@
|
||||
|
||||
use nalgebra::Vector3;
|
||||
|
||||
use crate::atmosphere::{self, AtmosphereTables};
|
||||
use crate::bodies::MOON_RADIUS;
|
||||
use crate::radiometry::SUN_RADIUS;
|
||||
use crate::sun::LimbDarkening;
|
||||
use crate::Scene;
|
||||
|
||||
/// 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.
|
||||
@@ -49,6 +52,57 @@ pub fn sun_visible_light(
|
||||
covered.map(|x| (1.0 - x).clamp(0.0, 1.0))
|
||||
}
|
||||
|
||||
/// Per channel, the fraction of the Sun's light that reaches inertial `point` past Earth
|
||||
/// and through its atmosphere (bent and dimmed, [`atmosphere::Grazing::sunlight`]), and
|
||||
/// past the Moon, with the solar disk's brightness spread as `limb` gives.
|
||||
///
|
||||
/// A point inside the atmosphere instead takes the transmittance of the straight path
|
||||
/// from it toward the Sun's center, times the part of the disk past Earth's hard edge.
|
||||
pub fn sun_visible_through_air(
|
||||
point: Vector3<f64>,
|
||||
scene: &Scene,
|
||||
tables: &AtmosphereTables,
|
||||
limb: &LimbDarkening,
|
||||
) -> [f64; 4] {
|
||||
let moon = sun_visible_light(point, scene.sun, &[(scene.moon, MOON_RADIUS)], limb);
|
||||
// Where the atmosphere is a sphere.
|
||||
let p = atmosphere::stretched(scene, point);
|
||||
let sun = atmosphere::stretched(scene, scene.sun);
|
||||
let earth = if p.norm() >= atmosphere::TOP_RADIUS {
|
||||
tables.grazing.sunlight(p, sun, limb)
|
||||
} else {
|
||||
let edge = sun_visible_light(
|
||||
p,
|
||||
sun,
|
||||
&[(Vector3::zeros(), atmosphere::BOTTOM_RADIUS)],
|
||||
limb,
|
||||
);
|
||||
let r = p.norm().max(atmosphere::BOTTOM_RADIUS);
|
||||
let air = tables.transmittance(r, p.dot(&(sun - p).normalize()) / p.norm());
|
||||
std::array::from_fn(|k| edge[k] * air[k])
|
||||
};
|
||||
// The Moon and Earth do not overlap on the sky (as `sun_visible_fraction`).
|
||||
std::array::from_fn(|k| (earth[k] + moon[k] - 1.0).clamp(0.0, 1.0))
|
||||
}
|
||||
|
||||
/// When some of the Moon may be in Earth's shadow, lit only through Earth's atmosphere
|
||||
/// or not at all: the distance of the Moon's center from Earth's, and the Sun's
|
||||
/// angular radius there, in the atmosphere's stretched frame
|
||||
/// ([`atmosphere::stretched`]), where [`atmosphere::EarthShadow`] takes them.
|
||||
pub fn moon_near_shadow(scene: &Scene) -> Option<(f64, f64)> {
|
||||
let moon = atmosphere::stretched(scene, scene.moon);
|
||||
let to_sun = atmosphere::stretched(scene, scene.sun) - moon;
|
||||
let distance = moon.norm();
|
||||
let a = (SUN_RADIUS / to_sun.norm()).min(1.0).asin();
|
||||
let c = to_sun.cross(&-moon).norm().atan2(to_sun.dot(&-moon));
|
||||
// Seen from anywhere on the Moon, Earth's center is less than the Moon's angular
|
||||
// radius from Earth (stretched by up to 0.34 %) from where it is seen from the
|
||||
// Moon's center.
|
||||
let near = distance - MOON_RADIUS;
|
||||
let reach = (atmosphere::TOP_RADIUS / near).asin() + a + 1.01 * MOON_RADIUS / near;
|
||||
(c < reach).then_some((distance, a))
|
||||
}
|
||||
|
||||
/// Seen from `point`: the Sun's angular radius, the occulter's, and the angle between
|
||||
/// their centers.
|
||||
fn angles(
|
||||
@@ -240,4 +294,46 @@ mod tests {
|
||||
max_relative = 1e-12
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn through_air_adds_the_atmosphere_to_the_eclipse_model() {
|
||||
let r = crate::radiometry::Radiometry::new(&crate::camera::CameraConfig::example());
|
||||
let tables = AtmosphereTables::new(&r);
|
||||
let limb = LimbDarkening::new(&r);
|
||||
let sun = Vector3::new(AU, 0.0, 0.0);
|
||||
let scene = |moon| Scene {
|
||||
earth_attitude: nalgebra::UnitQuaternion::identity(),
|
||||
sun,
|
||||
moon,
|
||||
moon_attitude: nalgebra::UnitQuaternion::identity(),
|
||||
targets: vec![],
|
||||
velocities: Default::default(),
|
||||
};
|
||||
// Over the day side, with the Moon partly over the Sun: the eclipse model.
|
||||
let day = Vector3::new(0.0, RE + 700e3, 0.0);
|
||||
let moon = day + Vector3::new(3.844e8, MOON_RADIUS, 0.0);
|
||||
let got = sun_visible_through_air(day, &scene(moon), &tables, &limb);
|
||||
let want = sun_visible_light(day, sun, &[(moon, MOON_RADIUS)], &limb);
|
||||
assert!(want[1] > 0.3 && want[1] < 0.9, "{want:?}");
|
||||
for k in 0..3 {
|
||||
assert_relative_eq!(got[k], want[k], max_relative = 1e-12);
|
||||
}
|
||||
// At sunset, with the Moon far from the Sun: the light through the limb.
|
||||
let far = Vector3::new(0.0, 0.0, -3.844e8);
|
||||
let alpha = std::f64::consts::PI - (RE / (RE + 700e3)).asin() - 0.01;
|
||||
let sunset = (RE + 700e3) * Vector3::new(alpha.cos(), alpha.sin(), 0.0);
|
||||
let got = sun_visible_through_air(sunset, &scene(far), &tables, &limb);
|
||||
let want = tables.grazing.sunlight(sunset, sun, &limb);
|
||||
for k in 0..3 {
|
||||
assert_relative_eq!(got[k], want[k], max_relative = 1e-12);
|
||||
}
|
||||
assert!(got[0] > 0.1 && got[0] < 0.9, "{got:?}");
|
||||
// 50 km up with the Sun overhead: the transmittance straight up.
|
||||
let inside = Vector3::new(RE + 50e3, 0.0, 0.0);
|
||||
let got = sun_visible_through_air(inside, &scene(far), &tables, &limb);
|
||||
let want = tables.transmittance(RE + 50e3, 1.0);
|
||||
for k in 0..3 {
|
||||
assert_relative_eq!(got[k], want[k], max_relative = 1e-12);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -12,7 +12,10 @@ use super::atmosphere::{self, AtmosphereTextures};
|
||||
use super::profile::PassTimestamps;
|
||||
use super::targets::{FrameTargets, DEPTH_FORMAT, RADIANCE_FORMAT};
|
||||
use super::GpuContext;
|
||||
use crate::atmosphere::AtmosphereTables;
|
||||
use crate::atmosphere::{
|
||||
AtmosphereTables, EarthShadow, BOTTOM_RADIUS, GRAZING_SIZE, MAX_BENDING, SHADOW_SIZE,
|
||||
TOP_RADIUS,
|
||||
};
|
||||
use crate::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS};
|
||||
use crate::camera::Intrinsics;
|
||||
use crate::earth::{DetailFormat, DetailLayout, EarthDetail, EarthTexture};
|
||||
@@ -202,6 +205,19 @@ pub struct BodiesUniforms {
|
||||
}
|
||||
|
||||
impl BodiesUniforms {
|
||||
/// Whether this frame needs the variant of the pass with Earth's atmosphere in the
|
||||
/// paths of light (`bending.wgsl`): the Moon or the Sun may be seen through it, or
|
||||
/// the Moon lit through it.
|
||||
fn through_air(&self) -> bool {
|
||||
self.earth_shape[1..].iter().any(|&flag| flag > 0.5)
|
||||
}
|
||||
|
||||
/// Lights the Moon this frame from the table of Earth's shadow last set with
|
||||
/// [`BodiesPass::set_earth_shadow`].
|
||||
pub fn light_moon_from_shadow(&mut self) {
|
||||
self.earth_shape[3] = 1.0;
|
||||
}
|
||||
|
||||
/// The uniforms for `scene` seen from `camera`. Geometry is differenced in f64
|
||||
/// before the cast, as in [`Scene::relative_to`]. With `apparent` (camera C1 and
|
||||
/// up), the bodies are drawn where light time and aberration make them appear
|
||||
@@ -245,8 +261,14 @@ impl BodiesUniforms {
|
||||
earth_axes: columns(rel.earth.rotation.to_rotation_matrix()),
|
||||
earth_shape: [
|
||||
(EARTH_EQUATORIAL_RADIUS / EARTH_POLAR_RADIUS) as f32,
|
||||
0.0,
|
||||
0.0,
|
||||
f32::from(u8::from(
|
||||
radiometry.atmosphere
|
||||
&& behind_air(scene.moon - camera.position, MOON_RADIUS, camera),
|
||||
)),
|
||||
f32::from(u8::from(
|
||||
radiometry.atmosphere
|
||||
&& behind_air(scene.sun - camera.position, SUN_RADIUS, camera),
|
||||
)),
|
||||
0.0,
|
||||
],
|
||||
earth_sun: vec4(sun_from(Vector3::zeros()), 0.0),
|
||||
@@ -288,6 +310,20 @@ impl BodiesUniforms {
|
||||
}
|
||||
}
|
||||
|
||||
/// Whether some of a body of `radius` at `offset` from `camera` may be seen through
|
||||
/// Earth's atmosphere: whether its disk comes within the atmosphere's bending of the
|
||||
/// cone of lines of sight that pass through it, with a margin for aberration. Rays come
|
||||
/// from no nearer Earth's center than the ground less the most any ray bends, so a body
|
||||
/// wholly behind that, as the Sun is from a camera in Earth's shadow, cannot be.
|
||||
fn behind_air(offset: Vector3<f64>, radius: f64, camera: &Pose) -> bool {
|
||||
let to_earth = -camera.position;
|
||||
let air = (TOP_RADIUS / to_earth.norm()).min(1.0).asin();
|
||||
let ground = (BOTTOM_RADIUS / to_earth.norm()).min(1.0).asin();
|
||||
let body = (radius / offset.norm()).min(1.0).asin();
|
||||
let apart = offset.cross(&to_earth).norm().atan2(offset.dot(&to_earth));
|
||||
apart < air + body + MAX_BENDING + 1e-3 && apart + body + 1e-3 > ground - MAX_BENDING
|
||||
}
|
||||
|
||||
/// The unit vector toward the Sun from Earth's center in Earth-fixed axes stretched
|
||||
/// along the polar axis, where Earth is a sphere: the frame of `atmosphere.wgsl`.
|
||||
pub(crate) fn sun_stretched(scene: &Scene) -> Vector3<f64> {
|
||||
@@ -304,6 +340,12 @@ pub(crate) fn sun_stretched(scene: &Scene) -> Vector3<f64> {
|
||||
#[derive(Debug)]
|
||||
pub struct BodiesPass {
|
||||
pipeline: wgpu::RenderPipeline,
|
||||
/// The pipeline with the Moon and the Sun seen through Earth's atmosphere
|
||||
/// (`bending.wgsl`), and its bind group, once the atmosphere is set: used in frames
|
||||
/// where either may be.
|
||||
bending: Option<(wgpu::RenderPipeline, wgpu::BindGroup)>,
|
||||
samples: u32,
|
||||
area: bool,
|
||||
uniforms: wgpu::Buffer,
|
||||
sun_disk: wgpu::Buffer,
|
||||
sampler: wgpu::Sampler,
|
||||
@@ -314,6 +356,12 @@ pub struct BodiesPass {
|
||||
atmosphere: AtmosphereTextures,
|
||||
detail: [wgpu::TextureView; DETAIL_GROUPS],
|
||||
detail_sampler: wgpu::Sampler,
|
||||
/// Rays through the atmosphere (`Grazing`), as `bodies.wgsl` declares them: the
|
||||
/// transmittance per entry, then the bending.
|
||||
grazing: wgpu::Buffer,
|
||||
/// Earth's shadow around the Moon ([`EarthShadow`]), for the atmosphere's variant:
|
||||
/// its range, then the light's log per entry, in f32.
|
||||
earth_shadow: wgpu::Buffer,
|
||||
bind_group: wgpu::BindGroup,
|
||||
}
|
||||
|
||||
@@ -340,63 +388,7 @@ impl BodiesPass {
|
||||
/// buffer.
|
||||
pub fn new(ctx: &GpuContext, samples: u32, area: bool) -> Self {
|
||||
let device = &ctx.device;
|
||||
let (radiance, depth) = if samples > 1 {
|
||||
(
|
||||
"texture_multisampled_2d<f32>",
|
||||
"texture_depth_multisampled_2d",
|
||||
)
|
||||
} else {
|
||||
("texture_2d<f32>", "texture_depth_2d")
|
||||
};
|
||||
let positions = sample_positions(samples)
|
||||
.iter()
|
||||
.map(|[x, y]| format!("vec2<f32>({x}, {y})"))
|
||||
.collect::<Vec<_>>()
|
||||
.join(", ");
|
||||
let prelude = format!(
|
||||
"const SAMPLES: u32 = {samples}u;\nconst AREA: bool = {area};\n\
|
||||
alias MeshRadiance = {radiance};\nalias MeshDepth = {depth};\n\
|
||||
const SAMPLE_POSITIONS = array<vec2<f32>, {samples}>({positions});\n"
|
||||
);
|
||||
let module = device.create_shader_module(wgpu::ShaderModuleDescriptor {
|
||||
label: Some("bodies"),
|
||||
source: wgpu::ShaderSource::Wgsl(
|
||||
(prelude
|
||||
+ include_str!("../shaders/common.wgsl")
|
||||
+ &atmosphere::wgsl(5)
|
||||
+ include_str!("../shaders/bodies.wgsl"))
|
||||
.into(),
|
||||
),
|
||||
});
|
||||
let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
|
||||
label: Some("bodies"),
|
||||
layout: None,
|
||||
vertex: wgpu::VertexState {
|
||||
module: &module,
|
||||
entry_point: Some("vs_main"),
|
||||
compilation_options: Default::default(),
|
||||
buffers: &[],
|
||||
},
|
||||
primitive: wgpu::PrimitiveState::default(),
|
||||
// Every pixel it does not discard is written: it has already sorted the
|
||||
// bodies against the meshes.
|
||||
depth_stencil: Some(wgpu::DepthStencilState {
|
||||
format: DEPTH_FORMAT,
|
||||
depth_write_enabled: Some(true),
|
||||
depth_compare: Some(wgpu::CompareFunction::Always),
|
||||
stencil: Default::default(),
|
||||
bias: Default::default(),
|
||||
}),
|
||||
multisample: Default::default(),
|
||||
fragment: Some(wgpu::FragmentState {
|
||||
module: &module,
|
||||
entry_point: Some("fs_main"),
|
||||
compilation_options: Default::default(),
|
||||
targets: &[Some(RADIANCE_FORMAT.into())],
|
||||
}),
|
||||
multiview_mask: None,
|
||||
cache: None,
|
||||
});
|
||||
let pipeline = compile(ctx, samples, area, false);
|
||||
let uniforms = device.create_buffer(&wgpu::BufferDescriptor {
|
||||
label: Some("bodies uniforms"),
|
||||
size: std::mem::size_of::<BodiesUniforms>() as u64,
|
||||
@@ -431,6 +423,19 @@ impl BodiesPass {
|
||||
..Default::default()
|
||||
});
|
||||
let detail = std::array::from_fn(|_| placeholder_detail(ctx));
|
||||
// Zeros until the atmosphere is set; the shader does not read them then.
|
||||
let grazing = device.create_buffer(&wgpu::BufferDescriptor {
|
||||
label: Some("grazing rays"),
|
||||
size: (GRAZING_SIZE * 5 * 4) as u64,
|
||||
usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST,
|
||||
mapped_at_creation: false,
|
||||
});
|
||||
let earth_shadow = device.create_buffer(&wgpu::BufferDescriptor {
|
||||
label: Some("earth shadow"),
|
||||
size: (16 * (1 + 2 * SHADOW_SIZE)) as u64,
|
||||
usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST,
|
||||
mapped_at_creation: false,
|
||||
});
|
||||
let bind_group = bind(
|
||||
ctx,
|
||||
&pipeline,
|
||||
@@ -441,9 +446,13 @@ impl BodiesPass {
|
||||
&moon_view,
|
||||
&atmosphere,
|
||||
(&detail, &detail_sampler),
|
||||
None,
|
||||
);
|
||||
Self {
|
||||
pipeline,
|
||||
bending: None,
|
||||
samples,
|
||||
area,
|
||||
uniforms,
|
||||
sun_disk,
|
||||
sampler,
|
||||
@@ -454,6 +463,8 @@ impl BodiesPass {
|
||||
atmosphere,
|
||||
detail,
|
||||
detail_sampler,
|
||||
grazing,
|
||||
earth_shadow,
|
||||
bind_group,
|
||||
}
|
||||
}
|
||||
@@ -576,18 +587,63 @@ impl BodiesPass {
|
||||
/// uniforms' atmosphere flag turns them on.
|
||||
pub fn set_atmosphere(&mut self, ctx: &GpuContext, tables: &AtmosphereTables) {
|
||||
self.atmosphere = AtmosphereTextures::new(ctx, tables);
|
||||
let grazing = &tables.grazing;
|
||||
let values: Vec<f32> = grazing
|
||||
.transmittance
|
||||
.iter()
|
||||
.flatten()
|
||||
.chain(&grazing.bending)
|
||||
.map(|&x| x as f32)
|
||||
.collect();
|
||||
ctx.queue
|
||||
.write_buffer(&self.grazing, 0, bytemuck::cast_slice(&values));
|
||||
let pipeline = compile(ctx, self.samples, self.area, true);
|
||||
let group = self.bind_to(ctx, &pipeline, true);
|
||||
self.bending = Some((pipeline, group));
|
||||
self.rebind(ctx);
|
||||
}
|
||||
|
||||
/// Uploads the table the Moon is lit from in frames that call
|
||||
/// [`BodiesUniforms::light_moon_from_shadow`].
|
||||
pub fn set_earth_shadow(&self, ctx: &GpuContext, shadow: &EarthShadow) {
|
||||
let range = [
|
||||
shadow.start,
|
||||
shadow.step,
|
||||
shadow.near / 1e3,
|
||||
shadow.far / 1e3,
|
||||
];
|
||||
let values: Vec<f32> = range
|
||||
.iter()
|
||||
.chain(shadow.light.iter().flatten())
|
||||
.map(|&x| x as f32)
|
||||
.collect();
|
||||
ctx.queue
|
||||
.write_buffer(&self.earth_shadow, 0, bytemuck::cast_slice(&values));
|
||||
}
|
||||
|
||||
/// The atmosphere's tables (placeholders until set), for the Earthshine pass.
|
||||
pub fn atmosphere(&self) -> &AtmosphereTextures {
|
||||
&self.atmosphere
|
||||
}
|
||||
|
||||
fn rebind(&mut self, ctx: &GpuContext) {
|
||||
self.bind_group = bind(
|
||||
self.bind_group = self.bind_to(ctx, &self.pipeline, false);
|
||||
if let Some((pipeline, _)) = &self.bending {
|
||||
let group = self.bind_to(ctx, pipeline, true);
|
||||
self.bending.as_mut().expect("set").1 = group;
|
||||
}
|
||||
}
|
||||
|
||||
/// The bind group for `pipeline`, which reads the grazing rays if `bending`.
|
||||
fn bind_to(
|
||||
&self,
|
||||
ctx: &GpuContext,
|
||||
pipeline: &wgpu::RenderPipeline,
|
||||
bending: bool,
|
||||
) -> wgpu::BindGroup {
|
||||
bind(
|
||||
ctx,
|
||||
&self.pipeline,
|
||||
pipeline,
|
||||
&self.uniforms,
|
||||
&self.sun_disk,
|
||||
&self.sampler,
|
||||
@@ -595,7 +651,8 @@ impl BodiesPass {
|
||||
&self.moon,
|
||||
&self.atmosphere,
|
||||
(&self.detail, &self.detail_sampler),
|
||||
);
|
||||
bending.then_some((&self.grazing, &self.earth_shadow)),
|
||||
)
|
||||
}
|
||||
|
||||
/// Draws the Sun with `limb`'s spread of brightness, and eclipses on Earth and the
|
||||
@@ -630,9 +687,13 @@ impl BodiesPass {
|
||||
) {
|
||||
ctx.queue
|
||||
.write_buffer(&self.uniforms, 0, bytemuck::bytes_of(uniforms));
|
||||
let (pipeline, group) = match &self.bending {
|
||||
Some((pipeline, group)) if uniforms.through_air() => (pipeline, group),
|
||||
_ => (&self.pipeline, &self.bind_group),
|
||||
};
|
||||
let meshes = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor {
|
||||
label: Some("bodies meshes"),
|
||||
layout: &self.pipeline.get_bind_group_layout(1),
|
||||
layout: &pipeline.get_bind_group_layout(1),
|
||||
entries: &[
|
||||
wgpu::BindGroupEntry {
|
||||
binding: 0,
|
||||
@@ -672,8 +733,8 @@ impl BodiesPass {
|
||||
timestamp_writes: timestamps.map(PassTimestamps::render),
|
||||
..Default::default()
|
||||
});
|
||||
pass.set_pipeline(&self.pipeline);
|
||||
pass.set_bind_group(0, &self.bind_group, &[]);
|
||||
pass.set_pipeline(pipeline);
|
||||
pass.set_bind_group(0, group, &[]);
|
||||
pass.set_bind_group(1, &meshes, &[]);
|
||||
pass.draw(0..3, 0..1);
|
||||
}
|
||||
@@ -778,6 +839,88 @@ fn upload<'a>(
|
||||
gpu.create_view(&Default::default())
|
||||
}
|
||||
|
||||
/// The bodies pipeline for mesh targets of `samples` samples per pixel, casting grids of
|
||||
/// rays over pixels on bodies' edges if `area`, and with the Moon and the Sun seen
|
||||
/// through Earth's atmosphere if `bending`.
|
||||
fn compile(ctx: &GpuContext, samples: u32, area: bool, bending: bool) -> wgpu::RenderPipeline {
|
||||
let device = &ctx.device;
|
||||
let (radiance, depth) = if samples > 1 {
|
||||
(
|
||||
"texture_multisampled_2d<f32>",
|
||||
"texture_depth_multisampled_2d",
|
||||
)
|
||||
} else {
|
||||
("texture_2d<f32>", "texture_depth_2d")
|
||||
};
|
||||
let positions = sample_positions(samples)
|
||||
.iter()
|
||||
.map(|[x, y]| format!("vec2<f32>({x}, {y})"))
|
||||
.collect::<Vec<_>>()
|
||||
.join(", ");
|
||||
let prelude = format!(
|
||||
"const SAMPLES: u32 = {samples}u;\nconst AREA: bool = {area};\n\
|
||||
alias MeshRadiance = {radiance};\nalias MeshDepth = {depth};\n\
|
||||
const SAMPLE_POSITIONS = array<vec2<f32>, {samples}>({positions});\n"
|
||||
);
|
||||
let module = device.create_shader_module(wgpu::ShaderModuleDescriptor {
|
||||
label: Some("bodies"),
|
||||
source: wgpu::ShaderSource::Wgsl(
|
||||
(prelude
|
||||
+ include_str!("../shaders/common.wgsl")
|
||||
+ &atmosphere::wgsl(5)
|
||||
+ include_str!("../shaders/bodies.wgsl")
|
||||
+ if bending {
|
||||
include_str!("../shaders/bending.wgsl")
|
||||
} else {
|
||||
STRAIGHT
|
||||
})
|
||||
.into(),
|
||||
),
|
||||
});
|
||||
device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
|
||||
label: Some("bodies"),
|
||||
layout: None,
|
||||
vertex: wgpu::VertexState {
|
||||
module: &module,
|
||||
entry_point: Some("vs_main"),
|
||||
compilation_options: Default::default(),
|
||||
buffers: &[],
|
||||
},
|
||||
primitive: wgpu::PrimitiveState::default(),
|
||||
// Every pixel it does not discard is written: it has already sorted the
|
||||
// bodies against the meshes.
|
||||
depth_stencil: Some(wgpu::DepthStencilState {
|
||||
format: DEPTH_FORMAT,
|
||||
depth_write_enabled: Some(true),
|
||||
depth_compare: Some(wgpu::CompareFunction::Always),
|
||||
stencil: Default::default(),
|
||||
bias: Default::default(),
|
||||
}),
|
||||
multisample: Default::default(),
|
||||
fragment: Some(wgpu::FragmentState {
|
||||
module: &module,
|
||||
entry_point: Some("fs_main"),
|
||||
compilation_options: Default::default(),
|
||||
targets: &[Some(RADIANCE_FORMAT.into())],
|
||||
}),
|
||||
multiview_mask: None,
|
||||
cache: None,
|
||||
})
|
||||
}
|
||||
|
||||
/// `through_air` and `earth_shadow` where no body is seen or lit through the
|
||||
/// atmosphere: every ray straight from the camera and whole, and Earth's shadow its hard
|
||||
/// edge's.
|
||||
const STRAIGHT: &str = "
|
||||
fn through_air(d: vec3<f32>, body: u32) -> Bent {
|
||||
return Bent(vec3<f32>(0.0), d, vec4<f32>(1.0));
|
||||
}
|
||||
|
||||
fn earth_shadow(p: vec3<f32>) -> vec4<f32> {
|
||||
return sun_visible(p, u.sun_center, u.earth_center);
|
||||
}
|
||||
";
|
||||
|
||||
/// The bodies pass's bind group.
|
||||
#[allow(clippy::too_many_arguments)]
|
||||
fn bind(
|
||||
@@ -790,6 +933,7 @@ fn bind(
|
||||
moon: &wgpu::TextureView,
|
||||
atmosphere: &AtmosphereTextures,
|
||||
detail: (&[wgpu::TextureView; DETAIL_GROUPS], &wgpu::Sampler),
|
||||
air: Option<(&wgpu::Buffer, &wgpu::Buffer)>,
|
||||
) -> wgpu::BindGroup {
|
||||
let entries = [
|
||||
wgpu::BindGroupEntry {
|
||||
@@ -830,6 +974,18 @@ fn bind(
|
||||
binding: 10 + DETAIL_GROUPS as u32,
|
||||
resource: wgpu::BindingResource::Sampler(detail.1),
|
||||
}])
|
||||
.chain(air.into_iter().flat_map(|(grazing, shadow)| {
|
||||
[
|
||||
wgpu::BindGroupEntry {
|
||||
binding: 15,
|
||||
resource: grazing.as_entire_binding(),
|
||||
},
|
||||
wgpu::BindGroupEntry {
|
||||
binding: 16,
|
||||
resource: shadow.as_entire_binding(),
|
||||
},
|
||||
]
|
||||
}))
|
||||
.collect();
|
||||
ctx.device.create_bind_group(&wgpu::BindGroupDescriptor {
|
||||
label: Some("bodies"),
|
||||
@@ -837,3 +993,33 @@ fn bind(
|
||||
entries: &entries,
|
||||
})
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::radiometry::{AU, SUN_RADIUS};
|
||||
|
||||
#[test]
|
||||
fn the_sun_is_behind_air_only_near_the_limb() {
|
||||
// 7000 km from Earth's center, `alpha` around from the Sun's side, in the plane
|
||||
// of the Sun: the Sun at the limb near alpha = 114.5 degrees.
|
||||
let sun = Vector3::new(AU, 0.0, 0.0);
|
||||
let seen = |alpha: f64| {
|
||||
let alpha = alpha.to_radians();
|
||||
let p = 7.0e6 * Vector3::new(alpha.cos(), alpha.sin(), 0.0);
|
||||
let camera = Pose::looking_at(p, Vector3::zeros(), Vector3::z()).unwrap();
|
||||
behind_air(sun - p, SUN_RADIUS, &camera)
|
||||
};
|
||||
assert!(!seen(0.0));
|
||||
assert!(!seen(100.0));
|
||||
assert!(seen(114.5));
|
||||
assert!(!seen(130.0));
|
||||
// Deep in Earth's shadow: hidden by the ground, however the air bends.
|
||||
assert!(!seen(180.0));
|
||||
// From the Moon's distance the lowest rays cross the shadow's axis, so the Sun
|
||||
// may be seen through the air from anywhere in the shadow.
|
||||
let p = Vector3::new(-3.8e8, 0.0, 0.0);
|
||||
let camera = Pose::looking_at(p, Vector3::zeros(), Vector3::z()).unwrap();
|
||||
assert!(behind_air(sun - p, SUN_RADIUS, &camera));
|
||||
}
|
||||
}
|
||||
@@ -15,9 +15,11 @@ use crate::bodies::{Appearance, MOON_RADIUS};
|
||||
use crate::camera::{Blooming, CameraConfig, Intrinsics, Psf};
|
||||
use std::sync::Arc;
|
||||
|
||||
use crate::atmosphere::AtmosphereTables;
|
||||
use crate::atmosphere::{AtmosphereTables, EarthShadow};
|
||||
use crate::earth::{EarthDetail, EarthTexture};
|
||||
use crate::eclipse::{sun_visible_fraction, sun_visible_light};
|
||||
use crate::eclipse::{
|
||||
moon_near_shadow, sun_visible_fraction, sun_visible_light, sun_visible_through_air,
|
||||
};
|
||||
use crate::fidelity::{CameraTier, Fidelity, LightingTier};
|
||||
use crate::fixed_pattern::FixedPattern;
|
||||
use crate::gpu::blur::BlurPass;
|
||||
@@ -174,6 +176,11 @@ pub struct Renderer {
|
||||
profiler: Option<Profiler>,
|
||||
/// How the Sun's brightness is spread across its disk: uniform below lighting L1.
|
||||
limb: LimbDarkening,
|
||||
/// Earth's atmosphere's tables, when targets' sunlight passes through it.
|
||||
atmosphere: Option<Arc<AtmosphereTables>>,
|
||||
/// The table of Earth's shadow the Moon was last lit from, kept while the Moon stays
|
||||
/// within it.
|
||||
earth_shadow: Option<EarthShadow>,
|
||||
}
|
||||
|
||||
impl Renderer {
|
||||
@@ -427,6 +434,8 @@ impl Renderer {
|
||||
radiance_staging,
|
||||
profiler,
|
||||
limb: LimbDarkening::uniform(),
|
||||
atmosphere: None,
|
||||
earth_shadow: None,
|
||||
};
|
||||
if options.fidelity.lighting >= LightingTier::L1 {
|
||||
renderer.limb = LimbDarkening::new(&renderer.radiometry);
|
||||
@@ -440,6 +449,7 @@ impl Renderer {
|
||||
let tables = AtmosphereTables::new(&renderer.radiometry);
|
||||
renderer.bodies.set_atmosphere(&renderer.ctx, &tables);
|
||||
renderer.body_radiometry.atmosphere = true;
|
||||
renderer.atmosphere = Some(tables);
|
||||
}
|
||||
}
|
||||
if let Some(texture) = &options.earth_texture {
|
||||
@@ -592,7 +602,7 @@ impl Renderer {
|
||||
(scene.moon, MOON_RADIUS),
|
||||
];
|
||||
|
||||
let bodies = BodiesUniforms::new(
|
||||
let mut bodies = BodiesUniforms::new(
|
||||
scene,
|
||||
&camera,
|
||||
&intrinsics,
|
||||
@@ -600,6 +610,30 @@ impl Renderer {
|
||||
&self.body_radiometry,
|
||||
self.fidelity.camera >= CameraTier::C1,
|
||||
);
|
||||
// The Moon in Earth's shadow, lit through the atmosphere (which comes with
|
||||
// lighting L1, and with it eclipses on the bodies): from a table made around the
|
||||
// Moon's distance and kept while the Moon stays within it.
|
||||
if let Some(tables) = &self.atmosphere {
|
||||
if let Some((distance, sun_radius)) = moon_near_shadow(scene) {
|
||||
let (near, far) = (distance - MOON_RADIUS, distance + MOON_RADIUS);
|
||||
if !self
|
||||
.earth_shadow
|
||||
.as_ref()
|
||||
.is_some_and(|s| s.covers(near, far, sun_radius))
|
||||
{
|
||||
let shadow = EarthShadow::new(
|
||||
&tables.grazing,
|
||||
&self.limb,
|
||||
distance - 1.5 * MOON_RADIUS,
|
||||
distance + 1.5 * MOON_RADIUS,
|
||||
sun_radius,
|
||||
);
|
||||
self.bodies.set_earth_shadow(&self.ctx, &shadow);
|
||||
self.earth_shadow = Some(shadow);
|
||||
}
|
||||
bodies.light_moon_from_shadow();
|
||||
}
|
||||
}
|
||||
// Shadow maps, from lighting L1: one per target, with the targets that can
|
||||
// shadow it drawn in.
|
||||
let bounded: Vec<Bounded> = scene
|
||||
@@ -645,11 +679,16 @@ impl Renderer {
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(i, t)| {
|
||||
// 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]
|
||||
// From lighting L1, through the limb-darkened disk, per channel, and
|
||||
// through Earth's atmosphere with it.
|
||||
let visible = match &self.atmosphere {
|
||||
Some(tables) => {
|
||||
sun_visible_through_air(t.pose.position, scene, tables, &self.limb)
|
||||
}
|
||||
None if l1 => {
|
||||
sun_visible_light(t.pose.position, scene.sun, &occulters, &self.limb)
|
||||
}
|
||||
None => [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]),
|
||||
|
||||
@@ -0,0 +1,96 @@
|
||||
// The Moon and the Sun seen through Earth's atmosphere, and the Moon lit through it:
|
||||
// `through_air` and `earth_shadow` for bodies.wgsl, in the variant of the bodies pass
|
||||
// used when any of these may be (bodies.rs). Kept out of the other, where the code
|
||||
// alone, never run, cost the RTX 3060 0.3 ms a frame.
|
||||
|
||||
// Rays through Earth's atmosphere from space to space (atmosphere/grazing.rs), by impact
|
||||
// parameter, ATMOSPHERE_GRAZING_SIZE entries evenly spaced from
|
||||
// ATMOSPHERE_GRAZING_BOTTOM to ATMOSPHERE_TOP: the transmittance along each per channel,
|
||||
// and its bending toward Earth in radians.
|
||||
struct Grazing {
|
||||
transmittance: array<vec4<f32>, ATMOSPHERE_GRAZING_SIZE>,
|
||||
bending: array<f32, ATMOSPHERE_GRAZING_SIZE>,
|
||||
}
|
||||
@group(0) @binding(15) var<storage, read> grazing: Grazing;
|
||||
|
||||
// The ray along `d` toward `body` (1 the Moon, 2 the Sun) through Earth's atmosphere,
|
||||
// if it passes through it: bent toward Earth and dimmed by the air along it, or stopped
|
||||
// if it meets the ground. Rays that miss the atmosphere, and all rays toward a body
|
||||
// that cannot be seen through it this frame, are `d` from the camera, whole.
|
||||
fn through_air(d: vec3<f32>, body: u32) -> Bent {
|
||||
var bent = Bent(vec3<f32>(0.0), d, vec4<f32>(1.0));
|
||||
if u.earth_shape[body] < 0.5 {
|
||||
return bent;
|
||||
}
|
||||
let camera = camera_stretched();
|
||||
let v = normalize(transpose(u.earth_axes) * d * vec3<f32>(1.0, 1.0, u.earth_shape.x));
|
||||
let along = dot(camera, v);
|
||||
// The impact parameter: the line's closest approach to Earth's center, in km.
|
||||
let b = length(camera - along * v);
|
||||
if along >= 0.0 || b >= ATMOSPHERE_TOP {
|
||||
return bent;
|
||||
}
|
||||
let last = f32(ATMOSPHERE_GRAZING_SIZE - 1u);
|
||||
let x = (b - ATMOSPHERE_GRAZING_BOTTOM) / (ATMOSPHERE_TOP - ATMOSPHERE_GRAZING_BOTTOM) * last;
|
||||
if x < 0.0 {
|
||||
bent.transmittance = vec4<f32>(0.0);
|
||||
return bent;
|
||||
}
|
||||
let i = min(u32(x), ATMOSPHERE_GRAZING_SIZE - 2u);
|
||||
let f = x - f32(i);
|
||||
bent.transmittance = mix(grazing.transmittance[i], grazing.transmittance[i + 1u], f);
|
||||
let bend = mix(grazing.bending[i], grazing.bending[i + 1u], f);
|
||||
// The path is symmetric about its lowest point, so the line it leaves on is the
|
||||
// line it came in on turned by `bend` about Earth's center, toward Earth.
|
||||
let e = u.earth_center.xyz;
|
||||
let axis = normalize(cross(d, e));
|
||||
bent.dir = rotate(d, axis, bend);
|
||||
bent.origin = e + rotate(-e, axis, bend);
|
||||
return bent;
|
||||
}
|
||||
|
||||
// `x` rotated by `angle` about the unit `axis`, right-handed.
|
||||
fn rotate(x: vec3<f32>, axis: vec3<f32>, angle: f32) -> vec3<f32> {
|
||||
let c = cos(angle);
|
||||
return x * c + cross(axis, x) * sin(angle) + axis * dot(axis, x) * (1.0 - c);
|
||||
}
|
||||
|
||||
// The sunlight that reaches points around the Moon past Earth and through its
|
||||
// atmosphere (atmosphere/grazing.rs, `EarthShadow`): two rows, at distances `range.z`
|
||||
// and `range.w` km from Earth's center, of ATMOSPHERE_SHADOW_SIZE entries, from
|
||||
// `range.x` radians in steps of `range.y`, of the angle between the Sun's center and
|
||||
// Earth's seen from the point, less the angular radius of the atmosphere's top. Per
|
||||
// channel, the natural log of the fraction of the Sun's light.
|
||||
struct EarthShadow {
|
||||
range: vec4<f32>,
|
||||
light: array<vec4<f32>, 2u * ATMOSPHERE_SHADOW_SIZE>,
|
||||
}
|
||||
@group(0) @binding(16) var<storage, read> shadow: EarthShadow;
|
||||
|
||||
// Per channel, the fraction of the Sun's light that reaches camera-frame point `p` on
|
||||
// the Moon past Earth: through the atmosphere from the table when it is set this frame,
|
||||
// else past Earth's hard edge.
|
||||
fn earth_shadow(p: vec3<f32>) -> vec4<f32> {
|
||||
if u.earth_shape.w < 0.5 {
|
||||
return sun_visible(p, u.sun_center, u.earth_center);
|
||||
}
|
||||
// In the atmosphere's frame, where Earth is a sphere.
|
||||
let to_earth = transpose(u.earth_axes) * (u.earth_center.xyz - p);
|
||||
let down = earth_stretched_km(to_earth, u.earth_shape.x);
|
||||
let stretch = vec3<f32>(1.0, 1.0, u.earth_shape.x);
|
||||
let to_sun = transpose(u.earth_axes) * (u.sun_center.xyz - p) * stretch;
|
||||
let distance = length(down);
|
||||
let c = atan2(length(cross(to_sun, down)), dot(to_sun, down));
|
||||
let last = f32(ATMOSPHERE_SHADOW_SIZE - 1u);
|
||||
let x = (c - asin(ATMOSPHERE_TOP / distance) - shadow.range.x) / shadow.range.y;
|
||||
if x >= last {
|
||||
return vec4<f32>(1.0);
|
||||
}
|
||||
let i = min(u32(max(x, 0.0)), ATMOSPHERE_SHADOW_SIZE - 2u);
|
||||
let f = max(x, 0.0) - f32(i);
|
||||
let w = clamp((distance - shadow.range.z) / (shadow.range.w - shadow.range.z), 0.0, 1.0);
|
||||
let far = ATMOSPHERE_SHADOW_SIZE + i;
|
||||
let near_light = mix(shadow.light[i], shadow.light[i + 1u], f);
|
||||
let far_light = mix(shadow.light[far], shadow.light[far + 1u], f);
|
||||
return exp(mix(near_light, far_light, w));
|
||||
}
|
||||
@@ -10,7 +10,8 @@
|
||||
//
|
||||
// The code that includes this file first declares SAMPLES, AREA, the types of the mesh
|
||||
// targets, MeshRadiance and MeshDepth (multisampled when SAMPLES > 1), and
|
||||
// SAMPLE_POSITIONS, where in a pixel each sample lies (bodies.rs).
|
||||
// SAMPLE_POSITIONS, where in a pixel each sample lies, and adds `through_air` and
|
||||
// `earth_shadow` (bodies.rs).
|
||||
|
||||
struct Uniforms {
|
||||
// fx, fy, cx, cy in pixels.
|
||||
@@ -22,7 +23,9 @@ struct Uniforms {
|
||||
earth_center: vec4<f32>,
|
||||
// Columns: Earth-fixed axes in the camera frame.
|
||||
earth_axes: mat3x3<f32>,
|
||||
// x: equatorial / polar radius.
|
||||
// x: equatorial / polar radius. y, z: 1 if the Moon, the Sun may be seen through
|
||||
// the atmosphere (bodies.rs, `behind_air`). w: 1 if the Moon is lit from the table
|
||||
// of Earth's shadow through the atmosphere (bending.wgsl, `earth_shadow`).
|
||||
earth_shape: vec4<f32>,
|
||||
// xyz: unit vector from Earth toward the Sun.
|
||||
earth_sun: vec4<f32>,
|
||||
@@ -274,17 +277,17 @@ fn moon_radiance(p: vec3<f32>, d: vec3<f32>, n: vec3<f32>, at: MapLookup) -> vec
|
||||
}
|
||||
// Earth's shadow: a lunar eclipse.
|
||||
if u.params.z > 0.5 && lit > 0.0 {
|
||||
radiance *= sun_visible(p, u.sun_center, u.earth_center);
|
||||
radiance *= earth_shadow(p);
|
||||
}
|
||||
return radiance;
|
||||
}
|
||||
|
||||
// The Sun's channel radiance along `d`, which meets it.
|
||||
fn sun_radiance(d: vec3<f32>) -> vec4<f32> {
|
||||
// The Sun's channel radiance along the ray from `o`, relative to the Sun's center, along
|
||||
// `d`, which meets it.
|
||||
fn sun_radiance(o: vec3<f32>, d: vec3<f32>) -> vec4<f32> {
|
||||
// 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;
|
||||
return u.sun_radiance * disk_profile(sqrt(max((r - miss) * (r + miss), 0.0)) / r);
|
||||
@@ -305,6 +308,26 @@ fn seen(d: vec3<f32>, b: u32) -> vec3<f32> {
|
||||
return d - beta + dot(beta, d) * d;
|
||||
}
|
||||
|
||||
// A ray from the camera toward the Moon or the Sun after Earth's atmosphere: the
|
||||
// straight line it leaves on, as an origin and a unit direction in the camera frame,
|
||||
// and per channel the light it lets through. `through_air(d, body)` gives it for the
|
||||
// ray along `d` toward body 1 (the Moon) or 2 (the Sun): bending.wgsl, or straight
|
||||
// and whole where no body is seen through the atmosphere (bodies.rs).
|
||||
struct Bent {
|
||||
origin: vec3<f32>,
|
||||
dir: vec3<f32>,
|
||||
transmittance: vec4<f32>,
|
||||
}
|
||||
|
||||
// The distance along `d` to the depth of the point `t` along the bent ray: `t` itself
|
||||
// for a ray that is not bent.
|
||||
fn unbent(bent: Bent, d: vec3<f32>, t: f32) -> f32 {
|
||||
if all(bent.origin == vec3<f32>(0.0)) || t <= 0.0 {
|
||||
return t;
|
||||
}
|
||||
return dot(bent.origin + t * bent.dir, d) / dot(d, d);
|
||||
}
|
||||
|
||||
// Earth's center and the ray `d`, both in Earth-fixed axes stretched along the polar
|
||||
// axis so that Earth is a sphere of its equatorial radius (`ellipsoid_hit`).
|
||||
fn earth_stretched(d: vec3<f32>) -> array<vec3<f32>, 2> {
|
||||
@@ -317,8 +340,10 @@ fn earth_stretched(d: vec3<f32>) -> array<vec3<f32>, 2> {
|
||||
fn nearest(d: vec3<f32>) -> f32 {
|
||||
let to_cam = transpose(u.earth_axes) * (-u.earth_center.xyz);
|
||||
let te = ellipsoid_hit(to_cam, transpose(u.earth_axes) * seen(d, 0u), u.earth_center.w, u.earth_shape.x);
|
||||
let tm = sphere_hit(-u.moon_center.xyz, seen(d, 1u), u.moon_center.w);
|
||||
let ts = sphere_hit(-u.sun_center.xyz, seen(d, 2u), u.sun_center.w);
|
||||
let bm = through_air(seen(d, 1u), 1u);
|
||||
let bs = through_air(seen(d, 2u), 2u);
|
||||
let tm = unbent(bm, d, sphere_hit(bm.origin - u.moon_center.xyz, bm.dir, u.moon_center.w));
|
||||
let ts = unbent(bs, d, sphere_hit(bs.origin - u.sun_center.xyz, bs.dir, u.sun_center.w));
|
||||
var best = FAR;
|
||||
if te > 0.0 {
|
||||
best = te;
|
||||
@@ -346,8 +371,12 @@ fn trace(d: vec3<f32>, spacing: f32) -> Hit {
|
||||
let to_cam = transpose(u.earth_axes) * (-u.earth_center.xyz);
|
||||
let dir = transpose(u.earth_axes) * de;
|
||||
let te = ellipsoid_hit(to_cam, dir, u.earth_center.w, u.earth_shape.x);
|
||||
let tm = sphere_hit(-u.moon_center.xyz, dm, u.moon_center.w);
|
||||
let ts = sphere_hit(-u.sun_center.xyz, ds, u.sun_center.w);
|
||||
let bm = through_air(dm, 1u);
|
||||
let bs = through_air(ds, 2u);
|
||||
let hm = sphere_hit(bm.origin - u.moon_center.xyz, bm.dir, u.moon_center.w);
|
||||
let hs = sphere_hit(bs.origin - u.sun_center.xyz, bs.dir, u.sun_center.w);
|
||||
let tm = unbent(bm, d, hm);
|
||||
let ts = unbent(bs, d, hs);
|
||||
var best = FAR;
|
||||
var radiance = vec4<f32>(0.0);
|
||||
if te > 0.0 && te < best {
|
||||
@@ -359,13 +388,14 @@ fn trace(d: vec3<f32>, spacing: f32) -> Hit {
|
||||
}
|
||||
if tm > 0.0 && tm < best {
|
||||
best = tm;
|
||||
let n = normalize(tm * dm - u.moon_center.xyz);
|
||||
let at = map_lookup_at(transpose(u.moon_axes) * n, tm, abs(dot(n, dm)), spacing, u.moon_center.w);
|
||||
radiance = moon_radiance(tm * dm, dm, n, at);
|
||||
let p = bm.origin + hm * bm.dir;
|
||||
let n = normalize(p - u.moon_center.xyz);
|
||||
let at = map_lookup_at(transpose(u.moon_axes) * n, tm, abs(dot(n, bm.dir)), spacing, u.moon_center.w);
|
||||
radiance = moon_radiance(p, bm.dir, n, at) * bm.transmittance;
|
||||
}
|
||||
if ts > 0.0 && ts < best {
|
||||
best = ts;
|
||||
radiance = sun_radiance(ds);
|
||||
radiance = sun_radiance(bs.origin - u.sun_center.xyz, bs.dir) * bs.transmittance;
|
||||
}
|
||||
if !(te > 0.0 && te == best) {
|
||||
radiance += earth_glow(d, best);
|
||||
@@ -433,10 +463,12 @@ fn limb(c: vec3<f32>, r: f32, d: vec3<f32>) -> vec2<f32> {
|
||||
// even, so that each mesh sample is nearest to the same number of rays.
|
||||
fn grid_size(d: vec3<f32>, pixel: f32) -> u32 {
|
||||
let earth = earth_stretched(d);
|
||||
let moon = through_air(seen(d, 1u), 1u);
|
||||
let sun = through_air(seen(d, 2u), 2u);
|
||||
var edges = array<vec2<f32>, 3>(
|
||||
limb(earth[0], u.earth_center.w, earth[1]),
|
||||
limb(u.moon_center.xyz, u.moon_center.w, seen(d, 1u)),
|
||||
limb(u.sun_center.xyz, u.sun_center.w, seen(d, 2u)),
|
||||
limb(u.moon_center.xyz - moon.origin, u.moon_center.w, moon.dir),
|
||||
limb(u.sun_center.xyz - sun.origin, u.sun_center.w, sun.dir),
|
||||
);
|
||||
var n = 1u;
|
||||
for (var i = 0; i < 3; i++) {
|
||||
@@ -474,19 +506,24 @@ fn fs_main(@builtin(position) pos: vec4<f32>) -> Out {
|
||||
radiance = earth_radiance(te * de, n_fixed, earth_at);
|
||||
}
|
||||
|
||||
let tm = sphere_hit(-u.moon_center.xyz, dm, u.moon_center.w);
|
||||
// The Moon and the Sun, each along its ray through Earth's atmosphere.
|
||||
let bm = through_air(dm, 1u);
|
||||
let hm = sphere_hit(bm.origin - u.moon_center.xyz, bm.dir, u.moon_center.w);
|
||||
let tm = unbent(bm, d, hm);
|
||||
// Where the Moon's map would be read, for every pixel, for the same reason as Earth's.
|
||||
let moon_n = normalize(max(tm, 0.0) * dm - u.moon_center.xyz);
|
||||
let moon_p = bm.origin + max(hm, 0.0) * bm.dir;
|
||||
let moon_n = normalize(moon_p - u.moon_center.xyz);
|
||||
let moon_at = map_lookup(transpose(u.moon_axes) * moon_n);
|
||||
if tm > 0.0 && tm < best {
|
||||
best = tm;
|
||||
radiance = moon_radiance(tm * dm, dm, moon_n, moon_at);
|
||||
radiance = moon_radiance(moon_p, bm.dir, moon_n, moon_at) * bm.transmittance;
|
||||
}
|
||||
|
||||
let ts = sphere_hit(-u.sun_center.xyz, ds, u.sun_center.w);
|
||||
let bs = through_air(ds, 2u);
|
||||
let ts = unbent(bs, d, sphere_hit(bs.origin - u.sun_center.xyz, bs.dir, u.sun_center.w));
|
||||
if ts > 0.0 && ts < best {
|
||||
best = ts;
|
||||
radiance = sun_radiance(ds);
|
||||
radiance = sun_radiance(bs.origin - u.sun_center.xyz, bs.dir) * bs.transmittance;
|
||||
}
|
||||
var glow = vec4<f32>(0.0);
|
||||
if !(te > 0.0 && te == best) {
|
||||
|
||||
@@ -86,12 +86,16 @@ pub const SAMPLES: usize = 64;
|
||||
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| {
|
||||
RULE.get_or_init(|| legendre_rule(RINGS).try_into().expect("one node per ring"))
|
||||
}
|
||||
|
||||
/// Gauss-Legendre nodes and weights for `n` points on [0, pi], as (angle, weight). Found
|
||||
/// by Newton's method on the Legendre polynomial's roots.
|
||||
pub(crate) fn legendre_rule(n: usize) -> Vec<(f64, f64)> {
|
||||
(0..n)
|
||||
.map(|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;
|
||||
@@ -114,7 +118,7 @@ pub fn ring_rule() -> &'static [(f64, f64); RINGS] {
|
||||
let half = std::f64::consts::PI / 2.0;
|
||||
(half * (x + 1.0), half * weight)
|
||||
})
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// How each camera channel's sunlight is spread over the solar disk.
|
||||
|
||||
@@ -4,10 +4,12 @@
|
||||
|
||||
use crate::common;
|
||||
|
||||
use resolved::atmosphere::{AtmosphereTables, BOTTOM_RADIUS, TOP_RADIUS};
|
||||
use resolved::atmosphere::{AtmosphereTables, BOTTOM_RADIUS, GRAZING_BOTTOM, TOP_RADIUS};
|
||||
use resolved::bodies::MOON_RADIUS;
|
||||
use resolved::gpu::AdapterChoice;
|
||||
use resolved::nalgebra::{UnitQuaternion, Vector3};
|
||||
use resolved::nalgebra::{Rotation3, Unit, UnitQuaternion, Vector2, Vector3};
|
||||
use resolved::radiometry::{Radiometry, AU};
|
||||
use resolved::sun::LimbDarkening;
|
||||
use std::f64::consts::PI;
|
||||
use std::sync::Arc;
|
||||
|
||||
@@ -225,3 +227,230 @@ fn a_camera_inside_the_atmosphere_is_refused() {
|
||||
"{error}"
|
||||
);
|
||||
}
|
||||
|
||||
/// Seen from 700 km above 0 N 90 E (+Y), the line of sight in the equatorial plane
|
||||
/// toward -X whose straight line passes `height` above the lowest ray that clears the
|
||||
/// ground, and a pose looking along it with the zenith up the image.
|
||||
fn toward_limb(height: f64) -> (Vector3<f64>, Pose) {
|
||||
let camera = Vector3::new(0.0, BOTTOM_RADIUS + 700e3, 0.0);
|
||||
let seen = ((GRAZING_BOTTOM + height) / camera.norm()).asin();
|
||||
let view = Vector3::new(-seen.sin(), -seen.cos(), 0.0);
|
||||
let pose = Pose::looking_at(camera, camera + view, Vector3::y()).unwrap();
|
||||
(view, pose)
|
||||
}
|
||||
|
||||
/// What `bodies.wgsl`'s `through_air` does to the line of sight through pixel (u, v), in
|
||||
/// f64: the line it leaves the atmosphere on (origin and direction, inertial) and its
|
||||
/// transmittance per channel, or `None` if it meets the ground.
|
||||
fn bent(
|
||||
tables: &AtmosphereTables,
|
||||
pose: &Pose,
|
||||
u: usize,
|
||||
v: usize,
|
||||
) -> Option<(Vector3<f64>, Vector3<f64>, [f64; 4])> {
|
||||
let k = common::example_intrinsics();
|
||||
let d = pose.attitude * k.ray(&Vector2::new(u as f64 + 0.5, v as f64 + 0.5));
|
||||
let c = pose.position;
|
||||
let b = (c - c.dot(&d) * d).norm();
|
||||
let (bend, transmittance) = tables.grazing.ray(b)?;
|
||||
let turn = Rotation3::from_axis_angle(&Unit::new_normalize(d.cross(&-c)), bend);
|
||||
Some((turn * c, turn * d, transmittance))
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_sun_through_the_limb_is_bent_and_dimmed() {
|
||||
// The Sun placed where the line of sight through the image center goes after the
|
||||
// atmosphere bends it, behind air 12, 20 and 30 km up. Down the image's middle
|
||||
// column, each pixel takes the Sun's radiance where its bent ray meets the disk,
|
||||
// times its transmittance. At 12 km the bending, 0.3 degrees, is more than the
|
||||
// Sun's radius: the Sun is seen where a straight line would miss it.
|
||||
let _devices = common::device_lock();
|
||||
let (r, _, _) = constants();
|
||||
let tables = AtmosphereTables::new(&r);
|
||||
let limb = LimbDarkening::new(&r);
|
||||
for adapter in &ADAPTERS {
|
||||
let mut renderer = renderer(adapter);
|
||||
let name = renderer.adapter_info().name.clone();
|
||||
for height in [12e3, 20e3, 30e3] {
|
||||
let (_, pose) = toward_limb(height);
|
||||
let (origin, dir, _) = bent(&tables, &pose, 959, 539).unwrap();
|
||||
let sun = origin + AU * dir;
|
||||
let (scene, _) = view(sun, pose.position, pose.position + dir);
|
||||
let frame = renderer.submit(&scene, &pose).unwrap().wait().unwrap();
|
||||
let radiance = frame.radiance.unwrap();
|
||||
let mut compared = 0;
|
||||
for v in 400..680 {
|
||||
let Some((origin, dir, t)) = bent(&tables, &pose, 959, v) else {
|
||||
continue;
|
||||
};
|
||||
let o = origin - sun;
|
||||
let miss = (o - o.dot(&dir) * dir).norm() / resolved::radiometry::SUN_RADIUS;
|
||||
// Inside the disk, away from its edge, where the air lets enough through
|
||||
// that the glow over the limb does not count.
|
||||
if miss > 0.97 || t[0] < 0.05 {
|
||||
continue;
|
||||
}
|
||||
let profile = limb.profile((1.0 - miss * miss).sqrt());
|
||||
let got = radiance[v * 1920 + 959];
|
||||
for k in 0..3 {
|
||||
let want = r.sun_radiance(k) * profile[k] * t[k];
|
||||
let error = f64::from(got[k]) / want - 1.0;
|
||||
// Measured: 7.3e-4 at most.
|
||||
assert!(
|
||||
error.abs() < 2e-3,
|
||||
"{name} {height} m, row {v}, channel {k}: {error}"
|
||||
);
|
||||
}
|
||||
compared += 1;
|
||||
}
|
||||
assert!(compared > 30, "{name} {height} m: {compared} rows");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_moon_through_the_limb_is_bent_and_dimmed() {
|
||||
// The Moon placed as the Sun is above, behind air 15 km up, uniform and lit from
|
||||
// beside the camera and above. Down the image's middle column the light falls from
|
||||
// the disk's center to its edge with the cosine of the Sun's incidence where each
|
||||
// bent ray meets it: each pixel's radiance over its transmittance, that cosine and
|
||||
// the phase curve is the same everywhere if the rays are bent where the shader
|
||||
// bends them. The air's glow, the same with or without the Moon behind it, is
|
||||
// taken out with an image of the Moon moved out of view.
|
||||
let _devices = common::device_lock();
|
||||
let (r, _, _) = constants();
|
||||
let tables = AtmosphereTables::new(&r);
|
||||
for adapter in &ADAPTERS {
|
||||
let mut renderer = renderer(adapter);
|
||||
let name = renderer.adapter_info().name.clone();
|
||||
let (view_dir, pose) = toward_limb(15e3);
|
||||
let (origin, dir, _) = bent(&tables, &pose, 959, 539).unwrap();
|
||||
let moon = origin + 3.844e8 * dir;
|
||||
let to_sun = (Vector3::z() * 2.0 - view_dir).normalize();
|
||||
let scene = Scene {
|
||||
earth_attitude: UnitQuaternion::identity(),
|
||||
sun: moon + AU * to_sun,
|
||||
moon,
|
||||
moon_attitude: UnitQuaternion::identity(),
|
||||
targets: vec![],
|
||||
velocities: Default::default(),
|
||||
};
|
||||
let frame = renderer.submit(&scene, &pose).unwrap().wait().unwrap();
|
||||
let radiance = frame.radiance.unwrap();
|
||||
let away = Scene {
|
||||
moon: -moon,
|
||||
..scene.clone()
|
||||
};
|
||||
let frame = renderer.submit(&away, &pose).unwrap().wait().unwrap();
|
||||
let glow = frame.radiance.unwrap();
|
||||
let mut ratios = vec![];
|
||||
for v in 300..780 {
|
||||
let Some((origin, dir, t)) = bent(&tables, &pose, 959, v) else {
|
||||
continue;
|
||||
};
|
||||
let o = origin - moon;
|
||||
let along = -o.dot(&dir);
|
||||
let miss = (o + along * dir).norm();
|
||||
// Away from the edge, where the cosine changes fastest with where the ray
|
||||
// meets the disk: a 3 m error in f32 in the impact parameter moves the
|
||||
// bending by 1.3e-6 rad, 500 m at the Moon.
|
||||
if miss > 0.9 * MOON_RADIUS || t[0] < 0.05 {
|
||||
continue;
|
||||
}
|
||||
let hit = o + (along - (MOON_RADIUS * MOON_RADIUS - miss * miss).sqrt()) * dir;
|
||||
let lit = (hit / MOON_RADIUS).dot(&to_sun);
|
||||
let phase = to_sun.cross(&-dir).norm().atan2(to_sun.dot(&-dir));
|
||||
let (got, glow) = (radiance[v * 1920 + 959], glow[v * 1920 + 959]);
|
||||
ratios.push(std::array::from_fn::<f64, 3, _>(|k| {
|
||||
f64::from(got[k] - glow[k]) / (t[k] * lit * resolved::moon::phase_curve(phase))
|
||||
}));
|
||||
}
|
||||
assert!(ratios.len() > 30, "{name}: {} rows", ratios.len());
|
||||
for k in 0..3 {
|
||||
let mean = ratios.iter().map(|q| q[k]).sum::<f64>() / ratios.len() as f64;
|
||||
let worst = ratios
|
||||
.iter()
|
||||
.map(|q| (q[k] / mean - 1.0).abs())
|
||||
.fold(0.0, f64::max);
|
||||
// Measured: 1.2e-3 at most, on lavapipe.
|
||||
assert!(worst < 3e-3, "{name} channel {k}: {worst}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_eclipsed_moon_is_lit_through_the_air() {
|
||||
// The uniform Moon in Earth's shadow, across the edge of the umbra (2300 to 5700 km
|
||||
// from its axis), seen from 20,000 km on its sunward side, and again out of the
|
||||
// shadow. Each pixel's radiance over the cosine of the Sun's incidence, the phase
|
||||
// curve and the sunlight that reaches its point on the Moon
|
||||
// (`Grazing::sunlight`, in f64) is the same everywhere, in and out of the shadow,
|
||||
// if the shader looks the light up where it should.
|
||||
let _devices = common::device_lock();
|
||||
let (r, _, _) = constants();
|
||||
let tables = AtmosphereTables::new(&r);
|
||||
let limb = LimbDarkening::new(&r);
|
||||
let k = common::example_intrinsics();
|
||||
let sun = Vector3::new(AU, 0.0, 0.0);
|
||||
for adapter in &ADAPTERS {
|
||||
let mut renderer = renderer(adapter);
|
||||
let name = renderer.adapter_info().name.clone();
|
||||
// Per channel: each sampled pixel's radiance over what it should be, but for a
|
||||
// constant.
|
||||
let mut ratios = |off: f64| {
|
||||
let moon = Vector3::new(-3.844e8, off, 0.0);
|
||||
let scene = Scene {
|
||||
earth_attitude: UnitQuaternion::identity(),
|
||||
sun,
|
||||
moon,
|
||||
moon_attitude: UnitQuaternion::identity(),
|
||||
targets: vec![],
|
||||
velocities: Default::default(),
|
||||
};
|
||||
let pose = Pose::looking_at(moon + Vector3::new(2e7, 0.0, 0.0), moon, Vector3::z());
|
||||
let pose = pose.unwrap();
|
||||
let frame = renderer.submit(&scene, &pose).unwrap().wait().unwrap();
|
||||
let radiance = frame.radiance.unwrap();
|
||||
let to_sun = (sun - moon).normalize();
|
||||
let mut ratios = vec![];
|
||||
for v in (0..1080).step_by(12) {
|
||||
for u in (0..1920).step_by(12) {
|
||||
let d = pose.attitude * k.ray(&Vector2::new(u as f64 + 0.5, v as f64 + 0.5));
|
||||
let o = pose.position - moon;
|
||||
let along = -o.dot(&d);
|
||||
let miss = (o + along * d).norm();
|
||||
if miss > 0.9 * MOON_RADIUS {
|
||||
continue;
|
||||
}
|
||||
let hit = o + (along - (MOON_RADIUS * MOON_RADIUS - miss * miss).sqrt()) * d;
|
||||
let lit = (hit / MOON_RADIUS).dot(&to_sun);
|
||||
let phase = to_sun.cross(&-d).norm().atan2(to_sun.dot(&-d));
|
||||
let light = tables.grazing.sunlight(
|
||||
resolved::atmosphere::stretched(&scene, moon + hit),
|
||||
resolved::atmosphere::stretched(&scene, sun),
|
||||
&limb,
|
||||
);
|
||||
let got = radiance[v * 1920 + u];
|
||||
ratios.push(std::array::from_fn::<f64, 3, _>(|c| {
|
||||
f64::from(got[c]) / (light[c] * lit * resolved::moon::phase_curve(phase))
|
||||
}));
|
||||
}
|
||||
}
|
||||
ratios
|
||||
};
|
||||
let outside = ratios(2e7);
|
||||
let inside = ratios(4e6);
|
||||
assert!(inside.len() > 3000, "{name}: {} pixels", inside.len());
|
||||
for c in 0..3 {
|
||||
let scale = outside.iter().map(|q| q[c]).sum::<f64>() / outside.len() as f64;
|
||||
let worst = inside
|
||||
.iter()
|
||||
.chain(&outside)
|
||||
.map(|q| (q[c] / scale - 1.0).abs())
|
||||
.fold(0.0, f64::max);
|
||||
// Measured: 4.7e-3 on the RTX 3060, 3.6e-3 on lavapipe, the table's
|
||||
// interpolation (`EarthShadow`).
|
||||
assert!(worst < 7e-3, "{name} channel {c}: {worst}");
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -42,9 +42,12 @@ fn renderer(adapter: &AdapterChoice, lighting: LightingTier) -> Renderer {
|
||||
moon_texture: Some(common::uniform_moon()),
|
||||
// The same Moon at L0 (flat albedo) and L1 (a uniform map with a flat color), so
|
||||
// the two images differ only by the eclipse and L1's phase factor (`with_phase`).
|
||||
// Without Earth's atmosphere, so Earth's shadow has its hard edge (the shadow
|
||||
// through the air is tested in atmosphere.rs).
|
||||
appearance: Appearance {
|
||||
moon_albedo: Spectrum::constant(0.125),
|
||||
moon_color: Spectrum::constant(1.0),
|
||||
earth_atmosphere: false,
|
||||
..Default::default()
|
||||
},
|
||||
read_radiance: true,
|
||||
|
||||
Reference in new issue
Block a user