1pub const FILM_SAMPLES: usize = 24;
2pub const FILM_NU: f32 = 0.0013513514;
3pub const FILM_STEP: f32 = 5.2728312e-05;
4pub const FILM_WEIGHTS: [[f32; 3]; 24] = [
5 [0.00049457996, 6.4814856e-05, 0.00010619715],
6 [0.002711246, -0.0001910701, -0.00037234355],
7 [0.01515967, 0.00028906547, 0.00012810956],
8 [0.06871473, -0.006184027, -0.00132862],
9 [0.19734336, -0.023679964, -0.0021775665],
10 [0.34738973, -0.02425651, -0.0060811765],
11 [0.3639468, 0.030559337, -0.014212847],
12 [0.23666404, 0.12368723, -0.022455476],
13 [0.08203017, 0.19909823, -0.027573477],
14 [-0.02977581, 0.22957213, -0.027457258],
15 [-0.091345824, 0.21851267, -0.020553917],
16 [-0.09862617, 0.16825487, -0.004000803],
17 [-0.071613036, 0.09727391, 0.021227753],
18 [-0.047891952, 0.05012592, 0.0553259],
19 [-0.024773184, 0.0206622, 0.112508],
20 [-0.01048964, -0.0007772626, 0.18015881],
21 [0.006274984, -0.01616651, 0.20803724],
22 [0.017522685, -0.023618592, 0.2038683],
23 [0.015057241, -0.021199455, 0.17673187],
24 [0.0156005435, -0.01456651, 0.10748986],
25 [0.0014033349, -0.00467564, 0.039753698],
26 [0.004500126, -0.0028086475, 0.014122494],
27 [-0.0012337026, 0.00018912496, 0.0036677718],
28 [0.00093538617, -0.00016521492, 0.0030877078],
29];
30pub const FILM_EDGE: f32 = 0.462_664_37;
31
32#[derive(Clone, Copy, Debug, PartialEq)]
33struct Stack {
34 lift: [f32; 2],
35 floor: [f32; 2],
36 cosine: [f32; 2],
37 sine: [f32; 2],
38}
39
40fn amplitudes(cos_i: f32, na: f32, nb: f32) -> [f32; 3] {
41 let sin2 = (na / nb) * (na / nb) * (1.0 - cos_i * cos_i);
42 if sin2 >= 1.0 {
43 return [1.0, 1.0, 0.0];
44 }
45 let cos_t = (1.0 - sin2).sqrt();
46 let rs = (na * cos_i - nb * cos_t) / (na * cos_i + nb * cos_t);
47 let rp = (nb * cos_i - na * cos_t) / (nb * cos_i + na * cos_t);
48 [rs, rp, cos_t]
49}
50
51fn stack(r12: [f32; 2], re: [f32; 2], im: [f32; 2]) -> Stack {
52 let r23 = [re[0] * re[0] + im[0] * im[0], re[1] * re[1] + im[1] * im[1]];
53 let floor = [
54 1.0 + r12[0] * r12[0] * r23[0],
55 1.0 + r12[1] * r12[1] * r23[1],
56 ];
57 Stack {
58 lift: [
59 0.5 * (r12[0] * r12[0] + r23[0] - floor[0]),
60 0.5 * (r12[1] * r12[1] + r23[1] - floor[1]),
61 ],
62 floor,
63 cosine: [2.0 * r12[0] * re[0], 2.0 * r12[1] * re[1]],
64 sine: [-2.0 * r12[0] * im[0], -2.0 * r12[1] * im[1]],
65 }
66}
67
68fn pair(z: [f32; 2], s: &Stack) -> f32 {
69 let x = [
70 s.cosine[0] * z[0] + s.sine[0] * z[1],
71 s.cosine[1] * z[0] + s.sine[1] * z[1],
72 ];
73 s.lift[0] / (s.floor[0] + x[0]) + s.lift[1] / (s.floor[1] + x[1])
74}
75
76fn sum(path: f32, stacks: [Stack; 3]) -> [f32; 3] {
77 let turn = [(path * FILM_STEP).cos(), (path * FILM_STEP).sin()];
78 let mut z = [(path * FILM_NU).cos(), (path * FILM_NU).sin()];
79 let mut total = [1.0_f32; 3];
80 for weight in FILM_WEIGHTS {
81 for channel in 0..3 {
82 total[channel] += weight[channel] * pair(z, &stacks[channel]);
83 }
84 z = [
85 z[0] * turn[0] - z[1] * turn[1],
86 z[0] * turn[1] + z[1] * turn[0],
87 ];
88 }
89 total.map(|v| v.clamp(0.0, 1.0))
90}
91
92pub fn film_wavelengths() -> [f32; FILM_SAMPLES] {
93 std::array::from_fn(|k| 1.0 / (FILM_NU + FILM_STEP * k as f32))
94}
95
96pub fn film_dielectric(cosine: f32, thickness: f32, film_ior: f32, substrate_ior: f32) -> [f32; 3] {
97 let nf = film_ior.max(1.01);
98 let first = amplitudes(cosine.clamp(0.0, 1.0), 1.0, nf);
99 let second = amplitudes(first[2], nf, substrate_ior.max(1.0));
100 if second[2] <= 0.0 {
101 return [1.0; 3];
102 }
103 let layer = stack([first[0], first[1]], [second[0], second[1]], [0.0, 0.0]);
104 sum(
105 4.0 * std::f32::consts::PI * nf * thickness * first[2],
106 [layer; 3],
107 )
108}
109
110pub fn conductor_ior(f0: f32) -> [f32; 2] {
111 let r = f0.clamp(0.0, 0.999);
112 let edge = f0.clamp(0.0, 1.0) + (1.0 - f0.clamp(0.0, 1.0)) * FILM_EDGE;
113 let root = r.sqrt();
114 let n = (1.0 + root) / (1.0 - root) * (1.0 - edge) + (1.0 - r) / (1.0 + r) * edge;
115 let k = ((r * (n + 1.0) * (n + 1.0) - (n - 1.0) * (n - 1.0)) / (1.0 - r))
116 .max(0.0)
117 .sqrt();
118 [n, k]
119}
120
121fn mul(a: [f32; 2], b: [f32; 2]) -> [f32; 2] {
122 [a[0] * b[0] - a[1] * b[1], a[0] * b[1] + a[1] * b[0]]
123}
124
125fn div(a: [f32; 2], b: [f32; 2]) -> [f32; 2] {
126 let size = b[0] * b[0] + b[1] * b[1];
127 [
128 (a[0] * b[0] + a[1] * b[1]) / size,
129 (a[1] * b[0] - a[0] * b[1]) / size,
130 ]
131}
132
133fn root(a: [f32; 2]) -> [f32; 2] {
134 let size = (a[0] * a[0] + a[1] * a[1]).sqrt();
135 let re = (0.5 * (size + a[0])).max(0.0).sqrt();
136 let im = (0.5 * (size - a[0])).max(0.0).sqrt();
137 [re, if a[1] < 0.0 { -im } else { im }]
138}
139
140fn metal_stack(r12: [f32; 2], cos_f: f32, nf: f32, f0: f32) -> Stack {
141 let metal = conductor_ior(f0);
142 let square = mul(metal, metal);
143 let u = root([square[0] - nf * nf * (1.0 - cos_f * cos_f), square[1]]);
144 let across = [nf * cos_f, 0.0];
145 let rs = div([across[0] - u[0], -u[1]], [across[0] + u[0], u[1]]);
146 let tilted = [square[0] * cos_f, square[1] * cos_f];
147 let rp = div(
148 [tilted[0] - nf * u[0], tilted[1] - nf * u[1]],
149 [tilted[0] + nf * u[0], tilted[1] + nf * u[1]],
150 );
151 stack(r12, [rs[0], rp[0]], [rs[1], rp[1]])
152}
153
154pub fn film_conductor(cosine: f32, thickness: f32, film_ior: f32, f0: [f32; 3]) -> [f32; 3] {
155 let nf = film_ior.max(1.01);
156 let first = amplitudes(cosine.clamp(0.0, 1.0), 1.0, nf);
157 let r12 = [first[0], first[1]];
158 sum(
159 4.0 * std::f32::consts::PI * nf * thickness * first[2],
160 f0.map(|f| metal_stack(r12, first[2], nf, f)),
161 )
162}
163
164#[cfg(test)]
165mod tests {
166 use super::*;
167 use crate::{FILM, fresnel_dielectric};
168
169 #[test]
170 fn a_flat_spectrum_stays_white() {
171 for channel in 0..3 {
172 let total: f64 = FILM_WEIGHTS.iter().map(|w| w[channel] as f64).sum();
173 assert!((total - 1.0).abs() < 1e-5, "{total}");
174 }
175 }
176
177 #[test]
178 fn the_shader_table_is_the_same_table() {
179 let start = FILM.find("const FILM_WEIGHTS").unwrap();
180 let end = start + FILM[start..].find(");").unwrap();
181 let numbers: Vec<f32> = FILM[start..end]
182 .split("vec3f(")
183 .skip(1)
184 .flat_map(|row| {
185 row.split(')')
186 .next()
187 .unwrap()
188 .split(',')
189 .map(|n| n.trim().parse::<f32>().unwrap())
190 .collect::<Vec<_>>()
191 })
192 .collect();
193 let ours: Vec<f32> = FILM_WEIGHTS.iter().flatten().copied().collect();
194 assert_eq!(numbers, ours);
195 assert!(FILM.contains(&format!("const FILM_SAMPLES: u32 = {FILM_SAMPLES}u;")));
196 for (name, value) in [
197 ("FILM_NU", FILM_NU),
198 ("FILM_STEP", FILM_STEP),
199 ("FILM_EDGE", FILM_EDGE),
200 ] {
201 let key = format!("const {name}: f32 = ");
202 let at = FILM.find(&key).unwrap() + key.len();
203 let text = &FILM[at..at + FILM[at..].find(';').unwrap()];
204 assert_eq!(text.parse::<f32>().unwrap(), value, "{name}");
205 }
206 }
207
208 #[test]
209 fn the_wavelengths_run_from_740_to_390_nm_evenly_in_wavenumber() {
210 let waves = film_wavelengths();
211 assert!((waves[0] - 740.0).abs() < 0.01);
212 assert!((waves[FILM_SAMPLES - 1] - 390.0).abs() < 0.01);
213 }
214
215 #[test]
216 fn no_film_is_the_plain_fresnel_of_the_substrate() {
217 for cosine in [1.0, 0.8, 0.5, 0.2] {
218 let bare = film_dielectric(cosine, 0.0, 1.9, 1.4);
219 for value in bare {
220 assert!(
221 (value - fresnel_dielectric(cosine, 1.4)).abs() < 1e-5,
222 "{value}"
223 );
224 }
225 }
226 }
227
228 #[test]
229 fn a_bare_conductor_keeps_its_colour_at_normal_incidence() {
230 let f0 = [0.95, 0.64, 0.54];
231 let bare = film_conductor(1.0, 0.0, 1.6, f0);
232 for channel in 0..3 {
233 assert!((bare[channel] - f0[channel]).abs() < 1e-4, "{bare:?}");
234 }
235 for cosine in [0.8, 0.4, 0.1] {
236 let grazing = film_conductor(cosine, 0.0, 1.6, f0);
237 for channel in 0..3 {
238 let [n, k] = conductor_ior(f0[channel]);
239 let exact = crate::fresnel_conductor(cosine, n, k);
240 assert!(
241 (grazing[channel] - exact).abs() < 1e-4,
242 "{grazing:?} {exact}"
243 );
244 }
245 }
246 }
247
248 #[test]
249 fn the_conductor_ior_is_cycles_f82_edge_on_gulbrandsen() {
250 for f0 in [0.04_f32, 0.3, 0.56, 0.95] {
251 let [n, k] = conductor_ior(f0);
252 let normal = ((n - 1.0) * (n - 1.0) + k * k) / ((n + 1.0) * (n + 1.0) + k * k);
253 assert!((normal - f0).abs() < 1e-4, "{f0} {n} {k}");
254 }
255 let [n, k] = conductor_ior(0.0);
256 assert!((n - 1.0).abs() < 1e-6 && k == 0.0);
257 }
258}