sparse_ir_minipole/minipole/
con_map.rs

1//! Holomorphic maps between the z plane and the unit disk (port of
2//! `mini_pole/con_map.py`).
3//!
4//! Ported from Green-Phys/MiniPole (commit 15e4a54, MIT License,
5//! Copyright (c) 2024 lzphy); see `LICENSE-THIRD-PARTY`.
6
7use num_complex::Complex;
8
9type C64 = Complex<f64>;
10
11/// A holomorphic map `z(w)` from the unit disk onto the z plane cut along
12/// the integration contour.
13pub trait ConMap {
14    /// `z(w)`.
15    fn z(&self, w: C64) -> C64;
16    /// The preimage `w(z)` inside the unit disk.
17    fn w(&self, z: C64) -> C64;
18    /// `dz/dw`.
19    fn dz(&self, w: C64) -> C64;
20}
21
22/// `z = (Δω_h/2)(w - 1/w) + iω_m`: the unit circle onto the segment
23/// `[i(ω_m - Δω_h), i(ω_m + Δω_h)]` (`ConMapGeneric`, `branch_in = True`).
24#[derive(Debug, Clone, Copy, PartialEq)]
25pub struct ConMapGeneric {
26    /// Midpoint `ω_m` of the segment.
27    pub w_m: f64,
28    /// Half-width `Δω_h` of the segment.
29    pub dw_h: f64,
30}
31
32impl ConMap for ConMapGeneric {
33    fn z(&self, w: C64) -> C64 {
34        if w == C64::new(0.0, 0.0) {
35            return C64::new(f64::INFINITY, 0.0);
36        }
37        (w - w.inv()) * (0.5 * self.dw_h) + C64::new(0.0, self.w_m)
38    }
39
40    fn w(&self, z: C64) -> C64 {
41        let x = (z - C64::new(0.0, self.w_m)) / self.dw_h;
42        let mut w = x - (x * x + 1.0).sqrt();
43        if w.norm() > 1.0 {
44            w = x * 2.0 - w;
45        }
46        w
47    }
48
49    fn dz(&self, w: C64) -> C64 {
50        if w == C64::new(0.0, 0.0) {
51            return C64::new(f64::INFINITY, 0.0);
52        }
53        (C64::new(1.0, 0.0) + (w * w).inv()) * (0.5 * self.dw_h)
54    }
55}
56
57/// `z = 2 ω_min w / (1 - w²)`: the unit circle onto `i(-∞, -ω_min] ∪
58/// i[ω_min, ∞)` (`ConMapGapless`), for data with up-down symmetry.
59#[derive(Debug, Clone, Copy, PartialEq)]
60pub struct ConMapGapless {
61    /// `ω_min > 0`.
62    pub w_min: f64,
63}
64
65impl ConMap for ConMapGapless {
66    fn z(&self, w: C64) -> C64 {
67        if w == C64::new(1.0, 0.0) || w == C64::new(-1.0, 0.0) {
68            return C64::new(f64::INFINITY, 0.0);
69        }
70        w * (2.0 * self.w_min) / (C64::new(1.0, 0.0) - w * w)
71    }
72
73    fn w(&self, z: C64) -> C64 {
74        if z == C64::new(0.0, 0.0) {
75            return C64::new(0.0, 0.0);
76        }
77        let wm = self.w_min;
78        let mut w = ((z * z).inv() + 1.0 / (wm * wm)).sqrt() - z.inv();
79        w *= wm;
80        if w.norm() > 1.0 {
81            w = -w - z.inv() * (2.0 * wm);
82        }
83        w
84    }
85
86    fn dz(&self, w: C64) -> C64 {
87        if w == C64::new(1.0, 0.0) || w == C64::new(-1.0, 0.0) {
88            return C64::new(f64::INFINITY, 0.0);
89        }
90        let one = C64::new(1.0, 0.0);
91        let d = one - w * w;
92        (one + w * w) * (2.0 * self.w_min) / (d * d)
93    }
94}
sparse_ir_minipole/minipole/con_map.rs

sparse_ir_minipole/minipole/
con_map.rs
