treams_core/lattice/
sheets.rs

1//! Continuation of the reciprocal integrals to the sheet of the outgoing sum.
2//!
3//! treams-rs extension: treams takes the principal branches.
4
5use std::f64::consts::{FRAC_PI_2, PI, TAU};
6
7use super::accuracy::below_automatic;
8use crate::{Complex, special};
9
10/// `z` with a zero imaginary part made `-0.0`, which puts it on the lower side of a
11/// branch cut along the real axis.
12pub(super) fn lower_side(mut z: Complex) -> Complex {
13    if z.im == 0.0 {
14        z.im = -0.0;
15    }
16    z
17}
18
19/// Angle beyond `|arg((k eta)^2)|` by which a diffraction order's `k_q^2` must stay off
20/// the positive real axis for [`split_sheet`] to take the principal argument of `v`
21/// unevaluated: far above the rounding of that angle (1e-2 next to a threshold).
22const SHEET_MARGIN: f64 = 0.5;
23
24/// The split of an Ewald sum as the sheets of its diffraction orders need it
25/// ([`split_sheet`]).
26#[derive(Clone, Copy, Debug)]
27pub(super) struct Split {
28    /// `k^2`.
29    square: Complex,
30    /// `arg((k eta)^2)`, which turns every sheet.
31    turn: f64,
32    /// `tan^2(|arg((k eta)^2)| + SHEET_MARGIN)` below `pi / 2`, within which a `k_q^2`
33    /// needs the argument of its sheet.
34    near: Option<f64>,
35    /// Whether the split lies below every automatic one ([`below_automatic`]).
36    pub(super) small: bool,
37}
38
39impl Split {
40    pub(super) fn new(k: Complex, eta: Complex) -> Self {
41        let split = k * eta;
42        let turn = (split * split).arg();
43        let angle = turn.abs() + SHEET_MARGIN;
44        Self {
45            square: k * k,
46            turn,
47            near: (angle < FRAC_PI_2).then(|| angle.tan().powi(2)),
48            small: below_automatic(eta),
49        }
50    }
51}
52
53/// The sheet of `v` of one diffraction order ([`split_sheet`]).
54#[derive(Clone, Copy, Debug)]
55pub(super) enum Sheet {
56    /// The principal argument of `v`.
57    Principal,
58    /// The argument of `v` on the sheet, which differs from the principal one by whole
59    /// turns up to rounding.
60    Argument(f64),
61}
62
63impl Sheet {
64    /// The whole turns from the principal argument of `value` (`= v` up to rounding) to
65    /// the sheet.
66    pub(super) fn turns(self, value: Complex) -> f64 {
67        match self {
68            Self::Principal => 0.0,
69            Self::Argument(sheet) => sheet_turns(sheet, value),
70        }
71    }
72
73    /// Whether `root`, a root of `v` up to rounding, is the opposite of the root
74    /// `|v|^(1/2) e^(i arg / 2)` of the sheet's argument `arg`.
75    pub(super) fn opposes(self, root: Complex, value: Complex) -> bool {
76        let direction = match self {
77            Self::Principal => crate::numerics::complex_sqrt(lower_side(value)).conj(),
78            Self::Argument(sheet) => Complex::from_polar(1.0, -0.5 * sheet),
79        };
80        (root * direction).re < 0.0
81    }
82}
83
84/// The sheet of `v = (q^2 - k^2) / (2 k^2 eta^2)` of a diffraction order of squared
85/// transverse modulus `q2` that continues the Ewald sum from the automatic split to
86/// `split`.
87///
88/// # Formula
89///
90/// The argument `arg(k_q^2) - pi - arg((k eta)^2)` with `arg(k_q^2)` in `[0, 2 pi)`.
91///
92/// # Branch
93///
94/// It takes the root `k_q = sqrt(k^2 - q^2)` with `Im k_q >= 0` of the spectral series
95/// ([`spectral_sw1d`]). Orders farther off the positive real axis than
96/// `|arg((k eta)^2)| +` [`SHEET_MARGIN`] take the principal argument.
97///
98/// # Why
99///
100/// The principal `ln v` misses the sheet where `k_q^2` lies between the positive real
101/// axis and `(k eta)^2` (propagating orders of real `k` at splits rotated to
102/// `Im(k eta) > 0`). It leaves the side of the cut to rounding where `v` lies on it
103/// (`q = 0` for complex `k` and a real split).
104///
105/// [`spectral_sw1d`]: super::spectral::spectral_sw1d
106pub(super) fn split_sheet(split: &Split, q2: Complex) -> Sheet {
107    let square = split.square - q2;
108    let near = |tan2: f64| square.re > 0.0 && square.im * square.im <= tan2 * square.re * square.re;
109    if split.near.is_some_and(|tan2| !near(tan2)) {
110        return Sheet::Principal;
111    }
112    let turn = if square.im == 0.0 {
113        if square.re > 0.0 { 0.0 } else { PI }
114    } else {
115        square.arg().rem_euclid(TAU)
116    };
117    Sheet::Argument(turn - PI - split.turn)
118}
119
120/// The whole turns from the principal argument of `value` to `sheet`, which differ by
121/// a multiple of `2 pi` up to rounding.
122pub(super) fn sheet_turns(sheet: f64, value: Complex) -> f64 {
123    ((sheet - value.arg()) / TAU).round()
124}
125
126/// The sheet `arg v = 2 arg k - pi - arg((k eta)^2)` of `v = -1 / (2 eta^2)` for
127/// [`self_term`]: `-pi` continues the outgoing side `Im v = -0` of real `k` and splits,
128/// and `k` takes `arg k = pi` on the negative real axis for either sign of zero, as the
129/// AMOS outgoing waves do.
130///
131/// [`self_term`]: super::reciprocal::self_term
132pub(super) fn self_sheet(k: Complex, eta: Complex) -> f64 {
133    let k = if k.im == 0.0 {
134        Complex::new(k.re, 0.0)
135    } else {
136        k
137    };
138    2.0 * k.arg() - PI - special::kambe_argument(k, eta).arg()
139}
treams_core/lattice/sheets.rs

treams_core/lattice/
sheets.rs
