1use crate::daylight::Daylight;
2use std::f32::consts::PI;
3
4pub const SUN_RADIUS: f32 = 0.004_675;
5
6pub fn adapt_exposure(previous: f32, luminances: &[f32], frame_seconds: f32) -> f32 {
7 let previous = if previous.is_finite() && previous > 0.0 {
8 previous
9 } else {
10 1.0
11 };
12 let mut log_sum = 0.0_f64;
13 let mut count = 0;
14 for &luminance in luminances {
15 if luminance.is_finite() && luminance >= 0.0 {
16 log_sum += luminance.max(0.001).ln() as f64;
17 count += 1;
18 }
19 }
20 if count == 0 || !frame_seconds.is_finite() || frame_seconds <= 0.0 {
21 return previous;
22 }
23 let average = (log_sum / count as f64).exp() as f32;
24 let target = (0.18 / average).clamp(0.25, 8.0);
25 let blend = 1.0 - (-frame_seconds / 0.5).exp();
26 previous + (target - previous) * blend
27}
28
29#[derive(Clone, Copy, Debug)]
30pub struct AnalyticSky {
31 pub sun: [f32; 3],
32 pub sun_colour: [f32; 3],
33 pub sun_intensity: f32,
34 pub ambient: f32,
35 pub turbidity: f32,
36 pub ground_albedo: [f32; 3],
37}
38
39impl AnalyticSky {
40 pub fn new(
41 daylight: Daylight,
42 reference_hour: f64,
43 turbidity: f32,
44 ground_albedo: [f32; 3],
45 ) -> Self {
46 let sun = daylight.sun();
47 let light = daylight.light(reference_hour);
48 Self {
49 sun: sun.y_up.map(|v| v as f32),
50 sun_colour: light.colour.map(|v| v as f32),
51 sun_intensity: light.intensity as f32,
52 ambient: light.ambient as f32,
53 turbidity: turbidity.clamp(1.0, 10.0),
54 ground_albedo: ground_albedo.map(|v| v.clamp(0.0, 1.0)),
55 }
56 }
57
58 pub fn diffuse(&self, direction: [f32; 3]) -> [f32; 3] {
59 let direction = unit(direction);
60 if direction[1] < 0.0 {
61 return self.ground_albedo.map(|v| v * self.ambient * 0.12);
62 }
63 let cos_theta = direction[1].max(0.001);
64 let cos_gamma = dot(direction, self.sun).clamp(-1.0, 1.0);
65 let gamma = cos_gamma.acos();
66 let t = self.turbidity;
67 let a = 0.1787 * t - 1.4630;
68 let b = -0.3554 * t + 0.4275;
69 let c = -0.0227 * t + 5.3251;
70 let d = 0.1206 * t - 2.5771;
71 let e = -0.0670 * t + 0.3703;
72 let shape = (1.0 + a * (b / cos_theta).exp())
73 * (1.0 + c * (d * gamma).exp() + e * cos_gamma * cos_gamma);
74 let zenith_gamma = self.sun[1].clamp(-1.0, 1.0).acos();
75 let zenith_cos = zenith_gamma.cos();
76 let norm = (1.0 + a * b.exp())
77 * (1.0 + c * (d * zenith_gamma).exp() + e * zenith_cos * zenith_cos);
78 let haze = (t - 1.0) / 9.0;
79 let warm = (1.0 - self.sun[1].max(0.0)).powi(2) * (-gamma * 2.0).exp();
80 let base = [0.34 + 0.24 * haze, 0.50 + 0.18 * haze, 1.05 - 0.30 * haze];
81 let tint = [1.0 + 0.75 * warm, 1.0 + 0.15 * warm, 1.0 - 0.45 * warm];
82 let rise = self.sun[1].max(0.0);
83 let twilight = ((rise - 0.2) / 0.4).clamp(0.0, 1.0);
84 let twilight = twilight * twilight * (3.0 - 2.0 * twilight);
85 let scale = 0.215 * rise.powf(2.5) + 0.05 * (1.0 - twilight);
86 let strength = (shape / norm.max(0.01)).max(0.0) * self.ambient * scale;
87 [0, 1, 2].map(|i| base[i] * tint[i] * strength)
88 }
89
90 pub fn radiance(&self, direction: [f32; 3]) -> [f32; 3] {
91 let direction = unit(direction);
92 let mut value = self.diffuse(direction);
93 let cosine = dot(direction, self.sun).clamp(-1.0, 1.0);
94 let edge = SUN_RADIUS.cos();
95 if self.sun[1] > 0.0 && cosine >= edge {
96 let radius = (2.0 * (1.0 - cosine)).sqrt() / SUN_RADIUS;
97 let limb = 0.4 + 0.6 * (1.0 - radius * radius).max(0.0).sqrt();
98 let radiance = self.sun_intensity * limb / (PI * SUN_RADIUS * SUN_RADIUS * 0.8);
99 for (channel, solar) in value.iter_mut().zip(self.sun_colour) {
100 *channel += solar * radiance;
101 }
102 }
103 value
104 }
105}
106
107fn dot(a: [f32; 3], b: [f32; 3]) -> f32 {
108 a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
109}
110
111fn unit(a: [f32; 3]) -> [f32; 3] {
112 let scale = 1.0 / dot(a, a).sqrt().max(1e-8);
113 a.map(|v| v * scale)
114}
115
116pub const SKY_MODEL_WGSL: &str = r#"
117struct SkyUniform {
118 sun: vec4f,
119 sun_colour_intensity: vec4f,
120 ambient_turbidity_albedo: vec4f,
121 albedo_mode: vec4f,
122 forward: vec4f,
123 right: vec4f,
124 up: vec4f,
125};
126@group(1) @binding(0) var<uniform> sky: SkyUniform;
127@group(1) @binding(1) var hdr: texture_2d<f32>;
128@group(1) @binding(2) var hdr_sampler: sampler;
129fn sky_diffuse(direction: vec3f) -> vec3f {
130 let d = normalize(direction);
131 if d.y < 0.0 {
132 return vec3f(sky.ambient_turbidity_albedo.zw, sky.albedo_mode.x) * sky.ambient_turbidity_albedo.x * 0.12;
133 }
134 let cos_theta = max(d.y, 0.001);
135 let cos_gamma = clamp(dot(d, sky.sun.xyz), -1.0, 1.0);
136 let gamma = acos(cos_gamma);
137 let t = sky.ambient_turbidity_albedo.y;
138 let a = 0.1787*t-1.4630;
139 let b = -0.3554*t+0.4275;
140 let c = -0.0227*t+5.3251;
141 let dd = 0.1206*t-2.5771;
142 let e = -0.0670*t+0.3703;
143 let shape = (1.0+a*exp(b/cos_theta))*(1.0+c*exp(dd*gamma)+e*cos_gamma*cos_gamma);
144 let zenith_gamma = acos(clamp(sky.sun.y, -1.0, 1.0));
145 let zenith_cos = cos(zenith_gamma);
146 let norm = (1.0+a*exp(b))*(1.0+c*exp(dd*zenith_gamma)+e*zenith_cos*zenith_cos);
147 let haze = (t-1.0)/9.0;
148 let warm = pow(1.0-max(sky.sun.y, 0.0), 2.0)*exp(-gamma*2.0);
149 let base = vec3f(0.34+0.24*haze, 0.50+0.18*haze, 1.05-0.30*haze);
150 let tint = vec3f(1.0+0.75*warm, 1.0+0.15*warm, 1.0-0.45*warm);
151 let rise = max(sky.sun.y, 0.0);
152 let scale = 0.215*pow(rise, 2.5)+0.05*(1.0-smoothstep(0.2, 0.6, rise));
153 return base*tint*max(shape/max(norm, 0.01), 0.0)*sky.ambient_turbidity_albedo.x*scale;
154}
155fn sky_analytic(direction: vec3f) -> vec3f {
156 let d = normalize(direction);
157 var value = sky_diffuse(d);
158 let cosine = clamp(dot(d, sky.sun.xyz), -1.0, 1.0);
159 let radius = 0.004675;
160 if sky.sun.y > 0.0 && cosine >= cos(radius) {
161 let radial = sqrt(2.0*(1.0-cosine))/radius;
162 let limb = 0.4+0.6*sqrt(max(1.0-radial*radial, 0.0));
163 value += sky.sun_colour_intensity.xyz*(sky.sun_colour_intensity.w*limb/(3.14159265*radius*radius*0.8));
164 }
165 return value;
166}
167fn sky_radiance(direction: vec3f) -> vec3f {
168 let d = normalize(direction);
169 if sky.albedo_mode.y < 0.5 {
170 return sky_analytic(d);
171 }
172 let u = fract(atan2(d.x, -d.z)/6.28318531+0.5);
173 let v = acos(clamp(d.y, -1.0, 1.0))/3.14159265;
174 return textureSampleLevel(hdr, hdr_sampler, vec2f(u, v), 0.0).rgb;
175}
176"#;
177
178#[cfg(test)]
179mod tests {
180 use super::*;
181 use crate::daylight::REFERENCE_HOUR;
182
183 #[test]
184 fn exposure_uses_the_frame_log_average_and_caller_time() {
185 let first = adapt_exposure(1.0, &[0.04, 0.16], 0.5);
186 let target = 0.18 / 0.08;
187 let expected = 1.0 + (target - 1.0) * (1.0 - (-1.0_f32).exp());
188 assert!((first - expected).abs() < 1e-6);
189 assert_eq!(adapt_exposure(first, &[0.04, 0.16], 0.0), first);
190 assert_eq!(adapt_exposure(first, &[f32::NAN], 1.0), first);
191 assert_eq!(adapt_exposure(1.0, &[0.04, 0.16], 0.5), first);
192 }
193
194 #[test]
195 fn model_keeps_the_solar_disc_and_calibrates_diffuse() {
196 let expected = [
197 ([0.0237, 0.0326, 0.0607], [0.1190, 0.1560, 0.2770]),
198 ([0.0684, 0.0952, 0.1792], [0.0506, 0.0706, 0.1330]),
199 ([0.0237, 0.0326, 0.0607], [0.0262, 0.0366, 0.0690]),
200 ];
201 for (hour, (zenith, side)) in [8.0, 12.0, 16.0].into_iter().zip(expected) {
202 let sky = AnalyticSky::new(
203 Daylight {
204 hour,
205 day: 172.0,
206 latitude: 45.0,
207 heading: 180.0,
208 },
209 REFERENCE_HOUR,
210 2.5,
211 [0.2, 0.25, 0.2],
212 );
213 for (actual, reference) in [
214 (sky.diffuse([0.0, 1.0, 0.0]), zenith),
215 (sky.diffuse([1.0, 0.5, 0.0]), side),
216 ] {
217 for (value, expected) in actual.into_iter().zip(reference) {
218 assert!(
219 (value - expected).abs() <= 0.001,
220 "{hour}: {value} != {expected}"
221 );
222 }
223 }
224 let direct = sky.radiance(sky.sun);
225 let diffuse = sky.diffuse(sky.sun);
226 for channel in 0..3 {
227 let expected = sky.sun_colour[channel] * sky.sun_intensity
228 / (PI * SUN_RADIUS * SUN_RADIUS * 0.8);
229 assert!((direct[channel] - diffuse[channel] - expected).abs() < 0.01);
230 }
231 }
232 }
233
234 #[test]
235 fn clear_afternoon_has_a_sun_dominated_horizontal_irradiance() {
236 fn luminance(rgb: [f32; 3]) -> f64 {
237 rgb.into_iter()
238 .zip([0.2126, 0.7152, 0.0722])
239 .map(|(channel, weight)| channel as f64 * weight)
240 .sum()
241 }
242
243 fn sky_share(hour: f64) -> f64 {
244 let sky = AnalyticSky::new(
245 Daylight {
246 hour,
247 day: 172.0,
248 latitude: 45.0,
249 heading: 180.0,
250 },
251 16.0,
252 2.5,
253 [0.2, 0.25, 0.2],
254 );
255 let mut diffuse = 0.0;
256 for row in 0..96 {
257 let theta = (row as f64 + 0.5) * std::f64::consts::FRAC_PI_2 / 96.0;
258 for column in 0..192 {
259 let phi = (column as f64 + 0.5) * std::f64::consts::TAU / 192.0;
260 let direction = [
261 (theta.sin() * phi.cos()) as f32,
262 theta.cos() as f32,
263 (theta.sin() * phi.sin()) as f32,
264 ];
265 let solid_angle = theta.sin() * std::f64::consts::FRAC_PI_2 / 96.0
266 * std::f64::consts::TAU
267 / 192.0;
268 diffuse += luminance(sky.diffuse(direction)) * theta.cos() * solid_angle;
269 }
270 }
271 let direct =
272 luminance(sky.sun_colour) * sky.sun_intensity as f64 * sky.sun[1].max(0.0) as f64;
273 diffuse / (diffuse + direct)
274 }
275
276 for hour in [10.0, 12.0, 14.0, 16.0] {
277 let share = sky_share(hour);
278 assert!((0.15..=0.25).contains(&share), "{hour}: {share}");
279 }
280 assert!(sky_share(6.0) > sky_share(8.0));
281 assert!(sky_share(19.0) > sky_share(16.0));
282 }
283}