nereids_physics/
flight_time_grid.rs

1//! The flight times whose neutrons can be counted in a set of detector time
2//! bins, and the counts a beam predicts there, on a uniform flight-time grid.
3//!
4//! A neutron of flight time `u` has energy `E = (TOF_FACTOR·L/u)²` and
5//! arrives at `t0 + u + delay`, the delay drawn from the Ikeda–Carpenter
6//! pulse at `E`.  The grid spans every flight time whose neutrons reach the
7//! bins with more than [`NEGLIGIBLE_ARRIVAL_PROBABILITY`](crate::ikeda_carpenter::NEGLIGIBLE_ARRIVAL_PROBABILITY) chance, and refuses
8//! a window that neutrons from outside the pulse's synthesis grid can reach.
9
10use std::fmt;
11use std::sync::Arc;
12
13use crate::counts_response::{CountsResponseError, DetectorBinResponseMatrix};
14use crate::ikeda_carpenter::{EnergyLaw, IkedaCarpenter};
15use crate::resolution::{ResolutionFunction, ResolutionParseError, TOF_FACTOR};
16
17/// Largest number of grid points a window may need.
18pub const MAX_POINTS: usize = 100_000;
19
20/// Why a window has no flight-time grid.
21#[derive(Debug)]
22pub enum FlightTimeGridError {
23    /// Fewer than two edges, or edges not finite and strictly ascending.
24    InvalidTimeEdges,
25    /// The timing offset is not finite.
26    InvalidTimingOffset(f64),
27    /// Not supported: the pulse's `parameter` law lengthens the pulse as
28    /// energy rises (α falling, β falling with a storage term, or R rising).
29    /// The range relies on
30    /// arrival times growing with flight time, which such a pulse does not
31    /// guarantee.
32    LengthensWithEnergy { parameter: &'static str },
33    /// The pulse cannot be evaluated.
34    Pulse(ResolutionParseError),
35    /// Neutrons with energies outside the pulse's synthesis grid,
36    /// `low_ev` to `high_ev`, can reach the window.
37    OutsideCalibration { low_ev: f64, high_ev: f64 },
38    /// A grid at `step_us` would need more than [`MAX_POINTS`] points.
39    TooManyPoints { step_us: f64 },
40    /// The bin probabilities could not be built.
41    Response(CountsResponseError),
42    /// `found` values were given for a grid of `expected` points.
43    ValuesLength { expected: usize, found: usize },
44}
45
46impl fmt::Display for FlightTimeGridError {
47    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
48        match self {
49            Self::InvalidTimeEdges => {
50                write!(
51                    f,
52                    "time edges must be at least two finite, strictly ascending values"
53                )
54            }
55            Self::InvalidTimingOffset(t0) => write!(f, "t0 = {t0} µs must be finite"),
56            Self::LengthensWithEnergy { parameter } => write!(
57                f,
58                "unsupported pulse: its {parameter} law lengthens the pulse as energy rises; \
59                 the flight-time grid needs α and β that do not fall and R that does not \
60                 rise with energy"
61            ),
62            Self::Pulse(e) => write!(f, "pulse: {e}"),
63            Self::OutsideCalibration { low_ev, high_ev } => write!(
64                f,
65                "neutrons from outside the pulse's synthesis grid, {low_ev} to \
66                 {high_ev} eV, can reach the window"
67            ),
68            Self::TooManyPoints { step_us } => write!(
69                f,
70                "a step of {step_us} µs needs more than {MAX_POINTS} grid points"
71            ),
72            Self::Response(e) => write!(f, "bin probabilities: {e}"),
73            Self::ValuesLength { expected, found } => write!(
74                f,
75                "{found} values were given for a grid of {expected} points"
76            ),
77        }
78    }
79}
80
81impl std::error::Error for FlightTimeGridError {}
82
83impl From<ResolutionParseError> for FlightTimeGridError {
84    fn from(e: ResolutionParseError) -> Self {
85        Self::Pulse(e)
86    }
87}
88
89/// Uniform flight-time grid over a window, with the probability that each
90/// grid point's neutrons are counted in each bin.
91#[derive(Debug, Clone)]
92pub struct FlightTimeGrid {
93    time_edges_us: Vec<f64>,
94    t0_us: f64,
95    pulse: Arc<IkedaCarpenter>,
96    range_us: (f64, f64),
97    step_us: f64,
98    flight_times_us: Vec<f64>,
99    response: DetectorBinResponseMatrix,
100}
101
102impl FlightTimeGrid {
103    /// The grid for bins `time_edges_us` (µs) with timing offset `t0_us`,
104    /// the flight path taken from `pulse`.  The step, half the pulse's
105    /// shorter rise at the two ends of the range, is a starting step: the
106    /// counts are converged once [`Self::halved`] no longer changes them.
107    ///
108    /// # Errors
109    /// See [`FlightTimeGridError`].
110    pub fn new(
111        time_edges_us: &[f64],
112        t0_us: f64,
113        pulse: &Arc<IkedaCarpenter>,
114    ) -> Result<Self, FlightTimeGridError> {
115        if time_edges_us.len() < 2
116            || !time_edges_us.iter().all(|t| t.is_finite())
117            || time_edges_us.windows(2).any(|w| w[0] >= w[1])
118        {
119            return Err(FlightTimeGridError::InvalidTimeEdges);
120        }
121        if !t0_us.is_finite() {
122            return Err(FlightTimeGridError::InvalidTimingOffset(t0_us));
123        }
124        let references = pulse.ref_energies();
125        let (e_min, e_max) = (references[0], references[references.len() - 1]);
126        let params = pulse.params();
127        let change = |law: &EnergyLaw| law.eval(e_max) - law.eval(e_min);
128        for (parameter, lengthens) in [
129            ("α", change(&params.alpha) < 0.0),
130            (
131                "β",
132                change(&params.beta) < 0.0 && params.r.eval(e_min) > 0.0,
133            ),
134            ("R", change(&params.r) > 0.0),
135        ] {
136            if lengthens {
137                return Err(FlightTimeGridError::LengthensWithEnergy { parameter });
138            }
139        }
140        let clock = TOF_FACTOR * pulse.flight_path_m();
141        let energy = |u: f64| (clock / u).powi(2);
142        let first_edge = time_edges_us[0] - t0_us;
143        let last_edge = time_edges_us[time_edges_us.len() - 1] - t0_us;
144        let latest = |u: f64| {
145            pulse
146                .delays_us(energy(u))
147                .map(|(_, last)| u + last - first_edge)
148        };
149        let earliest = |u: f64| {
150            pulse
151                .delays_us(energy(u))
152                .map(|(first, _)| u + first - last_edge)
153        };
154        let (fastest, slowest) = (clock / e_max.sqrt(), clock / e_min.sqrt());
155        let (Some(u_lo), Some(u_hi)) = (
156            crossing(&latest, fastest, slowest)?,
157            crossing(&earliest, fastest, slowest)?,
158        ) else {
159            return Err(FlightTimeGridError::OutsideCalibration {
160                low_ev: e_min,
161                high_ev: e_max,
162            });
163        };
164
165        let rise = pulse
166            .rise_us(energy(u_lo))?
167            .min(pulse.rise_us(energy(u_hi))?);
168        let intervals = ((u_hi - u_lo) / (0.5 * rise)).ceil();
169        Self::build(
170            time_edges_us,
171            t0_us,
172            pulse,
173            (u_lo, u_hi),
174            intervals as usize,
175        )
176    }
177
178    fn build(
179        time_edges_us: &[f64],
180        t0_us: f64,
181        pulse: &Arc<IkedaCarpenter>,
182        range_us: (f64, f64),
183        intervals: usize,
184    ) -> Result<Self, FlightTimeGridError> {
185        let (u_lo, u_hi) = range_us;
186        let step_us = (u_hi - u_lo) / intervals as f64;
187        if intervals >= MAX_POINTS {
188            return Err(FlightTimeGridError::TooManyPoints { step_us });
189        }
190        let flight_times_us: Vec<f64> =
191            (0..=intervals).map(|j| u_lo + step_us * j as f64).collect();
192        let clock = TOF_FACTOR * pulse.flight_path_m();
193        let energies: Vec<f64> = flight_times_us
194            .iter()
195            .map(|u| (clock / u).powi(2))
196            .collect();
197        let response = DetectorBinResponseMatrix::new(
198            &energies,
199            time_edges_us,
200            t0_us,
201            &ResolutionFunction::IkedaCarpenter(Arc::clone(pulse)),
202        )
203        .map_err(FlightTimeGridError::Response)?;
204        Ok(Self {
205            time_edges_us: time_edges_us.to_vec(),
206            t0_us,
207            pulse: Arc::clone(pulse),
208            range_us,
209            step_us,
210            flight_times_us,
211            response,
212        })
213    }
214
215    /// The same window at half the step.
216    ///
217    /// # Errors
218    /// As [`Self::new`], when the halved grid exceeds [`MAX_POINTS`].
219    pub fn halved(&self) -> Result<Self, FlightTimeGridError> {
220        Self::build(
221            &self.time_edges_us,
222            self.t0_us,
223            &self.pulse,
224            self.range_us,
225            2 * (self.flight_times_us.len() - 1),
226        )
227    }
228
229    /// `(u_lo, u_hi)`, the flight-time range in µs.
230    #[must_use]
231    pub fn range_us(&self) -> (f64, f64) {
232        self.range_us
233    }
234
235    /// The grid step in µs.
236    #[must_use]
237    pub fn step_us(&self) -> f64 {
238        self.step_us
239    }
240
241    /// The grid's flight times in µs, ascending.
242    #[must_use]
243    pub fn flight_times_us(&self) -> &[f64] {
244        &self.flight_times_us
245    }
246
247    /// The energy, in eV, of each grid point, descending.
248    #[must_use]
249    pub fn energies_ev(&self) -> Vec<f64> {
250        let clock = TOF_FACTOR * self.pulse.flight_path_m();
251        self.flight_times_us
252            .iter()
253            .map(|u| (clock / u).powi(2))
254            .collect()
255    }
256
257    /// `step · Σ_j v_j P_k(u_j)` for each bin `k`: the counts a beam of
258    /// `values[j]` neutrons per µs at each grid point predicts, or, for any
259    /// other per-point values, the same linear map.  It differs from the
260    /// trapezoid rule by half a step of the end points' terms, where `P_k`
261    /// is zero, or at most [`NEGLIGIBLE_ARRIVAL_PROBABILITY`](crate::ikeda_carpenter::NEGLIGIBLE_ARRIVAL_PROBABILITY) at the fast end
262    /// for an unfolded pulse.  Features of the values narrower than the step
263    /// are not resolved.
264    ///
265    /// # Errors
266    /// [`FlightTimeGridError::ValuesLength`] unless `values` has one entry
267    /// per grid point.
268    pub fn predict(&self, values: &[f64]) -> Result<Vec<f64>, FlightTimeGridError> {
269        if values.len() != self.flight_times_us.len() {
270            return Err(FlightTimeGridError::ValuesLength {
271                expected: self.flight_times_us.len(),
272                found: values.len(),
273            });
274        }
275        let mut counts = vec![0.0; self.response.n_detector_bins()];
276        for (j, &v) in values.iter().enumerate() {
277            for (k, p) in self.response.row_entries(j) {
278                counts[k] += self.step_us * v * p;
279            }
280        }
281        Ok(counts)
282    }
283}
284
285fn crossing(
286    signed: &dyn Fn(f64) -> Result<f64, ResolutionParseError>,
287    fast: f64,
288    slow: f64,
289) -> Result<Option<f64>, ResolutionParseError> {
290    if signed(fast)? >= 0.0 || signed(slow)? < 0.0 {
291        return Ok(None);
292    }
293    let (mut low, mut high) = (fast, slow);
294    while high - low > f64::EPSILON * high {
295        let middle = 0.5 * (low + high);
296        if signed(middle)? < 0.0 {
297            low = middle;
298        } else {
299            high = middle;
300        }
301    }
302    Ok(Some(high))
303}
nereids_physics/flight_time_grid.rs

nereids_physics/
flight_time_grid.rs
