From 970551b4e1dcb9d2322e88b27a5f466117dbe83c Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Fri, 16 Oct 2026 14:43:50 -0400 Subject: [PATCH] Added stars at L1 --- crates/resolved-cli/src/main.rs | 41 +- crates/resolved-cli/src/scene_file.rs | 18 + crates/resolved/src/fidelity.rs | 8 +- crates/resolved/src/gpu/mod.rs | 1 + crates/resolved/src/gpu/profile.rs | 17 +- crates/resolved/src/gpu/sensor.rs | 60 +- crates/resolved/src/gpu/stars.rs | 238 ++++++++ crates/resolved/src/gpu/targets.rs | 16 + crates/resolved/src/lib.rs | 2 + crates/resolved/src/radiometry.rs | 5 +- crates/resolved/src/render.rs | 59 +- crates/resolved/src/shaders/sensor.wgsl | 10 +- crates/resolved/src/shaders/stars.wgsl | 85 +++ crates/resolved/src/stars.rs | 547 ++++++++++++++++++ crates/resolved/tests/gpu/distortion.rs | 2 +- crates/resolved/tests/gpu/main.rs | 1 + crates/resolved/tests/gpu/sensor.rs | 4 +- crates/resolved/tests/gpu/stars.rs | 261 +++++++++ docs/optics.md | 98 ++++ examples/scenes/geo-telescope-1000km.toml | 37 ++ examples/scenes/geo-telescope-100km-side.toml | 35 ++ examples/scenes/geo-telescope-100km.toml | 36 ++ examples/scenes/stars-pleiades.toml | 24 + readme.md | 27 +- scripts/fetch-stars.py | 399 +++++++++++++ scripts/make-example-scenes.py | 142 ++++- 26 files changed, 2118 insertions(+), 55 deletions(-) create mode 100644 crates/resolved/src/gpu/stars.rs create mode 100644 crates/resolved/src/shaders/stars.wgsl create mode 100644 crates/resolved/src/stars.rs create mode 100644 crates/resolved/tests/gpu/stars.rs create mode 100644 examples/scenes/geo-telescope-1000km.toml create mode 100644 examples/scenes/geo-telescope-100km-side.toml create mode 100644 examples/scenes/geo-telescope-100km.toml create mode 100644 examples/scenes/stars-pleiades.toml create mode 100755 scripts/fetch-stars.py diff --git a/crates/resolved-cli/src/main.rs b/crates/resolved-cli/src/main.rs index 2e00fbe..54f6cc8 100644 --- a/crates/resolved-cli/src/main.rs +++ b/crates/resolved-cli/src/main.rs @@ -11,7 +11,7 @@ use clap::{Args, Parser, Subcommand}; use resolved::gpu::AdapterChoice; use resolved::{ CameraTier, EarthTexture, Fidelity, Frame, LightingTier, Look, MoonTexture, Preview, Renderer, - RendererOptions, + RendererOptions, StarCatalog, StarPhotometry, }; use crate::scene_file::SceneFile; @@ -52,13 +52,15 @@ struct RenderArgs { #[arg(long)] pretty: bool, /// Lighting fidelity tier: L0, or L1 (adds Earth's surface, shadows on targets, - /// light from Earth and the Moon, eclipses on them, and the Sun's limb darkening; - /// needs Blue Marble data). + /// light from Earth and the Moon, eclipses on them, the Sun's limb darkening, and + /// stars, which camera C1 shows; needs Blue Marble data, and with C1 the star + /// catalog). #[arg(long, default_value = "L0")] lighting: LightingTier, /// Camera fidelity tier: C0, or C1 (adds pixels that collect light over their area, /// for graded edges and continuous thin parts, lens distortion, lens blur, and the - /// Sun, Moon and Earth where light time and aberration make them appear). + /// Sun, Moon and Earth where light time and aberration make them appear, and at + /// lighting L1 stars). #[arg(long, default_value = "C0")] camera: CameraTier, /// From camera C1: render each pixel as N x N finer ones (1 to 4; default 2), which @@ -72,7 +74,8 @@ struct RenderArgs { #[arg(long, value_name = "FRAMES")] timings: Option, /// Directory of data the tiers need (earth/bmng-2004-MM.jpg from - /// scripts/fetch-earth.sh, moon/wac-hapke-643.png from scripts/fetch-moon.py). + /// scripts/fetch-earth.sh, moon/wac-hapke-643.png from scripts/fetch-moon.py, + /// stars/stars.bin from scripts/fetch-stars.py). /// Without it: RESOLVED_DATA, then ./data, then the repository's data/. #[arg(long, value_name = "DIR")] data: Option, @@ -124,12 +127,13 @@ fn earth_image(data: Option<&Path>, month: u32) -> Result { data, &format!("earth/bmng-2004-{month:02}.jpg"), "fetch-earth.sh", + "lighting L1", ) } /// `name` from the first data directory that has it; the error names `script`, which -/// makes it. -fn data_file(data: Option<&Path>, name: &str, script: &str) -> Result { +/// makes it, and `tier`, which needs it. +fn data_file(data: Option<&Path>, name: &str, script: &str, tier: &str) -> Result { let tried = data_directories(data); tried .iter() @@ -138,7 +142,7 @@ fn data_file(data: Option<&Path>, name: &str, script: &str) -> Result = tried.iter().map(|d| d.display().to_string()).collect(); format!( - "lighting L1 needs {name} (run scripts/{script}, or pass --data DIR); looked in {}", + "{tier} needs {name} (run scripts/{script}, or pass --data DIR); looked in {}", tried.join(", ") ) }) @@ -171,6 +175,7 @@ fn render(args: &RenderArgs) -> Result<(), String> { args.data.as_deref(), "moon/wac-hapke-643.png", "fetch-moon.py", + "lighting L1", )?; Some(Arc::new( MoonTexture::load(&image).map_err(|e| e.to_string())?, @@ -178,6 +183,22 @@ fn render(args: &RenderArgs) -> Result<(), String> { } else { None }; + let stars = if args.lighting >= LightingTier::L1 && args.camera >= CameraTier::C1 { + let catalog = data_file( + args.data.as_deref(), + "stars/stars.bin", + "fetch-stars.py", + "lighting L1 with camera C1", + )?; + // Only the stars this exposure can show. + let faint = StarPhotometry::new(&camera).faintest_drawn(camera.exposure); + Some(Arc::new( + StarCatalog::load(&catalog, file.years_since_j2000()?, faint) + .map_err(|e| e.to_string())?, + )) + } else { + None + }; let texturing = start.elapsed(); let options = RendererOptions { @@ -199,6 +220,7 @@ fn render(args: &RenderArgs) -> Result<(), String> { shadow_map_size: None, dn_in_place: None, supersampling: args.supersampling, + stars, }; let start = Instant::now(); let mut renderer = Renderer::new(camera, options).map_err(|e| e.to_string())?; @@ -260,12 +282,13 @@ fn render(args: &RenderArgs) -> Result<(), String> { v[v.len() / 2] }; println!( - "gpu shadow {:.3} ms, earthshine {:.3} ms, mesh {:.3} ms, bodies {:.3} ms, blur {:.3} ms, sensor {:.3} ms, first to last {:.3} ms (medians)", + "gpu shadow {:.3} ms, earthshine {:.3} ms, mesh {:.3} ms, bodies {:.3} ms, blur {:.3} ms, stars {:.3} ms, sensor {:.3} ms, first to last {:.3} ms (medians)", median(|t| t.shadow), median(|t| t.earthshine), median(|t| t.mesh), median(|t| t.bodies), median(|t| t.blur), + median(|t| t.stars), median(|t| t.sensor), median(|t| t.total) ); diff --git a/crates/resolved-cli/src/scene_file.rs b/crates/resolved-cli/src/scene_file.rs index 9727099..f5ae6e6 100644 --- a/crates/resolved-cli/src/scene_file.rs +++ b/crates/resolved-cli/src/scene_file.rs @@ -223,6 +223,24 @@ impl SceneFile { Ok(None) } + /// Julian years from J2000.0 to the scene's epoch, for the stars' proper motion; 0 + /// for a scene without one. + #[cfg(feature = "astrokit")] + pub fn years_since_j2000(&self) -> Result { + let Some(epoch) = &self.epoch else { + return Ok(0.0); + }; + let epoch: hifitime::Epoch = epoch.parse().map_err(|e| format!("epoch {epoch:?}: {e}"))?; + let j2000 = + hifitime::Epoch::from_gregorian_hms(2000, 1, 1, 12, 0, 0, hifitime::TimeScale::TT); + Ok((epoch - j2000).to_unit(hifitime::Unit::Day) / 365.25) + } + + #[cfg(not(feature = "astrokit"))] + pub fn years_since_j2000(&self) -> Result { + Ok(0.0) + } + /// The observer's pose. pub fn observer(&self) -> Result { let o = &self.observer; diff --git a/crates/resolved/src/fidelity.rs b/crates/resolved/src/fidelity.rs index 483f10f..614d47c 100644 --- a/crates/resolved/src/fidelity.rs +++ b/crates/resolved/src/fidelity.rs @@ -26,7 +26,10 @@ pub enum LightingTier { /// targets by Earth and the Moon (not shadowed by the targets); Earth and the Moon /// shadowing each other, lunar and solar eclipses; and a Sun brighter at the center /// of its disk than at the edge (limb darkening), in its image and in every - /// eclipse. + /// eclipse; and the stars of a catalog when one is given + /// ([`RendererOptions::stars`](crate::RendererOptions::stars)), seen only by camera + /// tier C1 and up, whose pixels collect light over their area (a C0 pixel samples + /// one point, which never meets a star). L1, } @@ -43,7 +46,8 @@ pub enum CameraTier { /// [`BrownConrady`](crate::BrownConrady)) and its blur (the camera's /// [`Psf`](crate::Psf)). The scene is rendered finer than the pixels, 2 x 2 by /// default, distorted, blurred, then averaged into pixels. Earth, the Moon and the - /// Sun are drawn at their apparent positions ([`Scene::velocities`](crate::Scene::velocities)). Target meshes are covered at several + /// Sun, and at lighting L1 the stars, are drawn at their apparent positions + /// ([`Scene::velocities`](crate::Scene::velocities)). Target meshes are covered at several /// points per rendered pixel, and rendered pixels on the edge of Earth, the Moon or /// the Sun cast a grid of rays, 4 x 4 and more for bodies only a few pixels across. /// More of the camera joins this tier as it is modeled. diff --git a/crates/resolved/src/gpu/mod.rs b/crates/resolved/src/gpu/mod.rs index 1632575..fb39241 100644 --- a/crates/resolved/src/gpu/mod.rs +++ b/crates/resolved/src/gpu/mod.rs @@ -8,6 +8,7 @@ pub mod mesh; pub mod profile; pub mod sensor; pub mod shadow; +pub mod stars; pub mod targets; pub mod warp; diff --git a/crates/resolved/src/gpu/profile.rs b/crates/resolved/src/gpu/profile.rs index f1c3f8e..23c272c 100644 --- a/crates/resolved/src/gpu/profile.rs +++ b/crates/resolved/src/gpu/profile.rs @@ -8,7 +8,15 @@ use super::GpuContext; /// The passes timed, in the order they run. -pub const PASSES: [&str; 6] = ["shadow", "earthshine", "mesh", "bodies", "blur", "sensor"]; +pub const PASSES: [&str; 7] = [ + "shadow", + "earthshine", + "mesh", + "bodies", + "blur", + "stars", + "sensor", +]; /// Where one pass writes its begin and end timestamps. #[derive(Clone, Copy, Debug)] @@ -74,6 +82,8 @@ pub struct PassTimings { /// Lens distortion, the point spread function and averaging finer rendering into /// pixels (zero work below camera C1). pub blur: f64, + /// Stars (zero work below camera C1 or without a catalog). + pub stars: f64, /// Radiance to DN. pub sensor: f64, /// From the start of the first pass to the end of the last, gaps included. @@ -155,8 +165,9 @@ impl Profiler { mesh: ms(ticks[4], ticks[5]), bodies: ms(ticks[6], ticks[7]), blur: ms(ticks[8], ticks[9]), - sensor: ms(ticks[10], ticks[11]), - total: ms(ticks[0], ticks[11]), + stars: ms(ticks[10], ticks[11]), + sensor: ms(ticks[12], ticks[13]), + total: ms(ticks[0], ticks[13]), } } } diff --git a/crates/resolved/src/gpu/sensor.rs b/crates/resolved/src/gpu/sensor.rs index b035367..72ae88e 100644 --- a/crates/resolved/src/gpu/sensor.rs +++ b/crates/resolved/src/gpu/sensor.rs @@ -55,17 +55,23 @@ const WORKGROUP: u32 = 256; pub struct SensorPass { pipeline: wgpu::ComputePipeline, uniforms: wgpu::Buffer, + stars: bool, + /// Bound for the stars' electrons without stars: the shader names the binding even + /// where it never reads it. + no_stars: wgpu::Buffer, } impl SensorPass { - /// Compiles the shader and allocates the uniform buffer. - pub fn new(ctx: &GpuContext) -> Self { + /// Compiles the shader, adding the stars' electrons if `stars`, and allocates the + /// uniform buffer. + pub fn new(ctx: &GpuContext, stars: bool) -> Self { let module = ctx .device .create_shader_module(wgpu::ShaderModuleDescriptor { label: Some("sensor"), source: wgpu::ShaderSource::Wgsl( [ + &format!("const STARS: bool = {stars};\n"), include_str!("../shaders/distortion.wgsl"), include_str!("../shaders/sensor.wgsl"), ] @@ -89,7 +95,17 @@ impl SensorPass { usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST, mapped_at_creation: false, }); - Self { pipeline, uniforms } + Self { + pipeline, + uniforms, + stars, + no_stars: ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("no stars"), + size: 4, + usage: wgpu::BufferUsages::STORAGE, + mapped_at_creation: false, + }), + } } /// Uploads `uniforms` and records the pass into `encoder`, reading `targets`' @@ -105,23 +121,33 @@ impl SensorPass { ctx.queue .write_buffer(&self.uniforms, 0, bytemuck::bytes_of(uniforms)); let radiance = targets.sensed().create_view(&Default::default()); + let stars = match (self.stars, &targets.stars) { + (true, Some(stars)) => stars, + (true, None) => panic!("targets allocated without stars"), + (false, _) => &self.no_stars, + }; + let entries = [ + wgpu::BindGroupEntry { + binding: 0, + resource: wgpu::BindingResource::TextureView(&radiance), + }, + wgpu::BindGroupEntry { + binding: 1, + resource: targets.dn.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 2, + resource: self.uniforms.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 3, + resource: stars.as_entire_binding(), + }, + ]; let bind_group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { label: Some("sensor"), layout: &self.pipeline.get_bind_group_layout(0), - entries: &[ - wgpu::BindGroupEntry { - binding: 0, - resource: wgpu::BindingResource::TextureView(&radiance), - }, - wgpu::BindGroupEntry { - binding: 1, - resource: targets.dn.as_entire_binding(), - }, - wgpu::BindGroupEntry { - binding: 2, - resource: self.uniforms.as_entire_binding(), - }, - ], + entries: &entries, }); // One invocation per u32 word, in rows of at most 65,535 workgroups. let groups = (targets.dn.size() / 4).div_ceil(WORKGROUP as u64) as u32; diff --git a/crates/resolved/src/gpu/stars.rs b/crates/resolved/src/gpu/stars.rs new file mode 100644 index 0000000..a7666ac --- /dev/null +++ b/crates/resolved/src/gpu/stars.rs @@ -0,0 +1,238 @@ +//! The pass that draws stars, and choosing which. +//! +//! Each frame [`draws`] finds the catalog's stars in the field on the CPU, puts each at +//! its apparent position (aberrated by the observer's motion, [`crate::scene::aberrate`]) +//! through the lens's distortion, and works out its electrons. The pass then spreads +//! them by the PSF over the pixels around each star, as electrons the sensor pass adds +//! before the full well (`stars.wgsl`). Stars are drawn this way, not rendered and +//! blurred, because a point of light rendered into a pixel loses where in the pixel it +//! fell (see [`super::blur`]); here the blurred spot is exact for any position. + +use nalgebra::{Vector2, Vector3}; + +use super::profile::PassTimestamps; +use super::targets::FrameTargets; +use super::GpuContext; +use crate::camera::{CameraConfig, Intrinsics, Psf}; +use crate::scene::{aberrate, Pose}; +use crate::stars::{StarCatalog, StarPhotometry}; + +/// Stars giving fewer electrons than this in every channel are not drawn: a tenth of +/// an electron, far below any detector's noise. In a wide field the many fainter stars +/// would cost more to place than they add. +pub const FAINTEST: f64 = 0.1; + +/// One star for the pass; the layout matches `StarDraw` in `stars.wgsl`. +#[repr(C)] +#[derive(Clone, Copy, Debug, PartialEq, bytemuck::Pod, bytemuck::Zeroable)] +pub struct StarDraw { + /// Position in image pixels, then in the rendered images' pixels. + pub at: [f32; 4], + /// Electrons per channel. + pub electrons: [f32; 4], +} + +/// The shader's uniform block; the layout matches `Params` in `stars.wgsl`. +#[repr(C)] +#[derive(Clone, Copy, Debug, PartialEq, bytemuck::Pod, bytemuck::Zeroable)] +struct Params { + sigma: [f32; 4], + dims: [u32; 4], + count: [u32; 4], +} + +/// The stars of `catalog` that appear to a camera at `camera` and give at least +/// [`FAINTEST`] electrons in some channel. +/// +/// `image` projects into the image (with distortion), `render` into the rendered +/// images (pinhole, for the depth test), `beta` is the observer's velocity relative to +/// the stars over c, in the inertial frame, and `reach` how far the PSF reaches in +/// pixels. +#[allow(clippy::too_many_arguments)] +pub fn draws( + catalog: &StarCatalog, + photometry: &StarPhotometry, + exposure: f64, + camera: &Pose, + image: &Intrinsics, + render: &Intrinsics, + beta: &Vector3, + reach: u32, +) -> Vec { + let to_camera = camera.attitude.inverse(); + let pad = reach as f64 + 1.0; + let (w, h) = (image.width as f64, image.height as f64); + // The field: as far from the boresight as the corners, with the PSF's reach, and + // with room for aberration, at most 1e-4 rad. + let radius = [ + (-pad, -pad), + (w + pad, -pad), + (-pad, h + pad), + (w + pad, h + pad), + ] + .iter() + .map(|&(u, v)| image.ray(&Vector2::new(u, v)).z.clamp(-1.0, 1.0).acos()) + .fold(0.0f64, f64::max) + + 2e-4; + let boresight = camera.attitude * Vector3::z(); + // No star fainter than this can reach FAINTEST: skipped before any work. + let faint = photometry.faintest_drawn(exposure) as f32; + catalog + .within(&boresight, radius.min(std::f64::consts::PI)) + .filter(|star| star.v <= faint) + .filter_map(|star| { + let d = to_camera * aberrate(&star.direction.cast::(), beta); + if d.z <= 0.0 { + return None; + } + let uv = image.project(&d)?; + if uv.x < -pad || uv.y < -pad || uv.x > w + pad || uv.y > h + pad { + return None; + } + // A point source's light is cut only by the aperture's foreshortening. + let collected = photometry.aperture * exposure * d.z / d.norm(); + let rates = photometry.rates(star.v, star.bv); + let electrons: [f32; 4] = std::array::from_fn(|k| (rates[k] * collected) as f32); + if electrons.iter().all(|&e| (e as f64) < FAINTEST) { + return None; + } + let r = render.project(&d)?; + Some(StarDraw { + at: [uv.x as f32, uv.y as f32, r.x as f32, r.y as f32], + electrons, + }) + }) + .collect() +} + +/// The compiled star pipeline, its uniforms and a buffer of stars that grows as needed. +#[derive(Debug)] +pub struct StarPass { + pipeline: wgpu::ComputePipeline, + params: wgpu::Buffer, + sigma: [f32; 4], + reach: u32, + stars: wgpu::Buffer, +} + +impl StarPass { + /// Compiles the shader for `camera`'s PSF. + pub fn new(ctx: &GpuContext, camera: &CameraConfig) -> Self { + let module = ctx + .device + .create_shader_module(wgpu::ShaderModuleDescriptor { + label: Some("stars"), + source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/stars.wgsl").into()), + }); + let pipeline = ctx + .device + .create_compute_pipeline(&wgpu::ComputePipelineDescriptor { + label: Some("stars"), + layout: None, + module: &module, + entry_point: Some("main"), + compilation_options: Default::default(), + cache: None, + }); + let sigma = match &camera.psf { + Psf::None => [0.0; 4], + Psf::Gaussian { sigma } => { + std::array::from_fn(|k| sigma.get(k).copied().unwrap_or(0.0) as f32) + } + }; + Self { + pipeline, + params: ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("star params"), + size: std::mem::size_of::() as u64, + usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }), + sigma, + reach: camera.psf.radius(1), + stars: Self::star_buffer(ctx, 1024), + } + } + + /// How far the PSF reaches, in image pixels. + pub fn reach(&self) -> u32 { + self.reach + } + + fn star_buffer(ctx: &GpuContext, capacity: u64) -> wgpu::Buffer { + ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("stars"), + size: capacity * std::mem::size_of::() as u64, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }) + } + + /// Uploads `stars` and records clearing `targets`' star electrons and drawing them. + pub fn encode( + &mut self, + ctx: &GpuContext, + encoder: &mut wgpu::CommandEncoder, + targets: &FrameTargets, + stars: &[StarDraw], + timestamps: Option, + ) { + let electrons = targets + .stars + .as_ref() + .expect("targets allocated with stars"); + encoder.clear_buffer(electrons, 0, None); + let needed = stars.len().max(1) as u64 * std::mem::size_of::() as u64; + if needed > self.stars.size() { + self.stars = Self::star_buffer(ctx, (stars.len() as u64).next_power_of_two()); + } + if !stars.is_empty() { + ctx.queue + .write_buffer(&self.stars, 0, bytemuck::cast_slice(stars)); + } + let params = Params { + sigma: self.sigma, + dims: [targets.width, targets.height, targets.channels, self.reach], + count: [stars.len() as u32, 0, 0, 0], + }; + ctx.queue + .write_buffer(&self.params, 0, bytemuck::bytes_of(¶ms)); + let bind_group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("stars"), + layout: &self.pipeline.get_bind_group_layout(0), + entries: &[ + wgpu::BindGroupEntry { + binding: 0, + resource: self.stars.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 1, + resource: self.params.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 2, + resource: wgpu::BindingResource::TextureView(&targets.depth.create_view( + &wgpu::TextureViewDescriptor { + aspect: wgpu::TextureAspect::DepthOnly, + ..Default::default() + }, + )), + }, + wgpu::BindGroupEntry { + binding: 3, + resource: electrons.as_entire_binding(), + }, + ], + }); + let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("stars"), + timestamp_writes: timestamps.map(PassTimestamps::compute), + }); + if !stars.is_empty() { + let groups = stars.len() as u32; + pass.set_pipeline(&self.pipeline); + pass.set_bind_group(0, &bind_group, &[]); + pass.dispatch_workgroups(groups.min(65_535), groups.div_ceil(65_535), 1); + } + } +} diff --git a/crates/resolved/src/gpu/targets.rs b/crates/resolved/src/gpu/targets.rs index 8f8dfc5..c896fe5 100644 --- a/crates/resolved/src/gpu/targets.rs +++ b/crates/resolved/src/gpu/targets.rs @@ -30,6 +30,8 @@ pub struct TargetLayout { /// With lens distortion, the size of the undistorted render that the focal plane is /// read from; `None` without, where the render is the focal plane. pub undistorted: Option<[u32; 2]>, + /// Whether stars are drawn, which needs a buffer of their electrons. + pub stars: bool, } impl TargetLayout { @@ -44,6 +46,7 @@ impl TargetLayout { scale: 1, margin: 0, undistorted: None, + stars: false, } } @@ -100,6 +103,9 @@ pub struct FrameTargets { pub blurred: Option<[wgpu::Texture; 2]>, /// Number of camera channels. pub channels: u32, + /// With stars: their electrons per pixel and channel, channels interleaved per pixel, + /// in fixed point, 1/64 of an electron per unit, which the sensor pass adds. + pub stars: Option, /// The sensor pass's output: u16 DN, channels interleaved per pixel, rows from the /// top, packed two to a u32 word. Mappable where [`FrameTargets::dn_mappable`]. pub dn: wgpu::Buffer, @@ -186,6 +192,16 @@ impl FrameTargets { ] }), channels, + stars: layout.stars.then(|| { + ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("star electrons"), + size: values * 4, + usage: wgpu::BufferUsages::STORAGE + | wgpu::BufferUsages::COPY_DST + | wgpu::BufferUsages::COPY_SRC, + mapped_at_creation: false, + }) + }), dn: ctx.device.create_buffer(&wgpu::BufferDescriptor { label: Some("dn"), size: values.div_ceil(2) * 4, diff --git a/crates/resolved/src/lib.rs b/crates/resolved/src/lib.rs index ddf04c7..5077726 100644 --- a/crates/resolved/src/lib.rs +++ b/crates/resolved/src/lib.rs @@ -35,6 +35,7 @@ pub mod scene; pub mod shadow; mod solar_data; pub mod spectrum; +pub mod stars; pub mod sun; pub use camera::{BrownConrady, CameraConfig, Channel, Detector, Intrinsics, Psf}; @@ -48,3 +49,4 @@ pub use radiometry::Radiometry; pub use render::{Frame, FrameView, Pending, Renderer, RendererOptions}; pub use scene::{CameraRelative, Placement, Pose, Scene, Target, Velocities}; pub use spectrum::Spectrum; +pub use stars::{StarCatalog, StarPhotometry}; diff --git a/crates/resolved/src/radiometry.rs b/crates/resolved/src/radiometry.rs index 2242340..fb6cad1 100644 --- a/crates/resolved/src/radiometry.rs +++ b/crates/resolved/src/radiometry.rs @@ -57,7 +57,10 @@ pub fn integrate(spectra: &[&Spectrum]) -> Option { integrate_weighted(spectra, |_| 1.0) } -fn integrate_weighted(spectra: &[&Spectrum], weight: impl Fn(f64) -> f64) -> Option { +pub(crate) fn integrate_weighted( + spectra: &[&Spectrum], + weight: impl Fn(f64) -> f64, +) -> Option { let (lo, hi) = spectra .iter() .filter_map(|s| s.support()) diff --git a/crates/resolved/src/render.rs b/crates/resolved/src/render.rs index bfcc589..be6969e 100644 --- a/crates/resolved/src/render.rs +++ b/crates/resolved/src/render.rs @@ -25,14 +25,17 @@ use crate::gpu::mesh::{GpuMesh, MeshDraw, MeshPass, ShadowDraw}; use crate::gpu::profile::{PassTimings, Profiler}; use crate::gpu::sensor::{SensorPass, SensorUniforms}; use crate::gpu::shadow::{CasterDraw, ShadowPass}; +use crate::gpu::stars::{draws as star_draws, StarPass}; use crate::gpu::targets::{FrameTargets, TargetLayout, DEPTH_FORMAT, RADIANCE_FORMAT}; use crate::gpu::warp::{WarpLayout, WarpUniforms}; use crate::gpu::{AdapterChoice, GpuContext}; use crate::model::{Mesh, ModelId}; use crate::moon::MoonTexture; use crate::radiometry::Radiometry; +use crate::scene::SPEED_OF_LIGHT; use crate::scene::{Pose, Scene}; use crate::shadow::{self, Bounded}; +use crate::stars::{StarCatalog, StarPhotometry}; use crate::sun::LimbDarkening; use crate::Error; @@ -72,6 +75,9 @@ pub struct RendererOptions { /// each pixel. `None` for the tier's choice: 2 at C1. 1 to 4; each step multiplies /// rendering's time and memory by about n^2. Must be 1 (or `None`) at C0. pub supersampling: Option, + /// The stars, drawn at lighting tier L1 and up when given, by camera tier C1 and up + /// (see [`StarCatalog::load`]); ignored below them. + pub stars: Option>, } /// One rendered image. @@ -125,6 +131,9 @@ pub struct Renderer { shadows: ShadowPass, earthshine: EarthshinePass, sensor: SensorPass, + /// From lighting L1 and camera C1 with a catalog: the stars, their photometry and + /// their pass. + stars: Option<(Arc, StarPhotometry, StarPass)>, sensor_uniforms: SensorUniforms, models: Vec, /// The buffer the DN image is read from: the sensor pass's own output, or a staging @@ -193,7 +202,14 @@ impl Renderer { [camera.width, camera.height], warp.as_ref(), ); + // Stars are scene content (lighting L1) that only a camera collecting light over + // each pixel's area (C1) can see: a single point per pixel never meets one. + let stars = options + .stars + .clone() + .filter(|_| area && options.fidelity.lighting >= LightingTier::L1); let layout = TargetLayout { + stars: stars.is_some(), samples, scale, margin, @@ -296,7 +312,14 @@ impl Renderer { blur, shadows: ShadowPass::new(&ctx, options.shadow_map_size.unwrap_or(4096)), earthshine: EarthshinePass::new(&ctx), - sensor: SensorPass::new(&ctx), + sensor: SensorPass::new(&ctx, stars.is_some()), + stars: stars.map(|catalog| { + ( + catalog, + StarPhotometry::new(&camera), + StarPass::new(&ctx, &camera), + ) + }), sensor_uniforms: SensorUniforms::new(&camera, &radiometry, warp.is_some()), distorted: warp.is_some(), render_intrinsics, @@ -540,20 +563,42 @@ impl Renderer { ); self.bodies .encode(ctx, &mut encoder, &self.targets, &bodies, stamps(3)); + // A pass with nothing to do is still opened, empty, to write its timestamps. + let empty = |encoder: &mut wgpu::CommandEncoder, label, index| { + drop(encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some(label), + timestamp_writes: stamps(index).map(|t| t.compute()), + })) + }; match &self.blur { Some(blur) => blur.encode(ctx, &mut encoder, &self.targets, stamps(4)), - // Nothing to blur: an empty pass still writes the timestamps. - None => drop(encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { - label: Some("blur (none)"), - timestamp_writes: stamps(4).map(|t| t.compute()), - })), + None => empty(&mut encoder, "blur (none)", 4), + } + match &mut self.stars { + Some((catalog, photometry, pass)) => { + // Stars are at rest relative to the Sun, to 12 m/s. + let beta = (scene.velocities.earth + scene.velocities.observer) / SPEED_OF_LIGHT; + let draws = star_draws( + catalog, + photometry, + self.camera.exposure, + &camera, + &self.camera.intrinsics(), + &self.render_intrinsics, + &beta, + pass.reach(), + ); + let timestamps = self.profiler.as_ref().map(|p| p.pass(5)); + pass.encode(ctx, &mut encoder, &self.targets, &draws, timestamps); + } + None => empty(&mut encoder, "stars (none)", 5), } self.sensor.encode( ctx, &mut encoder, &self.targets, &self.sensor_uniforms, - stamps(5), + stamps(6), ); if let Some(profiler) = &self.profiler { profiler.resolve(&mut encoder); diff --git a/crates/resolved/src/shaders/sensor.wgsl b/crates/resolved/src/shaders/sensor.wgsl index 7a18fc7..3c59274 100644 --- a/crates/resolved/src/shaders/sensor.wgsl +++ b/crates/resolved/src/shaders/sensor.wgsl @@ -18,6 +18,9 @@ struct Sensor { @group(0) @binding(0) var radiance: texture_2d; @group(0) @binding(1) var dn: array; @group(0) @binding(2) var s: Sensor; +// With STARS (declared by the code that includes this file): the stars' electrons per +// value, 1/64 of an electron per unit (stars.wgsl). +@group(0) @binding(3) var star_electrons: array; fn value(i: u32) -> u32 { let n = s.dims.z; @@ -41,7 +44,12 @@ fn value(i: u32) -> u32 { } let ray = vec3(sight, 1.0); let c2 = 1.0 / dot(ray, ray); - let electrons = clamp(l * s.detector.x * c2 * c2 / stretch, 0.0, s.detector.y); + var collected = l * s.detector.x * c2 * c2 / stretch; + if STARS { + // Light from the stars adds before the pixel fills. + collected += f32(star_electrons[i]) / 64.0; + } + let electrons = clamp(collected, 0.0, s.detector.y); let level = floor(electrons / s.detector.z + s.detector.w + 0.5); return u32(clamp(level, 0.0, f32(s.dims.w))); } diff --git a/crates/resolved/src/shaders/stars.wgsl b/crates/resolved/src/shaders/stars.wgsl new file mode 100644 index 0000000..3b7cf1c --- /dev/null +++ b/crates/resolved/src/shaders/stars.wgsl @@ -0,0 +1,85 @@ +// Stars: each star's electrons spread over the pixels around it by the point spread +// function, the Gaussian integrated over each pixel's area, which is exact for any +// position within a pixel. Stars lie at infinity, so anything drawn where a star is +// hides it. +// +// One workgroup per star. Electrons accumulate in fixed point, 1/SCALE of an electron, +// in a u32 per pixel and channel, channels interleaved per pixel like the DN image; WGSL +// has no floating-point atomics. + +struct StarDraw { + // xy: position in image pixels. zw: position in the rendered images' pixels, for + // the depth test. + at: vec4, + // Per channel, electrons. + electrons: vec4, +} + +struct Params { + // Per channel, the PSF's sigma in pixels; 0 for none. + sigma: vec4, + // Image width, height, channels, and how far the PSF reaches in pixels. + dims: vec4, + // x: number of stars. + count: vec4, +} + +const SCALE: f32 = 64.0; + +@group(0) @binding(0) var stars: array; +@group(0) @binding(1) var p: Params; +@group(0) @binding(2) var depth: texture_depth_2d; +@group(0) @binding(3) var electrons: array>; + +// The error function, to 1.5e-7 (Abramowitz and Stegun 7.1.26). +fn erf(x: f32) -> f32 { + let t = 1.0 / (1.0 + 0.3275911 * abs(x)); + let poly = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741 + t * (-1.453152027 + t * 1.061405429)))); + return sign(x) * (1.0 - poly * exp(-x * x)); +} + +// The fraction of a Gaussian of `sigma` centered at `c` that falls in [a, a + 1), or +// with no blur, 1 for the pixel holding `c`. +fn share(a: f32, c: f32, sigma: f32) -> f32 { + if sigma <= 0.0 { + return select(0.0, 1.0, c >= a && c < a + 1.0); + } + let s = 1.0 / (sigma * 1.41421356); + return 0.5 * (erf((a + 1.0 - c) * s) - erf((a - c) * s)); +} + +@compute @workgroup_size(64) +fn main(@builtin(workgroup_id) group: vec3, @builtin(local_invocation_index) i: u32) { + // Workgroups in rows of at most 65,535. + let index = group.x + group.y * 65535u; + if index >= p.count.x { + return; + } + let star = stars[index]; + let size = textureDimensions(depth); + let render = vec2(floor(star.at.zw)); + if all(render >= vec2(0)) && all(render < vec2(size)) { + if textureLoad(depth, render, 0) > 0.0 { + return; + } + } + let reach = i32(p.dims.w); + let side = 2 * reach + 1; + let center = vec2(floor(star.at.xy)); + for (var j = i32(i); j < side * side; j += 64) { + let pixel = center + vec2(j % side - reach, j / side - reach); + if any(pixel < vec2(0)) || pixel.x >= i32(p.dims.x) || pixel.y >= i32(p.dims.y) { + continue; + } + let corner = vec2(pixel); + let base = (u32(pixel.y) * p.dims.x + u32(pixel.x)) * p.dims.z; + for (var k = 0u; k < p.dims.z; k++) { + let sigma = p.sigma[k]; + let value = star.electrons[k] * share(corner.x, star.at.x, sigma) * share(corner.y, star.at.y, sigma); + let fixed = u32(value * SCALE + 0.5); + if fixed > 0u { + atomicAdd(&electrons[base + k], fixed); + } + } + } +} diff --git a/crates/resolved/src/stars.rs b/crates/resolved/src/stars.rs new file mode 100644 index 0000000..b9534c1 --- /dev/null +++ b/crates/resolved/src/stars.rs @@ -0,0 +1,547 @@ +//! Stars: the catalog, finding the stars in a field, and how many electrons each gives. +//! +//! The catalog is the merge of Hipparcos, Tycho-2 and Gaia DR3 that +//! `scripts/fetch-stars.py` writes: 37 million stars to about G 15 by default, each with +//! an ICRS position at epoch J2000, its proper motion, Johnson V and B - V. [`StarCatalog::load`] moves every star by its +//! proper motion to one epoch and indexes them in declination bands sorted by right +//! ascension, so the stars in a field are found without looking at the rest. +//! +//! A star's light: its B - V gives a temperature (Ballesteros 2012), and a blackbody at +//! that temperature, scaled so that its flux through the Johnson V passband matches the +//! star's V magnitude, is folded into the camera's channels ([`StarPhotometry`]). + +use std::path::Path; + +use nalgebra::Vector3; + +use crate::camera::CameraConfig; +use crate::radiometry::{PLANCK, SPEED_OF_LIGHT}; +use crate::spectrum::Spectrum; +use crate::Error; + +/// One star as the catalog file stores it. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct CatalogStar { + /// Right ascension, radians, ICRS, epoch J2000. + pub ra: f64, + /// Declination, radians, ICRS, epoch J2000. + pub dec: f64, + /// Proper motion in right ascension times cos(declination), mas per year. + pub pm_ra: f32, + /// Proper motion in declination, mas per year. + pub pm_dec: f32, + /// Johnson V magnitude. + pub v: f32, + /// Johnson B - V, NaN if unknown. + pub bv: f32, +} + +/// A star at the catalog's epoch. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Star { + /// Unit vector toward the star in the inertial frame (ICRS, which this crate's GCRF + /// matches to 0.1 urad). Single precision, good to 0.06 urad, keeps 37 million stars + /// in 0.9 GB. + pub direction: Vector3, + /// Johnson V magnitude. + pub v: f32, + /// Johnson B - V; [`DEFAULT_BV`] where the catalog has none. + pub bv: f32, +} + +/// The B - V a star without a measured color is given: about the Sun's, near the median +/// of the catalog's faint stars. 1,233 Tycho-2 stars have no BT magnitude. +pub const DEFAULT_BV: f32 = 0.6; + +/// Declination bands of the index, in radians. +const BAND: f64 = 0.25 * std::f64::consts::PI / 180.0; + +/// The stars, at one epoch, indexed for finding those in a field. +#[derive(Clone)] +pub struct StarCatalog { + /// Grouped by declination band, sorted by right ascension within each. + stars: Vec, + /// The right ascension of each star, in the same order. + ra: Vec, + /// Where each band starts in `stars`; one more entry than bands. + bands: Vec, +} + +impl std::fmt::Debug for StarCatalog { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("StarCatalog") + .field("stars", &self.stars.len()) + .finish_non_exhaustive() + } +} + +impl StarCatalog { + /// Reads a star file from `scripts/fetch-stars.py`, keeping the stars of V at most + /// `faint`, at `epoch`, in Julian years from J2000.0 (2026.0 is 26.0). + pub fn load(path: &Path, epoch: f64, faint: f64) -> Result { + use std::io::Read; + let fail = |what: String| Error::MissingInput(format!("{}: {what}", path.display())); + let file = std::fs::File::open(path).map_err(|e| fail(e.to_string()))?; + let mut reader = std::io::BufReader::with_capacity(1 << 20, file); + let mut header = [0u8; 16]; + reader + .read_exact(&mut header) + .map_err(|e| fail(e.to_string()))?; + if &header[..8] != b"RSTARS02" { + return Err(fail( + "not a star file from scripts/fetch-stars.py (run it again)".into(), + )); + } + let count = u32::from_le_bytes(header[8..12].try_into().expect("4 bytes")) as usize; + let at = |r: &[u8], i: usize| -> [u8; 4] { r[i..i + 4].try_into().expect("4 bytes") }; + let half = |r: &[u8], i: usize| -> [u8; 2] { r[i..i + 2].try_into().expect("2 bytes") }; + let (turn, half_turn) = ( + std::f64::consts::TAU / 2f64.powi(32), + std::f64::consts::PI / 2f64.powi(32), + ); + // Read in blocks, brightest first, stopping at the first star past the limit, + // straight into the index, so only the stars kept are ever in memory. + let mut block = vec![0u8; 16 * 65_536]; + let (mut left, mut filled, mut next) = (count, 0, 0); + let mut error = None; + let stars = std::iter::from_fn(|| loop { + if next < filled { + let r = &block[16 * next..16 * next + 16]; + next += 1; + let v = u16::from_le_bytes(half(r, 12)) as f32 / 1000.0 - 5.0; + if v as f64 > faint { + return None; + } + let bv = i16::from_le_bytes(half(r, 14)); + return Some(CatalogStar { + ra: u32::from_le_bytes(at(r, 0)) as f64 * turn, + dec: i32::from_le_bytes(at(r, 4)) as f64 * half_turn, + pm_ra: i16::from_le_bytes(half(r, 8)) as f32, + pm_dec: i16::from_le_bytes(half(r, 10)) as f32, + v, + bv: if bv == i16::MIN { + f32::NAN + } else { + bv as f32 / 1000.0 + }, + }); + } + if left == 0 { + return None; + } + let n = left.min(65_536); + if let Err(e) = reader.read_exact(&mut block[..16 * n]) { + error = Some(fail(format!("{count} stars but the file ends early: {e}"))); + return None; + } + (left, filled, next) = (left - n, n, 0); + }); + let catalog = Self::new(stars, epoch); + match error { + Some(e) => Err(e), + None => Ok(catalog), + } + } + + /// The catalog of `stars`, moved by their proper motion to `epoch`, in Julian years + /// from J2000.0. + pub fn new(stars: impl IntoIterator, epoch: f64) -> Self { + let bands = (std::f64::consts::PI / BAND).ceil() as usize; + let band = + |dec: f64| (((dec + std::f64::consts::FRAC_PI_2) / BAND) as usize).min(bands - 1); + let mut moved: Vec<(u32, f32, Star)> = stars + .into_iter() + .map(|s| { + let direction = moved_direction(&s, epoch); + let ra = direction + .y + .atan2(direction.x) + .rem_euclid(std::f64::consts::TAU); + let dec = direction.z.clamp(-1.0, 1.0).asin(); + let bv = if s.bv.is_nan() { DEFAULT_BV } else { s.bv }; + ( + band(dec) as u32, + ra as f32, + Star { + direction: direction.cast::(), + v: s.v, + bv, + }, + ) + }) + .collect(); + moved.sort_unstable_by(|a, b| a.0.cmp(&b.0).then(a.1.total_cmp(&b.1))); + let mut starts = vec![0; bands + 1]; + for (b, _, _) in &moved { + starts[*b as usize + 1] += 1; + } + for i in 0..bands { + starts[i + 1] += starts[i]; + } + Self { + ra: moved.iter().map(|m| m.1).collect(), + stars: moved.into_iter().map(|m| m.2).collect(), + bands: starts, + } + } + + /// Number of stars. + pub fn len(&self) -> usize { + self.stars.len() + } + + /// Whether there are no stars. + pub fn is_empty(&self) -> bool { + self.stars.is_empty() + } + + /// Every star, in the index's order. + pub fn stars(&self) -> &[Star] { + &self.stars + } + + /// The stars within `radius` radians of the unit vector `center`. + pub fn within(&self, center: &Vector3, radius: f64) -> impl Iterator { + use std::f64::consts::{FRAC_PI_2, PI, TAU}; + let dec = center.z.clamp(-1.0, 1.0).asin(); + let ra = center.y.atan2(center.x).rem_euclid(TAU); + let bands = self.bands.len() - 1; + let first = (((dec - radius + FRAC_PI_2) / BAND).floor().max(0.0) as usize).min(bands - 1); + let last = (((dec + radius + FRAC_PI_2) / BAND).floor().max(0.0) as usize).min(bands - 1); + let cos_radius = radius.cos(); + let mut ranges = Vec::new(); + for b in first..=last { + let (start, end) = (self.bands[b], self.bands[b + 1]); + // The band's declination farthest from the equator sets how wide in right + // ascension the field is there. + let edge = (b as f64 * BAND - FRAC_PI_2) + .abs() + .max(((b + 1) as f64 * BAND - FRAC_PI_2).abs()); + let width = if edge + radius >= FRAC_PI_2 - 1e-9 { + PI + } else { + // With room for right ascension held in single precision. + (radius.sin() / edge.cos()).min(1.0).asin() + 1e-6 + }; + if width >= PI { + ranges.push(start..end); + continue; + } + let search = |from: f64, to: f64| { + let slice = &self.ra[start..end]; + start + slice.partition_point(|&r| (r as f64) < from) + ..start + slice.partition_point(|&r| (r as f64) <= to) + }; + let (lo, hi) = (ra - width, ra + width); + if lo < 0.0 { + ranges.push(search(lo + TAU, TAU)); + ranges.push(search(0.0, hi)); + } else if hi > TAU { + ranges.push(search(lo, TAU)); + ranges.push(search(0.0, hi - TAU)); + } else { + ranges.push(search(lo, hi)); + } + } + let center = *center; + ranges + .into_iter() + .flat_map(move |r| self.stars[r].iter()) + .filter(move |s| s.direction.cast::().dot(¢er) >= cos_radius) + } +} + +/// The direction of `star` after `years` of its proper motion, moved along the sphere. +fn moved_direction(star: &CatalogStar, years: f64) -> Vector3 { + let (sa, ca) = star.ra.sin_cos(); + let (sd, cd) = star.dec.sin_cos(); + let p = Vector3::new(cd * ca, cd * sa, sd); + let east = Vector3::new(-sa, ca, 0.0); + let north = Vector3::new(-sd * ca, -sd * sa, cd); + let scale = (1.0 / 3.6e6f64).to_radians() * years; + (p + (east * star.pm_ra as f64 + north * star.pm_dec as f64) * scale).normalize() +} + +/// A star's effective temperature from its B - V, in kelvin: Ballesteros (2012, EPL 97, +/// 34008), from a blackbody fit to main-sequence stars. B - V is clamped to -0.4 to 2.0. +pub fn temperature(bv: f32) -> f64 { + let bv = (bv as f64).clamp(-0.4, 2.0); + 4600.0 * (1.0 / (0.92 * bv + 1.7) + 1.0 / (0.92 * bv + 0.62)) +} + +/// The Johnson V passband (Bessell 1990, PASP 102, 1181), 470 to 700 nm. +fn v_passband() -> Spectrum { + let response = [ + 0.0, 0.030, 0.163, 0.458, 0.780, 0.967, 1.0, 0.973, 0.898, 0.792, 0.684, 0.574, 0.461, + 0.359, 0.270, 0.197, 0.135, 0.081, 0.045, 0.025, 0.017, 0.013, 0.009, 0.0, + ]; + let wavelength = (0..response.len()) + .map(|i| (470.0 + 10.0 * i as f64) * 1e-9) + .collect(); + Spectrum::sampled(wavelength, response.to_vec()).expect("a valid table") +} + +/// The mean spectral irradiance through the V passband of a V = 0 star, in W/m^2 per +/// meter of wavelength: 3.631e-11 W/m^2/nm (Bessell, Castelli and Plez 1998, from Vega). +pub const V_ZERO_POINT: f64 = 3.631e-2; + +/// Planck's law, spectral radiance at wavelength `l` meters and temperature `t` kelvin, +/// up to a constant factor. +fn blackbody(l: f64, t: f64) -> f64 { + const BOLTZMANN: f64 = 1.380_649e-23; + let x = PLANCK * SPEED_OF_LIGHT / (l * BOLTZMANN * t); + 1.0 / (l.powi(5) * x.exp_m1()) +} + +/// Temperatures of the photometry table, log-spaced, kelvin. +const T_MIN: f64 = 2000.0; +const T_MAX: f64 = 50_000.0; +const T_STEPS: usize = 256; + +/// How many electrons a star gives in each of a camera's channels. +#[derive(Clone, Debug, PartialEq)] +pub struct StarPhotometry { + /// Per channel, per temperature step: the electron rate per square meter of + /// aperture from a V = 0 star on the boresight, through the optics' transmission + /// and the channel's QE, in electrons / s / m^2. + table: Vec>, + /// The aperture's area, in m^2. + pub aperture: f64, +} + +impl StarPhotometry { + /// The tables for `camera`. + pub fn new(camera: &CameraConfig) -> Self { + let passband = v_passband(); + let band_width = crate::radiometry::integrate(&[&passband]).expect("bounded"); + let table = camera + .channels + .iter() + .map(|channel| { + (0..T_STEPS) + .map(|i| { + let t = temperature_step(i); + let shape = blackbody_spectrum(t); + // Scale the blackbody so its mean through V is the zero point. + let through_v = + crate::radiometry::integrate(&[&shape, &passband]).expect("bounded"); + let scale = V_ZERO_POINT * band_width / through_v; + let photons = crate::radiometry::integrate_weighted( + &[&shape, &camera.transmission, &channel.quantum_efficiency], + |l| l / (PLANCK * SPEED_OF_LIGHT), + ) + .expect("bounded"); + scale * photons + }) + .collect() + }) + .collect(); + let diameter = camera.focal_length / camera.f_number; + Self { + table, + aperture: std::f64::consts::PI * diameter * diameter / 4.0, + } + } + + /// The faintest V magnitude whose star can give [`crate::gpu::stars::FAINTEST`] + /// electrons in some channel in an exposure of `exposure` seconds, at any color: the + /// limit to load a catalog to for that exposure. + pub fn faintest_drawn(&self, exposure: f64) -> f64 { + let brightest = self.table.iter().flatten().fold(0.0f64, |a, &b| a.max(b)); + 2.5 * (brightest * self.aperture * exposure / crate::gpu::stars::FAINTEST).log10() + } + + /// Number of channels. + pub fn channels(&self) -> usize { + self.table.len() + } + + /// The electron rate per square meter of aperture, in channel `channel`, from a star + /// of magnitude `v` and color `bv` on the boresight. + pub fn rate(&self, channel: usize, v: f32, bv: f32) -> f64 { + self.rates(v, bv)[channel] + } + + /// [`StarPhotometry::rate`] in every channel, 0 past the camera's channels. + pub fn rates(&self, v: f32, bv: f32) -> [f64; 4] { + let t = temperature(bv); + let x = ((t / T_MIN).ln() / (T_MAX / T_MIN).ln() * (T_STEPS - 1) as f64) + .clamp(0.0, (T_STEPS - 1) as f64); + let i = (x as usize).min(T_STEPS - 2); + let f = x - i as f64; + let scale = 10f64.powf(-0.4 * v as f64); + std::array::from_fn(|k| { + self.table.get(k).map_or(0.0, |table| { + (table[i] * (1.0 - f) + table[i + 1] * f) * scale + }) + }) + } +} + +fn temperature_step(i: usize) -> f64 { + T_MIN * (T_MAX / T_MIN).powf(i as f64 / (T_STEPS - 1) as f64) +} + +/// A blackbody's spectral shape at `t` kelvin, from 200 to 1200 nm every 2 nm. +fn blackbody_spectrum(t: f64) -> Spectrum { + let wavelength: Vec = (0..=500).map(|i| (200.0 + 2.0 * i as f64) * 1e-9).collect(); + let value = wavelength.iter().map(|&l| blackbody(l, t)).collect(); + Spectrum::sampled(wavelength, value).expect("a valid table") +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + fn random_stars(n: usize) -> Vec { + let mut state = 0x9e37_79b9_7f4a_7c15u64; + let mut random = move || { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + (state >> 11) as f64 / (1u64 << 53) as f64 + }; + (0..n) + .map(|_| CatalogStar { + ra: random() * std::f64::consts::TAU, + dec: (2.0 * random() - 1.0).asin(), + pm_ra: 0.0, + pm_dec: 0.0, + v: (random() * 12.0) as f32, + bv: 0.6, + }) + .collect() + } + + #[test] + fn a_field_holds_exactly_the_stars_within_it() { + // Against a brute-force search over 200,000 random stars, at the equator, near a + // pole, across right ascension 0, and for a wide and a narrow field. + let catalog = StarCatalog::new(random_stars(200_000), 0.0); + for (ra, dec, radius) in [ + (1.0f64, 0.1f64, 0.1f64), + (0.01, -0.3, 0.12), + (6.27, 0.5, 0.05), + (2.0, 1.53, 0.08), + (4.0, -1.5, 0.2), + (3.0, 0.0, 0.002), + ] { + let center = Vector3::new(dec.cos() * ra.cos(), dec.cos() * ra.sin(), dec.sin()); + let mut found: Vec<_> = catalog + .within(¢er, radius) + .map(|s| s.direction.cast::()) + .collect(); + let mut all: Vec<_> = catalog + .stars() + .iter() + .map(|s| s.direction.cast::()) + .filter(|d| d.dot(¢er) >= radius.cos()) + .collect(); + let key = |v: &Vector3| (v.x, v.y, v.z); + found.sort_by(|a, b| key(a).partial_cmp(&key(b)).unwrap()); + all.sort_by(|a, b| key(a).partial_cmp(&key(b)).unwrap()); + assert_eq!(found, all, "{ra} {dec} {radius}"); + } + } + + #[test] + fn the_star_file_reads_back() { + // Three stars in scripts/fetch-stars.py's format, brightest first: Sirius, a + // star with no color, and one past the faint limit. + let record = |ra_deg: f64, dec_deg: f64, pm: [i16; 2], v: f64, bv: Option| { + let mut r = Vec::new(); + r.extend(((ra_deg / 360.0 * 2f64.powi(32)).round() as u32).to_le_bytes()); + r.extend(((dec_deg / 180.0 * 2f64.powi(32)).round() as i32).to_le_bytes()); + r.extend(pm[0].to_le_bytes()); + r.extend(pm[1].to_le_bytes()); + r.extend((((v + 5.0) * 1000.0).round() as u16).to_le_bytes()); + r.extend( + bv.map_or(i16::MIN, |b| (b * 1000.0).round() as i16) + .to_le_bytes(), + ); + r + }; + let mut file = b"RSTARS02".to_vec(); + file.extend(3u32.to_le_bytes()); + file.extend(0u32.to_le_bytes()); + file.extend(record( + 101.2872, + -16.7161, + [-546, -1223], + -1.44, + Some(0.009), + )); + file.extend(record(10.0, 45.0, [0, 0], 9.5, None)); + file.extend(record(200.0, -60.0, [0, 0], 14.0, Some(1.2))); + let path = std::env::temp_dir().join(format!("resolved-stars-{}.bin", std::process::id())); + std::fs::write(&path, &file).unwrap(); + let catalog = StarCatalog::load(&path, 0.0, 12.0).unwrap(); + std::fs::remove_file(&path).unwrap(); + assert_eq!(catalog.len(), 2); + let sirius = catalog.stars().iter().find(|s| s.v < 0.0).unwrap(); + assert_relative_eq!(sirius.v, -1.44, epsilon = 1e-3); + assert_relative_eq!(sirius.bv, 0.009, epsilon = 1e-3); + let d = sirius.direction.cast::(); + assert_relative_eq!(d.y.atan2(d.x).to_degrees(), 101.2872, epsilon = 1e-6); + assert_relative_eq!(d.z.asin().to_degrees(), -16.7161, epsilon = 1e-6); + let colorless = catalog.stars().iter().find(|s| s.v > 0.0).unwrap(); + assert_eq!(colorless.bv, DEFAULT_BV); + } + + #[test] + fn proper_motion_moves_barnards_star() { + // Barnard's star: -798.6 mas/yr in RA times cos(dec), 10328.1 mas/yr in dec + // (Hipparcos). In 26 years, 268.9 arcseconds of motion. + let star = CatalogStar { + ra: 269.454_f64.to_radians(), + dec: 4.668_f64.to_radians(), + pm_ra: -798.58, + pm_dec: 10328.12, + v: 9.5, + bv: 1.57, + }; + let catalog = StarCatalog::new([star], 26.0); + let start = moved_direction(&star, 0.0); + let moved = catalog.stars()[0] + .direction + .cast::() + .angle(&start) + .to_degrees() + * 3600.0; + let expected = 26.0 * (798.58f64.hypot(10328.12)) / 1000.0; + // To the single precision directions are kept in, 0.012 arcseconds. + assert!((moved - expected).abs() < 0.02, "{moved} vs {expected}"); + } + + #[test] + fn color_gives_temperature() { + // The Sun: B - V 0.65 and 5772 K. For B - V 0 the formula gives 10,125 K, the + // temperature of the blackbody with that color; Vega's effective temperature, + // about 9,600 K, differs because an A star's spectrum is not a blackbody's. + assert_relative_eq!(temperature(0.65), 5772.0, max_relative = 0.01); + assert_relative_eq!(temperature(0.0), 10_125.2, max_relative = 1e-4); + } + + #[test] + fn a_zero_magnitude_star_in_the_example_camera() { + // An independent estimate: a V = 0 star gives about 1000 photons / s / cm^2 / A + // near 550 nm. Through the green channel (QE 0.6 peak, 80 nm wide: an equivalent + // width of 51.1 nm) and 0.9 transmission, onto the 12.5 mm aperture for 0.1 ms, + // that is about 56 electrons for a Sun-like star. + let camera = CameraConfig::example(); + let p = StarPhotometry::new(&camera); + let electrons = p.rate(1, 0.0, 0.65) * p.aperture * camera.exposure; + let estimate = 1000.0 * 1e4 * 1e10 * 0.6 * 1.0645 * 80e-9 * 0.9 * p.aperture * 1e-4; + assert_relative_eq!(electrons, estimate, max_relative = 0.05); + // Five magnitudes is a factor of 100. + assert_relative_eq!( + p.rate(1, 5.0, 0.65) * 100.0, + p.rate(1, 0.0, 0.65), + max_relative = 1e-9 + ); + // A hot star is bluer: more blue than red electrons; a cool one the reverse. + assert!(p.rate(2, 0.0, -0.2) > p.rate(0, 0.0, -0.2)); + assert!(p.rate(2, 0.0, 1.5) < p.rate(0, 0.0, 1.5)); + } +} diff --git a/crates/resolved/tests/gpu/distortion.rs b/crates/resolved/tests/gpu/distortion.rs index 3fbcf58..f943fb1 100644 --- a/crates/resolved/tests/gpu/distortion.rs +++ b/crates/resolved/tests/gpu/distortion.rs @@ -302,7 +302,7 @@ fn the_sensor_uses_the_undistorted_line_of_sight() { targets.radiance.size(), ); let mut encoder = ctx.device.create_command_encoder(&Default::default()); - SensorPass::new(ctx).encode( + SensorPass::new(ctx, false).encode( ctx, &mut encoder, &targets, diff --git a/crates/resolved/tests/gpu/main.rs b/crates/resolved/tests/gpu/main.rs index add40bf..e241924 100644 --- a/crates/resolved/tests/gpu/main.rs +++ b/crates/resolved/tests/gpu/main.rs @@ -19,3 +19,4 @@ mod psf; mod renderer; mod sensor; mod shadows; +mod stars; diff --git a/crates/resolved/tests/gpu/sensor.rs b/crates/resolved/tests/gpu/sensor.rs index 67bd2dc..5f0a5d9 100644 --- a/crates/resolved/tests/gpu/sensor.rs +++ b/crates/resolved/tests/gpu/sensor.rs @@ -32,7 +32,7 @@ fn sense(ctx: &GpuContext, name: &str, radiance: &[[f32; 4]]) -> Vec { }, targets.radiance.size(), ); - let pass = SensorPass::new(ctx); + let pass = SensorPass::new(ctx, false); let uniforms = SensorUniforms::new(&config, &Radiometry::new(&config), false); let mut encoder = ctx.device.create_command_encoder(&Default::default()); pass.encode(ctx, &mut encoder, &targets, &uniforms, None); @@ -164,7 +164,7 @@ fn sunlit_earth_through_the_whole_chain() { let targets = FrameTargets::new(ctx, &TargetLayout::new(config.width, config.height, 3)); let mut encoder = ctx.device.create_command_encoder(&Default::default()); BodiesPass::new(ctx, 1, false).encode(ctx, &mut encoder, &targets, &bodies, None); - SensorPass::new(ctx).encode(ctx, &mut encoder, &targets, &sensor, None); + SensorPass::new(ctx, false).encode(ctx, &mut encoder, &targets, &sensor, None); ctx.queue.submit([encoder.finish()]); let dn = targets.read_dn(ctx); common::save_dn("sensor-earth", &ctx.info, &dn); diff --git a/crates/resolved/tests/gpu/stars.rs b/crates/resolved/tests/gpu/stars.rs new file mode 100644 index 0000000..ab4e1ee --- /dev/null +++ b/crates/resolved/tests/gpu/stars.rs @@ -0,0 +1,261 @@ +//! Lighting tier L1's stars, seen by camera tier C1: where they land, how bright, and +//! what hides them, on the GPU and on the CPU (lavapipe). + +use crate::common; +use crate::pixel_area::plate; + +use std::sync::Arc; + +use resolved::camera::Detector; +use resolved::gpu::AdapterChoice; +use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3}; +use resolved::radiometry::AU; +use resolved::scene::{aberrate, Velocities, SPEED_OF_LIGHT}; +use resolved::stars::CatalogStar; +use resolved::{ + BrownConrady, CameraConfig, CameraTier, EarthTexture, Fidelity, LightingTier, Pose, Renderer, + RendererOptions, Scene, StarCatalog, StarPhotometry, Target, +}; + +const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu]; + +/// 256 x 256 pixels of the example camera's 110 urad, with a lens that distorts and a +/// detector of 1 electron per DN, so DN read back the electrons to 0.5. +fn camera(exposure: f64) -> CameraConfig { + CameraConfig { + width: 256, + height: 256, + exposure, + distortion: Some(BrownConrady { + k1: 0.3, + p1: 1e-3, + ..Default::default() + }), + detector: Detector { + gain: 1.0, + bias: 100.0, + full_well: 60_000.0, + adc_bits: 16, + }, + ..CameraConfig::example() + } +} + +/// Vega: V 0.03, B - V 0.00, at right ascension 279.2347 and declination 38.7837 +/// degrees, with no proper motion here. +fn vega() -> CatalogStar { + CatalogStar { + ra: 279.2347f64.to_radians(), + dec: 38.7837f64.to_radians(), + pm_ra: 0.0, + pm_dec: 0.0, + v: 0.03, + bv: 0.0, + } +} + +fn direction(s: &CatalogStar) -> Vector3 { + Vector3::new( + s.dec.cos() * s.ra.cos(), + s.dec.cos() * s.ra.sin(), + s.dec.sin(), + ) +} + +/// A renderer with `stars` for a catalog. At L1, Earth's surface is black. +fn renderer( + adapter: &AdapterChoice, + camera: CameraConfig, + lighting: LightingTier, + tier: CameraTier, + stars: &[CatalogStar], +) -> Renderer { + let black = EarthTexture::from_rgba8(1, 1, vec![0, 0, 0, 255]).unwrap(); + let options = RendererOptions { + adapter: adapter.clone(), + fidelity: Fidelity { + lighting, + camera: tier, + }, + earth_texture: Some(Arc::new(black)), + moon_texture: Some(common::uniform_moon()), + stars: Some(Arc::new(StarCatalog::new(stars.iter().copied(), 0.0))), + ..Default::default() + }; + Renderer::new(camera, options).unwrap() +} + +/// A scene with the Sun far behind `toward` and Earth and the Moon out of view of a +/// camera 4.2e7 m out on +X looking that way, with Earth's orbital velocity. +fn deep_space(toward: Vector3, targets: Vec) -> Scene { + Scene { + earth_attitude: UnitQuaternion::identity(), + sun: Vector3::new(4.2e7, 0.0, 0.0) - toward * AU, + moon: Vector3::new(0.0, -3.844e8, 0.0), + moon_attitude: UnitQuaternion::identity(), + targets, + velocities: Velocities { + earth: Vector3::new(0.0, 29.8e3, 0.0), + ..Default::default() + }, + } +} + +/// A camera at `position` looking toward `toward`, rolled so that its boresight's +/// image is `offset` pixels from the star's, with the image's up toward +Z. +fn pointing(position: Vector3, toward: Vector3) -> Pose { + Pose::looking_at(position, position + toward, Vector3::z()).unwrap() +} + +/// Electrons above bias per pixel, from a frame's DN in channel `k`. +fn electrons(frame: &resolved::Frame, k: u32) -> Vec { + (0..256 * 256) + .map(|i| frame.dn_at(i % 256, i / 256, k) as f64 - 100.0) + .collect() +} + +#[test] +fn a_star_lands_where_its_apparent_position_projects() { + // Vega swept across a pixel in tenths, behind a lens with distortion and seen with + // Earth's 30 km/s of aberration: the centroid of its electrons tracks where + // Intrinsics::project puts its apparent direction to 0.01 px, and they total what + // StarPhotometry gives, times the aperture, the exposure and cos(off axis), to + // 0.2 %. + let _devices = common::device_lock(); + let config = camera(0.02); + let k = config.intrinsics(); + let photometry = StarPhotometry::new(&config); + let star = direction(&vega()); + let beta = Vector3::new(0.0, 29.8e3, 0.0) / SPEED_OF_LIGHT; + let apparent = aberrate(&star, &beta); + let position = Vector3::new(4.2e7, 0.0, 0.0); + for adapter in &ADAPTERS { + let mut r = renderer( + adapter, + config.clone(), + LightingTier::L1, + CameraTier::C1, + &[vega()], + ); + let name = r.adapter_info().name.clone(); + for step in 0..10 { + // Look 40.3 px left of and 30.1 px above the star, plus a tenth of a pixel + // per step: off axis, where the distortion moves it about 0.3 px. + let look = Vector3::new((-40.3 - step as f64 / 10.0) / k.fx, 30.1 / k.fy, 1.0); + let base = pointing(position, star); + let toward = base.attitude * look.normalize(); + let camera = pointing(position, toward); + let scene = deep_space(star, Vec::new()); + let frame = r.submit(&scene, &camera).unwrap().wait().unwrap(); + let d = camera.attitude.inverse() * apparent; + let want = k.project(&d).unwrap(); + let image = electrons(&frame, 1); + let (u0, v0) = (want.x as usize, want.y as usize); + let (mut total, mut centroid) = (0.0, Vector2::zeros()); + for v in v0 - 5..=v0 + 5 { + for u in u0 - 5..=u0 + 5 { + let e = image[v * 256 + u]; + total += e; + centroid += Vector2::new(u as f64 + 0.5, v as f64 + 0.5) * e; + } + } + centroid /= total; + let off = (centroid - want).norm(); + assert!(off < 0.01, "{name}: step {step}: centroid off by {off} px"); + let expected = + photometry.rate(1, 0.03, 0.0) * photometry.aperture * 0.02 * d.normalize().z; + assert!( + (total / expected - 1.0).abs() < 2e-3, + "{name}: {total} electrons, expected {expected}" + ); + } + } +} + +#[test] +fn stars_hide_behind_what_is_nearer() { + // Vega in clear sky, behind Earth seen from low orbit, and behind a plate: the + // bias alone wherever something nearer is drawn. Lighting L0 has no stars, and a C0 + // camera, sampling a point per pixel, sees none. + let _devices = common::device_lock(); + let config = camera(0.02); + let star = direction(&vega()); + let peak = |frame: &resolved::Frame| electrons(frame, 1).into_iter().fold(0.0f64, f64::max); + for adapter in &ADAPTERS { + let mut r = renderer( + adapter, + config.clone(), + LightingTier::L1, + CameraTier::C1, + &[vega()], + ); + let name = r.adapter_info().name.clone(); + let model = r.load_model(&plate(5.0)).unwrap(); + let clear = Vector3::new(4.2e7, 0.0, 0.0); + let frame = r + .submit(&deep_space(star, Vec::new()), &pointing(clear, star)) + .unwrap() + .wait() + .unwrap(); + assert!(peak(&frame) > 1000.0, "{name}: {}", peak(&frame)); + + // 500 km above the point on Earth opposite the star, looking at the star + // through Earth's night side, with the Sun beyond Earth too. + let low = -star * 6.878e6; + let mut scene = deep_space(star, Vec::new()); + scene.sun = star * AU; + let frame = r + .submit(&scene, &pointing(low, star)) + .unwrap() + .wait() + .unwrap(); + assert_eq!(peak(&frame), 0.0, "{name}: behind Earth"); + + // A plate 50 m along the line of sight, facing the camera. + let camera = pointing(clear, star); + let target = Target { + model, + pose: Pose::new( + clear + star * 50.0 - camera.attitude * Vector3::new(5.0, 0.0, 0.0), + camera.attitude * UnitQuaternion::identity(), + ), + }; + let mut scene = deep_space(star, vec![target]); + // Light the plate from behind the star; Earth's light still reaches its face, so + // the frame is compared with the same one drawn with no stars. + scene.sun = clear + star * AU; + let frame = r.submit(&scene, &camera).unwrap().wait().unwrap(); + let mut none = renderer( + adapter, + config.clone(), + LightingTier::L1, + CameraTier::C1, + &[], + ); + let model = none.load_model(&plate(5.0)).unwrap(); + scene.targets[0].model = model; + let starless = none.submit(&scene, &camera).unwrap().wait().unwrap(); + assert_eq!( + electrons(&frame, 1), + electrons(&starless, 1), + "{name}: behind a plate" + ); + + for (lighting, tier) in [ + (LightingTier::L0, CameraTier::C1), + (LightingTier::L1, CameraTier::C0), + ] { + let mut r = renderer(adapter, config.clone(), lighting, tier, &[vega()]); + let frame = r + .submit(&deep_space(star, Vec::new()), &pointing(clear, star)) + .unwrap() + .wait() + .unwrap(); + assert_eq!( + peak(&frame), + 0.0, + "{name}: {lighting} {tier} draws no stars" + ); + } + } +} diff --git a/docs/optics.md b/docs/optics.md index 3669f43..d195a4e 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -1359,6 +1359,104 @@ under 0.01 microradians; the rotation and lighting of a body during its light time (the Moon turns 3 microradians in 1.3 s); and the light-time delay between the Sun lighting a body and the camera seeing it. +## Stars + +Code: `stars.rs` (`StarCatalog`, `StarPhotometry`); `gpu/stars.rs`, +`shaders/stars.wgsl`; `scripts/fetch-stars.py`. + +Below lighting tier L1 the sky is black. A real camera shows stars wherever +its exposure is long enough, and they matter: star trackers find their +attitude from them, a camera measuring the Moon or a target measures it +against them, and to an algorithm looking for a target a star is a false one. + +Stars are part of the scene, so they come with lighting tier L1, but only a +camera at tier C1 or above shows them. A C0 pixel takes its value from the +one point at its center, and a star is a point that almost never lies +exactly there; from C1 each pixel collects the light over its whole area, +which is what records a star. + +**The catalogs.** `scripts/fetch-stars.py` merges three of ESA's catalogs +into `data/stars/stars.bin`: Hipparcos, 118,000 stars complete to about +magnitude 7.3; Tycho-2, 2.5 million stars to about 11.5 (99 % complete to +11.0); and Gaia DR3, from which it takes every star brighter than G 15 by +default, 37 million. Each star comes from the first catalog that has it, so +the brightest, which saturate Gaia's detectors, keep their Hipparcos or Tycho +measurements. Gaia measures brightness in its own bands, G and the color +BP - RP: V comes from them by Gaia's published relation (to 0.03 mag), and +B - V by a fit to the Hipparcos stars Gaia matched. Each star +keeps its position on the sky (ICRS, at epoch J2000), its proper motion and +its brightness and color. The renderer moves each star by its proper motion +to the scene's epoch when the catalog is loaded. Barnard's star, the +fastest, moves 10 arcseconds a year; most move under 0.1. Parallax, under 1 +arcsecond for every star, is left out. + +**Magnitudes.** A star's brightness is given as a magnitude: each step of 5 +is a factor of 100, and larger is fainter. Vega is about 0, the faintest +stars the eye sees about 6, Tycho-2's limit about 11.5, 25,000 times fainter +than 0. The magnitude in the Johnson V band, a filter close to the eye's +green-yellow, fixes how much light the star sends through that band: V = 0 +is 3.631e-11 W/m^2/nm averaged over it (Bessell, Castelli and Plez 1998). +The color index B - V, the difference between the star's magnitudes in blue +and in V, fixes its temperature: hot stars are blue, B - V below 0, and cool +ones red, up to about 2. The renderer takes the temperature from B - V +(Ballesteros 2012: the Sun's 0.65 gives 5,780 K), takes a blackbody spectrum +at that temperature, scales it so its light through the V band matches the +star's magnitude, and folds it through the camera's optics and channels, as +for sunlight. A star's spectrum is not a blackbody's: absorption lines and, +in cool stars, molecular bands move a few percent of the light between +channels. For the example camera, a Sun-colored star of V = 0 on the +boresight gives 56 electrons in green in the 0.1 ms exposure, against an +independent estimate of 56.4. + +**Exposure.** That 56 electrons is why the example scenes show no stars: an +exposure set for a sunlit spacecraft or Earth, a fraction of a millisecond, +leaves even the brightest stars at a few hundred electrons and the rest +below one. Star trackers expose for 10 to 200 ms, and a telescope's larger +aperture collects more: a 1 m f/8 telescope, with a 12.5 cm aperture, gives +100 times what the example camera does. + +**Drawing.** Each frame the renderer finds the stars in the field on the +CPU, from bands of declination sorted by right ascension, so it looks only +at the stars near the field. It places each at its apparent position, the +annual aberration of the observer's motion (about 20 arcseconds, see the +previous section) through the lens's distortion. It works out its electrons +in each channel: a point source loses light off axis only to the aperture's +foreshortening, cos(theta), not the cos^4(theta) of an extended surface, +since all of a star's light that enters lands in its spot. A star is not +rendered and then blurred: a point of light rendered into a pixel loses +where in the pixel it fell (see the blur section). Instead its electrons are +spread over the pixels around it as the Gaussian PSF integrated over each +pixel's area, which is exact for any position within a pixel. The detector +adds them to the light from everything else before the full well, so a +bright star saturates like any other light. A star is hidden wherever +anything nearer is drawn, Earth, the Moon or a target: stars are at +infinity, so a star shows only where the depth image is empty. Stars giving +fewer than 0.1 electrons in every channel, far below any detector's noise, are +skipped, and the catalog is loaded only to the magnitude that can reach that. + +Measured: a star moved across a pixel in tenths, behind a lens with +distortion and with Earth's 30 km/s of aberration, lands where the +projection puts its apparent direction to 0.01 px, with its electrons +matching the photometry to 0.2 %. A star behind Earth's night side or +behind a target leaves the bias alone. + +Costs. A catalog is loaded only as deep as the camera's exposure can show: +for the example camera at 0.1 ms that is V 7.9, 42,000 stars, in 6 ms and 7 +MB; at 50 ms, V 14.7, 21 million stars, in 2.8 s and 670 MB; for a 1 m +telescope at 0.2 s, the whole catalog, 37 million stars, in 5 s and 1.2 GB. +Each frame, finding and placing the stars takes 23 microseconds for the +telescope's 0.6 degree field (490 stars toward the Milky Way in Monoceros) +and 4.8 ms for the example camera's 12 degree field at 50 ms (52,000 stars). +Drawing them takes 0.24 ms for 40,000 stars on an RTX 3060, and 5.6 ms on +the CPU through lavapipe. + +A catalog's depth matters most for a telescope's narrow field. Tycho-2 alone +averages 60 stars per square degree, so a 1 m telescope's 0.2 square degrees +would hold about a dozen; with no sky background in space such a telescope +records stars to magnitude 14 or 15 in half a second, and Gaia to G 15 +averages about 900 per square degree. Deeper still (G 16, 78 million stars) +is a run of the script with `--gaia 16`, at twice the size. + ## Not modeled yet - Earthshine on the Moon's night side diff --git a/examples/scenes/geo-telescope-1000km.toml b/examples/scenes/geo-telescope-1000km.toml new file mode 100644 index 0000000..fea63d0 --- /dev/null +++ b/examples/scenes/geo-telescope-1000km.toml @@ -0,0 +1,37 @@ +# GEO to GEO at 1000 km, exposed 0.2 s at f/4 as a surveillance camera would: the +# target a few pixels across and far past saturation, in a field of stars to the +# catalog's limit, about V 15, which gives a few DN. +# A 1 m, f/4 telescope (25 cm aperture, 5.5 urad per pixel, 0.6 by 0.34 +# degrees) on one geostationary satellite, imaging an SSL-1300 1000 km ahead along +# the arc, lit at a phase of 30 degrees, against the stars toward the Milky Way in +# Monoceros. Stars are drawn at lighting L1 with camera C1 (scripts/fetch-stars.py). +# Generated by scripts/make-example-scenes.py; the model comes from +# scripts/fetch-models.sh. +# +# resolved render --scene --lighting L1 --camera C1 --pretty --preview geo.png + +sun = [101746122115, -106547756265, 25977397630] +moon = [0, -384400000, 0] +earth_velocity = [-19407.62, -20723.58, -8984.77] + +[camera] +focal_length = 1.0 +f_number = 4.0 +exposure = 0.2 + +[observer] +position = [41083339.41, 9484836.26, 0.00] +look_at = [40846855.46, 10456447.53, 0.00] +up = [0, 0, 1] +velocity = [-691.66, 2995.90, 0.00] + +# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh), and lights +# the target's Earth-facing side. +[appearance] +earth_month = 7 + +[[target]] +model = "../models/nasa/ssl1300.glb" +scale = 0.1 +position = [40846855.0383, 10456447.7315, 0.2096] +attitude = [0.466967605, 0.081764177, 0.291661326, -0.830776472] diff --git a/examples/scenes/geo-telescope-100km-side.toml b/examples/scenes/geo-telescope-100km-side.toml new file mode 100644 index 0000000..d65f31b --- /dev/null +++ b/examples/scenes/geo-telescope-100km-side.toml @@ -0,0 +1,35 @@ +# GEO to GEO at 100 km, lit from the side: half the target in shadow. +# A 1 m, f/8 telescope (12.5 cm aperture, 5.5 urad per pixel, 0.6 by 0.34 +# degrees) on one geostationary satellite, imaging an SSL-1300 100 km ahead along +# the arc, lit at a phase of 90 degrees, against the stars toward the Milky Way in +# Monoceros. Stars are drawn at lighting L1 with camera C1 (scripts/fetch-stars.py). +# Generated by scripts/make-example-scenes.py; the model comes from +# scripts/fetch-models.sh. +# +# resolved render --scene --lighting L1 --camera C1 --pretty --preview geo.png + +sun = [143509808059, 33311150104, 25977397630] +moon = [0, -384400000, 0] +earth_velocity = [8161.41, -26276.53, -11392.27] + +[camera] +focal_length = 1.0 +f_number = 8.0 +exposure = 0.0006 + +[observer] +position = [41083339.41, 9484836.26, 0.00] +look_at = [41060728.78, 9582246.49, 0.00] +up = [0, 0, 1] +velocity = [-691.66, 2995.90, 0.00] + +# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh), and lights +# the target's Earth-facing side. +[appearance] +earth_month = 7 + +[[target]] +model = "../models/nasa/ssl1300.glb" +scale = 0.1 +position = [41060728.3613, 9582246.7015, 0.2096] +attitude = [0.462527700, 0.083319399, 0.291220856, -0.833256513] diff --git a/examples/scenes/geo-telescope-100km.toml b/examples/scenes/geo-telescope-100km.toml new file mode 100644 index 0000000..e37186b --- /dev/null +++ b/examples/scenes/geo-telescope-100km.toml @@ -0,0 +1,36 @@ +# GEO to GEO at 100 km: the target about 55 px across, lit nearly face on, exposed +# short of saturating it, with only the brighter stars behind it. +# A 1 m, f/8 telescope (12.5 cm aperture, 5.5 urad per pixel, 0.6 by 0.34 +# degrees) on one geostationary satellite, imaging an SSL-1300 100 km ahead along +# the arc, lit at a phase of 30 degrees, against the stars toward the Milky Way in +# Monoceros. Stars are drawn at lighting L1 with camera C1 (scripts/fetch-stars.py). +# Generated by scripts/make-example-scenes.py; the model comes from +# scripts/fetch-models.sh. +# +# resolved render --scene --lighting L1 --camera C1 --pretty --preview geo.png + +sun = [100603206248, -107627564420, 25977397630] +moon = [0, -384400000, 0] +earth_velocity = [-19658.75, -20523.36, -8897.97] + +[camera] +focal_length = 1.0 +f_number = 8.0 +exposure = 0.0006 + +[observer] +position = [41083339.41, 9484836.26, 0.00] +look_at = [41060728.78, 9582246.49, 0.00] +up = [0, 0, 1] +velocity = [-691.66, 2995.90, 0.00] + +# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh), and lights +# the target's Earth-facing side. +[appearance] +earth_month = 7 + +[[target]] +model = "../models/nasa/ssl1300.glb" +scale = 0.1 +position = [41060728.3613, 9582246.7015, 0.2096] +attitude = [0.462527700, 0.083319399, 0.291220856, -0.833256513] diff --git a/examples/scenes/stars-pleiades.toml b/examples/scenes/stars-pleiades.toml new file mode 100644 index 0000000..95c08eb --- /dev/null +++ b/examples/scenes/stars-pleiades.toml @@ -0,0 +1,24 @@ +# A star field: the Pleiades through the example camera (12 by 7 degrees), from +# geostationary orbit, long-exposed at 5 s. Stars are drawn at lighting L1 with camera +# C1, from the catalog scripts/fetch-stars.py writes. Alcyone, the brightest at V 2.9, +# gives about 200,000 electrons, ten times the full well; a V 8 star about 2,000 and a +# V 11 star about 130. Generated by scripts/make-example-scenes.py. +# +# resolved render --scene examples/scenes/stars-pleiades.toml --lighting L1 --camera C1 --pretty --preview pleiades.png + +sun = [-74862158105, -114183767096, -61133031444] +moon = [0, -384400000, 0] +earth_velocity = [-25760.94, 13707.76, 5943.04] + +[camera] +exposure = 5 + +[observer] +position = [42164000.00, 0.00, 0.00] +look_at = [542586618.01, 763271339.10, 408649074.74] +up = [0, 0, 1] +velocity = [0.00, 3074.70, 0.00] + +# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh). +[appearance] +earth_month = 7 diff --git a/readme.md b/readme.md index 5bd0730..da6524c 100644 --- a/readme.md +++ b/readme.md @@ -121,7 +121,14 @@ comes from the LRO Camera's normalized-reflectance mosaic, which `scripts/fetch-moon.py` downloads (1.2 GB) and assembles (with numpy and Pillow) into `data/moon/wac-hapke-643.png`. The Moon is turned by the scene's epoch, or by `moon_attitude` in a scene without one (by default facing Earth, which can -be off by up to about 10 degrees). +be off by up to about 10 degrees). Also at lighting L1, and shown by camera C1, +stars come from the Hipparcos, Tycho-2 and Gaia DR3 catalogs, which `scripts/fetch-stars.py` downloads (210 MB of +files and 2.3 GB of Gaia queries, about an hour) and merges into +`data/stars/stars.bin` (590 MB, 37 million stars to G 15); their proper motion +is taken to the scene's epoch, or J2000 in a scene without one. The CLI loads +only the stars bright enough to show at the camera's exposure. +`stars-pleiades.toml` is a star field, and `geo-telescope-*.toml` a telescope +on one geostationary satellite imaging another against the stars. ```sh cargo run --release -p resolved-cli -- render \ @@ -166,8 +173,8 @@ cargo run --release -p resolved-showcase -- --serve 0.0.0.0:8010 L0 and C0). `--supersampling N` renders each pixel as N x N finer ones from C1 (default 2). A tier the scene or the data cannot support stops with an error naming what is missing. -- `--data DIR` says where the tiers' data is (`earth/bmng-2004-MM.jpg` and - `moon/wac-hapke-643.png`). +- `--data DIR` says where the tiers' data is (`earth/bmng-2004-MM.jpg`, + `moon/wac-hapke-643.png` and `stars/stars.bin`). Without it the CLI tries `RESOLVED_DATA`, then `./data`, then the repository's `data/`. - `--adapter` takes `auto`, `gpu`, `cpu` or part of an adapter's name; @@ -191,14 +198,14 @@ Lighting L0 and L1 and camera C0 are in: the library renders Earth, Moon, Sun and targets to DN through `Renderer`, and the CLI renders scene files. Lighting L1 draws Earth's surface from Blue Marble and the Moon's from the LRO Camera's map, shadows on the target, light from Earth and the Moon on it, eclipses on -Earth and the Moon, and the Sun's limb darkening. Camera C1 is under way: so -far, pixels that collect light over their area, lens distortion -(Brown-Conrady), lens blur (a Gaussian PSF) applied to the scene rendered -2 x 2 finer than the pixels, and the Sun, Moon and Earth at their apparent -positions (light time and aberration). Planned, in order: +Earth and the Moon, the Sun's limb darkening, and stars from Hipparcos, +Tycho-2 and Gaia, which camera C1 shows. Camera C1 is under way: so far, +pixels that collect light over their area, lens distortion (Brown-Conrady), +lens blur (a Gaussian PSF) applied to the scene rendered 2 x 2 finer than the +pixels, and the Sun, Moon and Earth at their apparent positions (light time +and aberration). Planned, in order: -- Camera: stars, - detector noise, pixel-to-pixel variation, saturation and blooming. +- Camera: detector noise, pixel-to-pixel variation, saturation and blooming. - Lighting: BRDF and spectral materials, atmospheric scattering, clouds, lunar terrain. - Portable ray tracing, path tracing, and label outputs for ML. diff --git a/scripts/fetch-stars.py b/scripts/fetch-stars.py new file mode 100755 index 0000000..fbe92c5 --- /dev/null +++ b/scripts/fetch-stars.py @@ -0,0 +1,399 @@ +#!/usr/bin/env python3 +"""Downloads the Hipparcos, Tycho-2 and Gaia DR3 star catalogs and writes the renderer's +star file. + +Sources: + - Hipparcos main catalog, ESA 1997 (CDS I/239, hip_main.dat): 118,218 entries, + complete to about V 7.3, with Johnson V and B-V. Positions are ICRS at epoch + J1991.25. The 2007 new reduction (I/311) improves the astrometry by about a + milliarcsecond, far below any pixel here, but carries no V, so the original is used. + - Tycho-2, Hog et al. 2000 (CDS I/259, tyc2.dat.00-19): 2.54 million stars to about + VT 11.5 (99 % complete to V 11.0), with Tycho BT and VT magnitudes. Mean positions + are ICRS at epoch J2000. + - Gaia DR3, Gaia Collaboration 2023 (gaiadr3.gaia_source_lite, from the ESA Gaia + archive's TAP service, https://gea.esac.esa.int/archive/): every star brighter than + G 15 by default (--gaia), 36.9 million, with its G magnitude and BP - RP color. + Positions are ICRS at epoch J2016.0. Gaia's photometry saturates for the brightest + stars, which Hipparcos and Tycho-2 cover. + +Every star is taken from one catalog: Hipparcos first, then Tycho-2 for the stars it +does not have (by Hipparcos number), then Gaia for the stars neither has (by Gaia's own +cross-matches to them, gaiadr3.hipparcos2_best_neighbour and +gaiadr3.tycho2tdsc_merge_best_neighbour). Magnitudes are converted to Johnson V and +B - V: + - Tycho: V = VT - 0.090 (BT - VT), B - V = 0.850 (BT - VT) (ESA 1997, vol. 1, 1.3, + appendix 4). A star with no BT keeps V = VT and no color. + - Gaia: V from G and BP - RP by Riello et al. 2021 (A&A 649, A3, table C.2), + G - V = -0.02704 + 0.01424 c - 0.2156 c^2 + 0.01426 c^3 for c = BP - RP in -0.5 to + 2.75, 0.03 mag scatter. B - V from BP - RP by a cubic fitted here to the Hipparcos + stars Gaia matched (printed when the script runs). A star with no BP - RP gets + c = 0.82, the median of stars near G 15, and no B - V. +Positions are carried to J2000 by their proper motion. The few Tycho-2 stars with no +mean position (pflag X) use their observed position, at about epoch 1991.5, with no +proper motion. Parallax is not kept: it moves every star by under 1 arcsecond. + +Hipparcos and Tycho-2 come from the CDS archive (https://cdsarc.cds.unistra.fr) and are +checked by sha256. Gaia comes from about 2,000 queries of the archive, each a range of +source_id (a patch of sky), 4 at a time, split in half when a reply reaches the +service's 50,000-row limit; a query's results are kept, so a run can be resumed, but +they are not pinned by checksum. At G 15 that is about 2.3 GB of text and takes about +an hour. Everything goes to data/stars/source/; the output is data/stars/stars.bin, +about 600 MB at G 15, both gitignored. Needs numpy and pandas. + +The output's format, little-endian: + + 8 bytes b"RSTARS02" + u32 number of stars + u32 0 + then per star, brightest first, 16 bytes: + u32 right ascension, ICRS, epoch J2000, in units of 2 pi / 2^32 + i32 declination, ICRS, epoch J2000, in units of pi / 2^32 + i16 proper motion in right ascension times cos(declination), mas/yr + i16 proper motion in declination, mas/yr + u16 Johnson V magnitude plus 5, in thousandths + i16 Johnson B - V in thousandths, or -32768 if unknown + +Usage: scripts/fetch-stars.py [--gaia G] [--faint V] + --gaia G take Gaia stars brighter than G (default 15; 0 for none) + --faint V keep stars to V (default: all those the catalogs give) +""" + +import argparse +import concurrent.futures +import gzip +import hashlib +import io +import math +import os +import sys +import time +import urllib.parse +import urllib.request + +import numpy as np +import pandas as pd + +CDS = "https://cdsarc.cds.unistra.fr/ftp/I" +GAIA_TAP = "https://gea.esac.esa.int/tap-server/tap/sync" +# The service's limit on rows per synchronous query. +GAIA_ROWS = 50_000 +# Gaia source_id carries the HEALPix level-12 cell in its high bits. +CELLS = 12 * 4**12 +CELL = 2**35 + +FILES = { + "hip_main.dat": ("239/hip_main.dat", "58ceabb104d647160d9437ce6e513a02a036bb4ad9f8879a5a22fd52943616e0"), + "tyc2.dat.00.gz": ("259/tyc2.dat.00.gz", "1b224570fc6eb151984ce106fdf797728649d1ea77aaf3effd2d3444cfac6df6"), + "tyc2.dat.01.gz": ("259/tyc2.dat.01.gz", "36221f4a5cdd9c5009d64299d2e6409f8d4cdf924321768780165bd3e2a2a99a"), + "tyc2.dat.02.gz": ("259/tyc2.dat.02.gz", "e2ea4eeb7d204dde70ac54d172b1c429540c0d23d62ec0d7500561adc44cc57b"), + "tyc2.dat.03.gz": ("259/tyc2.dat.03.gz", "29e0cb57cba7e1651455efcb528276b0c28e1fe83f47f8ebe963e48d8b180afb"), + "tyc2.dat.04.gz": ("259/tyc2.dat.04.gz", "d04868329b470aa320f0fe5e9c525be619a7a228b364433e5f114f2178946d9b"), + "tyc2.dat.05.gz": ("259/tyc2.dat.05.gz", "51f08772197dcf0cd26ff098c866036340b797f71ee04be22712517c10c73c4e"), + "tyc2.dat.06.gz": ("259/tyc2.dat.06.gz", "f55503a66abce0e11d1a4ccbcde1f4c97a6e3e769de85e005f57ba963086ab03"), + "tyc2.dat.07.gz": ("259/tyc2.dat.07.gz", "bc15b5b1f6308477360fbb93f6a98d9029601cab68973202be0f8764c127f524"), + "tyc2.dat.08.gz": ("259/tyc2.dat.08.gz", "315ac34c678bc2e0f568b43b6b2d7c3d7c9cf9c089b51115d18b872864e50426"), + "tyc2.dat.09.gz": ("259/tyc2.dat.09.gz", "d66c671fa29aad10bc5fd697f68c918bb774c1ffc3c7db9c20234895268390d5"), + "tyc2.dat.10.gz": ("259/tyc2.dat.10.gz", "d86529df819ebdf3f4b8892510973ee10e1db06930814e014dd467193e740d52"), + "tyc2.dat.11.gz": ("259/tyc2.dat.11.gz", "fc0c3203b3d2787da54c43d9963e363c8a3926576c383b101a15350a7ca23e9c"), + "tyc2.dat.12.gz": ("259/tyc2.dat.12.gz", "68fd9d7e7353d52dea0de043cc62b45c04d234104eb05c063dc23ea7b6759576"), + "tyc2.dat.13.gz": ("259/tyc2.dat.13.gz", "a08ecc092e6d134742c0f594c56c5c2e02b8fd5e76cffb1779877f2fcefac3fc"), + "tyc2.dat.14.gz": ("259/tyc2.dat.14.gz", "2846e8cc489795a4d8160d2d5b64ff0e9a00598adba35a52a4f8ccf4c5e38b3e"), + "tyc2.dat.15.gz": ("259/tyc2.dat.15.gz", "74784a3656382090a70d778820335f4a781509dbe3fcb9e92d18cf96b2c46c72"), + "tyc2.dat.16.gz": ("259/tyc2.dat.16.gz", "320560e3d551cc40fa1f54ef7133709bf4bd45efa6688f808be7c09be5544a4d"), + "tyc2.dat.17.gz": ("259/tyc2.dat.17.gz", "a8154e940aa0a0e8f31d91b4bd9fd56e6849da42cddb13c7dda772439b06991a"), + "tyc2.dat.18.gz": ("259/tyc2.dat.18.gz", "5e951a1a51df7956f205b83cbfa5501d357c843640a253c1e6d3917bebe7d928"), + "tyc2.dat.19.gz": ("259/tyc2.dat.19.gz", "f59605a38116f517a31a7dbdee3469c077658f2f40b8afe5da2aeb832eaee3dd"), +} + +# A star's fields as arrays: right ascension and declination (radians, J2000), proper +# motion (mas/yr), V and B - V (NaN if unknown). +FIELDS = ("ra", "dec", "pm_ra", "pm_dec", "v", "bv") + + +def fetch(source): + os.makedirs(source, exist_ok=True) + for name, (path, digest) in FILES.items(): + target = os.path.join(source, name) + if not os.path.exists(target): + print(f"fetching {name}") + urllib.request.urlretrieve(f"{CDS}/{path}", target + ".part") + os.replace(target + ".part", target) + h = hashlib.sha256() + with open(target, "rb") as f: + for block in iter(lambda: f.read(1 << 20), b""): + h.update(block) + if h.hexdigest() != digest: + sys.exit(f"{target}: sha256 {h.hexdigest()}, expected {digest}") + + +def number(line, first, last): + """The value in 1-based columns first to last, or None if blank.""" + text = line[first - 1:last].strip() + return float(text) if text else None + + +def propagate(ra, dec, pm_ra, pm_dec, years): + """Moves positions (radians) by their proper motion (mas/yr, right ascension times + cos declination) over `years`, along the sphere.""" + scale = np.radians(1 / 3.6e6) * years + ca, sa, cd, sd = np.cos(ra), np.sin(ra), np.cos(dec), np.sin(dec) + x = cd * ca + scale * (pm_ra * -sa + pm_dec * -sd * ca) + y = cd * sa + scale * (pm_ra * ca + pm_dec * -sd * sa) + z = sd + scale * (pm_dec * cd) + norm = np.sqrt(x * x + y * y + z * z) + return np.mod(np.arctan2(y, x), 2 * np.pi), np.arcsin(z / norm) + + +def table(rows): + columns = list(zip(*rows)) if rows else [[]] * 6 + return {k: np.array(c, dtype=float) for k, c in zip(FIELDS, columns)} + + +def hipparcos(path): + rows, numbers = [], set() + with open(path, encoding="ascii") as f: + for line in f: + v, ra, dec = number(line, 42, 46), number(line, 52, 63), number(line, 65, 76) + if v is None or ra is None or dec is None: + continue + bv = number(line, 246, 251) + rows.append((math.radians(ra), math.radians(dec), number(line, 88, 95) or 0.0, + number(line, 97, 104) or 0.0, v, math.nan if bv is None else bv)) + numbers.add(int(line[8:14])) + stars = table(rows) + stars["ra"], stars["dec"] = propagate( + stars["ra"], stars["dec"], stars["pm_ra"], stars["pm_dec"], 2000.0 - 1991.25) + return stars, numbers + + +def hipparcos_colors(path): + """B - V by Hipparcos number.""" + colors = {} + with open(path, encoding="ascii") as f: + for line in f: + bv = number(line, 246, 251) + if bv is not None: + colors[int(line[8:14])] = bv + return colors + + +def tycho2(paths, hip_numbers): + rows, skipped = [], {"hipparcos": 0, "no VT": 0} + for path in paths: + with gzip.open(path, "rt", encoding="ascii") as f: + for line in f: + hip = line[142:148].strip() + if hip and int(hip) in hip_numbers: + skipped["hipparcos"] += 1 + continue + bt, vt = number(line, 111, 116), number(line, 124, 129) + if vt is None: + skipped["no VT"] += 1 + continue + if line[13] == "X": + ra, dec = number(line, 153, 164), number(line, 166, 177) + pm_ra = pm_dec = 0.0 + else: + ra, dec = number(line, 16, 27), number(line, 29, 40) + pm_ra, pm_dec = number(line, 42, 48) or 0.0, number(line, 50, 56) or 0.0 + if bt is None: + v, bv = vt, math.nan + else: + v, bv = vt - 0.090 * (bt - vt), 0.850 * (bt - vt) + rows.append((math.radians(ra), math.radians(dec), pm_ra, pm_dec, v, bv)) + return table(rows), skipped + + +def gaia_query(query): + """The rows of an ADQL query as a DataFrame, retrying a few times.""" + data = urllib.parse.urlencode({"REQUEST": "doQuery", "LANG": "ADQL", "FORMAT": "csv", + "QUERY": query}).encode() + for attempt in range(6): + try: + with urllib.request.urlopen(GAIA_TAP, data=data, timeout=600) as reply: + return pd.read_csv(io.BytesIO(reply.read())) + except Exception as e: # network or service errors: wait and retry + if attempt == 5: + raise + print(f" retrying after {e}") + time.sleep(10 * (attempt + 1)) + + +def gaia_range(directory, faint, lo, hi): + """Gaia stars brighter than G `faint` with source_id in [lo, hi), from the cache in + `directory` or the archive, split until each query is within the row limit.""" + path = os.path.join(directory, f"{lo}-{hi}.csv") + split = os.path.join(directory, f"{lo}-{hi}.split") + mid = (lo + hi) // 2 + if os.path.exists(path): + return [path] + if os.path.exists(split): + return gaia_range(directory, faint, lo, mid) + gaia_range(directory, faint, mid, hi) + query = ( + "SELECT g.ra, g.dec, g.pmra, g.pmdec, g.phot_g_mean_mag, g.bp_rp, " + "t.source_id AS tycho, h.source_id AS hip " + "FROM gaiadr3.gaia_source_lite AS g " + "LEFT OUTER JOIN gaiadr3.tycho2tdsc_merge_best_neighbour AS t ON g.source_id = t.source_id " + "LEFT OUTER JOIN gaiadr3.hipparcos2_best_neighbour AS h ON g.source_id = h.source_id " + f"WHERE g.phot_g_mean_mag < {faint} AND g.source_id >= {lo} AND g.source_id < {hi}" + ) + rows = gaia_query(query) + if len(rows) >= GAIA_ROWS: + # Too many for one reply: split, and remember it, so a rerun does not ask again. + open(split, "w").close() + return gaia_range(directory, faint, lo, mid) + gaia_range(directory, faint, mid, hi) + rows.to_csv(path + ".part", index=False) + os.replace(path + ".part", path) + return [path] + + +def gaia(source, faint): + directory = os.path.join(source, f"gaia-g{faint:g}") + os.makedirs(directory, exist_ok=True) + chunks = 2048 + edges = [CELLS * i // chunks * CELL for i in range(chunks + 1)] + paths = [] + with concurrent.futures.ThreadPoolExecutor(4) as pool: + jobs = [pool.submit(gaia_range, directory, faint, lo, hi) for lo, hi in zip(edges, edges[1:])] + for done, job in enumerate(concurrent.futures.as_completed(jobs), 1): + paths += job.result() + if done % 100 == 0: + print(f" Gaia: {done} of {chunks} ranges") + return pd.concat((pd.read_csv(p) for p in sorted(paths)), ignore_index=True) + + +def gaia_color_fit(source, hip_main): + """A cubic giving B - V from BP - RP, fitted to the Hipparcos stars Gaia matched, and + the scatter of the fit.""" + path = os.path.join(source, "gaia-hipparcos.csv") + if not os.path.exists(path): + parts = [] + edges = [CELLS * i // 8 * CELL for i in range(9)] + for lo, hi in zip(edges, edges[1:]): + parts.append(gaia_query( + "SELECT h.original_ext_source_id AS hip, g.bp_rp " + "FROM gaiadr3.hipparcos2_best_neighbour AS h " + "JOIN gaiadr3.gaia_source_lite AS g ON h.source_id = g.source_id " + f"WHERE g.source_id >= {lo} AND g.source_id < {hi}" + )) + pd.concat(parts).to_csv(path, index=False) + matched = pd.read_csv(path).dropna() + colors = hipparcos_colors(hip_main) + c = matched["bp_rp"].to_numpy() + bv = np.array([colors.get(int(h), np.nan) for h in matched["hip"]]) + ok = np.isfinite(bv) & (c > -0.5) & (c < 3.5) + coefficients = np.polyfit(c[ok], bv[ok], 3) + residual = bv[ok] - np.polyval(coefficients, c[ok]) + return coefficients, np.median(np.abs(residual)) * 1.4826, ok.sum() + + +def gaia_stars(rows, coefficients): + """Gaia rows that neither Hipparcos nor Tycho-2 has, as V and B - V at J2000.""" + rows = rows[rows["tycho"].isna() & rows["hip"].isna()] + c = rows["bp_rp"].to_numpy() + known = np.isfinite(c) + c = np.where(known, c, 0.82) + cc = np.clip(c, -0.5, 2.75) + g_minus_v = -0.02704 + 0.01424 * cc - 0.2156 * cc**2 + 0.01426 * cc**3 + pm_ra = np.nan_to_num(rows["pmra"].to_numpy()) + pm_dec = np.nan_to_num(rows["pmdec"].to_numpy()) + ra, dec = propagate(np.radians(rows["ra"].to_numpy()), np.radians(rows["dec"].to_numpy()), + pm_ra, pm_dec, 2000.0 - 2016.0) + return { + "ra": ra, "dec": dec, "pm_ra": pm_ra, "pm_dec": pm_dec, + "v": rows["phot_g_mean_mag"].to_numpy() - g_minus_v, + "bv": np.where(known, np.polyval(coefficients, np.clip(c, -0.5, 3.5)), np.nan), + } + + +def unit(stars): + cd = np.cos(stars["dec"]) + return np.stack([cd * np.cos(stars["ra"]), cd * np.sin(stars["ra"]), np.sin(stars["dec"])], 1) + + +def drop_near(stars, others, radius): + """`stars` without those within `radius` radians of any of `others`, and how many + were dropped: cells of the radius on a side, searched with their neighbors.""" + key = lambda c: ((c[:, 0] + (1 << 20)) << 42) | ((c[:, 1] + (1 << 20)) << 21) | (c[:, 2] + (1 << 20)) + points = unit(stars) + keys = key(np.floor(points / radius).astype(np.int64)) + order = np.argsort(keys) + # A key no cell has ends the list, so a walk past the last star stops there. + sorted_keys = np.append(keys[order], np.iinfo(np.int64).max) + near = unit(others) + cells = np.floor(near / radius).astype(np.int64) + drop = np.zeros(len(points), dtype=bool) + for offset in np.ndindex(3, 3, 3): + k = key(cells + np.array(offset) - 1) + i = np.searchsorted(sorted_keys, k) + # Every star in the cell: walk forward while the key repeats. + while True: + hit = sorted_keys[i] == k + if not hit.any(): + break + j = order[i[hit]] + close = np.linalg.norm(points[j] - near[hit], axis=1) < radius + drop[j[close]] = True + i = i + hit + return {f: stars[f][~drop] for f in FIELDS}, int(drop.sum()) + + +def write(path, stars, faint): + keep = stars["v"] <= faint + order = np.argsort(stars["v"][keep], kind="stable") + s = {k: stars[k][keep][order] for k in FIELDS} + record = np.dtype([("ra", " 0: + coefficients, scatter, used = gaia_color_fit(source, hip_main) + print(f"B - V from BP - RP: cubic {np.round(coefficients, 5).tolist()} " + f"(highest power first), {scatter:.3f} mag scatter over {used} Hipparcos stars") + rows = gaia(source, args.gaia) + extra = gaia_stars(rows, coefficients) + print(f"Gaia DR3 to G {args.gaia:g}: {len(rows)} stars, {len(extra['v'])} not in " + f"Hipparcos or Tycho-2 by Gaia's cross-match") + # Gaia's cross-match misses about 1 % of Tycho-2: drop Gaia stars within 1 + # arcsecond of one, which Tycho's and Hipparcos's photometry, unable to resolve + # pairs that close, already includes. + others = {k: np.concatenate([hip[k], tyc[k]]) for k in FIELDS} + extra, dropped = drop_near(extra, others, np.radians(1 / 3600)) + print(f" {dropped} more within 1 arcsecond of a Hipparcos or Tycho-2 star") + parts.append(extra) + stars = {k: np.concatenate([p[k] for p in parts]) for k in FIELDS} + out = os.path.join(root, "stars.bin") + count, faintest = write(out, stars, args.faint) + print(f"{count} stars to V {faintest:.2f} in {out}") + + +if __name__ == "__main__": + main() diff --git a/scripts/make-example-scenes.py b/scripts/make-example-scenes.py index d64602e..ba55bd7 100755 --- a/scripts/make-example-scenes.py +++ b/scripts/make-example-scenes.py @@ -1,8 +1,10 @@ #!/usr/bin/env python3 """Writes examples/scenes/.toml for the NASA models fetch-models.sh downloads, -and two scenes with Earth, the Moon and a spacecraft together (below): +two scenes with Earth, the Moon and a spacecraft together (below): examples/scenes/earthrise.toml from lunar orbit and examples/scenes/moonrise.toml -from geostationary height. +from geostationary height, and scenes for the stars of camera tier C1: a star field, +examples/scenes/stars-pleiades.toml, and a telescope on one geostationary satellite +imaging another along the arc, examples/scenes/geo-telescope-*.toml. Every scene has the same composition: the spacecraft in an isometric view (the camera looks along its body (+1, +1, +1) diagonal, body +Z toward image up) with @@ -447,6 +449,135 @@ attitude = {fmt(quaternion(attitude), 9)} """ +# Earth's orbital speed, m/s, and the geostationary orbital speed. +EARTH_SPEED = 29.78e3 +GEO_SPEED = 3074.7 +# The telescope's longitude on the geostationary arc, degrees: looking along the arc +# from here points at right ascension 103 degrees, where the Milky Way crosses the +# equator in Monoceros, the densest stretch of the equator (190 stars per square degree +# to V 12, against 50 where it crosses in Aquila, behind the dust of the Aquila Rift). +GEO_LONGITUDE = 13.0 +# The ecliptic's north pole in the inertial (equatorial J2000) frame. +ECLIPTIC_POLE = [0.0, -math.sin(math.radians(23.4393)), math.cos(math.radians(23.4393))] + + +def earth_velocity(sun): + """Earth's velocity around the Sun when the Sun is in direction `sun`: along the + ecliptic, a quarter turn ahead of the Sun's direction from Earth.""" + return [EARTH_SPEED * x for x in norm(cross(ECLIPTIC_POLE, [-x for x in sun]))] + + +def stars_pleiades(): + # The Pleiades, right ascension 56.75 and declination 24.12 degrees, from + # geostationary orbit, with the Sun behind the camera. + ra, dec = math.radians(56.75), math.radians(24.12) + look = [math.cos(dec) * math.cos(ra), math.cos(dec) * math.sin(ra), math.sin(dec)] + camera = [GEO, 0.0, 0.0] + sun = [-x for x in look] + fmt = lambda v, d=2: "[" + ", ".join(f"{x:.{d}f}" for x in v) + "]" + return f"""# A star field: the Pleiades through the example camera (12 by 7 degrees), from +# geostationary orbit, long-exposed at 5 s. Stars are drawn at lighting L1 with camera +# C1, from the catalog scripts/fetch-stars.py writes. Alcyone, the brightest at V 2.9, +# gives about 200,000 electrons, ten times the full well; a V 8 star about 2,000 and a +# V 11 star about 130. Generated by scripts/make-example-scenes.py. +# +# resolved render --scene examples/scenes/stars-pleiades.toml --lighting L1 --camera C1 --pretty --preview pleiades.png + +sun = {fmt([AU * x for x in sun], 0)} +moon = [0, -384400000, 0] +earth_velocity = {fmt(earth_velocity(sun))} + +[camera] +exposure = 5 + +[observer] +position = {fmt(camera)} +look_at = {fmt([camera[i] + 1e9 * look[i] for i in range(3)])} +up = [0, 0, 1] +velocity = [0.00, {GEO_SPEED:.2f}, 0.00] + +# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh). +[appearance] +earth_month = 7 +""" + + +def geo_telescope(range_m, phase_deg, f_number, exposure, description): + # The camera at GEO_LONGITUDE on the geostationary arc, the target ahead along it. + file, scale, center = "ssl1300.glb", 0.1, MODELS[3][3] + lon = math.radians(GEO_LONGITUDE) + camera = [GEO * math.cos(lon), GEO * math.sin(lon), 0.0] + angle = lon + range_m / GEO + spot = [GEO * math.cos(angle), GEO * math.sin(angle), 0.0] + velocity = [-GEO_SPEED * math.sin(lon), GEO_SPEED * math.cos(lon), 0.0] + look = norm([spot[i] - camera[i] for i in range(3)]) + image_up = [0.0, 0.0, 1.0] + attitude = triad([1, 1, 1], [0, 0, 1], [-x for x in look], image_up) + offset = mul(attitude, center) + position = [spot[i] - offset[i] for i in range(3)] + # The Sun `phase_deg` from the camera's direction as the target sees it, turned + # about the image's up, 10 degrees out of the equator so the target is not shadowed + # by Earth at any phase. + back = [-x for x in look] + side = cross(image_up, back) + p, t = math.radians(phase_deg), math.radians(10.0) + sun = norm([ + math.cos(t) * (math.cos(p) * back[i] + math.sin(p) * side[i]) + math.sin(t) * image_up[i] + for i in range(3) + ]) + fmt = lambda v, d=2: "[" + ", ".join(f"{x:.{d}f}" for x in v) + "]" + heading = textwrap.fill(description, 86, initial_indent="# ", subsequent_indent="# ") + return f"""{heading} +# A 1 m, f/{f_number:.0f} telescope ({100 / f_number:g} cm aperture, 5.5 urad per pixel, 0.6 by 0.34 +# degrees) on one geostationary satellite, imaging an SSL-1300 {range_m / 1000:.0f} km ahead along +# the arc, lit at a phase of {phase_deg} degrees, against the stars toward the Milky Way in +# Monoceros. Stars are drawn at lighting L1 with camera C1 (scripts/fetch-stars.py). +# Generated by scripts/make-example-scenes.py; the model comes from +# scripts/fetch-models.sh. +# +# resolved render --scene --lighting L1 --camera C1 --pretty --preview geo.png + +sun = {fmt([AU * x for x in sun], 0)} +moon = [0, -384400000, 0] +earth_velocity = {fmt(earth_velocity(sun))} + +[camera] +focal_length = 1.0 +f_number = {f_number:.1f} +exposure = {exposure} + +[observer] +position = {fmt(camera)} +look_at = {fmt(spot)} +up = [0, 0, 1] +velocity = {fmt(velocity)} + +# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh), and lights +# the target's Earth-facing side. +[appearance] +earth_month = 7 + +[[target]] +model = "../models/nasa/{file}" +scale = {scale} +position = {fmt(position, 4)} +attitude = {fmt(quaternion(attitude), 9)} +""" + + +GEO_TELESCOPE = [ + ("geo-telescope-100km", 100e3, 30, 8.0, 6e-4, + "GEO to GEO at 100 km: the target about 55 px across, lit nearly face on, exposed " + "short of saturating it, with only the brighter stars behind it."), + ("geo-telescope-100km-side", 100e3, 90, 8.0, 6e-4, + "GEO to GEO at 100 km, lit from the side: half the target in shadow."), + ("geo-telescope-1000km", 1000e3, 30, 4.0, 0.2, + "GEO to GEO at 1000 km, exposed 0.2 s at f/4 as a surveillance camera would: the " + "target a few pixels across and far past saturation, in a field of stars to the " + "catalog's limit, about V 15, which gives a few DN."), +] + + out = pathlib.Path(__file__).resolve().parent.parent / "examples" / "scenes" for model in MODELS: path = out / f"{model[0]}.toml" @@ -458,3 +589,10 @@ print(f"wrote {path}") path = out / "moonrise.toml" path.write_text(moonrise()) print(f"wrote {path}") +path = out / "stars-pleiades.toml" +path.write_text(stars_pleiades()) +print(f"wrote {path}") +for name, range_m, phase, f_number, exposure, description in GEO_TELESCOPE: + path = out / f"{name}.toml" + path.write_text(geo_telescope(range_m, phase, f_number, exposure, description)) + print(f"wrote {path}")