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` over the flight path `L` has energy
5//! `E = (TOF_FACTOR·L/u)²` and arrives at `t0 + u + delay`, the delay drawn
6//! from the Ikeda–Carpenter pulse at `E`.  The grid spans every flight time
7//! whose neutrons reach the 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 energy span can reach.
9//!
10//! Its points keep their energies at any other `t0'` and `L'`: point `j`,
11//! of flight time `u_j` at the grid's own `t0` and `L`, then arrives at
12//! `t0' + (L'/L)·u_j` and stands for `(L'/L)` times the step.
13
14use std::fmt;
15
16use rayon::prelude::*;
17
18use crate::ikeda_carpenter::{DetectorPulse, EnergyLaw, IkedaCarpenterParams};
19use crate::resolution::{ResolutionParseError, TOF_FACTOR};
20
21/// Largest number of grid points a window may need.
22pub const MAX_POINTS: usize = 100_000;
23
24/// Why a window has no flight-time grid.
25#[derive(Debug)]
26pub enum FlightTimeGridError {
27    /// Fewer than two edges, or edges not finite and strictly ascending.
28    InvalidTimeEdges,
29    /// The timing offset is not finite.
30    InvalidTimingOffset(f64),
31    /// The flight path is not finite and positive.
32    InvalidFlightPath(f64),
33    /// Not supported: the pulse's `parameter` law lengthens the pulse as
34    /// energy rises (α falling, or β falling with a storage term).  The range
35    /// relies on arrival times growing with flight time, which such a pulse
36    /// does not guarantee.
37    LengthensWithEnergy { parameter: &'static str },
38    /// The pulse cannot be evaluated.
39    Pulse(ResolutionParseError),
40    /// Neutrons with energies outside the pulse's energy span, `low_ev` to
41    /// `high_ev`, can reach the window.
42    OutsideCalibration { low_ev: f64, high_ev: f64 },
43    /// A grid at `step_us` would need more than [`MAX_POINTS`] points.
44    TooManyPoints { step_us: f64 },
45    /// `found` values were given for a grid of `expected` points.
46    ValuesLength { expected: usize, found: usize },
47}
48
49impl fmt::Display for FlightTimeGridError {
50    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
51        match self {
52            Self::InvalidTimeEdges => {
53                write!(
54                    f,
55                    "time edges must be at least two finite, strictly ascending values"
56                )
57            }
58            Self::InvalidTimingOffset(t0) => write!(f, "t0 = {t0} µs must be finite"),
59            Self::InvalidFlightPath(l) => {
60                write!(f, "the flight path, {l} m, must be finite and positive")
61            }
62            Self::LengthensWithEnergy { parameter } => write!(
63                f,
64                "unsupported pulse: its {parameter} law lengthens the pulse as energy rises; \
65                 the flight-time grid needs α and β that do not fall with energy"
66            ),
67            Self::Pulse(e) => write!(f, "pulse: {e}"),
68            Self::OutsideCalibration { low_ev, high_ev } => write!(
69                f,
70                "neutrons from outside the pulse's energy span, {low_ev} to \
71                 {high_ev} eV, can reach the window"
72            ),
73            Self::TooManyPoints { step_us } => write!(
74                f,
75                "a step of {step_us} µs needs more than {MAX_POINTS} grid points"
76            ),
77            Self::ValuesLength { expected, found } => write!(
78                f,
79                "{found} values were given for a grid of {expected} points"
80            ),
81        }
82    }
83}
84
85impl std::error::Error for FlightTimeGridError {}
86
87impl From<ResolutionParseError> for FlightTimeGridError {
88    fn from(e: ResolutionParseError) -> Self {
89        Self::Pulse(e)
90    }
91}
92
93/// One value per bin for each grid point, kept where non-zero, with the
94/// step each point stands for.
95#[derive(Debug, Clone)]
96pub struct Rows {
97    step_us: f64,
98    bins: usize,
99    row_offsets: Vec<usize>,
100    columns: Vec<usize>,
101    values: Vec<f64>,
102}
103
104impl Rows {
105    fn new<'a>(step_us: f64, bins: usize, rows: impl Iterator<Item = &'a [(usize, f64)]>) -> Self {
106        let mut row_offsets = vec![0];
107        let mut columns = Vec::new();
108        let mut values = Vec::new();
109        for row in rows {
110            for &(k, v) in row {
111                columns.push(k);
112                values.push(v);
113            }
114            row_offsets.push(values.len());
115        }
116        Self {
117            step_us,
118            bins,
119            row_offsets,
120            columns,
121            values,
122        }
123    }
124
125    /// `step · Σ_j v_j R_k(j)` for each bin `k`, with `values[j]` the value
126    /// at grid point `j`.
127    ///
128    /// # Errors
129    /// [`FlightTimeGridError::ValuesLength`] unless `values` has one entry
130    /// per grid point.
131    pub fn predict(&self, values: &[f64]) -> Result<Vec<f64>, FlightTimeGridError> {
132        let points = self.row_offsets.len() - 1;
133        if values.len() != points {
134            return Err(FlightTimeGridError::ValuesLength {
135                expected: points,
136                found: values.len(),
137            });
138        }
139        let mut counts = vec![0.0; self.bins];
140        for (&v, row) in values.iter().zip(self.row_offsets.windows(2)) {
141            for (&k, &p) in self.columns[row[0]..row[1]]
142                .iter()
143                .zip(&self.values[row[0]..row[1]])
144            {
145                counts[k] += self.step_us * v * p;
146            }
147        }
148        Ok(counts)
149    }
150}
151
152/// Uniform flight-time grid over a window, with the probability that each
153/// grid point's neutrons are counted in each bin.
154#[derive(Debug, Clone)]
155pub struct FlightTimeGrid {
156    time_edges_us: Vec<f64>,
157    t0_us: f64,
158    flight_path_m: f64,
159    pulse: DetectorPulse,
160    range_us: (f64, f64),
161    step_us: f64,
162    flight_times_us: Vec<f64>,
163    rows: Rows,
164}
165
166impl FlightTimeGrid {
167    /// The grid for bins `time_edges_us` (µs) with timing offset `t0_us`,
168    /// flight path `flight_path_m` and `pulse`.  The step, half the pulse's shorter rise at the two ends of the flight
169    /// times whose neutrons can reach the bins, is a starting step: the counts
170    /// are converged once [`Self::halved`] no longer changes them.  The grid
171    /// extends one step beyond those flight times at each end.  For a storage
172    /// fraction that rises with energy, the latest arrival that sets the
173    /// earliest flight time takes the fraction at its largest over the energy
174    /// span.
175    ///
176    /// # Errors
177    /// See [`FlightTimeGridError`].
178    pub fn new(
179        time_edges_us: &[f64],
180        t0_us: f64,
181        flight_path_m: f64,
182        pulse: &DetectorPulse,
183    ) -> Result<Self, FlightTimeGridError> {
184        Self::spanning(time_edges_us, (t0_us, t0_us), flight_path_m, pulse)
185    }
186
187    /// The grid of [`Self::new`] for every timing offset from `t0_us` to
188    /// `last_t0_us`: its flight times are those whose neutrons can reach the
189    /// bins at either, and [`Self::rows`] are at `t0_us`.
190    ///
191    /// # Errors
192    /// See [`FlightTimeGridError`].
193    pub fn spanning(
194        time_edges_us: &[f64],
195        (t0_us, last_t0_us): (f64, f64),
196        flight_path_m: f64,
197        pulse: &DetectorPulse,
198    ) -> Result<Self, FlightTimeGridError> {
199        if time_edges_us.len() < 2
200            || !time_edges_us.iter().all(|t| t.is_finite())
201            || time_edges_us.windows(2).any(|w| w[0] >= w[1])
202        {
203            return Err(FlightTimeGridError::InvalidTimeEdges);
204        }
205        let range = |t0_us: f64| {
206            range_us(
207                time_edges_us,
208                (t0_us, flight_path_m),
209                pulse.params(),
210                pulse.energy_span_ev(),
211            )
212        };
213        let ((first_lo, first_hi), (last_lo, last_hi)) = (range(t0_us)?, range(last_t0_us)?);
214        let (u_lo, u_hi) = (first_lo.min(last_lo), first_hi.max(last_hi));
215        let energy = |u: f64| (TOF_FACTOR * flight_path_m / u).powi(2);
216        let rise = pulse
217            .rise_us(energy(u_lo))?
218            .min(pulse.rise_us(energy(u_hi))?);
219        let intervals = ((u_hi - u_lo) / (0.5 * rise)).ceil();
220        let step_us = (u_hi - u_lo) / intervals;
221        Self::build(
222            time_edges_us,
223            (t0_us, flight_path_m),
224            pulse,
225            (u_lo - step_us, u_hi + step_us),
226            intervals as usize + 2,
227        )
228    }
229
230    /// Whether the grid's points, arriving as [`Self::rows_at`] places them,
231    /// cover every flight time whose neutrons can reach the bins when the
232    /// timing offset is `t0_us`, the flight path `flight_path_m` and the
233    /// pulse's laws `params`.
234    ///
235    /// # Errors
236    /// As [`Self::new`] for that timing offset, flight path and pulse.
237    pub fn covers(
238        &self,
239        t0_us: f64,
240        flight_path_m: f64,
241        params: &IkedaCarpenterParams,
242    ) -> Result<bool, FlightTimeGridError> {
243        let (u_lo, u_hi) = range_us(
244            &self.time_edges_us,
245            (t0_us, flight_path_m),
246            params,
247            self.pulse.energy_span_ev(),
248        )?;
249        let scale = self.flight_path_m / flight_path_m;
250        Ok(u_lo * scale >= self.range_us.0 && u_hi * scale <= self.range_us.1)
251    }
252}
253
254fn range_us(
255    time_edges_us: &[f64],
256    (t0_us, flight_path_m): (f64, f64),
257    params: &IkedaCarpenterParams,
258    (e_min, e_max): (f64, f64),
259) -> Result<(f64, f64), FlightTimeGridError> {
260    if !t0_us.is_finite() {
261        return Err(FlightTimeGridError::InvalidTimingOffset(t0_us));
262    }
263    if !(flight_path_m.is_finite() && flight_path_m > 0.0) {
264        return Err(FlightTimeGridError::InvalidFlightPath(flight_path_m));
265    }
266    let change = |law: &EnergyLaw| law.eval(e_max) - law.eval(e_min);
267    let storage = params.r.eval(e_min).max(params.r.eval(e_max));
268    for (parameter, lengthens) in [
269        ("α", change(&params.alpha) < 0.0),
270        ("β", change(&params.beta) < 0.0 && storage > 0.0),
271    ] {
272        if lengthens {
273            return Err(FlightTimeGridError::LengthensWithEnergy { parameter });
274        }
275    }
276    let longest = IkedaCarpenterParams {
277        r: if change(&params.r) > 0.0 {
278            EnergyLaw::Const(storage)
279        } else {
280            params.r.clone()
281        },
282        ..params.clone()
283    };
284    let clock = TOF_FACTOR * flight_path_m;
285    let energy = |u: f64| (clock / u).powi(2);
286    let first_edge = time_edges_us[0] - t0_us;
287    let last_edge = time_edges_us[time_edges_us.len() - 1] - t0_us;
288    let latest = |u: f64| {
289        longest
290            .delays_us(energy(u))
291            .map(|(_, last)| u + last - first_edge)
292    };
293    let earliest = |u: f64| {
294        params
295            .delays_us(energy(u))
296            .map(|(first, _)| u + first - last_edge)
297    };
298    let (fastest, slowest) = (clock / e_max.sqrt(), clock / e_min.sqrt());
299    let (Some(u_lo), Some(u_hi)) = (
300        crossing(&latest, fastest, slowest)?,
301        crossing(&earliest, fastest, slowest)?,
302    ) else {
303        return Err(FlightTimeGridError::OutsideCalibration {
304            low_ev: e_min,
305            high_ev: e_max,
306        });
307    };
308    Ok((u_lo, u_hi))
309}
310
311impl FlightTimeGrid {
312    fn build(
313        time_edges_us: &[f64],
314        (t0_us, flight_path_m): (f64, f64),
315        pulse: &DetectorPulse,
316        range_us: (f64, f64),
317        intervals: usize,
318    ) -> Result<Self, FlightTimeGridError> {
319        let (u_lo, u_hi) = range_us;
320        let step_us = (u_hi - u_lo) / intervals as f64;
321        if intervals >= MAX_POINTS {
322            return Err(FlightTimeGridError::TooManyPoints { step_us });
323        }
324        let flight_times_us: Vec<f64> =
325            (0..=intervals).map(|j| u_lo + step_us * j as f64).collect();
326        let mut grid = Self {
327            time_edges_us: time_edges_us.to_vec(),
328            t0_us,
329            flight_path_m,
330            pulse: pulse.clone(),
331            range_us,
332            step_us,
333            flight_times_us,
334            rows: Rows::new(step_us, 0, std::iter::empty()),
335        };
336        grid.rows = grid.rows_at(t0_us, flight_path_m, pulse.params())?;
337        Ok(grid)
338    }
339
340    /// The same window at half the step.
341    ///
342    /// # Errors
343    /// As [`Self::new`], when the halved grid exceeds [`MAX_POINTS`].
344    pub fn halved(&self) -> Result<Self, FlightTimeGridError> {
345        Self::build(
346            &self.time_edges_us,
347            (self.t0_us, self.flight_path_m),
348            &self.pulse,
349            self.range_us,
350            2 * (self.flight_times_us.len() - 1),
351        )
352    }
353
354    /// The probability that each grid point's neutrons are counted in each
355    /// bin when the timing offset is `t0_us`, the flight path `flight_path_m`
356    /// and the pulse's laws `params`, with the step each point then stands
357    /// for.  Neutrons from flight times the grid does not
358    /// [cover](Self::covers) there are not counted.
359    ///
360    /// # Errors
361    /// [`FlightTimeGridError::InvalidTimingOffset`] or
362    /// [`FlightTimeGridError::InvalidFlightPath`] unless `t0_us` is finite and
363    /// `flight_path_m` finite and positive; [`FlightTimeGridError::Pulse`] if
364    /// the pulse cannot be evaluated there.
365    pub fn rows_at(
366        &self,
367        t0_us: f64,
368        flight_path_m: f64,
369        params: &IkedaCarpenterParams,
370    ) -> Result<Rows, FlightTimeGridError> {
371        let [rows] = self.rows_with(t0_us, flight_path_m, |energy, arrival| {
372            Ok([params.bin_probabilities_at(energy, arrival, &self.time_edges_us)?])
373        })?;
374        Ok(rows)
375    }
376
377    /// The derivative of [`Self::rows_at`]'s probabilities with respect to
378    /// each point's arrival time, per µs, with the same step.
379    ///
380    /// # Errors
381    /// As [`Self::rows_at`].
382    pub fn arrival_slopes_at(
383        &self,
384        t0_us: f64,
385        flight_path_m: f64,
386        params: &IkedaCarpenterParams,
387    ) -> Result<Rows, FlightTimeGridError> {
388        let [rows] = self.rows_with(t0_us, flight_path_m, |energy, arrival| {
389            Ok([params.bin_arrival_slopes_at(energy, arrival, &self.time_edges_us)?])
390        })?;
391        Ok(rows)
392    }
393
394    /// The derivatives of [`Self::rows_at`]'s probabilities with respect to
395    /// `α`, `β`, `R` and the square of the triangle's FWHM, each at its point's
396    /// energy, with
397    /// the same step: [`IkedaCarpenterParams::bin_pulse_slopes_at`] for every
398    /// point.
399    ///
400    /// # Errors
401    /// As [`Self::rows_at`].
402    pub fn pulse_slopes_at(
403        &self,
404        t0_us: f64,
405        flight_path_m: f64,
406        params: &IkedaCarpenterParams,
407    ) -> Result<[Rows; 4], FlightTimeGridError> {
408        self.rows_with(t0_us, flight_path_m, |energy, arrival| {
409            params.bin_pulse_slopes_at(energy, arrival, &self.time_edges_us)
410        })
411    }
412
413    /// The derivatives of [`Self::rows_at`]'s probabilities with respect to
414    /// the detector tail's `f` and `τ_D`, with the same step:
415    /// [`IkedaCarpenterParams::bin_tail_slopes_at`] for every point.
416    ///
417    /// # Errors
418    /// As [`Self::rows_at`], and [`FlightTimeGridError::Pulse`] for a pulse
419    /// without a detector tail.
420    pub fn tail_slopes_at(
421        &self,
422        t0_us: f64,
423        flight_path_m: f64,
424        params: &IkedaCarpenterParams,
425    ) -> Result<[Rows; 2], FlightTimeGridError> {
426        self.rows_with(t0_us, flight_path_m, |energy, arrival| {
427            params.bin_tail_slopes_at(energy, arrival, &self.time_edges_us)
428        })
429    }
430
431    fn rows_with<const N: usize>(
432        &self,
433        t0_us: f64,
434        flight_path_m: f64,
435        row: impl Fn(f64, f64) -> Result<[Vec<f64>; N], ResolutionParseError> + Sync,
436    ) -> Result<[Rows; N], FlightTimeGridError> {
437        if !t0_us.is_finite() {
438            return Err(FlightTimeGridError::InvalidTimingOffset(t0_us));
439        }
440        if !(flight_path_m.is_finite() && flight_path_m > 0.0) {
441            return Err(FlightTimeGridError::InvalidFlightPath(flight_path_m));
442        }
443        let scale = flight_path_m / self.flight_path_m;
444        let clock = TOF_FACTOR * self.flight_path_m;
445        let points: Vec<[Vec<(usize, f64)>; N]> = self
446            .energies_ev()
447            .par_iter()
448            .map(|&energy| {
449                Ok(
450                    row(energy, t0_us + scale * (clock / energy.sqrt()))?.map(|values| {
451                        values
452                            .into_iter()
453                            .enumerate()
454                            .filter(|&(_, v)| v != 0.0)
455                            .collect()
456                    }),
457                )
458            })
459            .collect::<Result<_, ResolutionParseError>>()?;
460        let bins = self.time_edges_us.len() - 1;
461        Ok(std::array::from_fn(|n| {
462            Rows::new(
463                self.step_us * scale,
464                bins,
465                points.iter().map(|point| point[n].as_slice()),
466            )
467        }))
468    }
469
470    /// The probabilities at the grid's own timing offset, flight path and
471    /// pulse.
472    #[must_use]
473    pub fn rows(&self) -> &Rows {
474        &self.rows
475    }
476
477    /// The pulse the grid was built at.
478    #[must_use]
479    pub fn pulse(&self) -> &DetectorPulse {
480        &self.pulse
481    }
482
483    /// The timing offset in µs the grid was built at.
484    #[must_use]
485    pub fn t0_us(&self) -> f64 {
486        self.t0_us
487    }
488
489    /// The flight path in m the grid was built at.
490    #[must_use]
491    pub fn flight_path_m(&self) -> f64 {
492        self.flight_path_m
493    }
494
495    /// `(u_lo, u_hi)`, the flight-time range in µs at the timing offset and
496    /// flight path the grid was built at.
497    #[must_use]
498    pub fn range_us(&self) -> (f64, f64) {
499        self.range_us
500    }
501
502    /// The grid step in µs at the flight path the grid was built at.
503    #[must_use]
504    pub fn step_us(&self) -> f64 {
505        self.step_us
506    }
507
508    /// The grid's flight times in µs, ascending.
509    #[must_use]
510    pub fn flight_times_us(&self) -> &[f64] {
511        &self.flight_times_us
512    }
513
514    /// The time-bin edges in µs, ascending.
515    #[must_use]
516    pub fn time_edges_us(&self) -> &[f64] {
517        &self.time_edges_us
518    }
519
520    /// The energy, in eV, of each grid point, descending.
521    #[must_use]
522    pub fn energies_ev(&self) -> Vec<f64> {
523        let clock = TOF_FACTOR * self.flight_path_m;
524        self.flight_times_us
525            .iter()
526            .map(|u| (clock / u).powi(2))
527            .collect()
528    }
529
530    /// `step · Σ_j v_j P_k(u_j)` for each bin `k`: the counts a beam of
531    /// `values[j]` neutrons per µs at each grid point predicts, or, for any
532    /// other per-point values, the same linear map.  It differs from the
533    /// trapezoid rule by half a step of the end points' terms, where `P_k`
534    /// is zero, or at most [`NEGLIGIBLE_ARRIVAL_PROBABILITY`](crate::ikeda_carpenter::NEGLIGIBLE_ARRIVAL_PROBABILITY) at the fast end
535    /// for an unfolded pulse.  Features of the values narrower than the step
536    /// are not resolved.
537    ///
538    /// # Errors
539    /// [`FlightTimeGridError::ValuesLength`] unless `values` has one entry
540    /// per grid point.
541    pub fn predict(&self, values: &[f64]) -> Result<Vec<f64>, FlightTimeGridError> {
542        self.rows.predict(values)
543    }
544}
545
546fn crossing(
547    signed: &dyn Fn(f64) -> Result<f64, ResolutionParseError>,
548    fast: f64,
549    slow: f64,
550) -> Result<Option<f64>, ResolutionParseError> {
551    if signed(fast)? >= 0.0 || signed(slow)? < 0.0 {
552        return Ok(None);
553    }
554    let (mut low, mut high) = (fast, slow);
555    while high - low > f64::EPSILON * high {
556        let middle = 0.5 * (low + high);
557        if signed(middle)? < 0.0 {
558            low = middle;
559        } else {
560            high = middle;
561        }
562    }
563    Ok(Some(high))
564}
nereids_physics/flight_time_grid.rs

nereids_physics/
flight_time_grid.rs
