sva-samples 0.6.0

// Concern: the stiff-string FD grid both string models share, and its point coupling | Non-concern: what drives a string | IO: (Wire, f0, b, sr) -> StringGrid

#[derive(Clone)]
pub(crate) struct StringGrid {
    pub(crate) y_now: Vec<f64>,
    pub(crate) y_prev: Vec<f64>,
    pub(crate) y_next: Vec<f64>,
    pub(crate) n: usize,
    pub(crate) dx: f64,
    pub(crate) rho: f64,
    pub(crate) courant_sq: f64,
    pub(crate) stiff_sq: f64,
    pub(crate) damp_a: f64,
    pub(crate) damp_b: f64,
}

/// The ghost a biharmonic stencil reads past either end.
pub(crate) fn ghost_pinned(y: &[f64], n: usize, idx: isize, pin: f64) -> f64 {
    if idx < 0 {
        -y[(-idx) as usize]
    } else if idx as usize > n {
        2.0 * pin - y[(2 * n as isize - idx) as usize]
    } else {
        y[idx as usize]
    }
}

/// `pin` is `0.0` rigid, or a shared bridge's state.
pub(crate) fn stencil_update(grid: &StringGrid, j: usize, pin_now: f64, pin_prev: f64) -> f64 {
    let n = grid.n;
    let jj = j as isize;
    let lap_now = ghost_pinned(&grid.y_now, n, jj + 1, pin_now) - 2.0 * grid.y_now[j]
        + ghost_pinned(&grid.y_now, n, jj - 1, pin_now);
    let biharm = ghost_pinned(&grid.y_now, n, jj + 2, pin_now)
        - 4.0 * ghost_pinned(&grid.y_now, n, jj + 1, pin_now)
        + 6.0 * grid.y_now[j]
        - 4.0 * ghost_pinned(&grid.y_now, n, jj - 1, pin_now)
        + ghost_pinned(&grid.y_now, n, jj - 2, pin_now);
    let lap_prev = ghost_pinned(&grid.y_prev, n, jj + 1, pin_prev) - 2.0 * grid.y_prev[j]
        + ghost_pinned(&grid.y_prev, n, jj - 1, pin_prev);
    2.0 * grid.y_now[j] - grid.y_prev[j] + grid.courant_sq * lap_now
        - grid.stiff_sq * biharm
        - grid.damp_a * (grid.y_now[j] - grid.y_prev[j])
        + grid.damp_b * (lap_now - lap_prev)
}

/// The most intervals the lossless scheme keeps stable.
fn finest_stable_points(c: f64, length: f64, kappa: f64, dt: f64) -> usize {
    let (c2, dt2, k2) = (c * c, dt * dt, kappa * kappa);
    let dx_bound = ((c2 * dt2 + (c2 * c2 * dt2 * dt2 + 16.0 * k2 * dt2).sqrt()) / 2.0).sqrt();
    (length / dx_bound).floor() as usize
}

/// At rest, losses set.
fn lossy_grid(
    n: usize,
    rho: f64,
    length: f64,
    damp_dc: f64,
    damp_freq: f64,
    dt: f64,
) -> StringGrid {
    let dx = length / n as f64;
    StringGrid {
        y_now: vec![0.0; n + 1],
        y_prev: vec![0.0; n + 1],
        y_next: vec![0.0; n + 1],
        n,
        dx,
        rho,
        courant_sq: 0.0,
        stiff_sq: 0.0,
        damp_a: 2.0 * damp_dc * dt,
        damp_b: 2.0 * damp_freq * dt / (dx * dx),
    }
}

/// Pinned mode `m` is exactly `sin(m pi j/n)`.
pub(crate) fn mode_s(m: usize, n: usize) -> f64 {
    (m as f64 * std::f64::consts::PI / (2.0 * n as f64))
        .sin()
        .powi(2)
}

fn mode_sigma(grid: &StringGrid, s: f64) -> f64 {
    grid.damp_a + 4.0 * grid.damp_b * s
}

/// `2 - sigma - 2 sqrt(1 - sigma) cos(theta)`, without cancellation.
fn stiffness_for(theta: f64, sigma: f64) -> f64 {
    let r = (1.0 - sigma).sqrt();
    (sigma / (1.0 + r)).powi(2) + 4.0 * r * (theta / 2.0).sin().powi(2)
}

/// Jury's test on `z^2 + (D + sigma - 2) z + 1 - sigma`, every mode.
fn is_stable(grid: &StringGrid) -> bool {
    (1..grid.n).all(|m| {
        let s = mode_s(m, grid.n);
        let sigma = mode_sigma(grid, s);
        let d = 4.0 * grid.courant_sq * s + 16.0 * grid.stiff_sq * s * s;
        (0.0..2.0).contains(&sigma) && d > 0.0 && d < 4.0 - 2.0 * sigma
    })
}

pub(crate) struct Wire {
    pub(crate) rho: f64,
    pub(crate) c: f64,
    pub(crate) length: f64,
}

/// Partials 1 and 2 at `k f0 sqrt(1 + b k^2)`, on the finest grid that can.
pub(crate) fn dispersive_grid(
    wire: Wire,
    f0: f64,
    b: f64,
    damp_dc: f64,
    damp_freq: f64,
    sr: f64,
) -> Option<StringGrid> {
    let dt = 1.0 / sr;
    let Wire { rho, c, length } = wire;
    let kappa = c * length * b.sqrt() / std::f64::consts::PI;
    let theta = |k: f64| std::f64::consts::TAU * f0 * k * (1.0 + b * k * k).sqrt() * dt;
    let (theta_1, theta_2) = (theta(1.0), theta(2.0));
    if theta_2 >= std::f64::consts::PI {
        return None;
    }
    (3..=finest_stable_points(c, length, kappa, dt))
        .rev()
        .find_map(|n| {
            let mut grid = lossy_grid(n, rho, length, damp_dc, damp_freq, dt);
            let (s1, s2) = (mode_s(1, n), mode_s(2, n));
            let (sigma_1, sigma_2) = (mode_sigma(&grid, s1), mode_sigma(&grid, s2));
            if sigma_1 >= 1.0 || sigma_2 >= 1.0 {
                return None;
            }
            let (d1, d2) = (
                stiffness_for(theta_1, sigma_1),
                stiffness_for(theta_2, sigma_2),
            );
            // `D_m = 4 lambda^2 s_m + 16 mu^2 s_m^2`, m = 1, 2.
            let det = s1 * s2 * (s2 - s1);
            grid.courant_sq = (d1 * s2 * s2 - d2 * s1 * s1) / (4.0 * det);
            grid.stiff_sq = (s1 * d2 - s2 * d1) / (16.0 * det);
            (grid.courant_sq > 0.0 && is_stable(&grid)).then_some(grid)
        })
}

/// `rho (lambda dx/dt)^2`.
pub(crate) fn grid_tension(grid: &StringGrid, dt: f64) -> f64 {
    grid.rho * grid.courant_sq * grid.dx * grid.dx / (dt * dt)
}

/// Modal content `sin(m pi x)` of every grid mode; read and spread alike.
pub(crate) fn point_weights(n: usize, x: f64) -> Vec<f64> {
    let nf = n as f64;
    // sum_{m=1}^{n-1} cos(m theta)
    let cos_sum = |theta: f64| {
        let half = theta / 2.0;
        let sh = half.sin();
        if sh == 0.0 {
            nf - 1.0
        } else {
            (nf * half).sin() * ((nf - 1.0) * half).cos() / sh - 1.0
        }
    };
    (0..=n)
        .map(|j| match j {
            0 => 0.0,
            j if j == n => 0.0,
            j => {
                let xj = j as f64 / nf;
                let pi = std::f64::consts::PI;
                (cos_sum(pi * (x - xj)) - cos_sum(pi * (x + xj))) / nf
            }
        })
        .collect()
}

/// `spread`'s adjoint.
pub(crate) fn read_at(weights: &[f64], y: &[f64]) -> f64 {
    weights.iter().zip(y).map(|(w, y)| w * y).sum()
}

/// The readout's adjoint.
pub(crate) fn spread(grid: &mut StringGrid, weights: &[f64], force: f64, dt: f64) {
    if force == 0.0 {
        return;
    }
    let scale = (dt * dt / (grid.rho * grid.dx)) * force;
    for (y, w) in grid.y_next.iter_mut().zip(weights) {
        *y += scale * w;
    }
}