reaction_kinematics

 38class Reaction:
 39    """
 40    Defines a two-body nuclear reaction: projectile + target → ejectile + recoil.
 41
 42    All energy-dependent methods accept a ``beam_energy`` parameter (beam kinetic
 43    energy, MeV by default). Internal computation is cached per energy so repeated
 44    calls at the same energy are efficient.
 45
 46    Parameters
 47    ----------
 48    m1, m2, m3, m4 : str, MassInput, or float
 49        Masses of projectile, target, ejectile, and recoil.
 50        Strings like ``"p"``, ``"12C"``, ``"alpha"`` are looked up in the mass table.
 51        Floats require ``mass_unit``.
 52    mass_unit : str or EnergyUnit, optional
 53        Unit for numeric masses (e.g. ``"MeV"``, ``"keV"``).
 54
 55    Attributes
 56    ----------
 57    ncoscm : int
 58        Number of CM angle grid points (default 500). Changing this invalidates
 59        the cache.
 60
 61    Examples
 62    --------
 63    >>> rxn = Reaction("p", "3H", "n", "3He")
 64    >>> rxn.q_value
 65    -0.763...
 66    >>> data = rxn.kinematics_table_at_beam_energy(1.2)
 67    >>> result = rxn.kinematics_at_beam_energy_and_variable(1.2, "theta3", 0.5)
 68    >>> branches = rxn.kinematics_curve_at_angle(np.linspace(1.0, 5.0, 200), np.deg2rad(30))
 69    """
 70
 71    def __init__(
 72        self,
 73        m1: MassArg,
 74        m2: MassArg,
 75        m3: MassArg,
 76        m4: MassArg,
 77        *,
 78        mass_unit: str | EnergyUnit | None = None,
 79    ) -> None:
 80        self._m1 = _parse_mass(m1, mass_unit)
 81        self._m2 = _parse_mass(m2, mass_unit)
 82        self._m3 = _parse_mass(m3, mass_unit)
 83        self._m4 = _parse_mass(m4, mass_unit)
 84        self.__ncoscm: int = 500
 85        self._cached_ek: float | None = None
 86        self._table: dict[str, list[float]] | None = None
 87        # per-energy kinematic state, populated by _bind
 88        self._nogo: bool = False
 89        self._pcmp: float | None = None
 90        self._thesinh: float | None = None
 91        self._thecosh: float | None = None
 92        self._e03: float | None = None
 93        self._e04: float | None = None
 94        # lab-frame energy extrema
 95        self._emax3: float | None = None
 96        self._emin3: float | None = None
 97        self._emax4: float | None = None
 98        self._emin4: float | None = None
 99        # max lab angles and associated quantities (None = no forward maximum)
100        self._theta3max: float | None = None
101        self._theta4max: float | None = None
102        self._e3atmaxang: float | None = None
103        self._e4atmaxang: float | None = None
104        self._cmcos3max: float | None = None
105        self._cmcos4max: float | None = None
106
107    @property
108    def _ncoscm(self) -> int:
109        return self.__ncoscm
110
111    @_ncoscm.setter
112    def _ncoscm(self, val: int) -> None:
113        self.__ncoscm = val
114        self._cached_ek = None
115        self._table = None
116
117    @property
118    def q_value(self) -> float:
119        """Q-value of the reaction in MeV: Q = m1 + m2 - m3 - m4."""
120        return self._m1 + self._m2 - self._m3 - self._m4
121
122    def _bind(self, ek_mev: float) -> None:
123        """Compute and cache kinematic state for the given beam energy."""
124        if ek_mev == self._cached_ek:
125            return
126        self._cached_ek = ek_mev
127        self._table = None
128        self._compute(ek_mev)
129
130    def _compute(self, ek_mev: float) -> None:
131        self._nogo = False
132        self._pcmp = self._thesinh = self._thecosh = self._e03 = self._e04 = None
133
134        # Mandelstam s
135        s = (self._m1 + self._m2) ** 2 + 2.0 * self._m2 * ek_mev
136        if s <= 0.0:
137            self._nogo = True
138            return
139
140        # initial CM momentum
141        pcm2 = (s - self._m1**2 - self._m2**2) ** 2 - 4.0 * self._m1**2 * self._m2**2
142        if pcm2 < 0:
143            self._nogo = True
144            return
145        pcm = math.sqrt(pcm2 / (4.0 * s))
146
147        # CM rapidity → boost parameters
148        acmratio = (math.sqrt(self._m2**2 + pcm**2) + pcm) / self._m2
149        cmrap = math.log(acmratio)
150        self._thesinh = math.sinh(cmrap)
151        self._thecosh = math.cosh(cmrap)
152
153        # final-state CM momentum
154        pcmp2 = (s - self._m3**2 - self._m4**2) ** 2 - 4.0 * self._m3**2 * self._m4**2
155        if pcmp2 < 0:
156            self._nogo = True
157            return
158        self._pcmp = math.sqrt(pcmp2 / (4.0 * s))
159
160        # CM total energies of outgoing particles
161        self._e03 = math.sqrt(self._pcmp**2 + self._m3**2)
162        self._e04 = math.sqrt(self._pcmp**2 + self._m4**2)
163
164        # lab-frame energy extrema
165        self._emax3 = self._e03 * self._thecosh + self._pcmp * self._thesinh - self._m3
166        self._emin3 = self._e03 * self._thecosh - self._pcmp * self._thesinh - self._m3
167        self._emax4 = self._e04 * self._thecosh + self._pcmp * self._thesinh - self._m4
168        self._emin4 = self._e04 * self._thecosh - self._pcmp * self._thesinh - self._m4
169
170        # reset max-angle quantities to sentinel values
171        self._theta3max = None
172        self._theta4max = None
173        self._e3atmaxang = None
174        self._e4atmaxang = None
175        self._cmcos3max = None
176        self._cmcos4max = None
177
178        # max ejectile lab angle (only exists when pcmp < m3 * sinh(rapidity))
179        thetatest3: float | None = None
180        if self._m3 > 0.0:
181            thetatest3 = self._pcmp / (self._m3 * self._thesinh)
182            if thetatest3 < 1.0:
183                self._theta3max = math.asin(thetatest3)
184                patmax = (self._e03 * math.cos(self._theta3max) * self._thesinh) / (
185                    1.0 + thetatest3**2 * self._thesinh**2
186                )
187                eatmax = math.sqrt(patmax**2 + self._m3**2)
188                self._e3atmaxang = eatmax - self._m3
189                self._cmcos3max = (eatmax - self._e03 * self._thecosh) / (
190                    self._pcmp * self._thesinh
191                )
192
193        # elastic case: forward-angle symmetry forces theta3max = 90°
194        if (self._m1 + self._m2) == (self._m3 + self._m4) and thetatest3 is not None:
195            if abs(thetatest3 - 1.0) < 1e-3:
196                self._theta3max = math.pi / 2.0
197                self._cmcos3max = -1.0
198                self._e3atmaxang = (
199                    self._e03 * self._thecosh
200                    + self._cmcos3max * self._pcmp * self._thesinh
201                    - self._m3
202                )
203
204        # max recoil lab angle
205        thetatest4: float | None = None
206        if self._m4 > 0.0:
207            thetatest4 = self._pcmp / (self._m4 * self._thesinh)
208            if thetatest4 < 1.0:
209                self._theta4max = math.asin(thetatest4)
210                patmax = (self._e04 * math.cos(self._theta4max) * self._thesinh) / (
211                    1.0 + thetatest4**2 * self._thesinh**2
212                )
213                eatmax = math.sqrt(patmax**2 + self._m4**2)
214                self._e4atmaxang = eatmax - self._m4
215                self._cmcos4max = (eatmax - self._e04 * self._thecosh) / (
216                    self._pcmp * self._thesinh
217                )
218
219        # elastic case: theta4max = 90°
220        if (self._m1 + self._m2) == (self._m3 + self._m4) and thetatest4 is not None:
221            if abs(thetatest4 - 1.0) < 1e-3:
222                self._theta4max = math.pi / 2.0
223                self._cmcos4max = 1.0
224                self._e4atmaxang = (
225                    self._e04 * self._thecosh
226                    - self._cmcos4max * self._pcmp * self._thesinh
227                    - self._m4
228                )
229
230    def _kinematics_at_coscm(self, coscm: float) -> dict[str, float]:
231        if (
232            self._pcmp is None
233            or self._thecosh is None
234            or self._e03 is None
235            or self._thesinh is None
236            or self._e04 is None
237        ):
238            raise ValueError("Kinematic state not computed — call _bind first")
239
240        sincm = math.sqrt(max(0.0, 1.0 - coscm**2))
241
242        ppar3 = self._pcmp * self._thecosh * coscm + self._e03 * self._thesinh
243        pperp3 = self._pcmp * sincm
244        ptot3 = math.hypot(ppar3, pperp3)
245
246        ppar4 = -self._pcmp * self._thecosh * coscm + self._e04 * self._thesinh
247        pperp4 = self._pcmp * sincm
248        ptot4 = math.hypot(ppar4, pperp4)
249
250        e3 = self._e03 * self._thecosh + coscm * self._pcmp * self._thesinh - self._m3
251        e4 = self._e04 * self._thecosh - coscm * self._pcmp * self._thesinh - self._m4
252
253        return {
254            "coscm": coscm,
255            "theta_cm": math.acos(coscm),
256            "theta3": math.acos(ppar3 / ptot3) if ptot3 > 0 else 0.0,
257            "theta4": math.acos(ppar4 / ptot4) if ptot4 > 0 else 0.0,
258            "e3": e3,
259            "e4": e4,
260            "v3": ptot3 / (e3 + self._m3),
261            "v4": ptot4 / (e4 + self._m4),
262            "p3": ptot3,
263            "p4": ptot4,
264        }
265
266    def _build_table(self) -> None:
267        """Build the interpolation table over the full CM angle grid."""
268        keys = ["coscm", "theta_cm", "theta3", "theta4", "e3", "e4", "v3", "v4", "p3", "p4"]
269        table: dict[str, list[float]] = {k: [] for k in keys}
270        for i in range(-self._ncoscm, self._ncoscm + 1):
271            row = self._kinematics_at_coscm(i / self._ncoscm)
272            for k in keys:
273                table[k].append(row[k])
274        self._table = table
275
276    def kinematics_table_at_beam_energy(
277        self, beam_energy: float, *, energy_unit: EnergyUnit = EnergyUnit.MeV
278    ) -> dict[str, npt.NDArray[np.float64]]:
279        """
280        Compute full kinematics over a CM angle grid.
281
282        Parameters
283        ----------
284        beam_energy : float
285            Beam kinetic energy.
286        energy_unit : EnergyUnit, optional
287            Unit of ``beam_energy`` (default MeV).
288
289        Returns
290        -------
291        dict[str, np.ndarray]
292            Keys: ``"coscm"``, ``"theta_cm"``, ``"theta3"``, ``"theta4"``,
293            ``"e3"``, ``"e4"``, ``"v3"``, ``"v4"``.
294
295        Raises
296        ------
297        ValueError
298            If the reaction is kinematically forbidden at this energy.
299        """
300        ek_mev = _parse_energy(beam_energy, energy_unit)
301        self._bind(ek_mev)
302        if self._nogo:
303            raise ValueError(f"Reaction kinematically forbidden at beam_energy={ek_mev} MeV")
304        keys = ["coscm", "theta_cm", "theta3", "theta4", "e3", "e4", "v3", "v4"]
305        rows = [
306            self._kinematics_at_coscm(i / self._ncoscm)
307            for i in range(-self._ncoscm, self._ncoscm + 1)
308        ]
309        return {k: np.array([row[k] for row in rows]) for k in keys}
310
311    def kinematics_at_beam_energy_and_variable(
312        self,
313        beam_energy: float,
314        var_name: str,
315        var_value: float,
316        *,
317        energy_unit: EnergyUnit = EnergyUnit.MeV,
318        return_variables: str | list[str] | None = None,
319        duplicate_tol: float = 1e-6,
320    ) -> dict[str, list[float]]:
321        """
322        Interpolate kinematic quantities at a fixed beam energy and kinematic variable value.
323
324        Always returns lists to handle multi-valued cases (e.g. two ejectile
325        energies at the same lab angle).
326
327        Parameters
328        ----------
329        beam_energy : float
330            Beam kinetic energy.
331        var_name : str
332            Independent variable name, e.g. ``"theta3"``, ``"theta_cm"``,
333            ``"coscm"``, ``"theta4"``.
334        var_value : float
335            Value to evaluate at (radians for angles).
336        energy_unit : EnergyUnit, optional
337            Unit of ``beam_energy`` (default MeV).
338        return_variables : str or list of str, optional
339            Dependent variables to return. ``None`` returns all.
340        duplicate_tol : float, optional
341            Tolerance for merging near-duplicate solutions (default 1e-6).
342
343        Returns
344        -------
345        dict of str -> list
346            Each value is a list of solutions, sorted descending by ``e3``.
347
348        Raises
349        ------
350        ValueError
351            If ``var_value`` is outside the physical range.
352
353        Examples
354        --------
355        >>> rxn = Reaction("p", "3H", "n", "3He")
356        >>> rxn.kinematics_at_beam_energy_and_variable(1.2, "theta3", 0.5, return_variables=["e3", "v3"])
357        {'e3': [...], 'v3': [...]}
358        """
359        ek_mev = _parse_energy(beam_energy, energy_unit)
360        self._bind(ek_mev)
361
362        if self._table is None:
363            self._build_table()
364        assert self._table is not None
365
366        xs = self._table[var_name]
367
368        if isinstance(return_variables, str):
369            return_variables = [return_variables]
370        if return_variables is None:
371            return_variables = list(self._table.keys())
372
373        solutions = []
374
375        exact_idx = np.where(np.isclose(xs, var_value, atol=1e-12))[0]
376        if len(exact_idx) > 0:
377            for i in exact_idx:
378                solutions.append({k: self._table[k][i] for k in return_variables})
379        else:
380            found = False
381            for i in range(len(xs) - 1):
382                x0, x1 = xs[i], xs[i + 1]
383                if (x0 - var_value) * (x1 - var_value) <= 0 and x0 != x1:
384                    found = True
385                    t = (var_value - x0) / (x1 - x0)
386                    solutions.append(
387                        {
388                            k: self._table[k][i] + t * (self._table[k][i + 1] - self._table[k][i])
389                            for k in return_variables
390                        }
391                    )
392            if not found:
393                raise ValueError(f"{var_name}={var_value} outside physical range")
394
395        e_key = "e3" if "e3" in return_variables else return_variables[0]
396        unique: list = []
397        for sol in solutions:
398            if not any(abs(sol[e_key] - u[e_key]) < duplicate_tol for u in unique):
399                unique.append(sol)
400        unique.sort(key=lambda s: s[e_key], reverse=True)
401
402        return {k: [s[k] for s in unique] for k in return_variables}
403
404    def kinematics_curve_at_angle(
405        self,
406        beam_energy_array: Iterable[float],
407        theta: float,
408        *,
409        angle_unit: AngleUnit = AngleUnit.rad,
410        energy_unit: EnergyUnit = EnergyUnit.MeV,
411    ) -> list[dict[str, npt.NDArray[np.float64]]]:
412        """
413        Compute ejectile kinematics at a fixed lab angle over a range of beam energies.
414
415        Returns two branches (high- and low-energy) as a list of two dicts. Each
416        dict contains arrays indexed by beam energy, with ``NaN`` where that branch
417        does not exist. Branch 0 is always the higher-energy solution.
418
419        Parameters
420        ----------
421        beam_energy_array : array-like
422            Beam energies to sweep.
423        theta : float
424            Fixed lab angle of the ejectile.
425        angle_unit : AngleUnit, optional
426            Unit of ``theta`` (default radians).
427        energy_unit : EnergyUnit, optional
428            Unit of ``beam_energy_array`` values (default MeV).
429
430        Returns
431        -------
432        list of two dicts, each with keys:
433            ``"ek"``, ``"e3"``, ``"e4"``, ``"theta4"``, ``"v3"``, ``"v4"``.
434
435        Examples
436        --------
437        >>> rxn = Reaction("p", "3H", "n", "3He")
438        >>> branches = rxn.kinematics_curve_at_angle(np.linspace(1.0, 5.0, 200), np.deg2rad(30))
439        >>> for b in branches:
440        ...     plt.plot(b["ek"], b["e3"])
441        """
442        if isinstance(angle_unit, str):
443            angle_unit = AngleUnit[angle_unit]
444        theta_rad = theta * angle_unit.value
445
446        keys = ["e3", "e4", "theta4", "v3", "v4"]
447        branches = [
448            {"ek": [], **{k: [] for k in keys}},
449            {"ek": [], **{k: [] for k in keys}},
450        ]
451
452        for ek in beam_energy_array:
453            ek_mev = _parse_energy(ek, energy_unit)
454            try:
455                row = self.kinematics_at_beam_energy_and_variable(
456                    ek_mev, "theta3", theta_rad, return_variables=keys
457                )
458            except ValueError:
459                solutions = []
460            else:
461                n = len(row["e3"])
462                solutions = [{k: row[k][i] for k in keys} for i in range(n)]
463
464            for i, branch in enumerate(branches):
465                branch["ek"].append(ek_mev)
466                sol = solutions[i] if i < len(solutions) else None
467                for k in keys:
468                    branch[k].append(sol[k] if sol is not None else float("nan"))
469
470        return [{k: np.array(v) for k, v in branch.items()} for branch in branches]