proof_engine/relativistic/
lensing.rs

1//! Gravitational lensing.
2
3use glam::{Vec2, Vec3, Vec4};
4
5/// Deflection angle for a light ray passing a mass M at impact parameter b.
6/// alpha = 4GM / (b * c^2)
7#[allow(non_snake_case)]
8pub fn deflection_angle(mass: f64, impact_param: f64, c: f64, G: f64) -> f64 {
9    if impact_param.abs() < 1e-30 {
10        return std::f64::consts::PI; // captured
11    }
12    4.0 * G * mass / (impact_param * c * c)
13}
14
15/// Einstein radius: the angular radius of a perfect ring image.
16/// theta_E = sqrt(4GM / c^2 * D_ls / (D_l * D_s))
17#[allow(non_snake_case)]
18pub fn einstein_radius(mass: f64, d_lens: f64, d_source: f64, c: f64, G: f64) -> f64 {
19    let d_ls = d_source - d_lens;
20    if d_lens <= 0.0 || d_source <= 0.0 || d_ls <= 0.0 {
21        return 0.0;
22    }
23    (4.0 * G * mass * d_ls / (c * c * d_lens * d_source)).sqrt()
24}
25
26/// Lens equation: theta - beta = theta_E^2 / theta
27/// Returns theta - beta - theta_E^2/theta (should be zero for image positions).
28#[allow(non_snake_case)]
29pub fn lens_equation(theta: f64, beta: f64, theta_E: f64) -> f64 {
30    if theta.abs() < 1e-30 {
31        return f64::INFINITY;
32    }
33    theta - beta - theta_E * theta_E / theta
34}
35
36/// Image positions for a point lens. Returns two image angles (theta_+, theta_-).
37/// theta_+/- = (beta +/- sqrt(beta^2 + 4*theta_E^2)) / 2
38#[allow(non_snake_case)]
39pub fn image_positions(beta: f64, theta_E: f64) -> (f64, f64) {
40    let disc = (beta * beta + 4.0 * theta_E * theta_E).sqrt();
41    let theta_plus = (beta + disc) / 2.0;
42    let theta_minus = (beta - disc) / 2.0;
43    (theta_plus, theta_minus)
44}
45
46/// Magnification of an image at angle theta for Einstein radius theta_E.
47/// mu = (theta / theta_E)^4 / ((theta/theta_E)^4 - 1)
48/// Equivalently: mu = |theta / beta * d(theta)/d(beta)|
49#[allow(non_snake_case)]
50pub fn magnification(theta: f64, theta_E: f64) -> f64 {
51    if theta_E.abs() < 1e-30 {
52        return 1.0;
53    }
54    let u = theta / theta_E;
55    let u4 = u * u * u * u;
56    if (u4 - 1.0).abs() < 1e-15 {
57        return f64::INFINITY; // on the Einstein ring (caustic)
58    }
59    (u4 / (u4 - 1.0)).abs()
60}
61
62/// 2D grid of deflection vectors around a point mass.
63#[derive(Debug, Clone)]
64pub struct LensingField {
65    pub width: usize,
66    pub height: usize,
67    pub center: Vec2,
68    pub mass: f64,
69    pub c: f64,
70    pub g_const: f64,
71    /// Deflection vectors at each grid point.
72    pub deflections: Vec<Vec2>,
73    pub cell_size: f32,
74}
75
76impl LensingField {
77    #[allow(non_snake_case)]
78    pub fn new(width: usize, height: usize, center: Vec2, mass: f64, c: f64, G: f64, cell_size: f32) -> Self {
79        let mut deflections = Vec::with_capacity(width * height);
80        for iy in 0..height {
81            for ix in 0..width {
82                let x = (ix as f32 - width as f32 / 2.0) * cell_size + center.x;
83                let y = (iy as f32 - height as f32 / 2.0) * cell_size + center.y;
84                let pos = Vec2::new(x, y);
85                let r = pos - center;
86                let dist = r.length() as f64;
87                if dist < 1e-10 {
88                    deflections.push(Vec2::ZERO);
89                    continue;
90                }
91                let alpha = deflection_angle(mass, dist, c, G);
92                let dir = r.normalize();
93                deflections.push(dir * alpha as f32);
94            }
95        }
96        Self {
97            width,
98            height,
99            center,
100            mass,
101            c,
102            g_const: G,
103            deflections,
104            cell_size,
105        }
106    }
107
108    /// Get deflection at a grid index.
109    pub fn get(&self, ix: usize, iy: usize) -> Vec2 {
110        if ix < self.width && iy < self.height {
111            self.deflections[iy * self.width + ix]
112        } else {
113            Vec2::ZERO
114        }
115    }
116
117    /// Interpolate deflection at an arbitrary position.
118    pub fn sample(&self, pos: Vec2) -> Vec2 {
119        let r = pos - self.center;
120        let dist = r.length() as f64;
121        if dist < 1e-10 {
122            return Vec2::ZERO;
123        }
124        let alpha = deflection_angle(self.mass, dist, self.c, self.g_const);
125        let dir = r.normalize();
126        dir * alpha as f32
127    }
128
129    /// Total magnification at a position.
130    pub fn magnification_at(&self, pos: Vec2) -> f64 {
131        let r = (pos - self.center).length() as f64;
132        let rs = 4.0 * self.g_const * self.mass / (self.c * self.c);
133        let theta_E = (rs).sqrt(); // simplified
134        if theta_E < 1e-30 {
135            return 1.0;
136        }
137        magnification(r, theta_E)
138    }
139}
140
141/// Apply gravitational lensing to glyph positions.
142/// Deflects each position away from the lens center.
143#[allow(non_snake_case)]
144pub fn apply_lensing(
145    glyph_positions: &[Vec2],
146    lens_center: Vec2,
147    lens_mass: f64,
148    c: f64,
149    G: f64,
150) -> Vec<Vec2> {
151    glyph_positions.iter().map(|pos| {
152        let r = *pos - lens_center;
153        let dist = r.length() as f64;
154        if dist < 1e-10 {
155            return *pos;
156        }
157        let alpha = deflection_angle(lens_mass, dist, c, G);
158        let dir = r.normalize();
159        *pos + dir * alpha as f32
160    }).collect()
161}
162
163/// Renderer for gravitational lensing effects.
164#[derive(Debug, Clone)]
165pub struct LensingRenderer {
166    pub lens_center: Vec2,
167    pub lens_mass: f64,
168    pub c: f64,
169    pub g_const: f64,
170    pub einstein_ring_visible: bool,
171}
172
173impl LensingRenderer {
174    #[allow(non_snake_case)]
175    pub fn new(lens_center: Vec2, lens_mass: f64, c: f64, G: f64) -> Self {
176        Self {
177            lens_center,
178            lens_mass,
179            c,
180            g_const: G,
181            einstein_ring_visible: true,
182        }
183    }
184
185    /// Compute lensed positions for background glyphs.
186    pub fn lens_positions(&self, positions: &[Vec2]) -> Vec<Vec2> {
187        apply_lensing(positions, self.lens_center, self.lens_mass, self.c, self.g_const)
188    }
189
190    /// Compute magnification for each position.
191    pub fn magnifications(&self, positions: &[Vec2]) -> Vec<f64> {
192        let d_lens = 1.0; // normalized
193        let d_source = 2.0;
194        let theta_E = einstein_radius(self.lens_mass, d_lens, d_source, self.c, self.g_const);
195
196        positions.iter().map(|pos| {
197            let theta = (*pos - self.lens_center).length() as f64;
198            magnification(theta, theta_E)
199        }).collect()
200    }
201
202    /// Render: returns (lensed_positions, magnification_factors).
203    pub fn render(&self, positions: &[Vec2]) -> (Vec<Vec2>, Vec<f64>) {
204        let lensed = self.lens_positions(positions);
205        let mags = self.magnifications(positions);
206        (lensed, mags)
207    }
208
209    /// Generate Einstein ring points for visualization.
210    pub fn einstein_ring_points(&self, d_lens: f64, d_source: f64, n_points: usize) -> Vec<Vec2> {
211        let theta_E = einstein_radius(self.lens_mass, d_lens, d_source, self.c, self.g_const);
212        let mut points = Vec::with_capacity(n_points);
213        for i in 0..n_points {
214            let angle = (i as f64 / n_points as f64) * std::f64::consts::TAU;
215            let p = self.lens_center + Vec2::new(
216                (theta_E * angle.cos()) as f32,
217                (theta_E * angle.sin()) as f32,
218            );
219            points.push(p);
220        }
221        points
222    }
223}
224
225/// Microlensing light curve: amplification as a function of time.
226/// Paczynski formula: A(u) = (u^2 + 2) / (u * sqrt(u^2 + 4))
227/// where u = sqrt(u_min^2 + ((t - t0)/t_E)^2)
228pub fn microlensing_lightcurve(u_min: f64, t_E: f64, times: &[f64]) -> Vec<f64> {
229    let t0 = if times.is_empty() {
230        0.0
231    } else {
232        (times[0] + times[times.len() - 1]) / 2.0
233    };
234
235    times.iter().map(|&t| {
236        let tau = (t - t0) / t_E;
237        let u = (u_min * u_min + tau * tau).sqrt();
238        if u < 1e-15 {
239            return f64::INFINITY;
240        }
241        (u * u + 2.0) / (u * (u * u + 4.0).sqrt())
242    }).collect()
243}
244
245/// Total magnification from both images of a point lens.
246/// A_total = (u^2 + 2) / (u * sqrt(u^2 + 4))
247pub fn total_magnification(u: f64) -> f64 {
248    if u < 1e-15 {
249        return f64::INFINITY;
250    }
251    (u * u + 2.0) / (u * (u * u + 4.0).sqrt())
252}
253
254/// Compute the critical curve radius for a singular isothermal sphere lens.
255pub fn sis_critical_radius(velocity_dispersion: f64, c: f64, d_lens: f64, d_source: f64) -> f64 {
256    let d_ls = d_source - d_lens;
257    if d_source <= 0.0 || d_ls <= 0.0 {
258        return 0.0;
259    }
260    4.0 * std::f64::consts::PI * velocity_dispersion * velocity_dispersion * d_ls / (c * c * d_source)
261}
262
263#[cfg(test)]
264mod tests {
265    use super::*;
266
267    const C: f64 = 299_792_458.0;
268    const G: f64 = 6.674e-11;
269
270    #[test]
271    fn test_deflection_angle_sun() {
272        // Sun deflection of starlight: ~1.75 arcseconds
273        let m_sun = 1.989e30;
274        let r_sun = 6.96e8; // solar radius as impact parameter
275        let alpha = deflection_angle(m_sun, r_sun, C, G);
276        let arcsec = alpha * 206265.0; // convert radians to arcseconds
277        assert!(
278            (arcsec - 1.75).abs() < 0.1,
279            "Solar deflection: {} arcsec, expected ~1.75",
280            arcsec
281        );
282    }
283
284    #[test]
285    fn test_deflection_decreases_with_distance() {
286        let m = 1e30;
287        let a1 = deflection_angle(m, 1e8, C, G);
288        let a2 = deflection_angle(m, 1e9, C, G);
289        assert!(a2 < a1, "Deflection should decrease with impact param");
290    }
291
292    #[test]
293    fn test_deflection_approaches_zero() {
294        let m = 1e30;
295        let a = deflection_angle(m, 1e20, C, G);
296        assert!(a < 1e-10, "Deflection at large b should be tiny: {}", a);
297    }
298
299    #[test]
300    fn test_einstein_radius() {
301        let m = 1e30;
302        let d_l = 1e20;
303        let d_s = 2e20;
304        let theta_E = einstein_radius(m, d_l, d_s, C, G);
305        assert!(theta_E > 0.0);
306    }
307
308    #[test]
309    fn test_image_positions() {
310        let theta_E = 1.0;
311        let beta = 0.5;
312        let (tp, tm) = image_positions(beta, theta_E);
313        // theta+ should be positive and > theta_E
314        assert!(tp > 0.0);
315        // theta- should be negative
316        assert!(tm < 0.0);
317        // Verify lens equation
318        let err_p = lens_equation(tp, beta, theta_E);
319        let err_m = lens_equation(tm, beta, theta_E);
320        assert!(err_p.abs() < 1e-10, "Lens equation error +: {}", err_p);
321        assert!(err_m.abs() < 1e-10, "Lens equation error -: {}", err_m);
322    }
323
324    #[test]
325    fn test_einstein_ring() {
326        // When beta = 0 (source directly behind lens), both images merge at theta_E
327        let theta_E = 1.0;
328        let (tp, tm) = image_positions(0.0, theta_E);
329        assert!((tp - theta_E).abs() < 1e-10);
330        assert!((tm + theta_E).abs() < 1e-10); // negative side
331    }
332
333    #[test]
334    fn test_magnification_at_einstein_ring() {
335        // Magnification diverges at the Einstein ring
336        let theta_E = 1.0;
337        let mag = magnification(theta_E, theta_E);
338        assert!(mag.is_infinite() || mag > 1e10, "Magnification at caustic: {}", mag);
339    }
340
341    #[test]
342    fn test_magnification_far_from_lens() {
343        let theta_E = 1.0;
344        let mag = magnification(100.0 * theta_E, theta_E);
345        assert!((mag - 1.0).abs() < 0.01, "Far magnification should be ~1: {}", mag);
346    }
347
348    #[test]
349    fn test_microlensing_lightcurve_peak() {
350        let u_min = 0.5;
351        let t_E = 10.0;
352        let times: Vec<f64> = (-50..=50).map(|i| i as f64).collect();
353        let curve = microlensing_lightcurve(u_min, t_E, &times);
354
355        // Peak should be at the middle (t=0)
356        let mid = curve.len() / 2;
357        assert!(curve[mid] > curve[0], "Peak should be at center");
358        assert!(curve[mid] > curve[curve.len() - 1], "Peak should be at center");
359
360        // All values should be >= 1 (lensing always brightens)
361        for &a in &curve {
362            assert!(a >= 1.0, "Amplification should be >= 1: {}", a);
363        }
364    }
365
366    #[test]
367    fn test_apply_lensing() {
368        let positions = vec![Vec2::new(10.0, 0.0), Vec2::new(0.0, 10.0)];
369        let lensed = apply_lensing(&positions, Vec2::ZERO, 1e30, C, G);
370        // Lensed positions should be deflected outward from center
371        assert!(lensed[0].x > positions[0].x || (lensed[0] - positions[0]).length() > 0.0);
372    }
373
374    #[test]
375    fn test_total_magnification() {
376        // At u=1 (source at Einstein radius): A = 3/sqrt(5) ~ 1.342
377        let a = total_magnification(1.0);
378        let expected = 3.0 / 5.0_f64.sqrt();
379        assert!((a - expected).abs() < 1e-10);
380    }
381
382    #[test]
383    fn test_lensing_field() {
384        let field = LensingField::new(10, 10, Vec2::ZERO, 1e30, C, G, 1.0);
385        assert_eq!(field.deflections.len(), 100);
386        // Center deflection should be zero
387        let center_defl = field.get(5, 5);
388        // Corner should have some deflection
389        let corner_defl = field.get(0, 0);
390        assert!(corner_defl.length() > 0.0 || center_defl.length() >= 0.0);
391    }
392}
proof_engine/relativistic/lensing.rs

proof_engine/relativistic/
lensing.rs
