694 lines
25 KiB
Rust
694 lines
25 KiB
Rust
//! Light from Earth and the Moon on targets at lighting L1 (Earthshine and Moonshine),
|
|
//! on the GPU and on the CPU (lavapipe), against a quadrature over each body in f64.
|
|
|
|
mod common;
|
|
|
|
use std::f64::consts::{PI, TAU};
|
|
use std::io::Cursor;
|
|
use std::sync::Arc;
|
|
|
|
use resolved::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS};
|
|
use resolved::eclipse::sun_visible_fraction;
|
|
use resolved::gpu::bodies::BodyRadiometry;
|
|
use resolved::gpu::AdapterChoice;
|
|
use resolved::nalgebra::{Matrix3, UnitQuaternion, Vector3};
|
|
use resolved::radiometry::{Radiometry, AU};
|
|
use resolved::{
|
|
CameraConfig, EarthTexture, Fidelity, LightingTier, LoadOptions, Material, Mesh, Pose,
|
|
Renderer, RendererOptions, Scene, Spectrum, Target,
|
|
};
|
|
|
|
const ALBEDO: f64 = 0.5;
|
|
const MOON_DISTANCE: f64 = 3.844e8;
|
|
|
|
/// An 8 m square plate facing body +X, gray.
|
|
fn plate() -> Mesh {
|
|
let mut s = String::from("solid plate\n");
|
|
let c = [
|
|
[0.0, -4.0, -4.0],
|
|
[0.0, 4.0, -4.0],
|
|
[0.0, 4.0, 4.0],
|
|
[0.0, -4.0, 4.0],
|
|
];
|
|
for tri in [[c[0], c[1], c[2]], [c[0], c[2], c[3]]] {
|
|
s += "facet normal 0 0 0\nouter loop\n";
|
|
for [x, y, z] in tri {
|
|
s += &format!("vertex {x} {y} {z}\n");
|
|
}
|
|
s += "endloop\nendfacet\n";
|
|
}
|
|
s += "endsolid plate\n";
|
|
let options = LoadOptions {
|
|
default_material: Material::gray("gray", ALBEDO),
|
|
..Default::default()
|
|
};
|
|
Mesh::from_stl(&mut Cursor::new(s), "plate", &options).unwrap()
|
|
}
|
|
|
|
/// What a test varies about the bodies: Earth's texture (one sRGB gray level, no haze)
|
|
/// and the Moon's albedo.
|
|
#[derive(Clone, Copy)]
|
|
struct Surfaces {
|
|
earth_level: u8,
|
|
moon_albedo: f64,
|
|
}
|
|
|
|
impl Surfaces {
|
|
const WHITE_EARTH: Self = Self {
|
|
earth_level: 255,
|
|
moon_albedo: 0.12,
|
|
};
|
|
|
|
fn appearance(&self) -> Appearance {
|
|
Appearance {
|
|
earth_haze: Spectrum::constant(0.0),
|
|
moon_albedo: Spectrum::constant(self.moon_albedo),
|
|
..Default::default()
|
|
}
|
|
}
|
|
|
|
/// Per channel: Earth's albedo and the Moon's, as the renderer folds them.
|
|
fn albedo(&self) -> ([f64; 3], [f64; 3]) {
|
|
let r = Radiometry::new(&CameraConfig::example());
|
|
let b = BodyRadiometry::new(&r, &self.appearance()).unwrap();
|
|
let level = srgb_to_linear(self.earth_level);
|
|
let earth = [0, 1, 2].map(|k| (0..3).map(|j| b.earth_rgb[j][k] as f64 * level).sum());
|
|
(earth, [self.moon_albedo; 3])
|
|
}
|
|
|
|
fn renderer(&self, adapter: &AdapterChoice, lighting: LightingTier) -> Renderer {
|
|
let v = self.earth_level;
|
|
let texture = EarthTexture::from_rgba8(1, 1, vec![v, v, v, 255]).unwrap();
|
|
let options = RendererOptions {
|
|
adapter: adapter.clone(),
|
|
fidelity: Fidelity {
|
|
lighting,
|
|
..Default::default()
|
|
},
|
|
earth_texture: Some(Arc::new(texture)),
|
|
appearance: self.appearance(),
|
|
read_radiance: true,
|
|
..Default::default()
|
|
};
|
|
Renderer::new(CameraConfig::example(), options).unwrap()
|
|
}
|
|
}
|
|
|
|
fn srgb_to_linear(level: u8) -> f64 {
|
|
let c = level as f64 / 255.0;
|
|
if c <= 0.04045 {
|
|
c / 12.92
|
|
} else {
|
|
((c + 0.055) / 1.055).powf(2.4)
|
|
}
|
|
}
|
|
|
|
/// One plate, its normal and the bodies around it.
|
|
#[derive(Clone, Copy, Debug)]
|
|
struct Case {
|
|
name: &'static str,
|
|
target: Vector3<f64>,
|
|
normal: Vector3<f64>,
|
|
sun: Vector3<f64>,
|
|
moon: Vector3<f64>,
|
|
earth_attitude: UnitQuaternion<f64>,
|
|
/// None: the camera looks straight at the plate. Some((psi, off)): from `off`
|
|
/// degrees off its normal, placed so the normal points `psi` degrees around the
|
|
/// image center from the image's +X axis (toward +Y).
|
|
view: Option<(f64, f64)>,
|
|
}
|
|
|
|
impl Case {
|
|
fn new(name: &'static str, target: Vector3<f64>, normal: Vector3<f64>) -> Self {
|
|
Self {
|
|
name,
|
|
target,
|
|
normal: normal.normalize(),
|
|
sun: Vector3::new(AU, 0.0, 0.0),
|
|
// Well away from every plate normal used, and visible from most places.
|
|
moon: MOON_DISTANCE * Vector3::new(-0.5, -0.5, -0.7).normalize(),
|
|
earth_attitude: UnitQuaternion::identity(),
|
|
view: None,
|
|
}
|
|
}
|
|
|
|
fn sun_at(self, degrees_from_x: f64) -> Self {
|
|
let a = degrees_from_x.to_radians();
|
|
Self {
|
|
sun: AU * Vector3::new(a.cos(), a.sin(), 0.0),
|
|
..self
|
|
}
|
|
}
|
|
|
|
fn scene(&self, model: resolved::ModelId) -> Scene {
|
|
let attitude = UnitQuaternion::rotation_between(&Vector3::x(), &self.normal)
|
|
.unwrap_or_else(|| UnitQuaternion::from_axis_angle(&Vector3::z_axis(), PI));
|
|
Scene {
|
|
earth_attitude: self.earth_attitude,
|
|
sun: self.sun,
|
|
moon: self.moon,
|
|
targets: vec![Target {
|
|
model,
|
|
pose: Pose::new(self.target, attitude),
|
|
}],
|
|
}
|
|
}
|
|
|
|
/// 20 m from the plate's center, along its normal or off it, looking at the center.
|
|
fn observer(&self) -> Pose {
|
|
let up = if self.normal.z.abs() < 0.9 {
|
|
Vector3::z()
|
|
} else {
|
|
Vector3::x()
|
|
};
|
|
let straight = Pose::looking_at(self.target + 20.0 * self.normal, self.target, up).unwrap();
|
|
let Some((psi, off)) = self.view else {
|
|
return straight;
|
|
};
|
|
// The camera frame normal wanted, turned onto the head-on view's (0, 0, -1).
|
|
let (psi, off) = (psi.to_radians(), off.to_radians());
|
|
let normal = Vector3::new(off.sin() * psi.cos(), off.sin() * psi.sin(), -off.cos());
|
|
let turn = UnitQuaternion::rotation_between(&normal, &-Vector3::z()).unwrap();
|
|
let attitude = straight.attitude * turn;
|
|
Pose::new(self.target - 20.0 * (attitude * Vector3::z()), attitude)
|
|
}
|
|
|
|
/// Direct sunlight on the plate, as a fraction of the Sun's irradiance there.
|
|
fn direct(&self) -> f64 {
|
|
let visible = sun_visible_fraction(
|
|
self.target,
|
|
self.sun,
|
|
&[
|
|
(Vector3::zeros(), EARTH_EQUATORIAL_RADIUS),
|
|
(self.moon, MOON_RADIUS),
|
|
],
|
|
);
|
|
self.normal
|
|
.dot(&(self.sun - self.target).normalize())
|
|
.max(0.0)
|
|
* visible
|
|
}
|
|
}
|
|
|
|
/// A body for the quadrature: an ellipsoid with axes `attitude` (fixed to inertial).
|
|
struct Body {
|
|
center: Vector3<f64>,
|
|
radius: f64,
|
|
/// Body-fixed axes in the inertial frame, as columns, with Z stretched by the ratio
|
|
/// of equatorial to polar radius: inertial offsets times its transpose land in a
|
|
/// frame where the body is a sphere of `radius`.
|
|
to_sphere: Matrix3<f64>,
|
|
albedo: [f64; 3],
|
|
}
|
|
|
|
impl Body {
|
|
fn new(
|
|
center: Vector3<f64>,
|
|
radius: f64,
|
|
flattening: f64,
|
|
attitude: UnitQuaternion<f64>,
|
|
albedo: [f64; 3],
|
|
) -> Self {
|
|
let stretch = Matrix3::from_diagonal(&Vector3::new(1.0, 1.0, flattening));
|
|
let to_sphere = stretch * attitude.to_rotation_matrix().matrix().transpose();
|
|
Self {
|
|
center,
|
|
radius,
|
|
to_sphere,
|
|
albedo,
|
|
}
|
|
}
|
|
|
|
/// Distance along unit `d` from `o` to the body, if the ray meets it.
|
|
fn distance(&self, o: Vector3<f64>, d: Vector3<f64>) -> Option<f64> {
|
|
let os = self.to_sphere * (o - self.center);
|
|
let ds = self.to_sphere * d;
|
|
// |os + s ds|^2 = r^2, s the distance along the unstretched ray.
|
|
let a = ds.norm_squared();
|
|
let b = os.dot(&ds);
|
|
let c = os.norm_squared() - self.radius * self.radius;
|
|
let disc = b * b - a * c;
|
|
let s = (-b - disc.max(0.0).sqrt()) / a;
|
|
(disc >= 0.0 && s > 0.0).then_some(s)
|
|
}
|
|
|
|
/// The outward unit normal at the inertial point `p` on the surface: the gradient
|
|
/// of |to_sphere (p - center)|^2.
|
|
fn normal(&self, p: Vector3<f64>) -> Vector3<f64> {
|
|
(self.to_sphere.transpose() * (self.to_sphere * (p - self.center))).normalize()
|
|
}
|
|
}
|
|
|
|
/// The light Earth and the Moon send onto a surface at `p` facing `n`, per channel, in
|
|
/// units of the Sun's irradiance at `p`. A midpoint rule over each body's cone, with
|
|
/// `rings` x `sectors` cells of equal solid angle.
|
|
fn quadrature(case: &Case, surfaces: &Surfaces, rings: usize, sectors: usize) -> [f64; 3] {
|
|
let (earth_albedo, moon_albedo) = surfaces.albedo();
|
|
let earth = Body::new(
|
|
Vector3::zeros(),
|
|
EARTH_EQUATORIAL_RADIUS,
|
|
EARTH_EQUATORIAL_RADIUS / EARTH_POLAR_RADIUS,
|
|
case.earth_attitude,
|
|
earth_albedo,
|
|
);
|
|
let moon = Body::new(
|
|
case.moon,
|
|
MOON_RADIUS,
|
|
1.0,
|
|
UnitQuaternion::identity(),
|
|
moon_albedo,
|
|
);
|
|
let (p, n) = (case.target, case.normal);
|
|
let mut total = [0.0; 3];
|
|
for (body, other) in [(&earth, &moon), (&moon, &earth)] {
|
|
let to = body.center - p;
|
|
let sin = (body.radius / to.norm()).min(1.0);
|
|
let vers = sin * sin / (1.0 + (1.0 - sin * sin).sqrt());
|
|
let cell = TAU * vers / (rings * sectors) as f64;
|
|
let axis = to.normalize();
|
|
let u = axis
|
|
.cross(&Vector3::z())
|
|
.try_normalize(1e-6)
|
|
.unwrap_or(Vector3::x());
|
|
let v = axis.cross(&u);
|
|
let to_sun = (case.sun - body.center).normalize();
|
|
let ratio = ((case.sun - p).norm() / (case.sun - body.center).norm()).powi(2);
|
|
// Rings in parallel: the test builds without optimization.
|
|
let ring = |i: usize| {
|
|
let mut sum = [0.0; 3];
|
|
let vi = (i as f64 + 0.5) / rings as f64 * vers;
|
|
let sin_t = (vi * (2.0 - vi)).sqrt();
|
|
for j in 0..sectors {
|
|
let phi = TAU * (j as f64 + 0.5) / sectors as f64;
|
|
let d = (1.0 - vi) * axis + sin_t * (phi.cos() * u + phi.sin() * v);
|
|
let facing = n.dot(&d);
|
|
if facing <= 0.0 {
|
|
continue;
|
|
}
|
|
let Some(t) = body.distance(p, d) else {
|
|
continue;
|
|
};
|
|
if other.distance(p, d).is_some_and(|s| s < t) {
|
|
continue;
|
|
}
|
|
// The other body's shadow on this one: a lunar or a solar eclipse.
|
|
let point = p + t * d;
|
|
let shadow = sun_visible_fraction(point, case.sun, &[(other.center, other.radius)]);
|
|
let lit = body.normal(point).dot(&to_sun).max(0.0) * shadow;
|
|
for (s, albedo) in sum.iter_mut().zip(body.albedo) {
|
|
*s += albedo * lit / PI * ratio * cell * facing;
|
|
}
|
|
}
|
|
sum
|
|
};
|
|
let threads = std::thread::available_parallelism().map_or(1, |n| n.get());
|
|
let sums: Vec<[f64; 3]> = std::thread::scope(|scope| {
|
|
let workers: Vec<_> = (0..threads)
|
|
.map(|w| {
|
|
scope.spawn(move || (w..rings).step_by(threads).map(ring).collect::<Vec<_>>())
|
|
})
|
|
.collect();
|
|
workers
|
|
.into_iter()
|
|
.flat_map(|h| h.join().unwrap())
|
|
.collect()
|
|
});
|
|
for sum in sums {
|
|
for (t, s) in total.iter_mut().zip(sum) {
|
|
*t += s;
|
|
}
|
|
}
|
|
}
|
|
total
|
|
}
|
|
|
|
/// The light from the bodies on each case's plate as rendered, per channel, in units of
|
|
/// the Sun's irradiance there: the center pixel's radiance over a Sun-facing plate's,
|
|
/// less the direct sunlight.
|
|
fn rendered(renderer: &mut Renderer, cases: &[Case]) -> Vec<[f64; 3]> {
|
|
let model = renderer.load_model(&plate()).unwrap();
|
|
let radiometry = Radiometry::new(&CameraConfig::example());
|
|
cases
|
|
.iter()
|
|
.map(|case| {
|
|
let frame = renderer
|
|
.submit(&case.scene(model), &case.observer())
|
|
.unwrap()
|
|
.wait()
|
|
.unwrap();
|
|
let radiance = frame.radiance.unwrap();
|
|
common::save(
|
|
&format!("earthshine-{}", case.name),
|
|
renderer.adapter_info(),
|
|
&radiance,
|
|
);
|
|
let pixel = radiance[540 * 1920 + 960];
|
|
let distance = (case.sun - case.target).norm();
|
|
[0, 1, 2].map(|k| {
|
|
pixel[k] as f64 / radiometry.lambert_radiance(k, ALBEDO, 1.0, distance)
|
|
- case.direct()
|
|
})
|
|
})
|
|
.collect()
|
|
}
|
|
|
|
/// Renders `cases` on both devices and compares each channel with the quadrature:
|
|
/// within 0.5 % plus `floor` (in units of the Sun's irradiance), and below 1e-7 where
|
|
/// the quadrature finds nothing at all.
|
|
///
|
|
/// Measured: 0.43 % at most wherever the bodies send more than 1e-3 of sunlight. The
|
|
/// largest errors are at the peaks (a plate facing a full Earth or Moon), where the
|
|
/// cube's bilinear interpolation between texels 5.6 degrees apart rounds the top off;
|
|
/// 32-texel faces cut those to 0.13 % but left the cone sampling's 0.44 % on a half-lit
|
|
/// GEO Earth. A thin crescent at Earth's limb from low orbit is off by 4 %, 7e-5 of
|
|
/// sunlight, which is what `floor` allows for there (about 0.1 DN on a gray plate in
|
|
/// the example camera); 64 rings of 16 cells halve that but put the half-lit GEO Earth
|
|
/// at 0.6 %.
|
|
fn compare(surfaces: Surfaces, floor: f64, cases: &[Case]) {
|
|
let want: Vec<_> = cases
|
|
.iter()
|
|
.map(|c| quadrature(c, &surfaces, 300, 300))
|
|
.collect();
|
|
let _devices = common::device_lock();
|
|
for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
|
|
let mut renderer = surfaces.renderer(&adapter, LightingTier::L1);
|
|
let name = renderer.adapter_info().name.clone();
|
|
let got = rendered(&mut renderer, cases);
|
|
for ((case, got), want) in cases.iter().zip(&got).zip(&want) {
|
|
for k in 0..3 {
|
|
let (got, want) = (got[k], want[k]);
|
|
let ok = if want == 0.0 {
|
|
got.abs() < 1e-7
|
|
} else {
|
|
(got - want).abs() <= 0.005 * want + floor
|
|
};
|
|
assert!(
|
|
ok,
|
|
"{name}, {}, channel {k}: {got:e}, want {want:e}",
|
|
case.name
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
const LEO: Vector3<f64> = Vector3::new(6_878_000.0, 0.0, 0.0);
|
|
const GEO: Vector3<f64> = Vector3::new(42_164_000.0, 0.0, 0.0);
|
|
|
|
#[test]
|
|
fn low_orbit() {
|
|
let down = -Vector3::x();
|
|
let tilted = UnitQuaternion::from_axis_angle(&Vector3::y_axis(), 0.5);
|
|
compare(
|
|
Surfaces::WHITE_EARTH,
|
|
1e-4,
|
|
&[
|
|
// Over the subsolar point, looking down, sideways and up.
|
|
Case::new("leo-down", LEO, down),
|
|
Case::new("leo-side", LEO, Vector3::y()),
|
|
Case::new("leo-up", LEO, Vector3::x()),
|
|
// Toward the terminator: less lit Earth in view, then a thin crescent.
|
|
Case::new("leo-down-sun60", LEO, down).sun_at(60.0),
|
|
Case::new("leo-down-sun100", LEO, down).sun_at(100.0),
|
|
Case::new("leo-side-sun100", LEO, Vector3::new(-1.0, 1.0, 0.0)).sun_at(100.0),
|
|
// In Earth's shadow, 500 km up only the night side is in view.
|
|
Case::new("leo-down-umbra", LEO, down).sun_at(150.0),
|
|
// Earth's axes turned: the flattened ellipsoid moves with them.
|
|
Case {
|
|
earth_attitude: tilted,
|
|
..Case::new("leo-down-tilted", LEO, Vector3::new(-1.0, 0.0, 0.3))
|
|
},
|
|
],
|
|
);
|
|
}
|
|
|
|
#[test]
|
|
fn geostationary() {
|
|
compare(
|
|
Surfaces::WHITE_EARTH,
|
|
1e-8,
|
|
&[
|
|
Case::new("geo-down", GEO, -Vector3::x()),
|
|
Case::new("geo-down-sun90", GEO, -Vector3::x()).sun_at(90.0),
|
|
Case::new("geo-side", GEO, Vector3::new(-0.2, 1.0, 0.0)),
|
|
// 1000 Earth radii sunward of Earth, where the target gets 8 % more sunlight
|
|
// than Earth does.
|
|
Case::new(
|
|
"far-down",
|
|
Vector3::new(1000.0 * EARTH_EQUATORIAL_RADIUS, 0.0, 0.0),
|
|
-Vector3::x(),
|
|
),
|
|
],
|
|
);
|
|
}
|
|
|
|
#[test]
|
|
fn every_cube_face_gives_the_same_light() {
|
|
// The irradiance cube is laid out in the camera's axes, and a surface seen head-on
|
|
// reads the face behind the camera (-Z) at its center, where the four texels around
|
|
// weigh equally. Seen 60 degrees off its normal from four sides, the same plate
|
|
// reads +X, +Y, -X and -Y, and 25 degrees off, -Z away from its center; each off
|
|
// center in both texel directions. (+Z faces away from the camera, and only a
|
|
// surface seen edge-on at the image's edge reads it.) The Sun is out of every
|
|
// symmetry plane, so Earth's light differs on each side of the plate's normal.
|
|
let base = Case {
|
|
sun: AU * Vector3::new(0.5, 0.7, 0.5).normalize(),
|
|
..Case::new("faces", LEO, Vector3::new(-1.0, 0.3, 0.2))
|
|
};
|
|
let views = [
|
|
(20.0, 60.0),
|
|
(110.0, 60.0),
|
|
(200.0, 60.0),
|
|
(290.0, 60.0),
|
|
(35.0, 25.0),
|
|
];
|
|
let cases: Vec<Case> = views
|
|
.into_iter()
|
|
.map(|v| Case {
|
|
view: Some(v),
|
|
..base
|
|
})
|
|
.collect();
|
|
compare(Surfaces::WHITE_EARTH, 1e-4, &cases);
|
|
}
|
|
|
|
#[test]
|
|
fn lunar_orbit() {
|
|
let moon = Vector3::new(MOON_DISTANCE, 0.0, 0.0);
|
|
let target = moon + Vector3::new(0.0, MOON_RADIUS + 100e3, 0.0);
|
|
let case = |name, sun| Case {
|
|
moon,
|
|
..Case::new(name, target, -Vector3::y()).sun_at(sun)
|
|
};
|
|
compare(
|
|
Surfaces::WHITE_EARTH,
|
|
1e-8,
|
|
&[case("lunar-full", 90.0), case("lunar-quarter", 0.0)],
|
|
);
|
|
}
|
|
|
|
#[test]
|
|
fn bodies_hide_each_other() {
|
|
let moon = Vector3::new(MOON_DISTANCE, 0.0, 0.0);
|
|
let case = |name, target, normal| Case {
|
|
moon,
|
|
..Case::new(name, target, normal).sun_at(90.0)
|
|
};
|
|
let above = MOON_RADIUS + 100e3;
|
|
// A black Moon in front of Earth, and beside it: Earth's light is gone, then not.
|
|
compare(
|
|
Surfaces {
|
|
earth_level: 255,
|
|
moon_albedo: 0.0,
|
|
},
|
|
1e-8,
|
|
&[
|
|
case(
|
|
"moon-hides-earth",
|
|
moon + Vector3::new(above, 0.0, 0.0),
|
|
-Vector3::x(),
|
|
),
|
|
case(
|
|
"moon-beside-earth",
|
|
moon + Vector3::new(0.0, above, 0.0),
|
|
-Vector3::x(),
|
|
),
|
|
],
|
|
);
|
|
// A black Earth in front of a white Moon, and beside it.
|
|
compare(
|
|
Surfaces {
|
|
earth_level: 0,
|
|
moon_albedo: 1.0,
|
|
},
|
|
1e-8,
|
|
&[
|
|
case(
|
|
"earth-hides-moon",
|
|
Vector3::new(-7.0e6, 0.0, 0.0),
|
|
Vector3::x(),
|
|
),
|
|
case(
|
|
"earth-beside-moon",
|
|
Vector3::new(0.0, 7.0e6, 0.0),
|
|
Vector3::x(),
|
|
),
|
|
],
|
|
);
|
|
}
|
|
|
|
#[test]
|
|
fn eclipses_dim_the_light_from_the_bodies() {
|
|
// The Moon in Earth's umbra sheds nothing; across the umbra's edge, part. Earth
|
|
// with the Moon's shadow on it, seen from inside that shadow.
|
|
let eclipsed_moon = |name, offset| {
|
|
let moon = Vector3::new(-MOON_DISTANCE, offset, 0.0);
|
|
Case {
|
|
moon,
|
|
..Case::new(
|
|
name,
|
|
moon + Vector3::new(MOON_RADIUS + 100e3, 0.0, 0.0),
|
|
-Vector3::x(),
|
|
)
|
|
}
|
|
};
|
|
compare(
|
|
Surfaces::WHITE_EARTH,
|
|
1e-8,
|
|
&[
|
|
eclipsed_moon("moon-in-umbra", 0.0),
|
|
eclipsed_moon("moon-in-penumbra", 5.5e6),
|
|
],
|
|
);
|
|
let solar = Case {
|
|
moon: Vector3::new(3.6e8, 0.0, 0.0),
|
|
..Case::new("earth-under-the-moon", LEO, -Vector3::x())
|
|
};
|
|
compare(Surfaces::WHITE_EARTH, 1e-4, &[solar]);
|
|
}
|
|
|
|
#[test]
|
|
fn lighting_l0_has_no_earthshine() {
|
|
let _devices = common::device_lock();
|
|
let case = Case::new("l0-leo-down", LEO, -Vector3::x());
|
|
for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
|
|
let mut renderer = Surfaces::WHITE_EARTH.renderer(&adapter, LightingTier::L0);
|
|
assert_eq!(rendered(&mut renderer, &[case])[0], [0.0; 3]);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn quadrature_has_converged() {
|
|
// Halving the cells in each direction moves the result by less than 5e-4 (measured
|
|
// 2e-4), against the 0.5 % the renders are held to.
|
|
for case in [
|
|
Case::new("leo-down-sun60", LEO, -Vector3::x()).sun_at(60.0),
|
|
Case::new("geo-down-sun90", GEO, -Vector3::x()).sun_at(90.0),
|
|
] {
|
|
let fine = quadrature(&case, &Surfaces::WHITE_EARTH, 300, 300);
|
|
let coarse = quadrature(&case, &Surfaces::WHITE_EARTH, 150, 150);
|
|
for k in 0..3 {
|
|
assert!(
|
|
(fine[k] - coarse[k]).abs() < 5e-4 * fine[k],
|
|
"{}: {fine:?} {coarse:?}",
|
|
case.name
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn quadrature_matches_the_far_field_phase_law() {
|
|
// A Lambertian sphere of albedo rho and radius R at distance D sends a surface facing
|
|
// it (2/3) rho (R/D)^2 Phi(alpha) of the irradiance the Sun gives the sphere, at
|
|
// phase angle alpha, with Phi(alpha) = (sin alpha + (pi - alpha) cos alpha) / pi.
|
|
// Earth's flattening shrinks the disk seen from the equator by b / a. That is the
|
|
// limit for large D. At 1000 Earth radii the quadrature is within 0.008 %, 0.035 %
|
|
// and 0.24 % of it at 0, 60 and 120 degrees; the last is the near-field term, where
|
|
// the lit crescent sits at the limb (0.08 % at ten times the distance). The Moon is
|
|
// black here: at this distance it would add about 1 %.
|
|
let surfaces = Surfaces {
|
|
earth_level: 255,
|
|
moon_albedo: 0.0,
|
|
};
|
|
let far = Vector3::new(1000.0 * EARTH_EQUATORIAL_RADIUS, 0.0, 0.0);
|
|
let (albedo, _) = surfaces.albedo();
|
|
for degrees in [0.0, 60.0, 120.0] {
|
|
let case = Case::new("far", far, -Vector3::x()).sun_at(degrees);
|
|
let alpha = f64::to_radians(degrees);
|
|
let phi = (alpha.sin() + (PI - alpha) * alpha.cos()) / PI;
|
|
// In units of the Sun's irradiance at the target, which is nearer or farther
|
|
// from the Sun than Earth by a few percent.
|
|
let sunlight = ((case.sun - far).norm() / case.sun.norm()).powi(2);
|
|
let got = quadrature(&case, &surfaces, 300, 300);
|
|
for k in 0..3 {
|
|
let want = 2.0 / 3.0 * albedo[k] * 1e-6 * phi * sunlight * EARTH_POLAR_RADIUS
|
|
/ EARTH_EQUATORIAL_RADIUS;
|
|
assert!(
|
|
(got[k] - want).abs() < 3e-3 * want,
|
|
"{degrees} deg, channel {k}: {:e} against {want:e}",
|
|
got[k]
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn earth_texture_is_read_in_earth_fixed_axes() {
|
|
// Earth white over longitudes 0 to 180 degrees (the Earth-fixed +Y side) and black
|
|
// over the rest. 500 km up, a target sees 22 degrees of longitude either side of
|
|
// the point below it, so over 90 E it sees only white and over 90 W only black,
|
|
// until Earth is turned half a revolution.
|
|
let (width, height) = (64, 32);
|
|
let pixels = (0..width * height)
|
|
.flat_map(|i| {
|
|
let v = if i % width >= width / 2 { 255 } else { 0 };
|
|
[v, v, v, 255]
|
|
})
|
|
.collect();
|
|
let texture = EarthTexture::from_rgba8(width, height, pixels).unwrap();
|
|
let east = Vector3::new(0.0, 6_878_000.0, 0.0);
|
|
let over = |name, target: Vector3<f64>, turned: bool| Case {
|
|
earth_attitude: if turned {
|
|
UnitQuaternion::from_axis_angle(&Vector3::z_axis(), PI)
|
|
} else {
|
|
UnitQuaternion::identity()
|
|
},
|
|
sun: AU * target.normalize(),
|
|
..Case::new(name, target, -target)
|
|
};
|
|
let cases = [
|
|
over("texture-east", east, false),
|
|
over("texture-west", -east, false),
|
|
over("texture-east-turned", east, true),
|
|
];
|
|
let white = quadrature(&cases[0], &Surfaces::WHITE_EARTH, 300, 300);
|
|
let _devices = common::device_lock();
|
|
for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
|
|
let options = RendererOptions {
|
|
adapter,
|
|
fidelity: Fidelity {
|
|
lighting: LightingTier::L1,
|
|
..Default::default()
|
|
},
|
|
earth_texture: Some(Arc::new(texture.clone())),
|
|
appearance: Surfaces::WHITE_EARTH.appearance(),
|
|
read_radiance: true,
|
|
..Default::default()
|
|
};
|
|
let mut renderer = Renderer::new(CameraConfig::example(), options).unwrap();
|
|
let name = renderer.adapter_info().name.clone();
|
|
let got = rendered(&mut renderer, &cases);
|
|
for k in 0..3 {
|
|
assert!(
|
|
(got[0][k] - white[k]).abs() < 0.005 * white[k],
|
|
"{name}: east {:?}, all white {white:?}",
|
|
got[0]
|
|
);
|
|
// Black ground, and the Moon's 1e-6 or so.
|
|
assert!(got[1][k] < 1e-5, "{name}: west {:?}", got[1]);
|
|
assert!(got[2][k] < 1e-5, "{name}: east, turned {:?}", got[2]);
|
|
}
|
|
}
|
|
}
|