phasesmith-core 0.3.0

//! Neutron time-of-flight calibration, profile parameters, and accumulation.

use std::error::Error;
use std::fmt::{Display, Formatter};
use std::sync::Arc;

use crate::fcj::{QUADRATURE_NODES, QUADRATURE_ORDER, QUADRATURE_WEIGHTS};
use crate::profile::{
    Accumulation, DenseJacobian, GridView, PatternDerivatives, ProfileError, SupportJacobian,
    SupportRange, zeroed_f64_vec,
};
use crate::tch::{TchError, TchShape, TchWidths};
use phasesmith_execution::ExecutionContext;

const GAUSSIAN_FWHM_PER_SIGMA: f64 = 2.354_820_045_030_949_3;
/// Public order: zero, difC, difA, difB, alpha, beta0, beta1, betaq,
/// sigma0, sigma1, sigma2, sigmaq, X, Y, Z.
pub const TOF_GLOBAL_PARAMETER_COUNT: usize = 15;
/// Stable names matching the dense global derivative-row order.
pub const TOF_GLOBAL_PARAMETER_NAMES: [&str; TOF_GLOBAL_PARAMETER_COUNT] = [
    "zero", "difc", "difa", "difb", "alpha", "beta0", "beta1", "betaq", "sigma0", "sigma1",
    "sigma2", "sigmaq", "x", "y", "z",
];
/// Number of coefficients in the Maxwellian-plus-Chebyshev incident spectrum.
pub const TOF_INCIDENT_SPECTRUM_COEFFICIENT_COUNT: usize = 12;

/// One selectable TOF calibration/profile coefficient.
#[repr(usize)]
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub enum TofInstrumentParameter {
    /// Additive time zero.
    Zero,
    /// Linear d-to-TOF calibration.
    Difc,
    /// Quadratic d-to-TOF calibration.
    Difa,
    /// Reciprocal d-to-TOF calibration.
    Difb,
    /// Leading-edge rate numerator.
    Alpha,
    /// Constant trailing-edge rate.
    Beta0,
    /// Inverse-fourth-power trailing-edge rate.
    Beta1,
    /// Inverse-square trailing-edge rate.
    Betaq,
    /// Constant Gaussian variance.
    Sigma0,
    /// Quadratic-d Gaussian variance.
    Sigma1,
    /// Quartic-d Gaussian variance.
    Sigma2,
    /// Linear-d Gaussian variance.
    Sigmaq,
    /// Linear-d Lorentzian width.
    X,
    /// Quadratic-d Lorentzian width.
    Y,
    /// Constant Lorentzian width.
    Z,
}

impl TofInstrumentParameter {
    /// Every parameter in the fused dense-Jacobian order.
    pub const ALL: [Self; TOF_GLOBAL_PARAMETER_COUNT] = [
        Self::Zero,
        Self::Difc,
        Self::Difa,
        Self::Difb,
        Self::Alpha,
        Self::Beta0,
        Self::Beta1,
        Self::Betaq,
        Self::Sigma0,
        Self::Sigma1,
        Self::Sigma2,
        Self::Sigmaq,
        Self::X,
        Self::Y,
        Self::Z,
    ];

    /// Dense global derivative-row index.
    #[must_use]
    pub const fn index(self) -> usize {
        self as usize
    }

    /// Stable short name.
    #[must_use]
    pub const fn name(self) -> &'static str {
        TOF_GLOBAL_PARAMETER_NAMES[self.index()]
    }
}
const LOCAL_PARAMETER_COUNT: usize = 2;
const TOF_QUADRATURE_PANELS: usize = 8;
const TOF_QUADRATURE_PANELS_F64: f64 = 8.0;
const TOF_SUPPORT_QUADRATURE_PANELS: usize = 4;
const TOF_SUPPORT_QUADRATURE_PANELS_F64: f64 = 4.0;
const TOF_QUADRATURE_COUNT: usize = TOF_QUADRATURE_PANELS * QUADRATURE_ORDER;

/// One calibrated TOF incident-spectrum value and coordinate derivative.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct TofIncidentSpectrumPoint {
    /// Positive incident intensity at the requested TOF.
    pub value: f64,
    /// Analytical derivative with respect to TOF in microseconds.
    pub d_value_d_tof_us: f64,
}

/// Facility-neutral Maxwellian-plus-Chebyshev TOF incident spectrum.
///
/// The twelve coefficients are `P1..P12`. With TOF `t` in milliseconds and
/// `x = 2/t - 1`, the calibrated intensity is
///
/// `P1 + P2 t^-5 exp(-P3/t^2) + sum(Pj T_(j-3)(x), j=4..12)`.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct TofIncidentSpectrum {
    /// Inclusive lower validity bound in microseconds.
    pub min_tof_us: f64,
    /// Inclusive upper validity bound in microseconds.
    pub max_tof_us: f64,
    /// Coefficients `P1..P12` in the documented function order.
    pub coefficients: [f64; TOF_INCIDENT_SPECTRUM_COEFFICIENT_COUNT],
}

impl TofIncidentSpectrum {
    /// Construct and validate a fixed incident-spectrum calibration.
    ///
    /// # Errors
    ///
    /// Returns [`TofIncidentSpectrumError`] for an invalid range or
    /// non-finite coefficient.
    pub fn new(
        min_tof_us: f64,
        max_tof_us: f64,
        coefficients: [f64; TOF_INCIDENT_SPECTRUM_COEFFICIENT_COUNT],
    ) -> Result<Self, TofIncidentSpectrumError> {
        if !min_tof_us.is_finite()
            || !max_tof_us.is_finite()
            || min_tof_us <= 0.0
            || max_tof_us <= min_tof_us
        {
            return Err(TofIncidentSpectrumError::InvalidRange);
        }
        if coefficients.iter().any(|value| !value.is_finite()) {
            return Err(TofIncidentSpectrumError::NonFiniteCoefficient);
        }
        Ok(Self {
            min_tof_us,
            max_tof_us,
            coefficients,
        })
    }

    /// Evaluate the calibrated intensity and analytical TOF derivative.
    ///
    /// The validity interval is inclusive at both ends.
    ///
    /// # Errors
    ///
    /// Returns [`TofIncidentSpectrumError`] when TOF lies outside the
    /// calibration interval or the fitted function is not positive and finite.
    pub fn evaluate(
        self,
        tof_us: f64,
    ) -> Result<TofIncidentSpectrumPoint, TofIncidentSpectrumError> {
        if !tof_us.is_finite() || tof_us < self.min_tof_us || tof_us > self.max_tof_us {
            return Err(TofIncidentSpectrumError::TofOutsideRange);
        }
        let time_milliseconds = tof_us / 1_000.0;
        let inverse_t = time_milliseconds.recip();
        let inverse_t2 = inverse_t * inverse_t;
        let x = 2.0 * inverse_t - 1.0;
        let d_x_d_t_ms = -2.0 * inverse_t2;
        let maxwell =
            self.coefficients[1] * inverse_t.powi(5) * (-self.coefficients[2] * inverse_t2).exp();
        let mut value = self.coefficients[0] + maxwell;
        let mut d_value_d_t_ms =
            maxwell * (-5.0 * inverse_t + 2.0 * self.coefficients[2] * inverse_t.powi(3));

        let mut previous = 1.0;
        let mut d_previous = 0.0;
        let mut current = x;
        let mut d_current = d_x_d_t_ms;
        for (index, &coefficient) in self.coefficients[3..].iter().enumerate() {
            if index > 0 {
                let next = 2.0 * x * current - previous;
                let d_next = 2.0 * (d_x_d_t_ms * current + x * d_current) - d_previous;
                previous = current;
                d_previous = d_current;
                current = next;
                d_current = d_next;
            }
            value += coefficient * current;
            d_value_d_t_ms += coefficient * d_current;
        }
        let d_value_d_tof_us = d_value_d_t_ms / 1_000.0;
        if !value.is_finite() || value <= 0.0 || !d_value_d_tof_us.is_finite() {
            return Err(TofIncidentSpectrumError::NonPositiveIntensity);
        }
        Ok(TofIncidentSpectrumPoint {
            value,
            d_value_d_tof_us,
        })
    }
}

/// Invalid TOF incident-spectrum calibration or evaluation.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum TofIncidentSpectrumError {
    /// The inclusive validity interval is not finite, positive, and increasing.
    InvalidRange,
    /// At least one fitted coefficient is non-finite.
    NonFiniteCoefficient,
    /// Requested TOF is non-finite or outside the inclusive validity interval.
    TofOutsideRange,
    /// Evaluation produced an incident intensity that is not positive and finite.
    NonPositiveIntensity,
}

impl Display for TofIncidentSpectrumError {
    fn fmt(&self, formatter: &mut Formatter<'_>) -> std::fmt::Result {
        formatter.write_str(match self {
            Self::InvalidRange => {
                "TOF incident-spectrum range must be finite, positive, and increasing"
            }
            Self::NonFiniteCoefficient => "TOF incident-spectrum coefficients must be finite",
            Self::TofOutsideRange => "TOF lies outside the incident-spectrum validity interval",
            Self::NonPositiveIntensity => {
                "TOF incident-spectrum intensity must be positive and finite"
            }
        })
    }
}

impl Error for TofIncidentSpectrumError {}

/// Facility-neutral fixed geometry for one focused TOF detector bank.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct TofBankGeometry {
    /// Nominal bank scattering angle in degrees `2theta`.
    pub two_theta_deg: f64,
}

impl TofBankGeometry {
    /// Validate the strict physical scattering-angle domain.
    ///
    /// # Errors
    ///
    /// Returns [`TofError::InvalidBankTwoTheta`] unless the angle is finite and
    /// strictly within `0 < 2theta < 180°`.
    pub fn validate(self) -> Result<(), TofError> {
        if !self.two_theta_deg.is_finite()
            || self.two_theta_deg <= 0.0
            || self.two_theta_deg >= 180.0
        {
            return Err(TofError::InvalidBankTwoTheta);
        }
        Ok(())
    }

    /// Return the half scattering angle in radians after validation.
    ///
    /// # Errors
    ///
    /// Returns [`TofError::InvalidBankTwoTheta`] for invalid geometry.
    pub fn theta_radians(self) -> Result<f64, TofError> {
        self.validate()?;
        Ok((0.5 * self.two_theta_deg).to_radians())
    }
}

/// TOF calibration and d-dependent profile coefficients in public units.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct TofInstrument {
    /// Additive time zero in microseconds.
    pub zero_us: f64,
    /// Linear calibration coefficient in microseconds per ångström.
    pub difc_us_per_angstrom: f64,
    /// Quadratic calibration coefficient in microseconds per ångström squared.
    pub difa_us_per_angstrom2: f64,
    /// Reciprocal calibration coefficient in microsecond ångströms.
    pub difb_us_angstrom: f64,
    /// Exponential-rise numerator; `alpha = alpha_coefficient / d`.
    pub alpha_coefficient: f64,
    /// Constant exponential-decay rate term in inverse microseconds.
    pub beta0_per_us: f64,
    /// `d^-4` exponential-decay coefficient.
    pub beta1_angstrom4_per_us: f64,
    /// `d^-2` exponential-decay coefficient.
    pub betaq_angstrom2_per_us: f64,
    /// Constant Gaussian variance term in microseconds squared.
    pub sigma0_us2: f64,
    /// `d^2` Gaussian variance coefficient.
    pub sigma1_us2_per_angstrom2: f64,
    /// `d^4` Gaussian variance coefficient.
    pub sigma2_us2_per_angstrom4: f64,
    /// Linear-d Gaussian variance coefficient.
    pub sigmaq_us2_per_angstrom: f64,
    /// Linear-d Lorentzian FWHM coefficient.
    pub x_us_per_angstrom: f64,
    /// Quadratic-d Lorentzian FWHM coefficient.
    pub y_us_per_angstrom2: f64,
    /// Constant Lorentzian FWHM in microseconds.
    pub z_us: f64,
}

impl TofInstrument {
    /// Coefficients in the stable dense global derivative-row order.
    #[must_use]
    pub const fn values(self) -> [f64; TOF_GLOBAL_PARAMETER_COUNT] {
        [
            self.zero_us,
            self.difc_us_per_angstrom,
            self.difa_us_per_angstrom2,
            self.difb_us_angstrom,
            self.alpha_coefficient,
            self.beta0_per_us,
            self.beta1_angstrom4_per_us,
            self.betaq_angstrom2_per_us,
            self.sigma0_us2,
            self.sigma1_us2_per_angstrom2,
            self.sigma2_us2_per_angstrom4,
            self.sigmaq_us2_per_angstrom,
            self.x_us_per_angstrom,
            self.y_us_per_angstrom2,
            self.z_us,
        ]
    }

    /// Construct from coefficients in the stable dense global-row order.
    ///
    /// # Errors
    ///
    /// Returns [`TofError`] when the complete instrument is invalid.
    pub fn from_values(values: [f64; TOF_GLOBAL_PARAMETER_COUNT]) -> Result<Self, TofError> {
        let result = Self {
            zero_us: values[0],
            difc_us_per_angstrom: values[1],
            difa_us_per_angstrom2: values[2],
            difb_us_angstrom: values[3],
            alpha_coefficient: values[4],
            beta0_per_us: values[5],
            beta1_angstrom4_per_us: values[6],
            betaq_angstrom2_per_us: values[7],
            sigma0_us2: values[8],
            sigma1_us2_per_angstrom2: values[9],
            sigma2_us2_per_angstrom4: values[10],
            sigmaq_us2_per_angstrom: values[11],
            x_us_per_angstrom: values[12],
            y_us_per_angstrom2: values[13],
            z_us: values[14],
        };
        result.validate()?;
        Ok(result)
    }

    /// Replace one coefficient without bypassing instrument validation.
    ///
    /// # Errors
    ///
    /// Returns [`TofError`] when the replacement makes the instrument invalid.
    pub fn with_parameter(
        self,
        parameter: TofInstrumentParameter,
        value: f64,
    ) -> Result<Self, TofError> {
        let mut values = self.values();
        values[parameter.index()] = value;
        Self::from_values(values)
    }

    /// Validate finite coefficients and a positive linear calibration scale.
    ///
    /// # Errors
    ///
    /// Returns [`TofError`] for non-finite coefficients or non-positive
    /// linear calibration.
    pub fn validate(self) -> Result<(), TofError> {
        let values = [
            self.zero_us,
            self.difc_us_per_angstrom,
            self.difa_us_per_angstrom2,
            self.difb_us_angstrom,
            self.alpha_coefficient,
            self.beta0_per_us,
            self.beta1_angstrom4_per_us,
            self.betaq_angstrom2_per_us,
            self.sigma0_us2,
            self.sigma1_us2_per_angstrom2,
            self.sigma2_us2_per_angstrom4,
            self.sigmaq_us2_per_angstrom,
            self.x_us_per_angstrom,
            self.y_us_per_angstrom2,
            self.z_us,
        ];
        if values.iter().any(|value| !value.is_finite()) {
            return Err(TofError::NonFiniteInstrumentParameter);
        }
        if self.difc_us_per_angstrom <= 0.0 {
            return Err(TofError::NonPositiveDifc);
        }
        Ok(())
    }
}

/// Derived TOF position, widths, exponential rates, and analytical chains.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct TofProfileParameters {
    /// Calibrated position in microseconds.
    pub position_us: f64,
    /// Leading-edge exponential rate in inverse microseconds.
    pub alpha_per_us: f64,
    /// Trailing-edge exponential rate in inverse microseconds.
    pub beta_per_us: f64,
    /// Gaussian variance in microseconds squared.
    pub gaussian_variance_us2: f64,
    /// Gaussian component FWHM in microseconds.
    pub gaussian_fwhm_us: f64,
    /// Lorentzian component FWHM in microseconds.
    pub lorentzian_fwhm_us: f64,
    /// TCH transform of the component widths.
    pub tch: TchShape,
    /// Direct-input derivatives with respect to d-spacing.
    pub d_position_d_d: f64,
    /// Alpha-rate derivative with respect to d-spacing.
    pub d_alpha_d_d: f64,
    /// Beta-rate derivative with respect to d-spacing.
    pub d_beta_d_d: f64,
    /// Gaussian-FWHM derivative with respect to d-spacing.
    pub d_gaussian_fwhm_d_d: f64,
    /// Lorentzian-FWHM derivative with respect to d-spacing.
    pub d_lorentzian_fwhm_d_d: f64,
    /// Parameter-major chains in [`TOF_GLOBAL_PARAMETER_COUNT`] order.
    pub d_position_d_instrument: [f64; TOF_GLOBAL_PARAMETER_COUNT],
    /// Alpha-rate chains in global instrument order.
    pub d_alpha_d_instrument: [f64; TOF_GLOBAL_PARAMETER_COUNT],
    /// Beta-rate chains in global instrument order.
    pub d_beta_d_instrument: [f64; TOF_GLOBAL_PARAMETER_COUNT],
    /// Gaussian-FWHM chains in global instrument order.
    pub d_gaussian_fwhm_d_instrument: [f64; TOF_GLOBAL_PARAMETER_COUNT],
    /// Lorentzian-FWHM chains in global instrument order.
    pub d_lorentzian_fwhm_d_instrument: [f64; TOF_GLOBAL_PARAMETER_COUNT],
}

impl TofProfileParameters {
    /// Derive all profile quantities for one positive d-spacing.
    ///
    /// # Errors
    ///
    /// Returns [`TofError`] for invalid inputs or nonphysical derived profile
    /// parameters.
    pub fn from_instrument(d: f64, instrument: TofInstrument) -> Result<Self, TofError> {
        instrument.validate()?;
        if !d.is_finite() || d <= 0.0 {
            return Err(TofError::InvalidDSpacing);
        }
        let d2 = d * d;
        let d3 = d2 * d;
        let d4 = d2 * d2;
        let inverse_d = d.recip();
        let inverse_d2 = inverse_d * inverse_d;
        let inverse_d3 = inverse_d2 * inverse_d;
        let inverse_d4 = inverse_d2 * inverse_d2;
        let inverse_d5 = inverse_d4 * inverse_d;
        let position_us = instrument.zero_us
            + instrument.difc_us_per_angstrom * d
            + instrument.difa_us_per_angstrom2 * d2
            + instrument.difb_us_angstrom * inverse_d;
        let alpha_per_us = instrument.alpha_coefficient * inverse_d;
        let beta_per_us = instrument.beta0_per_us
            + instrument.beta1_angstrom4_per_us * inverse_d4
            + instrument.betaq_angstrom2_per_us * inverse_d2;
        let gaussian_variance_us2 = instrument.sigma0_us2
            + instrument.sigma1_us2_per_angstrom2 * d2
            + instrument.sigma2_us2_per_angstrom4 * d4
            + instrument.sigmaq_us2_per_angstrom * d;
        let lorentzian_fwhm_us =
            instrument.z_us + instrument.x_us_per_angstrom * d + instrument.y_us_per_angstrom2 * d2;
        if !position_us.is_finite() {
            return Err(TofError::InvalidPosition);
        }
        if !alpha_per_us.is_finite() || alpha_per_us <= 0.0 {
            return Err(TofError::NonPositiveAlpha);
        }
        if !beta_per_us.is_finite() || beta_per_us <= 0.0 {
            return Err(TofError::NonPositiveBeta);
        }
        if !gaussian_variance_us2.is_finite() || gaussian_variance_us2 <= 0.0 {
            return Err(TofError::NonPositiveGaussianVariance);
        }
        if !lorentzian_fwhm_us.is_finite() || lorentzian_fwhm_us < 0.0 {
            return Err(TofError::NegativeLorentzianFwhm);
        }
        let sigma = gaussian_variance_us2.sqrt();
        let gaussian_fwhm_us = GAUSSIAN_FWHM_PER_SIGMA * sigma;
        let tch = TchShape::from_component_fwhm(TchWidths {
            gaussian_fwhm: gaussian_fwhm_us,
            lorentzian_fwhm: lorentzian_fwhm_us,
        })
        .map_err(|reason| TofError::InvalidTch { reason })?;
        let d_gaussian_d_variance = GAUSSIAN_FWHM_PER_SIGMA / (2.0 * sigma);
        let d_variance_d_d = 2.0 * instrument.sigma1_us2_per_angstrom2 * d
            + 4.0 * instrument.sigma2_us2_per_angstrom4 * d3
            + instrument.sigmaq_us2_per_angstrom;

        let mut d_position_d_instrument = [0.0; TOF_GLOBAL_PARAMETER_COUNT];
        d_position_d_instrument[..4].copy_from_slice(&[1.0, d, d2, inverse_d]);
        let mut d_alpha_d_instrument = [0.0; TOF_GLOBAL_PARAMETER_COUNT];
        d_alpha_d_instrument[4] = inverse_d;
        let mut d_beta_d_instrument = [0.0; TOF_GLOBAL_PARAMETER_COUNT];
        d_beta_d_instrument[5..8].copy_from_slice(&[1.0, inverse_d4, inverse_d2]);
        let mut d_gaussian_fwhm_d_instrument = [0.0; TOF_GLOBAL_PARAMETER_COUNT];
        d_gaussian_fwhm_d_instrument[8..12].copy_from_slice(&[
            d_gaussian_d_variance,
            d_gaussian_d_variance * d2,
            d_gaussian_d_variance * d4,
            d_gaussian_d_variance * d,
        ]);
        let mut d_lorentzian_fwhm_d_instrument = [0.0; TOF_GLOBAL_PARAMETER_COUNT];
        d_lorentzian_fwhm_d_instrument[12..15].copy_from_slice(&[d, d2, 1.0]);
        Ok(Self {
            position_us,
            alpha_per_us,
            beta_per_us,
            gaussian_variance_us2,
            gaussian_fwhm_us,
            lorentzian_fwhm_us,
            tch,
            d_position_d_d: instrument.difc_us_per_angstrom
                + 2.0 * instrument.difa_us_per_angstrom2 * d
                - instrument.difb_us_angstrom * inverse_d2,
            d_alpha_d_d: -instrument.alpha_coefficient * inverse_d2,
            d_beta_d_d: -4.0 * instrument.beta1_angstrom4_per_us * inverse_d5
                - 2.0 * instrument.betaq_angstrom2_per_us * inverse_d3,
            d_gaussian_fwhm_d_d: d_gaussian_d_variance * d_variance_d_d,
            d_lorentzian_fwhm_d_d: instrument.x_us_per_angstrom
                + 2.0 * instrument.y_us_per_angstrom2 * d,
            d_position_d_instrument,
            d_alpha_d_instrument,
            d_beta_d_instrument,
            d_gaussian_fwhm_d_instrument,
            d_lorentzian_fwhm_d_instrument,
        })
    }
}

/// One asymmetric TOF profile value and direct-input derivatives.
#[derive(Clone, Copy, Debug, Default, PartialEq)]
pub struct TofProfilePoint {
    /// Unit-area profile density in inverse microseconds.
    pub value: f64,
    /// Derivative with respect to ideal position in microseconds.
    pub d_position: f64,
    /// Derivative with respect to the leading-edge alpha rate.
    pub d_alpha: f64,
    /// Derivative with respect to the trailing-edge beta rate.
    pub d_beta: f64,
    /// Derivative with respect to Gaussian component FWHM.
    pub d_gaussian_fwhm: f64,
    /// Derivative with respect to Lorentzian component FWHM.
    pub d_lorentzian_fwhm: f64,
}

const SUPPORTED_DELTA: usize = 0;
const SUPPORTED_ALPHA: usize = 1;
const SUPPORTED_BETA: usize = 2;
const SUPPORTED_GAUSSIAN: usize = 3;
const SUPPORTED_LORENTZIAN: usize = 4;
const SUPPORTED_VARIABLE_COUNT: usize = 5;

#[derive(Clone, Copy, Default)]
struct SupportedScalar {
    value: f64,
    derivative: [f64; SUPPORTED_VARIABLE_COUNT],
}

impl SupportedScalar {
    fn clamped(self, lower: f64, upper: f64) -> Self {
        if lower < self.value && self.value < upper {
            self
        } else {
            Self {
                value: self.value.clamp(lower, upper),
                derivative: [0.0; SUPPORTED_VARIABLE_COUNT],
            }
        }
    }
}

/// Prepared truncated double-exponential convolution of a TCH profile.
#[derive(Clone, Debug)]
pub struct TofProfile {
    shape: TchShape,
    alpha: f64,
    beta: f64,
    quadrature: Arc<TofQuadrature>,
}

#[derive(Debug)]
struct TofQuadrature {
    tail_log: f64,
    nodes: [f64; TOF_QUADRATURE_COUNT],
    weights: [f64; TOF_QUADRATURE_COUNT],
}

impl TofQuadrature {
    fn new(tail_log: f64) -> Result<Self, TofError> {
        if !tail_log.is_finite() || tail_log <= 0.0 {
            return Err(TofError::InvalidTailLog);
        }
        let mut nodes = [0.0; TOF_QUADRATURE_COUNT];
        let mut weights = [0.0; TOF_QUADRATURE_COUNT];
        let mut normalization = 0.0;
        let panel_scale = TOF_QUADRATURE_PANELS_F64.recip();
        let mut panel_offset = 0.0;
        for panel in 0..TOF_QUADRATURE_PANELS {
            for quadrature in 0..QUADRATURE_ORDER {
                let index = panel * QUADRATURE_ORDER + quadrature;
                let unit_node = panel_offset + panel_scale * QUADRATURE_NODES[quadrature];
                nodes[index] = tail_log * unit_node;
                weights[index] =
                    tail_log * panel_scale * QUADRATURE_WEIGHTS[quadrature] * (-nodes[index]).exp();
                normalization += weights[index];
            }
            panel_offset += panel_scale;
        }
        if !normalization.is_finite() || normalization <= 0.0 {
            return Err(TofError::InvalidQuadrature);
        }
        for weight in &mut weights {
            *weight /= normalization;
        }
        Ok(Self {
            tail_log,
            nodes,
            weights,
        })
    }
}

impl TofProfile {
    /// Prepare a unit-area profile. Both exponential tails are truncated at
    /// `exp(-tail_log)` and renormalized before convolution.
    ///
    /// # Errors
    ///
    /// Returns [`TofError`] for invalid rates, widths, tail cutoff, or
    /// quadrature normalization.
    pub fn new(
        alpha_per_us: f64,
        beta_per_us: f64,
        widths: TchWidths,
        tail_log: f64,
    ) -> Result<Self, TofError> {
        Self::validate_rates(alpha_per_us, beta_per_us)?;
        let quadrature = Arc::new(TofQuadrature::new(tail_log)?);
        Self::from_validated_rates(alpha_per_us, beta_per_us, widths, quadrature)
    }

    fn validate_rates(alpha_per_us: f64, beta_per_us: f64) -> Result<(), TofError> {
        if !alpha_per_us.is_finite() || alpha_per_us <= 0.0 {
            return Err(TofError::NonPositiveAlpha);
        }
        if !beta_per_us.is_finite() || beta_per_us <= 0.0 {
            return Err(TofError::NonPositiveBeta);
        }
        Ok(())
    }

    fn from_validated_rates(
        alpha_per_us: f64,
        beta_per_us: f64,
        widths: TchWidths,
        quadrature: Arc<TofQuadrature>,
    ) -> Result<Self, TofError> {
        let shape = TchShape::from_component_fwhm(widths)
            .map_err(|reason| TofError::InvalidTch { reason })?;
        Ok(Self {
            shape,
            alpha: alpha_per_us,
            beta: beta_per_us,
            quadrature,
        })
    }

    /// Evaluate without finite TCH support truncation.
    #[must_use]
    pub fn evaluate(&self, x_minus_position_us: f64) -> TofProfilePoint {
        self.evaluate_with_radius(x_minus_position_us, f64::INFINITY)
    }

    fn evaluate_with_radius(&self, delta: f64, base_radius: f64) -> TofProfilePoint {
        if base_radius.is_finite() {
            return self.evaluate_supported(delta, base_radius);
        }
        let sum = self.alpha + self.beta;
        let left_fraction = self.beta / sum;
        let right_fraction = self.alpha / sum;
        let d_left_d_alpha = -self.beta / (sum * sum);
        let d_left_d_beta = self.alpha / (sum * sum);
        let mut left = TofProfilePoint::default();
        let mut right = TofProfilePoint::default();
        let mut left_alpha_shift = 0.0;
        let mut right_beta_shift = 0.0;
        for index in 0..TOF_QUADRATURE_COUNT {
            let node = self.quadrature.nodes[index];
            let weight = self.quadrature.weights[index];
            let left_delta = delta + node / self.alpha;
            if left_delta.abs() <= base_radius {
                let point = self.shape.evaluate(left_delta);
                left.value += weight * point.value;
                left.d_position += weight * point.d_delta;
                left.d_gaussian_fwhm += weight * point.d_gaussian_fwhm;
                left.d_lorentzian_fwhm += weight * point.d_lorentzian_fwhm;
                left_alpha_shift += weight * point.d_delta * (-node / self.alpha.powi(2));
            }
            let right_delta = delta - node / self.beta;
            if right_delta.abs() <= base_radius {
                let point = self.shape.evaluate(right_delta);
                right.value += weight * point.value;
                right.d_position += weight * point.d_delta;
                right.d_gaussian_fwhm += weight * point.d_gaussian_fwhm;
                right.d_lorentzian_fwhm += weight * point.d_lorentzian_fwhm;
                right_beta_shift += weight * point.d_delta * (node / self.beta.powi(2));
            }
        }
        TofProfilePoint {
            value: left_fraction * left.value + right_fraction * right.value,
            d_position: -(left_fraction * left.d_position + right_fraction * right.d_position),
            d_alpha: d_left_d_alpha * left.value + left_fraction * left_alpha_shift
                - d_left_d_alpha * right.value,
            d_beta: d_left_d_beta * left.value + right_fraction * right_beta_shift
                - d_left_d_beta * right.value,
            d_gaussian_fwhm: left_fraction * left.d_gaussian_fwhm
                + right_fraction * right.d_gaussian_fwhm,
            d_lorentzian_fwhm: left_fraction * left.d_lorentzian_fwhm
                + right_fraction * right.d_lorentzian_fwhm,
        }
    }

    #[allow(clippy::too_many_lines)]
    fn evaluate_supported(&self, delta: f64, base_radius: f64) -> TofProfilePoint {
        let sum = self.alpha + self.beta;
        let left_fraction = self.beta / sum;
        let right_fraction = self.alpha / sum;
        let d_left_d_alpha = -self.beta / (sum * sum);
        let d_left_d_beta = self.alpha / (sum * sum);
        let tail_log = self.quadrature.tail_log;
        let normalization = 1.0 - (-tail_log).exp();
        let support_multiple = base_radius / self.shape.total_fwhm;
        let mut d_radius = [0.0; SUPPORTED_VARIABLE_COUNT];
        d_radius[SUPPORTED_GAUSSIAN] = support_multiple * self.shape.d_total_fwhm_d_gaussian_fwhm;
        d_radius[SUPPORTED_LORENTZIAN] =
            support_multiple * self.shape.d_total_fwhm_d_lorentzian_fwhm;
        let (left_low, left_high) = self.supported_bounds(delta, base_radius, d_radius, true);
        let (right_low, right_high) = self.supported_bounds(delta, base_radius, d_radius, false);
        let mut left = TofProfilePoint::default();
        let mut right = TofProfilePoint::default();
        let mut left_alpha_shift = 0.0;
        let mut right_beta_shift = 0.0;

        if left_low.value < left_high.value {
            let panel_width =
                (left_high.value - left_low.value) / TOF_SUPPORT_QUADRATURE_PANELS_F64;
            let mut panel_left = left_low.value;
            for _ in 0..TOF_SUPPORT_QUADRATURE_PANELS {
                for quadrature in 0..QUADRATURE_ORDER {
                    let node = panel_left + panel_width * QUADRATURE_NODES[quadrature];
                    let weight = panel_width * QUADRATURE_WEIGHTS[quadrature] * (-node).exp()
                        / normalization;
                    let point = self.shape.evaluate(delta + node / self.alpha);
                    left.value += weight * point.value;
                    left.d_position += weight * point.d_delta;
                    left.d_gaussian_fwhm += weight * point.d_gaussian_fwhm;
                    left.d_lorentzian_fwhm += weight * point.d_lorentzian_fwhm;
                    left_alpha_shift += weight * point.d_delta * (-node / self.alpha.powi(2));
                }
                panel_left += panel_width;
            }
        }
        if right_low.value < right_high.value {
            let panel_width =
                (right_high.value - right_low.value) / TOF_SUPPORT_QUADRATURE_PANELS_F64;
            let mut panel_left = right_low.value;
            for _ in 0..TOF_SUPPORT_QUADRATURE_PANELS {
                for quadrature in 0..QUADRATURE_ORDER {
                    let node = panel_left + panel_width * QUADRATURE_NODES[quadrature];
                    let weight = panel_width * QUADRATURE_WEIGHTS[quadrature] * (-node).exp()
                        / normalization;
                    let point = self.shape.evaluate(delta - node / self.beta);
                    right.value += weight * point.value;
                    right.d_position += weight * point.d_delta;
                    right.d_gaussian_fwhm += weight * point.d_gaussian_fwhm;
                    right.d_lorentzian_fwhm += weight * point.d_lorentzian_fwhm;
                    right_beta_shift += weight * point.d_delta * (node / self.beta.powi(2));
                }
                panel_left += panel_width;
            }
        }
        let left_boundary =
            self.supported_boundary_chain(delta, left_low, left_high, true, normalization);
        let right_boundary =
            self.supported_boundary_chain(delta, right_low, right_high, false, normalization);
        TofProfilePoint {
            value: left_fraction * left.value + right_fraction * right.value,
            d_position: -(left_fraction * (left.d_position + left_boundary[SUPPORTED_DELTA])
                + right_fraction * (right.d_position + right_boundary[SUPPORTED_DELTA])),
            d_alpha: d_left_d_alpha * left.value
                + left_fraction * (left_alpha_shift + left_boundary[SUPPORTED_ALPHA])
                - d_left_d_alpha * right.value
                + right_fraction * right_boundary[SUPPORTED_ALPHA],
            d_beta: d_left_d_beta * left.value
                + left_fraction * left_boundary[SUPPORTED_BETA]
                + right_fraction * (right_beta_shift + right_boundary[SUPPORTED_BETA])
                - d_left_d_beta * right.value,
            d_gaussian_fwhm: left_fraction
                * (left.d_gaussian_fwhm + left_boundary[SUPPORTED_GAUSSIAN])
                + right_fraction * (right.d_gaussian_fwhm + right_boundary[SUPPORTED_GAUSSIAN]),
            d_lorentzian_fwhm: left_fraction
                * (left.d_lorentzian_fwhm + left_boundary[SUPPORTED_LORENTZIAN])
                + right_fraction * (right.d_lorentzian_fwhm + right_boundary[SUPPORTED_LORENTZIAN]),
        }
    }

    fn supported_bounds(
        &self,
        delta: f64,
        base_radius: f64,
        d_radius: [f64; SUPPORTED_VARIABLE_COUNT],
        left_side: bool,
    ) -> (SupportedScalar, SupportedScalar) {
        let tail_log = self.quadrature.tail_log;
        let rate = if left_side { self.alpha } else { self.beta };
        let (low_sign, high_sign) = if left_side {
            (-base_radius - delta, base_radius - delta)
        } else {
            (delta - base_radius, delta + base_radius)
        };
        let mut low_derivative = [0.0; SUPPORTED_VARIABLE_COUNT];
        let mut high_derivative = [0.0; SUPPORTED_VARIABLE_COUNT];
        if left_side {
            low_derivative[SUPPORTED_DELTA] = -rate;
            high_derivative[SUPPORTED_DELTA] = -rate;
            low_derivative[SUPPORTED_ALPHA] = low_sign;
            high_derivative[SUPPORTED_ALPHA] = high_sign;
            for parameter in [SUPPORTED_GAUSSIAN, SUPPORTED_LORENTZIAN] {
                low_derivative[parameter] = -rate * d_radius[parameter];
                high_derivative[parameter] = rate * d_radius[parameter];
            }
        } else {
            low_derivative[SUPPORTED_DELTA] = rate;
            high_derivative[SUPPORTED_DELTA] = rate;
            low_derivative[SUPPORTED_BETA] = low_sign;
            high_derivative[SUPPORTED_BETA] = high_sign;
            for parameter in [SUPPORTED_GAUSSIAN, SUPPORTED_LORENTZIAN] {
                low_derivative[parameter] = -rate * d_radius[parameter];
                high_derivative[parameter] = rate * d_radius[parameter];
            }
        }
        (
            SupportedScalar {
                value: rate * low_sign,
                derivative: low_derivative,
            }
            .clamped(0.0, tail_log),
            SupportedScalar {
                value: rate * high_sign,
                derivative: high_derivative,
            }
            .clamped(0.0, tail_log),
        )
    }

    fn supported_boundary_chain(
        &self,
        delta: f64,
        low: SupportedScalar,
        high: SupportedScalar,
        left_side: bool,
        normalization: f64,
    ) -> [f64; SUPPORTED_VARIABLE_COUNT] {
        let rate = if left_side { self.alpha } else { self.beta };
        let direction = if left_side { 1.0 } else { -1.0 };
        let integrand = |node: f64| {
            (-node).exp() / normalization
                * self.shape.evaluate(delta + direction * node / rate).value
        };
        let low_value = integrand(low.value);
        let high_value = integrand(high.value);
        let mut derivative = [0.0; SUPPORTED_VARIABLE_COUNT];
        for (parameter, value) in derivative.iter_mut().enumerate() {
            *value =
                high_value * high.derivative[parameter] - low_value * low.derivative[parameter];
        }
        derivative
    }

    fn support_range(&self, position: f64, base_radius: f64) -> SupportRange {
        SupportRange {
            left: position - base_radius - self.quadrature.tail_log / self.alpha,
            right: position + base_radius + self.quadrature.tail_log / self.beta,
        }
    }
}

/// TOF domain or accumulation error.
#[derive(Clone, Debug, PartialEq)]
pub enum TofError {
    /// Bank scattering angle is not finite or outside `0 < 2theta < 180°`.
    InvalidBankTwoTheta,
    /// At least one instrument coefficient is non-finite.
    NonFiniteInstrumentParameter,
    /// The linear calibration coefficient is not positive.
    NonPositiveDifc,
    /// A reflection d-spacing is not positive and finite.
    InvalidDSpacing,
    /// Calibration produced a non-finite position.
    InvalidPosition,
    /// The derived leading-edge exponential rate is invalid.
    NonPositiveAlpha,
    /// The derived trailing-edge exponential rate is invalid.
    NonPositiveBeta,
    /// The derived Gaussian variance is invalid.
    NonPositiveGaussianVariance,
    /// The derived Lorentzian width is invalid.
    NegativeLorentzianFwhm,
    /// The requested exponential tail cutoff is invalid.
    InvalidTailLog,
    /// Quadrature normalization failed.
    InvalidQuadrature,
    /// Component widths could not be transformed into a TCH shape.
    InvalidTch {
        /// Underlying component-width transformation failure.
        reason: TchError,
    },
    /// Reflection d-spacing and intensity arrays have different lengths.
    LengthMismatch,
    /// One reflection intensity is non-finite.
    NonFiniteIntensity {
        /// Index of the invalid reflection.
        reflection: usize,
    },
    /// Allocation-size arithmetic overflowed or allocation failed.
    AllocationOverflow,
    /// Grid or support validation failed.
    Profile {
        /// Underlying grid or support failure.
        reason: ProfileError,
    },
}

impl Display for TofError {
    fn fmt(&self, formatter: &mut Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::InvalidBankTwoTheta => formatter
                .write_str("TOF bank two_theta_deg must be finite and strictly within (0, 180)"),
            Self::NonFiniteInstrumentParameter => {
                write!(formatter, "TOF coefficients must be finite")
            }
            Self::NonPositiveDifc => write!(formatter, "difC must be positive"),
            Self::InvalidDSpacing => write!(formatter, "d-spacing must be positive and finite"),
            Self::InvalidPosition => write!(formatter, "derived TOF position must be finite"),
            Self::NonPositiveAlpha => write!(formatter, "TOF alpha must be positive and finite"),
            Self::NonPositiveBeta => write!(formatter, "TOF beta must be positive and finite"),
            Self::NonPositiveGaussianVariance => write!(
                formatter,
                "derived TOF Gaussian variance must be positive and finite"
            ),
            Self::NegativeLorentzianFwhm => write!(
                formatter,
                "derived TOF Lorentzian FWHM must be non-negative and finite"
            ),
            Self::InvalidTailLog => write!(formatter, "tail_log must be positive and finite"),
            Self::InvalidQuadrature => write!(formatter, "TOF quadrature normalization is invalid"),
            Self::InvalidTch { reason } => write!(formatter, "invalid TOF TCH widths: {reason}"),
            Self::LengthMismatch => {
                write!(formatter, "TOF reflection arrays must have equal length")
            }
            Self::NonFiniteIntensity { reflection } => write!(
                formatter,
                "reflection {reflection} intensity must be finite"
            ),
            Self::AllocationOverflow => write!(formatter, "TOF allocation size overflow"),
            Self::Profile { reason } => Display::fmt(reason, formatter),
        }
    }
}

impl Error for TofError {}

impl From<ProfileError> for TofError {
    fn from(reason: ProfileError) -> Self {
        Self::Profile { reason }
    }
}

struct PreparedReflection {
    parameters: TofProfileParameters,
    profile: TofProfile,
}

struct TofReflectionBlock {
    start: usize,
    y: Vec<f64>,
    local: Vec<f64>,
    global: Vec<f64>,
}

/// Fused finite-support TOF accumulation with local intensity/d-spacing rows.
///
/// # Errors
///
/// Returns [`TofError`] for invalid grids, reflection arrays, coefficients,
/// derived profiles, supports, or checked allocation failures.
#[allow(clippy::too_many_lines)]
pub fn accumulate_tof_batch(
    grid: GridView<'_>,
    d_spacings: &[f64],
    intensities: &[f64],
    instrument: TofInstrument,
    support_fwhm: f64,
    tail_log: f64,
) -> Result<Accumulation, TofError> {
    accumulate_tof_batch_with_context(
        grid,
        d_spacings,
        intensities,
        instrument,
        support_fwhm,
        tail_log,
        &ExecutionContext::serial(),
    )
}

/// Accumulate TOF reflections with a bounded execution context.
///
/// # Errors
///
/// Returns [`TofError`] for invalid inputs, profiles, supports, or allocations.
#[allow(clippy::too_many_lines)]
pub fn accumulate_tof_batch_with_context(
    grid: GridView<'_>,
    d_spacings: &[f64],
    intensities: &[f64],
    instrument: TofInstrument,
    support_fwhm: f64,
    tail_log: f64,
    execution: &ExecutionContext,
) -> Result<Accumulation, TofError> {
    if d_spacings.len() != intensities.len() {
        return Err(TofError::LengthMismatch);
    }
    if !support_fwhm.is_finite() || support_fwhm <= 0.0 {
        return Err(TofError::Profile {
            reason: ProfileError::InvalidSupport,
        });
    }
    instrument.validate()?;
    let x = grid.as_slice();
    let count = d_spacings.len();
    let mut prepared = Vec::new();
    let mut starts = Vec::new();
    let mut offsets = Vec::new();
    prepared
        .try_reserve_exact(count)
        .map_err(|_| TofError::AllocationOverflow)?;
    starts
        .try_reserve_exact(count)
        .map_err(|_| TofError::AllocationOverflow)?;
    let offset_count = count.checked_add(1).ok_or(TofError::AllocationOverflow)?;
    offsets
        .try_reserve_exact(offset_count)
        .map_err(|_| TofError::AllocationOverflow)?;
    offsets.push(0usize);
    let mut shared_quadrature = None;
    for reflection in 0..count {
        if !intensities[reflection].is_finite() {
            return Err(TofError::NonFiniteIntensity { reflection });
        }
        let parameters = TofProfileParameters::from_instrument(d_spacings[reflection], instrument)?;
        TofProfile::validate_rates(parameters.alpha_per_us, parameters.beta_per_us)?;
        let quadrature = if let Some(quadrature) = &shared_quadrature {
            Arc::clone(quadrature)
        } else {
            let quadrature = Arc::new(TofQuadrature::new(tail_log)?);
            shared_quadrature = Some(Arc::clone(&quadrature));
            quadrature
        };
        let profile = TofProfile::from_validated_rates(
            parameters.alpha_per_us,
            parameters.beta_per_us,
            TchWidths {
                gaussian_fwhm: parameters.gaussian_fwhm_us,
                lorentzian_fwhm: parameters.lorentzian_fwhm_us,
            },
            quadrature,
        )?;
        let base_radius = support_fwhm * parameters.tch.total_fwhm;
        let range = profile.support_range(parameters.position_us, base_radius);
        let lower = x.partition_point(|value| *value < range.left);
        let upper = x.partition_point(|value| *value <= range.right);
        offsets.push(
            offsets[reflection]
                .checked_add(upper - lower)
                .ok_or(TofError::AllocationOverflow)?,
        );
        starts.push(lower);
        prepared.push(PreparedReflection {
            parameters,
            profile,
        });
    }
    let active = offsets.last().copied().unwrap_or(0);
    let mut y = zeroed_f64_vec(x.len())?;
    let mut local = zeroed_f64_vec(
        active
            .checked_mul(LOCAL_PARAMETER_COUNT)
            .ok_or(TofError::AllocationOverflow)?,
    )?;
    let mut global = zeroed_f64_vec(
        TOF_GLOBAL_PARAMETER_COUNT
            .checked_mul(x.len())
            .ok_or(TofError::AllocationOverflow)?,
    )?;
    if execution.threads() == 1 || count < 16 {
        for reflection in 0..count {
            let item = &prepared[reflection];
            let intensity = intensities[reflection];
            let base_radius = support_fwhm * item.parameters.tch.total_fwhm;
            let begin = offsets[reflection];
            let end = offsets[reflection + 1];
            for active_index in begin..end {
                let sample = starts[reflection] + active_index - begin;
                let point = item
                    .profile
                    .evaluate_with_radius(x[sample] - item.parameters.position_us, base_radius);
                y[sample] += intensity * point.value;
                local[active_index * LOCAL_PARAMETER_COUNT] = point.value;
                local[active_index * LOCAL_PARAMETER_COUNT + 1] = intensity
                    * (point.d_position * item.parameters.d_position_d_d
                        + point.d_alpha * item.parameters.d_alpha_d_d
                        + point.d_beta * item.parameters.d_beta_d_d
                        + point.d_gaussian_fwhm * item.parameters.d_gaussian_fwhm_d_d
                        + point.d_lorentzian_fwhm * item.parameters.d_lorentzian_fwhm_d_d);
                for parameter in 0..TOF_GLOBAL_PARAMETER_COUNT {
                    let derivative = point.d_position
                        * item.parameters.d_position_d_instrument[parameter]
                        + point.d_alpha * item.parameters.d_alpha_d_instrument[parameter]
                        + point.d_beta * item.parameters.d_beta_d_instrument[parameter]
                        + point.d_gaussian_fwhm
                            * item.parameters.d_gaussian_fwhm_d_instrument[parameter]
                        + point.d_lorentzian_fwhm
                            * item.parameters.d_lorentzian_fwhm_d_instrument[parameter];
                    global[parameter * x.len() + sample] += intensity * derivative;
                }
            }
        }
    } else {
        let blocks = execution.map_ordered(count, 16, |reflection| {
            let item = &prepared[reflection];
            let intensity = intensities[reflection];
            let base_radius = support_fwhm * item.parameters.tch.total_fwhm;
            let begin = offsets[reflection];
            let end = offsets[reflection + 1];
            let support_count = end - begin;
            let mut block = TofReflectionBlock {
                start: starts[reflection],
                y: vec![0.0; support_count],
                local: vec![0.0; support_count * LOCAL_PARAMETER_COUNT],
                global: vec![0.0; support_count * TOF_GLOBAL_PARAMETER_COUNT],
            };
            for support_index in 0..support_count {
                let sample = block.start + support_index;
                let point = item
                    .profile
                    .evaluate_with_radius(x[sample] - item.parameters.position_us, base_radius);
                block.y[support_index] = intensity * point.value;
                block.local[support_index * LOCAL_PARAMETER_COUNT] = point.value;
                block.local[support_index * LOCAL_PARAMETER_COUNT + 1] = intensity
                    * (point.d_position * item.parameters.d_position_d_d
                        + point.d_alpha * item.parameters.d_alpha_d_d
                        + point.d_beta * item.parameters.d_beta_d_d
                        + point.d_gaussian_fwhm * item.parameters.d_gaussian_fwhm_d_d
                        + point.d_lorentzian_fwhm * item.parameters.d_lorentzian_fwhm_d_d);
                for parameter in 0..TOF_GLOBAL_PARAMETER_COUNT {
                    let derivative = point.d_position
                        * item.parameters.d_position_d_instrument[parameter]
                        + point.d_alpha * item.parameters.d_alpha_d_instrument[parameter]
                        + point.d_beta * item.parameters.d_beta_d_instrument[parameter]
                        + point.d_gaussian_fwhm
                            * item.parameters.d_gaussian_fwhm_d_instrument[parameter]
                        + point.d_lorentzian_fwhm
                            * item.parameters.d_lorentzian_fwhm_d_instrument[parameter];
                    block.global[parameter * support_count + support_index] =
                        intensity * derivative;
                }
            }
            block
        });
        for (reflection, block) in blocks.into_iter().enumerate() {
            let begin = offsets[reflection];
            let support_count = block.y.len();
            let local_begin = begin * LOCAL_PARAMETER_COUNT;
            let local_end = local_begin + block.local.len();
            local[local_begin..local_end].copy_from_slice(&block.local);
            for support_index in 0..support_count {
                let sample = block.start + support_index;
                y[sample] += block.y[support_index];
                for parameter in 0..TOF_GLOBAL_PARAMETER_COUNT {
                    global[parameter * x.len() + sample] +=
                        block.global[parameter * support_count + support_index];
                }
            }
        }
    }
    Ok(Accumulation {
        y,
        derivatives: PatternDerivatives {
            local: SupportJacobian {
                starts,
                offsets,
                values: local,
                parameter_count: LOCAL_PARAMETER_COUNT,
            },
            global: Some(DenseJacobian {
                values: global,
                parameter_count: TOF_GLOBAL_PARAMETER_COUNT,
                sample_count: x.len(),
            }),
        },
        sample_count: x.len(),
    })
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn incident_spectrum_matches_closed_form_and_centered_difference() {
        let spectrum = TofIncidentSpectrum::new(
            500.0,
            10_000.0,
            [
                12.0, 40_000.0, 3.0, 2.0, -0.5, 0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            ],
        )
        .expect("spectrum");
        let tof_us: f64 = 2_500.0;
        let t = tof_us / 1_000.0;
        let x = 2.0 / t - 1.0;
        let expected = 12.0 + 40_000.0 / t.powi(5) * (-3.0 / t.powi(2)).exp() + 2.0 * x
            - 0.5 * (2.0 * x * x - 1.0)
            + 0.25 * (4.0 * x.powi(3) - 3.0 * x);
        let actual = spectrum.evaluate(tof_us).expect("evaluation");
        assert!((actual.value - expected).abs() < 1.0e-12 * expected.abs());

        let step_us = 1.0e-3;
        let plus = spectrum.evaluate(tof_us + step_us).unwrap().value;
        let minus = spectrum.evaluate(tof_us - step_us).unwrap().value;
        let finite = (plus - minus) / (2.0 * step_us);
        assert!((actual.d_value_d_tof_us - finite).abs() < 1.0e-9);
    }

    #[test]
    fn incident_spectrum_enforces_constructor_range_and_positive_evaluation() {
        assert_eq!(
            TofIncidentSpectrum::new(1_000.0, 1_000.0, [1.0; 12]),
            Err(TofIncidentSpectrumError::InvalidRange)
        );
        assert_eq!(
            TofIncidentSpectrum::new(1_000.0, 2_000.0, [f64::NAN; 12]),
            Err(TofIncidentSpectrumError::NonFiniteCoefficient)
        );
        let positive = TofIncidentSpectrum::new(
            1_000.0,
            2_000.0,
            [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        )
        .unwrap();
        assert!(positive.evaluate(1_000.0).is_ok());
        assert!(positive.evaluate(2_000.0).is_ok());
        assert_eq!(
            positive.evaluate(999.0),
            Err(TofIncidentSpectrumError::TofOutsideRange)
        );
        let zero = TofIncidentSpectrum::new(1_000.0, 2_000.0, [0.0; 12]).unwrap();
        assert_eq!(
            zero.evaluate(1_500.0),
            Err(TofIncidentSpectrumError::NonPositiveIntensity)
        );
    }

    #[test]
    fn bank_geometry_requires_a_strict_physical_scattering_angle() {
        let geometry = TofBankGeometry {
            two_theta_deg: 88.05,
        };
        geometry.validate().expect("valid bank geometry");
        assert!((geometry.theta_radians().unwrap() - 44.025_f64.to_radians()).abs() < 1.0e-15);
        for two_theta_deg in [f64::NAN, 0.0, 180.0, f64::INFINITY] {
            assert_eq!(
                TofBankGeometry { two_theta_deg }.validate(),
                Err(TofError::InvalidBankTwoTheta)
            );
        }
    }

    fn instrument() -> TofInstrument {
        TofInstrument {
            zero_us: -0.773_346_536_757,
            difc_us_per_angstrom: 5_084.827_630_65,
            difa_us_per_angstrom2: -2.630_417_748_6,
            difb_us_angstrom: 0.0,
            alpha_coefficient: 5.0,
            beta0_per_us: 0.033_276_398_966_5,
            beta1_angstrom4_per_us: 0.000_964_057_827_372,
            betaq_angstrom2_per_us: 0.0,
            sigma0_us2: 0.0,
            sigma1_us2_per_angstrom2: 15.140_286_726_8,
            sigma2_us2_per_angstrom4: 0.0,
            sigmaq_us2_per_angstrom: 0.0,
            x_us_per_angstrom: 0.0,
            y_us_per_angstrom2: 0.0,
            z_us: 0.0,
        }
    }

    #[test]
    fn parameter_chains_match_centered_differences() {
        let d = 1.7;
        let step = 1.0e-6;
        let actual = TofProfileParameters::from_instrument(d, instrument()).expect("parameters");
        let plus = TofProfileParameters::from_instrument(d + step, instrument()).expect("plus");
        let minus = TofProfileParameters::from_instrument(d - step, instrument()).expect("minus");
        let finite = |high: f64, low: f64| (high - low) / (2.0 * step);
        assert!((actual.d_position_d_d - finite(plus.position_us, minus.position_us)).abs() < 1e-6);
        assert!((actual.d_alpha_d_d - finite(plus.alpha_per_us, minus.alpha_per_us)).abs() < 1e-9);
        assert!((actual.d_beta_d_d - finite(plus.beta_per_us, minus.beta_per_us)).abs() < 1e-9);
        assert!(
            (actual.d_gaussian_fwhm_d_d - finite(plus.gaussian_fwhm_us, minus.gaussian_fwhm_us))
                .abs()
                < 1e-8
        );
    }

    #[test]
    fn selectable_instrument_parameters_follow_dense_row_order() {
        let original = instrument();
        let values = original.values();
        for parameter in TofInstrumentParameter::ALL {
            assert_eq!(
                parameter.name(),
                TOF_GLOBAL_PARAMETER_NAMES[parameter.index()]
            );
            let replacement = values[parameter.index()] + 1.0e-6;
            let updated = original
                .with_parameter(parameter, replacement)
                .expect("valid replacement");
            for (index, value) in updated.values().iter().copied().enumerate() {
                let expected = if index == parameter.index() {
                    replacement
                } else {
                    values[index]
                };
                assert_eq!(value.to_bits(), expected.to_bits());
            }
        }
        let mut invalid = values;
        invalid[TofInstrumentParameter::Difc.index()] = 0.0;
        assert_eq!(
            TofInstrument::from_values(invalid),
            Err(TofError::NonPositiveDifc)
        );
        invalid = values;
        invalid[TofInstrumentParameter::Zero.index()] = f64::NAN;
        assert_eq!(
            TofInstrument::from_values(invalid),
            Err(TofError::NonFiniteInstrumentParameter)
        );
    }

    #[test]
    fn profiles_can_share_quadrature_storage() {
        let quadrature = Arc::new(TofQuadrature::new(20.0).expect("quadrature"));
        let widths = TchWidths {
            gaussian_fwhm: 22.0,
            lorentzian_fwhm: 4.0,
        };
        let first = TofProfile::from_validated_rates(0.08, 0.03, widths, Arc::clone(&quadrature))
            .expect("first profile");
        let second = TofProfile::from_validated_rates(0.09, 0.04, widths, Arc::clone(&quadrature))
            .expect("second profile");

        assert!(Arc::ptr_eq(&first.quadrature, &second.quadrature));
        assert!(std::mem::size_of::<TofProfile>() < 128);
    }

    #[test]
    fn supported_bound_derivative_is_zero_at_exact_clamp() {
        let derivative = [1.0, 2.0, 3.0, 4.0, 5.0];
        for value in [0.0, 20.0] {
            let bounded = SupportedScalar { value, derivative }.clamped(0.0, 20.0);
            assert_eq!(bounded.value.to_bits(), value.to_bits());
            assert!(bounded.derivative.iter().all(|value| value.to_bits() == 0));
        }
        let interior = SupportedScalar {
            value: 10.0,
            derivative,
        }
        .clamped(0.0, 20.0);
        assert!(
            interior
                .derivative
                .iter()
                .zip(derivative)
                .all(|(actual, expected)| actual.to_bits() == expected.to_bits())
        );
    }

    #[test]
    fn direct_profile_is_numerically_unit_area() {
        let profile = TofProfile::new(
            0.08,
            0.03,
            TchWidths {
                gaussian_fwhm: 22.0,
                lorentzian_fwhm: 4.0,
            },
            20.0,
        )
        .expect("profile");
        let step = 0.5;
        let radius = 6_000.0;
        let sample_count = 24_000_u32;
        let mut area = 0.0;
        let mut previous = profile.evaluate(-radius).value;
        for index in 1..=sample_count {
            let x = -radius + f64::from(index) * step;
            let value = profile.evaluate(x).value;
            area += 0.5 * step * (previous + value);
            previous = value;
        }
        // The remaining discrepancy is the analytically infinite Lorentzian
        // tail outside this deliberately finite integration interval.
        assert!((area - 1.0).abs() < 3.0e-4, "integrated area={area:.12}");
    }

    #[test]
    fn accumulation_blocks_are_bitwise_identical_across_worker_counts() {
        let x = (0..=8_000)
            .map(|index| 1_000.0 + 2.0 * f64::from(index))
            .collect::<Vec<_>>();
        let d_spacings = (0..36)
            .map(|index| 0.5 + 0.065 * f64::from(index))
            .collect::<Vec<_>>();
        let intensities = (0..36)
            .map(|index| 3.0 + 0.3 * f64::from(index))
            .collect::<Vec<_>>();
        let grid = GridView::new(&x).expect("grid");
        let serial = ExecutionContext::serial();
        let expected = accumulate_tof_batch_with_context(
            grid,
            &d_spacings,
            &intensities,
            instrument(),
            20.0,
            20.0,
            &serial,
        )
        .expect("serial TOF");
        for threads in [2, 3] {
            let context = ExecutionContext::new(threads).expect("parallel context");
            assert_eq!(
                accumulate_tof_batch_with_context(
                    grid,
                    &d_spacings,
                    &intensities,
                    instrument(),
                    20.0,
                    20.0,
                    &context,
                )
                .expect("parallel TOF"),
                expected
            );
        }
    }

    #[test]
    fn direct_profile_derivatives_match_centered_differences() {
        let alpha = 0.08;
        let beta = 0.03;
        let gaussian = 22.0;
        let lorentzian = 4.0;
        let delta = 5.0;
        let tail = 20.0;
        let point = TofProfile::new(
            alpha,
            beta,
            TchWidths {
                gaussian_fwhm: gaussian,
                lorentzian_fwhm: lorentzian,
            },
            tail,
        )
        .expect("profile")
        .evaluate(delta);
        let step = 1.0e-6;
        let value = |a, b, g, l, x| {
            TofProfile::new(
                a,
                b,
                TchWidths {
                    gaussian_fwhm: g,
                    lorentzian_fwhm: l,
                },
                tail,
            )
            .expect("profile")
            .evaluate(x)
            .value
        };
        let fd = |plus, minus| (plus - minus) / (2.0 * step);
        assert!(
            (point.d_position
                - fd(
                    value(alpha, beta, gaussian, lorentzian, delta - step),
                    value(alpha, beta, gaussian, lorentzian, delta + step)
                ))
            .abs()
                < 1e-8
        );
        assert!(
            (point.d_alpha
                - fd(
                    value(alpha + step, beta, gaussian, lorentzian, delta),
                    value(alpha - step, beta, gaussian, lorentzian, delta)
                ))
            .abs()
                < 1e-7
        );
        assert!(
            (point.d_beta
                - fd(
                    value(alpha, beta + step, gaussian, lorentzian, delta),
                    value(alpha, beta - step, gaussian, lorentzian, delta)
                ))
            .abs()
                < 1e-7
        );
        assert!(
            (point.d_gaussian_fwhm
                - fd(
                    value(alpha, beta, gaussian + step, lorentzian, delta),
                    value(alpha, beta, gaussian - step, lorentzian, delta)
                ))
            .abs()
                < 1e-8
        );
        assert!(
            (point.d_lorentzian_fwhm
                - fd(
                    value(alpha, beta, gaussian, lorentzian + step, delta),
                    value(alpha, beta, gaussian, lorentzian - step, delta)
                ))
            .abs()
                < 1e-8
        );
    }

    #[test]
    fn supported_profile_derivatives_include_moving_integration_bounds() {
        let alpha = 0.08;
        let beta = 0.03;
        let gaussian = 22.0;
        let lorentzian = 4.0;
        let delta = 35.0;
        let tail = 20.0;
        let support_multiple = 1.25;
        let profile = TofProfile::new(
            alpha,
            beta,
            TchWidths {
                gaussian_fwhm: gaussian,
                lorentzian_fwhm: lorentzian,
            },
            tail,
        )
        .expect("profile");
        let point =
            profile.evaluate_with_radius(delta, support_multiple * profile.shape.total_fwhm);
        let step = 1.0e-6;
        let value = |a, b, g, l, x| {
            let profile = TofProfile::new(
                a,
                b,
                TchWidths {
                    gaussian_fwhm: g,
                    lorentzian_fwhm: l,
                },
                tail,
            )
            .expect("profile");
            profile
                .evaluate_with_radius(x, support_multiple * profile.shape.total_fwhm)
                .value
        };
        let fd = |plus, minus| (plus - minus) / (2.0 * step);
        assert!(
            (point.d_position
                - fd(
                    value(alpha, beta, gaussian, lorentzian, delta - step),
                    value(alpha, beta, gaussian, lorentzian, delta + step),
                ))
            .abs()
                < 2.0e-9
        );
        assert!(
            (point.d_alpha
                - fd(
                    value(alpha + step, beta, gaussian, lorentzian, delta),
                    value(alpha - step, beta, gaussian, lorentzian, delta),
                ))
            .abs()
                < 2.0e-8
        );
        assert!(
            (point.d_beta
                - fd(
                    value(alpha, beta + step, gaussian, lorentzian, delta),
                    value(alpha, beta - step, gaussian, lorentzian, delta),
                ))
            .abs()
                < 2.0e-8
        );
        assert!(
            (point.d_gaussian_fwhm
                - fd(
                    value(alpha, beta, gaussian + step, lorentzian, delta),
                    value(alpha, beta, gaussian - step, lorentzian, delta),
                ))
            .abs()
                < 2.0e-9
        );
        assert!(
            (point.d_lorentzian_fwhm
                - fd(
                    value(alpha, beta, gaussian, lorentzian + step, delta),
                    value(alpha, beta, gaussian, lorentzian - step, delta),
                ))
            .abs()
                < 2.0e-9
        );
    }
}