kshana 0.37.0

// SPDX-License-Identifier: AGPL-3.0-only
//! **FFT parallel code-phase acquisition for any [`SpreadingCode`].**
//!
//! For every Doppler bin the carrier (intermediate frequency plus the bin's Doppler) is
//! wiped off the samples, each coherent block of `N` code periods is folded onto one
//! period (the sum of its periods, which is the coherent integration because the replica
//! repeats every period), and the circular correlation against one sampled period of the
//! code is formed for every code phase at once with a fast Fourier transform:
//! `IFFT(FFT(fold) · conj(FFT(replica)))`. The squared magnitudes of `M` consecutive
//! blocks are summed (non-coherent integration).
//!
//! The FFT is the crate's existing mixed-radix transform ([`crate::acquisition::fft_forward`]
//! and its plan in `crate::portable_math`), so a search returns the same bits on every
//! platform. The detector statistics are the existing square-law / Marcum-Q ones in
//! [`crate::acquisition`]: each cell is normalised by the measured sample power `σ²` to
//! `2·G/(L·σ²)`, which under noise alone is chi-square with `2M` degrees of freedom and with
//! a signal is non-central chi-square with non-centrality `2M·ρ`, `ρ = (C/N0)·N·T_code`.
//! The threshold is [`crate::acquisition::threshold_for_pfa`] at the per-cell false-alarm
//! probability that gives the requested probability over the whole grid,
//! `1 − (1 − P_fa)^(1/N_cells)`, and [`predicted_pd`] is
//! [`crate::acquisition::pd_square_law`] at that threshold.
//!
//! The code is evaluated through [`SpreadingCode::value_at`] at the sample instants
//! (point sampling, no front-end filter), so BPSK and BOC codes both work. A sample rate
//! must hold a whole number of samples per code period.
//!
//! Scope, stated: no code-Doppler compensation across non-coherent blocks, no data-bit
//! wipe-off (a bit edge inside a coherent block loses correlation), no straddle-loss
//! recovery beyond choosing a fine Doppler step. MODELLED: the
//! detector identities are checked against the analytic Marcum-Q values by seeded Monte
//! Carlo in `tests/iq_receiver.rs`, not against an external dataset.
//!
//! An opt-in cell-averaging CFAR detector ([`ca_cfar`], [`CfarConfig`]) decides on the same
//! grid by comparing a cell with the mean of reference cells in its row, so its false-alarm
//! probability does not depend on the sample-power normalisation (quantised or band-limited
//! input). Its closed form is [`ca_cfar_pfa`]; it is checked in `tests/iq_acq_cfar.rs`.

use super::{Cf64, IqError, IqSource, SampleSpec, SpreadingCode};
use crate::acquisition::{pd_square_law, threshold_for_pfa};
use crate::portable_math::{FftPlan, PortableFloat};

/// The parameters of one acquisition search.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct AcqConfig {
    /// Coherent integration in whole code periods (`N`).
    pub coherent_periods: usize,
    /// Number of consecutive coherent blocks whose powers are summed (`M`).
    pub noncoherent: usize,
    /// The Doppler grid is `j · doppler_step_hz` for every integer `j` with
    /// `|j · doppler_step_hz| ≤ doppler_max_hz` (symmetric, always containing 0).
    pub doppler_max_hz: f64,
    /// Doppler bin spacing (Hz). A common choice is `2 / (3 · N · T_code)`.
    pub doppler_step_hz: f64,
    /// False-alarm probability for the whole search grid.
    pub pfa: f64,
}

/// Coherent integration time (s) the tracking hand-off searches with by default: about
/// 4 ms. See [`auto_coherent_periods`].
pub const AUTO_COHERENT_S: f64 = 4.0e-3;

/// The default coherent length, in whole code periods, of the acquisition that initialises
/// tracking: `ceil(4 ms / T_code)`, at least one period.
///
/// `T_code` is the code's FULL period, [`SpreadingCode::period_s`] (primary length times
/// secondary length for a tiered code), because that is the unit [`acquire`] integrates
/// over: its replica is the whole tiered code. So every signal gets at least ~4 ms of
/// coherent integration. The untiered 1 ms codes (GPS L1 C/A, BeiDou B1I, GLONASS L1OF) get
/// 4 periods (4 ms, ~167 Hz bins with the default Doppler step `2 / (3 · N · T_code)`).
/// Every code whose full period is already 4 ms or longer gets 1: Galileo E1-B (4 ms),
/// BeiDou B1C (10 ms), GPS L5-I/L5-Q (10/20 ms with their Neuman-Hofman overlay), Galileo
/// E5a-I/E5a-Q (20/100 ms), Galileo E1-C (100 ms with CS25) and GPS L2C. That is the same
/// one tiered period those codes searched before, and finer bins than 167 Hz.
///
/// Why not one period: a 1 ms search has ~667 Hz Doppler bins, and a hand-off up to
/// ~333 Hz off leaves the FLL outside its pull-in. On a seeded sweep of 180 GPS L1 C/A
/// channels (±5 kHz, 38 to 47 dB-Hz; `docs/design/evidence/iq-track-acq-default/`) one
/// period false-locked 27 channels ~500 Hz off and missed 102; this default false-locked 1
/// (36 Hz off at 38 dB-Hz) and missed 11, for ~40 ms more search per PRN.
pub fn auto_coherent_periods(code_period_s: f64) -> usize {
    if !(code_period_s.is_finite() && code_period_s > 0.0) {
        return 1;
    }
    // The small tolerance keeps 4 ms / 1 ms (4.000…01 in binary) at 4, not 5.
    ((AUTO_COHERENT_S / code_period_s - 1e-9).ceil() as usize).max(1)
}

/// The default Doppler step is capped at `PULL_IN_STEP_MARGIN / T_track`: `0.2 / T_track`,
/// which is 0.8 of the default Atan2 FLL's pull-in `1 / (4 T_track)`. The cap is a step, not a
/// half step, because the acquisition's winning bin is not always the nearest one: with a
/// short coherent time the main lobe (`1 / (N T_code)` wide) is much wider than a bin, and
/// noise can make the neighbouring bin win, leaving a residual of up to a whole step.
/// 0.4 (half-step residual at 80 % of pull-in) was tried first; on Galileo E1-B it false-locked
/// in 1 of 42 noisy runs, and with 1 ms coherent GPS L1 C/A a neighbouring bin 350 Hz off won
/// and tracking locked at +500 Hz. 0.25 and 0.2 gave no false lock in the same 42 runs.
pub const PULL_IN_STEP_MARGIN: f64 = 0.2;

/// The default Doppler step (Hz) for the acquisition that hands a channel to tracking.
///
/// The textbook step `2 / (3 · N · T_code)` leaves a worst-case residual of half a step,
/// `1 / (3 · N · T_code)`, and up to a whole step when the neighbouring bin wins. The default
/// carrier FLL (two-quadrant `atan2`) pulls in only `±1 / (4 T_track)`, so with `N = 1`
/// (every code of 4 ms or longer) the residual is outside the pull-in and the loop can lock
/// `1 / T_track` away. The step is therefore also capped at [`PULL_IN_STEP_MARGIN`]
/// `/ T_track`. For the 1 ms codes (`N = 4`, 166.7 Hz against a 200 Hz cap) the textbook step
/// is already the smaller and is returned unchanged.
///
/// `t_track_s` is the tracking loop's integration time (code periods per loop update times
/// the full code period); a non-positive or non-finite value gives the textbook step.
pub fn default_step_hz(code_period_s: f64, n_coh: usize, t_track_s: f64) -> f64 {
    let textbook = 2.0 / (3.0 * n_coh.max(1) as f64 * code_period_s);
    if t_track_s.is_finite() && t_track_s > 0.0 {
        textbook.min(PULL_IN_STEP_MARGIN / t_track_s)
    } else {
        textbook
    }
}

impl AcqConfig {
    /// The Doppler grid (Hz), ascending.
    pub fn doppler_bins(&self) -> Vec<f64> {
        let n = (self.doppler_max_hz / self.doppler_step_hz + 1e-9).floor() as i64;
        (-n..=n).map(|j| j as f64 * self.doppler_step_hz).collect()
    }
}

/// Samples in one period of `code` at the sample rate of `spec`, or an error when that is
/// not a whole number (within 1e-6 of a sample).
pub fn samples_per_period(spec: &SampleSpec, code: &dyn SpreadingCode) -> Result<usize, String> {
    let spc = spec.fs_hz * code.period_s();
    if !(spc.is_finite() && spc >= 2.0) || (spc - spc.round()).abs() > 1e-6 {
        return Err(format!(
            "{} Hz does not hold a whole number of samples per {} period ({spc} samples)",
            spec.fs_hz,
            code.name()
        ));
    }
    Ok(spc.round() as usize)
}

/// Number of samples one search of `cfg` consumes: `N · M` code periods.
pub fn samples_needed(
    spec: &SampleSpec,
    code: &dyn SpreadingCode,
    cfg: &AcqConfig,
) -> Result<usize, String> {
    Ok(samples_per_period(spec, code)? * cfg.coherent_periods * cfg.noncoherent)
}

/// The outcome of one search.
#[derive(Clone, Debug, PartialEq)]
pub struct AcqResult {
    /// Name of the code searched.
    pub code_name: String,
    /// Doppler of the peak bin (Hz), relative to the intermediate frequency.
    pub doppler_hz: f64,
    /// Index of the peak Doppler bin in [`AcqConfig::doppler_bins`].
    pub doppler_index: usize,
    /// Lag of the peak (samples, `0 .. samples_per_period`): the sample, counted from the
    /// first sample searched, at which a code period starts.
    pub delay_samples: usize,
    /// Code phase of the signal at the first sample searched (chips, `0 .. len`).
    pub code_phase_chips: f64,
    /// Normalised peak `2·G_max/(L·σ²)` (chi-square with `2M` degrees of freedom under noise).
    pub statistic: f64,
    /// Largest normalised cell in the peak's Doppler bin more than one chip (plus one
    /// sample) from the peak.
    pub second_peak: f64,
    /// `statistic / second_peak`: the peak-to-second-peak ratio.
    pub peak_ratio: f64,
    /// Per-cell false-alarm probability the threshold was set for.
    pub pfa_cell: f64,
    /// The normalised threshold `statistic` is compared with.
    pub threshold: f64,
    /// Whether `statistic` exceeded `threshold`.
    pub acquired: bool,
    /// Samples per code period (`L / N`).
    pub samples_per_period: usize,
    /// Number of Doppler bins searched.
    pub n_doppler_bins: usize,
    /// Mean sample power `σ²` used to normalise.
    pub sample_power: f64,
}

/// A search together with its full normalised grid `grid[doppler_index][lag]`.
#[derive(Clone, Debug, PartialEq)]
pub struct AcqGrid {
    /// The search outcome.
    pub result: AcqResult,
    /// Normalised cell values `2·G/(L·σ²)`, one row per Doppler bin.
    pub grid: Vec<Vec<f64>>,
}

/// The per-search state the correlation rows share: the FFT plan, the conjugated code
/// spectrum and the normalisation. One [`RowEngine::row`] is one Doppler bin's lag
/// profile, so a caller can keep every row (the grid) or only a running best.
struct RowEngine<'a> {
    samples: &'a [Cf64],
    fs: f64,
    base_hz: f64,
    spc: usize,
    block: usize,
    noncoherent: usize,
    plan: FftPlan,
    code_fft: Vec<(f64, f64)>,
    norm: f64,
    sigma2: f64,
}

impl<'a> RowEngine<'a> {
    /// `samples` must hold [`samples_needed`] samples.
    fn new(
        samples: &'a [Cf64],
        spec: &SampleSpec,
        code: &dyn SpreadingCode,
        cfg: &AcqConfig,
    ) -> Result<Self, String> {
        let spc = samples_per_period(spec, code)?;
        let block = spc * cfg.coherent_periods;
        let needed = block * cfg.noncoherent;
        if samples.len() < needed {
            return Err(format!(
                "{} samples supplied, {needed} needed",
                samples.len()
            ));
        }
        let fs = spec.fs_hz;
        // Where the code's carrier sits in this baseband (an FDMA channel's offset included).
        let base_hz = spec.baseband_hz(code.carrier_hz());
        let chips_per_sample = code.chip_rate_hz() / fs;

        let plan = FftPlan::new(spc);
        let code_fft: Vec<(f64, f64)> = {
            let rep: Vec<(f64, f64)> = (0..spc)
                .map(|i| (code.value_at(i as f64 * chips_per_sample), 0.0))
                .collect();
            plan.forward(&rep)
                .into_iter()
                .map(|(r, i)| (r, -i))
                .collect()
        };

        let power: f64 = samples[..needed]
            .iter()
            .map(|x| x.re * x.re + x.im * x.im)
            .sum();
        let sigma2 = power / needed as f64;
        if sigma2.is_nan() || sigma2 <= 0.0 {
            return Err("the samples carry no power".into());
        }
        Ok(Self {
            samples,
            fs,
            base_hz,
            spc,
            block,
            noncoherent: cfg.noncoherent,
            plan,
            code_fft,
            norm: 2.0 / (block as f64 * sigma2),
            sigma2,
        })
    }

    /// The normalised power at every lag for Doppler `d` (Hz).
    fn row(&self, d: f64) -> Vec<f64> {
        let (spc, block, fs) = (self.spc, self.block, self.fs);
        let f = self.base_hz + d;
        let mut row = vec![0.0_f64; spc];
        for m in 0..self.noncoherent {
            let start = m * block;
            let mut fold = vec![(0.0_f64, 0.0_f64); spc];
            for i in 0..block {
                let n = start + i;
                let cyc = (f * n as f64 / fs).fract();
                let (sn, cs) = (-core::f64::consts::TAU * cyc).psin_cos();
                let s = self.samples[n];
                let r = &mut fold[i % spc];
                r.0 += s.re * cs - s.im * sn;
                r.1 += s.re * sn + s.im * cs;
            }
            let mut y = self.plan.forward(&fold);
            drop(fold);
            for (v, &(c, e)) in y.iter_mut().zip(&self.code_fft) {
                let (a, b) = *v;
                *v = (a * c - b * e, a * e + b * c);
            }
            for (cell, (re, im)) in row.iter_mut().zip(self.plan.inverse_owned(y)) {
                *cell += re * re + im * im;
            }
        }
        for cell in &mut row {
            *cell *= self.norm;
        }
        row
    }
}

/// The normalised correlation power `2·G/(L·σ²)` of `samples` against `code` at every
/// Doppler in `dopplers` (Hz, relative to the carrier's place in the baseband) and every
/// lag: `rows[doppler][lag]`, with the mean sample power `σ²`. [`acquire`] searches
/// `cfg.doppler_bins()` with it; the surface export evaluates finer Dopplers with the same
/// arithmetic. `samples` must hold [`samples_needed`] samples.
pub(crate) fn power_rows(
    samples: &[Cf64],
    spec: &SampleSpec,
    code: &dyn SpreadingCode,
    cfg: &AcqConfig,
    dopplers: &[f64],
) -> Result<(Vec<Vec<f64>>, f64), String> {
    let eng = RowEngine::new(samples, spec, code, cfg)?;
    let grid = dopplers.iter().map(|&d| eng.row(d)).collect();
    Ok((grid, eng.sigma2))
}

fn check_cfg(cfg: &AcqConfig) -> Result<(), String> {
    if cfg.coherent_periods == 0 || cfg.noncoherent == 0 {
        return Err("coherent_periods and noncoherent must be positive".into());
    }
    if !(cfg.doppler_step_hz > 0.0 && cfg.doppler_max_hz >= 0.0) {
        return Err("the Doppler grid needs a positive step and a non-negative maximum".into());
    }
    if !(cfg.pfa > 0.0 && cfg.pfa < 1.0) {
        return Err("pfa must lie in (0, 1)".into());
    }
    Ok(())
}

/// The result from the peak `(power, doppler index, lag)` and the peak's row.
#[allow(clippy::too_many_arguments)]
fn result_from_peak(
    code: &dyn SpreadingCode,
    cfg: &AcqConfig,
    bins: &[f64],
    spc: usize,
    sigma2: f64,
    chips_per_sample: f64,
    best: (f64, usize, usize),
    peak_row: &[f64],
) -> AcqResult {
    let (peak, jb, tb) = best;
    let guard = (1.0 / chips_per_sample).ceil() as usize + 1;
    let second_peak = peak_row
        .iter()
        .enumerate()
        .filter(|&(t, _)| {
            let d = t.abs_diff(tb);
            d.min(spc - d) > guard
        })
        .map(|(_, &v)| v)
        .fold(0.0_f64, f64::max);

    let n_cells = (spc * bins.len()) as f64;
    // Per-cell false-alarm probability 1 − (1 − P_fa)^(1/N), formed without cancellation.
    let pfa_cell = -((-cfg.pfa).pln_1p() / n_cells).pexp_m1();
    let threshold = threshold_for_pfa(pfa_cell, cfg.noncoherent as f64);
    let len = code.len_chips() as f64;
    let code_phase_chips = (-(tb as f64) * chips_per_sample).rem_euclid(len);
    AcqResult {
        code_name: code.name(),
        doppler_hz: bins[jb],
        doppler_index: jb,
        delay_samples: tb,
        code_phase_chips,
        statistic: peak,
        second_peak,
        peak_ratio: if second_peak > 0.0 {
            peak / second_peak
        } else {
            f64::INFINITY
        },
        pfa_cell,
        threshold,
        acquired: peak > threshold,
        samples_per_period: spc,
        n_doppler_bins: bins.len(),
        sample_power: sigma2,
    }
}

/// Search `samples` (at least [`samples_needed`] long; the first that many are used),
/// sampled as `spec`, for `code`, keeping the whole grid. The grid is
/// `Doppler bins × samples per period` f64 cells, which is gigabytes for the long tiered
/// codes at high rates: use [`acquire_peak`] when only the result is needed.
pub fn acquire(
    samples: &[Cf64],
    spec: &SampleSpec,
    code: &dyn SpreadingCode,
    cfg: &AcqConfig,
) -> Result<AcqGrid, String> {
    check_cfg(cfg)?;
    let bins = cfg.doppler_bins();
    let chips_per_sample = code.chip_rate_hz() / spec.fs_hz;
    let (grid, sigma2) = power_rows(samples, spec, code, cfg, &bins)?;
    let spc = samples_per_period(spec, code)?;

    let mut best = (f64::NEG_INFINITY, 0usize, 0usize);
    for (j, row) in grid.iter().enumerate() {
        for (t, &cell) in row.iter().enumerate() {
            if cell > best.0 {
                best = (cell, j, t);
            }
        }
    }
    let result = result_from_peak(
        code,
        cfg,
        &bins,
        spc,
        sigma2,
        chips_per_sample,
        best,
        &grid[best.1],
    );
    Ok(AcqGrid { result, grid })
}

/// The same search as [`acquire`] without keeping the grid: a running best and the row it
/// lies in (recomputed at the end), so memory is O(samples + samples per period), not
/// O(bins × samples). The result
/// is bit-identical to `acquire(..).result` (the cells come from the same arithmetic and ties
/// resolve to the first cell in the same order).
pub fn acquire_peak(
    samples: &[Cf64],
    spec: &SampleSpec,
    code: &dyn SpreadingCode,
    cfg: &AcqConfig,
) -> Result<AcqResult, String> {
    check_cfg(cfg)?;
    let bins = cfg.doppler_bins();
    let chips_per_sample = code.chip_rate_hz() / spec.fs_hz;
    let eng = RowEngine::new(samples, spec, code, cfg)?;
    let mut best = (f64::NEG_INFINITY, 0usize, 0usize);
    for (j, &d) in bins.iter().enumerate() {
        let row = eng.row(d);
        for (t, &cell) in row.iter().enumerate() {
            if cell > best.0 {
                best = (cell, j, t);
            }
        }
    }
    // The peak's row is recomputed (the arithmetic is deterministic) rather than held for
    // the whole search: one more row of work instead of one more row of memory.
    let best_row = eng.row(bins[best.1]);
    Ok(result_from_peak(
        code,
        cfg,
        &bins,
        eng.spc,
        eng.sigma2,
        chips_per_sample,
        best,
        &best_row,
    ))
}

/// Read [`samples_needed`] samples from `src` and search them for `code`. The search is
/// relative to the first sample read.
pub fn acquire_source(
    src: &mut dyn IqSource,
    code: &dyn SpreadingCode,
    cfg: &AcqConfig,
) -> Result<AcqGrid, IqError> {
    let spec = src.spec();
    let needed = samples_needed(&spec, code, cfg).map_err(IqError::Format)?;
    let mut buf = vec![Cf64::default(); needed];
    let mut got = 0;
    while got < needed {
        let n = src.read(&mut buf[got..])?;
        if n == 0 {
            break;
        }
        got += n;
    }
    if got < needed {
        return Err(IqError::Format(format!(
            "stream ended after {got} of {needed} samples"
        )));
    }
    acquire(&buf, &spec, code, cfg).map_err(IqError::Format)
}

/// Analytic detection probability of one correctly aligned cell of a search with `cfg`
/// on a code of period `period_s` at carrier-to-noise density `cn0_dbhz`, at per-cell
/// false-alarm probability `pfa_cell`: [`crate::acquisition::pd_square_law`] with
/// `ρ = (C/N0)·N·T_code` and `M` non-coherent blocks. No straddle or bit-edge loss.
pub fn predicted_pd(cfg: &AcqConfig, period_s: f64, cn0_dbhz: f64, pfa_cell: f64) -> f64 {
    let m = cfg.noncoherent as f64;
    let rho = 10f64.ppowf(cn0_dbhz / 10.0) * cfg.coherent_periods as f64 * period_s;
    pd_square_law(threshold_for_pfa(pfa_cell, m), m, rho)
}

// --- Cell-averaging CFAR (opt-in) ----------------------------------------------------

/// Settings of the opt-in cell-averaging CFAR detector ([`ca_cfar`]), which decides on the
/// normalised grid an [`acquire`] search returns instead of on the sample-power threshold.
///
/// A cell `X` at Doppler bin `j` and lag `ℓ` is compared with the mean of `reference_cells`
/// cells of the same row, half on each side at lags `ℓ ± (guard_lags + k·stride_lags)`,
/// `k = 1 ..= reference_cells/2`, wrapping round the period: it is declared when
/// `X > α·mean`, with `α` from [`ca_cfar_alpha`] for `pfa_cell`. The stride keeps the
/// reference cells and the cell close to independent under noise alone (it should exceed the
/// code's autocorrelation main lobe; [`CfarConfig::new`] uses `ceil(samples per chip) + 1`).
/// MODELLED; pre-registered in `docs/design/evidence/acq-cfar/`.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct CfarConfig {
    /// Number of reference cells `N` (even, at least 2).
    pub reference_cells: usize,
    /// Lags between the cell and the first reference cell on each side, before the stride.
    pub guard_lags: usize,
    /// Lags between consecutive reference cells (at least 1).
    pub stride_lags: usize,
    /// Per-cell false-alarm probability under noise alone.
    pub pfa_cell: f64,
}

impl CfarConfig {
    /// `reference_cells` cells, `guard_lags` guard, and the stride `ceil(fs/chip rate) + 1`
    /// for `code` sampled as `spec`.
    pub fn new(
        spec: &SampleSpec,
        code: &dyn SpreadingCode,
        reference_cells: usize,
        guard_lags: usize,
        pfa_cell: f64,
    ) -> Self {
        Self {
            reference_cells,
            guard_lags,
            stride_lags: (spec.fs_hz / code.chip_rate_hz()).ceil() as usize + 1,
            pfa_cell,
        }
    }

    fn check(&self, samples_per_period: usize) -> Result<(), String> {
        let n = self.reference_cells;
        if n < 2 || n % 2 != 0 {
            return Err(format!(
                "CFAR reference cells must be even and at least 2 (got {n})"
            ));
        }
        if self.stride_lags == 0 {
            return Err("CFAR stride must be at least 1 lag".into());
        }
        if !(self.pfa_cell > 0.0 && self.pfa_cell < 1.0) {
            return Err(format!(
                "CFAR per-cell Pfa must lie in (0, 1) (got {})",
                self.pfa_cell
            ));
        }
        let reach = self.guard_lags + (n / 2) * self.stride_lags;
        if 2 * reach >= samples_per_period {
            return Err(format!(
                "CFAR window reaches {reach} lags each side, too wide for a {samples_per_period}-lag period"
            ));
        }
        Ok(())
    }

    /// The mean of the reference cells of lag `lag` in `row` (one period of lags).
    fn reference_mean(&self, row: &[f64], lag: usize) -> f64 {
        let len = row.len();
        let half = self.reference_cells / 2;
        let mut sum = 0.0;
        for k in 1..=half {
            let d = self.guard_lags + k * self.stride_lags;
            sum += row[(lag + d) % len] + row[(lag + len - d % len) % len];
        }
        sum / self.reference_cells as f64
    }
}

/// The closed-form per-cell false-alarm probability of the cell-averaging CFAR test
/// `X > α·Z/N` when the cell `X` is `Gamma(M, θ)` and the `N` reference cells are independent
/// `Gamma(M, θ)` (square-law cells after `M` non-coherent blocks; `θ` cancels):
/// `Σ_{k=0}^{M−1} C(NM + k − 1, k)·t^k/(1 + t)^(NM + k)` with `t = α/N`, which for `M = 1` is
/// `(1 + α/N)^(−N)`.
pub fn ca_cfar_pfa(alpha: f64, n_ref: usize, m: usize) -> f64 {
    let nm = (n_ref * m) as f64;
    let t = alpha / n_ref as f64;
    // term_k = C(NM + k − 1, k)·t^k/(1 + t)^(NM + k), built by ratio from term_0.
    let mut term = (-nm * t.pln_1p()).pexp();
    let mut sum = term;
    for k in 1..m {
        term *= (nm + k as f64 - 1.0) / k as f64 * t / (1.0 + t);
        sum += term;
    }
    sum
}

/// The CFAR factor `α` whose [`ca_cfar_pfa`] equals `pfa` (bisection to 1e-12 relative).
pub fn ca_cfar_alpha(pfa: f64, n_ref: usize, m: usize) -> f64 {
    let (mut lo, mut hi) = (0.0f64, 1.0f64);
    while ca_cfar_pfa(hi, n_ref, m) > pfa {
        hi *= 2.0;
    }
    for _ in 0..200 {
        let mid = 0.5 * (lo + hi);
        if ca_cfar_pfa(mid, n_ref, m) > pfa {
            lo = mid;
        } else {
            hi = mid;
        }
        if hi - lo <= 1e-12 * hi {
            break;
        }
    }
    0.5 * (lo + hi)
}

/// The cell-averaging CFAR decision on a search's peak cell.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct CfarResult {
    /// The factor `α` the cells are compared with.
    pub alpha: f64,
    /// The peak cell over the mean of its reference cells.
    pub peak_ratio: f64,
    /// Whether `peak_ratio > alpha`.
    pub detected: bool,
    /// The peak's Doppler bin and lag (the same cell as the search result's).
    pub doppler_index: usize,
    /// The peak's lag (samples).
    pub delay_samples: usize,
}

/// The cell-averaging CFAR decision on the peak of `grid`, a search made with `cfg`.
pub fn ca_cfar(grid: &AcqGrid, cfg: &AcqConfig, cfar: &CfarConfig) -> Result<CfarResult, String> {
    let r = &grid.result;
    cfar.check(r.samples_per_period)?;
    let alpha = ca_cfar_alpha(cfar.pfa_cell, cfar.reference_cells, cfg.noncoherent);
    let row = &grid.grid[r.doppler_index];
    let peak_ratio = row[r.delay_samples] / cfar.reference_mean(row, r.delay_samples);
    Ok(CfarResult {
        alpha,
        peak_ratio,
        detected: peak_ratio > alpha,
        doppler_index: r.doppler_index,
        delay_samples: r.delay_samples,
    })
}

/// The cell-averaging CFAR decision on every cell of `grid` (a search made with `cfg`):
/// `decisions[j][ℓ]` is true when cell `(j, ℓ)` exceeds `α` times its reference mean.
pub fn ca_cfar_cells(
    grid: &AcqGrid,
    cfg: &AcqConfig,
    cfar: &CfarConfig,
) -> Result<Vec<Vec<bool>>, String> {
    cfar.check(grid.result.samples_per_period)?;
    let alpha = ca_cfar_alpha(cfar.pfa_cell, cfar.reference_cells, cfg.noncoherent);
    Ok(grid
        .grid
        .iter()
        .map(|row| {
            (0..row.len())
                .map(|l| row[l] > alpha * cfar.reference_mean(row, l))
                .collect()
        })
        .collect())
}

#[cfg(test)]
mod tests {
    use super::{
        auto_coherent_periods, ca_cfar_alpha, ca_cfar_pfa, default_step_hz, PULL_IN_STEP_MARGIN,
    };
    use crate::portable_math::PortableFloat;

    #[test]
    fn ca_cfar_closed_form_reduces_to_the_textbook_case_and_inverts() {
        // M = 1: (1 + α/N)^(−N).
        for n in [8usize, 16, 32] {
            for a in [1.0, 5.0, 12.5] {
                let want = (1.0 + a / n as f64).ppowi(-(n as i32));
                assert!((ca_cfar_pfa(a, n, 1) - want).abs() <= 1e-15 * want.max(1e-300));
            }
        }
        // α = 0: the cell exceeds zero with probability one, for every M.
        assert!((ca_cfar_pfa(0.0, 16, 4) - 1.0).abs() < 1e-12);
        // The inversion recovers the requested Pfa.
        for m in [1usize, 2, 4] {
            for pfa in [1e-2, 1e-3, 1e-6] {
                let a = ca_cfar_alpha(pfa, 16, m);
                assert!(
                    (ca_cfar_pfa(a, 16, m) / pfa - 1.0).abs() < 1e-9,
                    "{m} {pfa}"
                );
            }
        }
    }

    #[test]
    fn the_auto_coherent_length_is_about_four_milliseconds() {
        assert_eq!(auto_coherent_periods(1.0e-3), 4); // GPS L1 C/A, B1I, L1OF
        assert_eq!(auto_coherent_periods(4.0e-3), 1); // Galileo E1-B
        assert_eq!(auto_coherent_periods(10.0e-3), 1); // BeiDou B1C, GPS L5-I (tiered)
        assert_eq!(auto_coherent_periods(100.0e-3), 1); // Galileo E1-C, E5a-Q (tiered)
        assert_eq!(auto_coherent_periods(1.5e-3), 3);
        assert_eq!(auto_coherent_periods(0.0), 1);
    }

    /// Pre-registered bars (D8): for every signal, the hand-off residual is at most 0.8 of the
    /// default Atan2 FLL's pull-in, 1 / (4 T), both when the nearest bin wins (half the default
    /// step) and when the neighbouring bin wins (a whole step), for
    /// the full-period code and for the primary-only replica, with the default one-period
    /// loop update and with five periods per update.
    #[test]
    fn the_default_step_leaves_a_residual_inside_the_fll_pull_in() {
        use crate::iq::signals::{beidou, galileo, glonass, gps};
        use crate::iq::SpreadingCode;
        let mut codes = vec![
            gps::l1ca(3).unwrap(),
            gps::l5_i5(3).unwrap(),
            gps::l5_q5(3).unwrap(),
            gps::l2c_cm(3).unwrap(),
            galileo::e5a_i(3).unwrap(),
            galileo::e5a_q(3).unwrap(),
            beidou::b1i(6).unwrap(),
            beidou::b1c_data(19).unwrap(),
            glonass::l1of(1).unwrap(),
        ];
        // Galileo E1 only when the user has supplied the ICD tables (they are not bundled).
        codes.extend(galileo::e1b(3));
        codes.extend(galileo::e1c(3));
        let mut worst: f64 = 0.0;
        for code in &codes {
            for rx in [code.clone(), code.primary_only()] {
                let t_code = rx.period_s();
                let n = auto_coherent_periods(t_code);
                for periods in [1usize, 5] {
                    let t_track = periods as f64 * t_code;
                    let step = default_step_hz(t_code, n, t_track);
                    let pull_in = 1.0 / (4.0 * t_track);
                    // Bar 1: the nearest bin wins, residual half a step.
                    // Bar 2: the neighbouring bin wins, residual a whole step.
                    worst = worst.max(step / pull_in);
                    for (what, residual) in [("half-step", step / 2.0), ("neighbour", step)] {
                        assert!(
                            residual <= 0.8 * pull_in + 1e-9,
                            "{} ({periods} period(s)): {what} residual {residual} Hz vs pull-in {pull_in} Hz",
                            rx.name()
                        );
                    }
                }
            }
        }
        assert!(worst > 0.0);
    }

    #[test]
    fn the_cap_changes_nothing_for_the_one_millisecond_codes() {
        // 1 ms codes: N = 4, the textbook 166.7 Hz step is already below 0.2 / T = 200 Hz.
        let t = 1.0e-3;
        assert_eq!(default_step_hz(t, 4, t), 2.0 / (3.0 * 4.0 * t));
        // E1-B: N = 1, T = 4 ms: the cap (50 Hz) replaces the textbook 166.7 Hz.
        assert!((default_step_hz(4.0e-3, 1, 4.0e-3) - 50.0).abs() < 1e-9);
        assert_eq!(PULL_IN_STEP_MARGIN, 0.2);
        // No usable tracking time: the textbook step.
        assert_eq!(default_step_hz(4.0e-3, 1, 0.0), 2.0 / (3.0 * 4.0e-3));
    }
}