photonoxide 0.3.2

//! The modes of a bent slab waveguide, exactly: a radial shooting method matched to the
//! outgoing cylinder wave.
//!
//! A slab bent around an axis y, its layers on cylinders of radius ρ around it, carries modes
//! that travel around the axis as e^(i(νφ − ωt)) and radiate outwards: ν is complex, and the
//! mode loses power along the bend. Outside the last layer the field is exactly the outgoing
//! Hankel function H⁽¹⁾_ν(k n ρ) (D. Marcuse, Bell Syst. Tech. J. 50, 2551 (1971),
//! [doi:10.1002/j.1538-7305.1971.tb02620.x](https://doi.org/10.1002/j.1538-7305.1971.tb02620.x),
//! Eq. 10, there H⁽²⁾ for e^(+iωt)). Marcuse approximates the field inside by the straight
//! guide's; here it isn't approximated:
//!
//! - in each layer of index n the field ψ (E_y for TE, E along the axis; H_y for TM) satisfies
//!   ψ'' + ψ'/ρ + (k²n² − ν²/ρ²) ψ = 0, integrated by fourth-order Runge–Kutta in ρ;
//! - at an interface ψ and ψ' (TE) or ψ'/n² (TM) are continuous;
//! - inside, the integration starts in the inner cladding, below the inner turning point, with
//!   the solution growing outwards (Bessel's J_ν); the other solution dies out on the way;
//! - outside, beyond the outer turning point (k n ρ = 2.5 Re ν), ψ'/ψ must be that of H⁽¹⁾_ν,
//!   from Debye's expansion (DLMF 10.19.6), its polynomials generated by DLMF 10.41.9 and
//!   10.41.11 and checked against the large-argument expansion (DLMF 10.17.5);
//! - ν is a root of the mismatch, found by secant iteration.
//!
//! The effective index is ν/(kR) at a reference radius R: the propagation constant along the
//! arc there, so it can be compared with a straight guide's, and with photonoxide's bent
//! cross-sections ([`crate::mode::vector::CrossSection::bent`]).

use num_complex::Complex64 as c64;

use super::Polarization;
use crate::units::{Length, Wavelength};
use crate::{Error, Result};

/// A bent slab: layers between an inner and an outer cladding, at distances x from a reference
/// radius R (ρ = R + x, positive x towards the outside of the bend).
#[derive(Clone, Debug, PartialEq)]
pub struct SlabBend {
    radius: f64,
    inner: f64,
    outer: f64,
    /// The interfaces' x, increasing, and the index outside each: layer j lies between
    /// interface j and j + 1, with index `layers[j]`.
    interfaces: Vec<f64>,
    layers: Vec<f64>,
}

impl SlabBend {
    /// Films (index and thickness) starting at x = `start` from the reference radius `radius`,
    /// between the `inner` cladding (smaller radii) and the `outer` one. Lossless indices.
    ///
    /// # Errors
    ///
    /// [`Error::InvalidValue`] for a radius that isn't positive and finite, no films, a
    /// thickness that isn't positive, an index below 1, or a stack reaching the centre
    /// (R + start ≤ 0).
    pub fn new(
        radius: Length,
        inner: f64,
        films: &[(f64, Length)],
        outer: f64,
        start: Length,
    ) -> Result<SlabBend> {
        let r = radius.to_um();
        if !(r.is_finite() && r > 0.0) {
            return Err(Error::invalid(
                "bent slab",
                format!("the radius must be positive, got {radius}"),
            ));
        }
        if films.is_empty() {
            return Err(Error::invalid("bent slab", "needs at least one film"));
        }
        let ok = |n: f64| n.is_finite() && n >= 1.0;
        if !ok(inner) || !ok(outer) || films.iter().any(|&(n, _)| !ok(n)) {
            return Err(Error::invalid(
                "bent slab",
                "indices must be finite and at least 1",
            ));
        }
        let mut x = start.to_um();
        if !(x.is_finite() && r + x > 0.0) {
            return Err(Error::invalid(
                "bent slab",
                "the stack must lie at a positive radius",
            ));
        }
        let mut interfaces = vec![x];
        for &(_, d) in films {
            let t = d.to_um();
            if !(t.is_finite() && t > 0.0) {
                return Err(Error::invalid(
                    "bent slab",
                    format!("film thicknesses must be positive, got {d}"),
                ));
            }
            x += t;
            interfaces.push(x);
        }
        Ok(SlabBend {
            radius: r,
            inner,
            outer,
            interfaces,
            layers: films.iter().map(|&(n, _)| n).collect(),
        })
    }

    /// The index at radius ρ (an interface belongs to the layer outside it).
    fn index_at(&self, rho: f64) -> f64 {
        let x = rho - self.radius;
        if x < self.interfaces[0] {
            return self.inner;
        }
        for (j, &n) in self.layers.iter().enumerate() {
            if x < self.interfaces[j + 1] {
                return n;
            }
        }
        self.outer
    }

    /// Integrates ψ'' + ψ'/ρ + (k²n² − ν²/ρ²)ψ = 0 from `from` to `to` (either way) through the
    /// layers, with ψ' scaled by the ratio of n² at each interface for TM; returns ψ'/ψ at
    /// `to`, starting from ψ'/ψ = `ratio` at `from`.
    fn integrate(
        &self,
        polarization: Polarization,
        k: f64,
        nu: c64,
        from: f64,
        to: f64,
        ratio: c64,
    ) -> c64 {
        let r = self.radius;
        let mut stops: Vec<f64> = self
            .interfaces
            .iter()
            .map(|&x| r + x)
            .filter(|&p| (p - from) * (p - to) < 0.0)
            .collect();
        if to < from {
            stops.reverse();
        }
        stops.push(to);
        let (mut psi, mut dpsi) = (c64::new(1.0, 0.0), ratio);
        let mut rho = from;
        let mut n = self.index_at(if to > from { from } else { from - 1e-12 });
        for stop in stops {
            let k2n2 = c64::new(k * k * n * n, 0.0);
            let span = stop - rho;
            let steps = (span.abs() / (0.03 / (k * n)).min(0.01)).ceil().max(1.0) as usize;
            let h = span / steps as f64;
            let f = |p: f64, y: c64, dy: c64| -> (c64, c64) {
                (dy, -dy / p - (k2n2 - nu * nu / (p * p)) * y)
            };
            for _ in 0..steps {
                let (k1y, k1d) = f(rho, psi, dpsi);
                let (k2y, k2d) = f(rho + 0.5 * h, psi + 0.5 * h * k1y, dpsi + 0.5 * h * k1d);
                let (k3y, k3d) = f(rho + 0.5 * h, psi + 0.5 * h * k2y, dpsi + 0.5 * h * k2d);
                let (k4y, k4d) = f(rho + h, psi + h * k3y, dpsi + h * k3d);
                psi += h / 6.0 * (k1y + 2.0 * k2y + 2.0 * k3y + k4y);
                dpsi += h / 6.0 * (k1d + 2.0 * k2d + 2.0 * k3d + k4d);
                rho += h;
                // only ψ'/ψ matters: keep the numbers in range
                let scale = psi.norm();
                if scale > 1e100 || scale < 1e-100 {
                    psi /= scale;
                    dpsi /= scale;
                }
            }
            rho = stop;
            if stop != to {
                // across an interface: the layer on the far side
                let next = self.index_at(if to > from {
                    stop + 1e-12
                } else {
                    stop - 1e-12
                });
                if polarization == Polarization::Tm {
                    dpsi *= (next * next) / (n * n);
                }
                n = next;
            }
        }
        dpsi / psi
    }

    /// ψ'/ψ of the inner solution minus that of the outgoing one, at the middle of the stack:
    /// zero at a mode.
    fn mismatch(&self, polarization: Polarization, k: f64, nu: c64) -> c64 {
        let r = self.radius;
        let first = r + self.interfaces[0];
        let last = r + self.interfaces[self.interfaces.len() - 1];
        let middle = 0.5 * (first + last);
        // inside: a few decay lengths below the stack, the solution growing outwards
        let kappa = |rho: f64| {
            (nu * nu / (rho * rho) - c64::new(k * k * self.inner * self.inner, 0.0)).sqrt()
        };
        let start = (first - 6.0 / kappa(first).re.max(1e-3)).max(0.5 * first);
        let inside = self.integrate(
            polarization,
            k,
            nu,
            start,
            middle,
            kappa(start) - 1.0 / (2.0 * start),
        );
        // outside: beyond the outer turning point, the outgoing cylinder wave
        let rho_out = (2.5 * nu.re / (k * self.outer)).max(last + 1.0);
        let wave = k * self.outer * debye_log_derivative(nu, k * self.outer * rho_out);
        let outside = self.integrate(polarization, k, nu, rho_out, middle, wave);
        inside - outside
    }

    /// The fundamental mode of one polarization: from the straight guide's effective index at
    /// eight times the radius, where the bend barely changes it, followed in steps of √2 down
    /// to the radius, each step's mode the next one's guess.
    ///
    /// # Errors
    ///
    /// [`Error::InvalidValue`] if the straight guide (the bend's layers between its claddings)
    /// guides no mode, or the iteration fails on the way.
    pub fn fundamental(&self, polarization: Polarization, wavelength: Wavelength) -> Result<c64> {
        let films: Vec<(c64, Length)> = self
            .layers
            .iter()
            .zip(self.interfaces.windows(2))
            .map(|(&n, x)| (c64::new(n, 0.0), Length::um(x[1] - x[0])))
            .collect();
        let straight = crate::mode::multilayer::Multilayer::new(
            c64::new(self.inner, 0.0),
            &films,
            c64::new(self.outer, 0.0),
        )?
        .bound_modes(polarization, wavelength)?
        .first()
        .map(|m| m.effective_index())
        .ok_or_else(|| Error::invalid("bent slab", "the straight guide guides no mode"))?;
        let mut guess = straight;
        for step in (0..=6).rev() {
            let radius = self.radius * 2f64.powf(step as f64 / 2.0);
            let bend = SlabBend {
                radius,
                ..self.clone()
            };
            guess = bend.mode_near(polarization, wavelength, guess)?;
        }
        Ok(guess)
    }

    /// The mode whose angular order the secant iteration reaches from the effective index
    /// `guess` (at the reference radius): Re n_eff, and Im n_eff > 0, the radiation loss.
    ///
    /// # Errors
    ///
    /// [`Error::InvalidValue`] if the iteration doesn't converge in 60 steps.
    pub fn mode_near(
        &self,
        polarization: Polarization,
        wavelength: Wavelength,
        guess: c64,
    ) -> Result<c64> {
        let k = wavelength.wavenumber();
        let kr = k * self.radius;
        let f = |n: c64| self.mismatch(polarization, k, n * kr);
        let (mut x0, mut x1) = (guess, guess * (1.0 + 1e-7) + c64::new(0.0, 1e-9));
        let (mut f0, mut f1) = (f(x0), f(x1));
        // the integration's round-off makes f noisy near 1e-13: keep the best iterate, and
        // accept it once the steps have come down to 1e-9
        let mut best = (f1.norm(), x1, f64::INFINITY);
        for _ in 0..60 {
            if f1 == f0 {
                break;
            }
            let x2 = x1 - f1 * (x1 - x0) / (f1 - f0);
            (x0, f0) = (x1, f1);
            x1 = x2;
            f1 = f(x1);
            let step = (x1 - x0).norm() / x1.norm();
            if f1.norm() < best.0 {
                best = (f1.norm(), x1, step);
            }
            if step < 1e-13 {
                return Ok(x1);
            }
        }
        if best.2 < 1e-9 {
            return Ok(best.1);
        }
        Err(Error::invalid(
            "bent slab",
            format!("no mode found near {guess}: the iteration didn't converge"),
        ))
    }
}

/// The polynomials u_k and v_k of Debye's expansions, k = 0 … `count` − 1, as coefficients of
/// powers of p: u_{k+1}(p) = ½p²(1 − p²)u_k'(p) + ⅛∫₀^p (1 − 5t²)u_k(t) dt (DLMF 10.41.9) and
/// v_k(p) = u_k(p) + p(p² − 1)(½u_{k−1}(p) + p u_{k−1}'(p)) (DLMF 10.41.11).
fn debye_polynomials(count: usize) -> (Vec<Vec<f64>>, Vec<Vec<f64>>) {
    let deriv = |c: &[f64]| -> Vec<f64> {
        c.iter()
            .enumerate()
            .skip(1)
            .map(|(i, &a)| i as f64 * a)
            .collect()
    };
    let mul = |a: &[f64], b: &[f64]| -> Vec<f64> {
        let mut out = vec![0.0; a.len() + b.len() - 1];
        for (i, &x) in a.iter().enumerate() {
            for (j, &y) in b.iter().enumerate() {
                out[i + j] += x * y;
            }
        }
        out
    };
    let add = |a: &[f64], b: &[f64]| -> Vec<f64> {
        (0..a.len().max(b.len()))
            .map(|i| a.get(i).unwrap_or(&0.0) + b.get(i).unwrap_or(&0.0))
            .collect()
    };
    let mut u = vec![vec![1.0]];
    for k in 0..count.saturating_sub(1) {
        let a = mul(&[0.0, 0.0, 0.5, 0.0, -0.5], &deriv(&u[k]));
        let integrand = mul(&[1.0, 0.0, -5.0], &u[k]);
        let mut integral = vec![0.0];
        integral.extend(
            integrand
                .iter()
                .enumerate()
                .map(|(i, &c)| c / (i + 1) as f64 / 8.0),
        );
        u.push(add(&a, &integral));
    }
    let mut v = vec![vec![1.0]];
    for k in 1..count {
        let inner = add(
            &mul(&[0.5], &u[k - 1]),
            &mul(&[0.0, 1.0], &deriv(&u[k - 1])),
        );
        v.push(add(&u[k], &mul(&[0.0, -1.0, 0.0, 1.0], &inner)));
    }
    (u, v)
}

/// H⁽¹⁾_ν'(z) / H⁽¹⁾_ν(z) beyond the turning point (z > Re ν, ν large), from Debye's
/// expansions with z = ν sec β: H⁽¹⁾_ν ∝ Σ u_k(−i cot β)/ν^k and H⁽¹⁾_ν' ∝ i sin β Σ
/// v_k(−i cot β)/ν^k, with the same prefactor once (sin 2β tan β / 2)^½ = sin β is taken out
/// (DLMF 10.19.6 writes H⁽²⁾'s, whose conjugate under real ν this is). The sign of the
/// polynomials' argument is fixed by the large-argument limit, H'/H → i − 1/(2z), and every
/// order is checked against that expansion in the tests.
pub(crate) fn debye_log_derivative(nu: c64, z: f64) -> c64 {
    let (u, v) = debye_polynomials(10);
    let cos = nu / z;
    let sin = (1.0 - cos * cos).sqrt();
    let p = c64::new(0.0, -1.0) * cos / sin;
    let eval = |c: &[f64]| {
        c.iter()
            .rev()
            .fold(c64::new(0.0, 0.0), |acc, &a| acc * p + a)
    };
    let (mut su, mut sv) = (c64::new(0.0, 0.0), c64::new(0.0, 0.0));
    let mut factor = c64::new(1.0, 0.0);
    for k in 0..u.len() {
        su += factor * eval(&u[k]);
        sv += factor * eval(&v[k]);
        factor /= nu;
    }
    c64::new(0.0, 1.0) * sin * sv / su
}

#[cfg(test)]
/// H⁽¹⁾_ν'(z) / H⁽¹⁾_ν(z) for complex ν and large real z (z ≥ 2|ν|²), from the large-argument
/// expansion H⁽¹⁾_ν(z) ~ (2/πz)^½ e^(i(z − νπ/2 − π/4)) Σ_k i^k a_k(ν)/z^k, a_k(ν) =
/// Π_{j ≤ k} (4ν² − (2j − 1)²) / (k! 8^k) (DLMF 10.17.5), and H_ν' = H_{ν−1} − (ν/z) H_ν
/// (DLMF 10.6.2).
pub(crate) fn hankel1_log_derivative(nu: c64, z: f64) -> c64 {
    let series = |mu: c64| -> c64 {
        let mut sum = c64::new(1.0, 0.0);
        let mut term = c64::new(1.0, 0.0);
        let four_mu2 = 4.0 * mu * mu;
        for k in 1..400 {
            let odd = (2 * k - 1) as f64;
            term *= c64::new(0.0, 1.0) * (four_mu2 - odd * odd) / (8.0 * k as f64 * z);
            sum += term;
            if term.norm() < 1e-18 * sum.norm() {
                break;
            }
        }
        sum
    };
    // H_{ν−1}/H_ν = e^(iπ/2) S(ν−1)/S(ν)
    c64::new(0.0, 1.0) * series(nu - 1.0) / series(nu) - nu / z
}

/// Marcuse's bending loss of a symmetric slab's fundamental TE mode, as an imaginary effective
/// index: his Eqs. (32)–(33) with the exact U, for a core of index `core` and thickness 2d in
/// `cladding`, bent at radius R (to the slab's centre), from the straight slab's effective
/// index `n_eff`. An approximation for large R (his inequalities 35–36).
pub fn marcuse_loss(
    core: f64,
    cladding: f64,
    thickness: Length,
    radius: Length,
    wavelength: Wavelength,
    n_eff: f64,
) -> f64 {
    let k = wavelength.wavenumber();
    let d = thickness.to_um() / 2.0;
    let r = radius.to_um();
    let beta = k * n_eff;
    let gamma = (beta * beta - k * k * cladding * cladding).sqrt();
    let kappa = (k * k * core * core - beta * beta).sqrt();
    let u = (beta / gamma * ((1.0 + gamma / beta) / (1.0 - gamma / beta)).ln() - 2.0) * gamma * r;
    // 2α, power loss per unit length; symmetric: θ = γ
    let two_alpha = 2.0 * gamma * kappa * kappa * (2.0 * gamma * d).exp() * (-u).exp()
        / ((core * core - cladding * cladding) * k * k * beta * (2.0 * d + 2.0 / gamma));
    // the amplitude decays as e^(−αz): Im β = α
    two_alpha / 2.0 / k
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::mode::slab::Slab;

    fn lam(um: f64) -> Wavelength {
        Wavelength::um(um).unwrap()
    }

    #[test]
    fn the_hankel_ratio_matches_its_closed_form_at_half_order() {
        // H⁽¹⁾_{1/2}(z) = −i (2/πz)^½ e^(iz): H'/H = i − 1/(2z), exactly
        for z in [30.0, 200.0, 1000.0] {
            let got = hankel1_log_derivative(c64::new(0.5, 0.0), z);
            let want = c64::new(-0.5 / z, 1.0);
            assert!((got - want).norm() < 1e-14, "{z}: {got} vs {want}");
        }
    }

    #[test]
    fn debyes_expansion_agrees_with_the_large_argument_one() {
        // u_1, v_1 as DLMF prints them, and the two expansions where both hold (z ≥ 2|ν|²)
        let (u, v) = debye_polynomials(3);
        let close = |a: &[f64], b: &[f64]| a.iter().zip(b).all(|(x, y)| (x - y).abs() < 1e-15);
        assert!(
            close(&u[1], &[0.0, 3.0 / 24.0, 0.0, -5.0 / 24.0]),
            "{:?}",
            u[1]
        );
        assert!(
            close(&v[1], &[0.0, -9.0 / 24.0, 0.0, 7.0 / 24.0]),
            "{:?}",
            v[1]
        );
        for nu in [
            c64::new(8.0, 0.0),
            c64::new(12.5, 0.01),
            c64::new(20.0, 0.3),
        ] {
            let z = 2.0 * nu.norm_sqr() + 60.0;
            let (a, b) = (debye_log_derivative(nu, z), hankel1_log_derivative(nu, z));
            assert!((a - b).norm() < 1e-12, "{nu}: {a} vs {b}");
        }
    }

    #[test]
    fn debyes_ratio_is_right_near_the_turning_point() {
        // integrate Bessel's equation from z = 2.5ν, started with Debye's ratio, out to where
        // the large-argument expansion holds: the ratios must agree there
        for nu in [c64::new(12.0, 0.0), c64::new(30.0, 0.05)] {
            let z0 = 2.5 * nu.re;
            let z1 = 2.0 * nu.norm_sqr() + 60.0;
            let (mut y, mut dy) = (c64::new(1.0, 0.0), debye_log_derivative(nu, z0));
            let steps = ((z1 - z0) / 0.005).ceil() as usize;
            let h = (z1 - z0) / steps as f64;
            let mut z = z0;
            let f = |z: f64, y: c64, dy: c64| (dy, -dy / z - (1.0 - nu * nu / (z * z)) * y);
            for _ in 0..steps {
                let (a1, b1) = f(z, y, dy);
                let (a2, b2) = f(z + h / 2.0, y + h / 2.0 * a1, dy + h / 2.0 * b1);
                let (a3, b3) = f(z + h / 2.0, y + h / 2.0 * a2, dy + h / 2.0 * b2);
                let (a4, b4) = f(z + h, y + h * a3, dy + h * b3);
                y += h / 6.0 * (a1 + 2.0 * a2 + 2.0 * a3 + a4);
                dy += h / 6.0 * (b1 + 2.0 * b2 + 2.0 * b3 + b4);
                z += h;
            }
            let (got, want) = (dy / y, hankel1_log_derivative(nu, z1));
            assert!((got - want).norm() < 1e-7, "{nu}: {got} vs {want}");
        }
    }

    #[test]
    fn a_large_radius_gives_the_straight_slab() {
        // the bend's index at the slab's centre tends to the straight one's as 1/R²
        let (core, clad, t) = (2.0, 1.5, 0.6);
        let straight = Slab::new(clad, core, clad, Length::um(t))
            .unwrap()
            .modes(Polarization::Te, lam(1.0))[0]
            .effective_index();
        let mut last = f64::NAN;
        for r in [50.0, 100.0, 200.0] {
            let bend = SlabBend::new(
                Length::um(r),
                clad,
                &[(core, Length::um(t))],
                clad,
                Length::um(-t / 2.0),
            )
            .unwrap();
            let n = bend
                .mode_near(Polarization::Te, lam(1.0), c64::new(straight, 0.0))
                .unwrap();
            let err = (n.re - straight).abs();
            if last.is_finite() {
                assert!((last / err - 4.0).abs() < 0.5, "R {r}: {last} then {err}");
            }
            last = err;
            // the loss is below double precision here: zero to round-off
            assert!(n.im > -1e-15 && n.im < 1e-6, "{n}");
        }
    }
}