Complete general maths functionality for absorbers in i19 transmission/attenuation

src/dodal/common/general_maths/absorber_geometry.py

@@ -0,0 +1,130 @@
+from collections.abc import Callable
+from typing import Annotated, Literal, Self
+
+from pydantic import (
+    BaseModel,
+    Field,
+    PrivateAttr,
+    StrictFloat,
+    StrictInt,
+    field_validator,
+    model_validator,
+)
+
+from .arithmetic_conversions import (
+    convert_cm_to_mm,
+    convert_microns_to_cm,
+    convert_mm_to_cm,
+    get_straight_line_y,
+)
+
+SupportedThicknessUnits = Annotated[
+    Literal["cm", "mm", "um", "micron"], "Supported thickness units"
+]
+
+
+class FoilGeometry(BaseModel):
+    """Represents flat shape of a specific foil absorber.
+
+    Initialised with a foil thickness quantity in a supported thickness unit and,
+    where necessary, converting this to a cm thickness.
+    N.B. The missing "air" slot in a foil wheel is represented as a canonical / system "OUT" position,
+    of the absorber - rather than as a zero thickness foil.
+
+    Args:
+        unit (str): one element from the set of supported choices - "cm", "mm", "um" or "micron"
+        numerical_value (float): non-zero positive - physical thickness of the foil in the given units
+    """
+
+    # Class-level private dictionary - defined once only when the class is loaded
+    _convertors: dict[str, Callable[[float], float]] = {
+        "cm": lambda t: t,
+        "mm": convert_mm_to_cm,
+        "um": convert_microns_to_cm,
+        "micron": convert_microns_to_cm,
+    }
+
+    unit: str
+    numerical_value: Annotated[float, Field(gt=0.0)]
+    _foil_thickness: float = PrivateAttr(default=0.0)
+
+    @model_validator(mode="after")
+    def _calculate_thickness(self) -> "FoilGeometry":
+        # Look up directly from the conversion dictionary
+        self._foil_thickness = self._convertors[self.unit](self.numerical_value)
+        return self
+
+    def get_thickness_cm(self) -> float:
+        return self._foil_thickness
+
+
+class WedgeGeometry(BaseModel):
+    """Represents a semi infinite triangular prism wedge shape, growing out indefinitely from a tip.
+    Initialised with cotangent of the wedge angle and offset between motor axis zero and wedge tip.
+    N.B. absorber thickness are best worked out in cm (science), motor positions prefer mm (controls).
+    Values used for inputs come from beamline scientist calibrations on real absorber.
+
+    Args:
+        tip_mm (float): unrestricted value, point on motor axis scale where modelled taper reaches zero thickness.
+        taper_cotangent (float): taper angle cotangent, sign depends on motor direction convention, bounded to keep wedge angle thin.
+    """
+
+    tip_mm: StrictFloat | StrictInt
+    taper_cotangent: StrictFloat | StrictInt
+    _taper_gradient: float = PrivateAttr(default=0.0)
+
+    @field_validator("taper_cotangent")
+    @classmethod
+    def validate_taper_is_thin(cls, value: float) -> float:
+        bound = 5.0  # round number cotangent limit for thin wedge
+        if -bound < value < bound:
+            msg: str = f"taper_contangent must correspond to thin wedge, of magnitude {bound} or larger"
+            raise ValueError(msg)
+        return value
+
+    @model_validator(mode="after")
+    def _calculate_gradient(self) -> Self:
+        self._taper_gradient = 1.0 / self.taper_cotangent
+        return self
+
+    def motor_position_mm_for_thickness_cm(
+        self, thickness_cm: Annotated[float, Field(ge=0.0)]
+    ) -> float:
+        """Hypothetical motor position at which the wedge thickness would take on the requested value.
+        N.B. absorber thickness are best worked out in cm (science), motor positions prefer mm (controls).
+
+        Args:
+            thickness_cm: (float) non-negative float, zero permitted - the required target thickness in cm
+
+        Returns:
+            motor position in mm where a mathematical model wedge puts the target thickness in the x-ray beam
+        """
+        _line_offset = (
+            self.tip_mm
+        )  # Given that thickness (line x) is zero at this motor position (line y)
+        _line_gradient = (
+            self.taper_cotangent
+        )  # small changes in thickness lead to large changes in position
+        _x = convert_cm_to_mm(thickness_cm)  # do the line calculations in mm
+        return get_straight_line_y(_line_offset, _line_gradient, _x)
+
+    def thickness_cm_at_motor_position_mm(self, motor_position_mm: float) -> float:
+        """Thickness of modelled wedge a given motor position in most useful units.
+        N.B. The wedge geometry is not bounded by lateral position limits, but does reject,
+        negative thickness,
+        unlikely large taper angles (beyond about 11.31 degrees).
+
+        Args:
+            motor_position_mm: (float) unrestricted float, zero permitted - the required target thickness
+
+        Returns:
+            position in mm on the motor axis scale where a mathematical wedge would have the target thickness
+        """
+        _line_offset = (
+            -self.tip_mm * self._taper_gradient
+        )  # not immediately obvious - but after some algebra
+        _line_gradient = self._taper_gradient
+        _thickness_mm = get_straight_line_y(
+            _line_offset, _line_gradient, motor_position_mm
+        )
+        return convert_mm_to_cm(_thickness_mm)

src/dodal/common/general_maths/arithmetic_conversions.py

@@ -1,5 +1,6 @@
 """Provides functions to convert between common units for attenuation."""
 
+from numpy.polynomial import polynomial
 from pydantic import StrictFloat, validate_call
 
 
@@ -105,3 +106,21 @@ def convert_ev_to_kev(energy_ev: StrictFloat) -> float:
         energy_ev (float): a value of kilo-electron volts
     """
     return energy_ev * 1e-3
+
+
+@validate_call
+def get_straight_line_y(
+    line_offset: StrictFloat, line_gradient: StrictFloat, x: StrictFloat
+) -> float:
+    """Convenient internal conversion method which delegates to numpy for a 1st order polynomial (line).
+
+    The traditional y = mx + c to compute y.
+
+    Args:
+        line_offset: the offset, where the modelled line intercepts the abscissa.
+        line_gradient: the gradient of the modelled line
+        x: the input coordinate of the point on the line, from which y is computed.
+    """
+    _coefficients = [line_offset, line_gradient]
+    _y = polynomial.polyval(x, _coefficients)
+    return float(_y)

src/dodal/common/general_maths/material_absorption_maths.py

@@ -1,8 +1,11 @@
+from collections.abc import Callable
 from typing import Annotated
 
+from numpy.polynomial import polynomial
 from pydantic import (
     Field,
     StrictFloat,
+    StrictInt,
     validate_call,
 )
 
@@ -12,18 +15,108 @@
 )
 
 
+class AbsorptionCalculator:
+    """Interface for calculating absorption (due to mass attenuation) per cm as a function of x-ray energy in keV.
+    Reports back a natural log based absorption factor as function of x-ray energy.
+    """
+
+    def __init__(self, _calculation: Callable[[float], float]):
+        self._calculate = _calculation
+
+    def absorption_coefficient_per_cm(
+        self, energy_kev: Annotated[StrictFloat | StrictInt, Field(gt=0)]
+    ) -> float:
+        """Logarithmic contribution to x-ray absorption per cm (depth into absorber) at a particular photon energy.
+
+        Notes:
+        1) The cm as a typical unit length scale is a de facto standard in the field.
+        2a) The absorption coefficient is in factors of e.
+        2b) Barnett absorption units are only used once this coefficient is combined with depth in cm.
+        3a) Negative contributions to a sum have to be permitted so the output,
+            of a specific calculator instance, can be a float of either sign.
+        3b) Sub-classes may implicitly restrict this expectation.
+            Real materials the final absorption will be positive.
+
+        Args:
+            energy_kev (StrictFloat|StrictInt): The individual energy per photon in kilo-electronvolts
+
+        Returns:
+            (float): Model adjustment of attenuation per cm of absorber material depth.
+        """
+        return self._calculate(energy_kev)
+
+
+class CompoundAbsorptionCalculator(AbsorptionCalculator):
+    """Advanced physics model for mass attenuation per cm as a function of x-ray energy in keV.
+    Applies effects from several absorption term calculations.
+
+    Attributes:
+        contributers (list[AbsorptionCalculator])
+    """
+
+    def __init__(self, contributions: list[AbsorptionCalculator]):
+        super().__init__(
+            # lambda k is energy in keV
+            lambda k: sum(c.absorption_coefficient_per_cm(k) for c in contributions)
+        )
+
+
+class PolynomialAbsorptionCorrection(AbsorptionCalculator):
+    """Corrective model for mass attenuation per cm as a function of x-ray energy in keV.
+    Provides terms for correcting a baseline mass attenuation.
+
+    Attributes:
+        corrective_terms (float): Polynomial coefficients to correct the baseline modelled absorption per cm
+    """
+
+    def __init__(self, coefficients_per_cm: list[float]):
+
+        def _calculate_correction(energy_kev: float) -> float:
+            correction = polynomial.polyval(energy_kev, coefficients_per_cm)
+            return float(correction)  # numpy did not specify float as the return type
+
+        super().__init__(_calculate_correction)
+
+
+class SingleRollOffAbsorptionCalculator(AbsorptionCalculator):
+    """Simplest physics model for mass attenuation per cm as a function of x-ray energy in keV.
+    Typically appropriate where one element dominates the absorption,
+    and energies passed in as calculation arguments are above the elements absorption resonance edge.
+
+    Attributes:
+        material_factor_per_cm (StrictFloat): Positive non-zero material constant
+        (hypothetical trend offset equivalent to an absorption per cm at 1 keV)
+        roll_off: negative exponent of energy dependence above the resonant edge.
+    """
+
+    def __init__(
+        self,
+        material_factor_per_cm: Annotated[StrictFloat | StrictInt, Field(gt=0)],
+        roll_off: Annotated[StrictFloat | StrictInt, Field(lt=0)],
+    ):
+        # lambda k is energy in keV
+        super().__init__(
+            lambda k: photon_mass_attenuation_per_unit_length(
+                k, material_factor_per_cm, roll_off
+            )
+        )
+
+
 @validate_call
 def photon_mass_attenuation_per_unit_length(
-    energy_kev: Annotated[StrictFloat, Field(gt=0)],
-    photon_absorption_factor_per_unit_length: Annotated[StrictFloat, Field(gt=0)],
-    energy_dependence_exponent: Annotated[StrictFloat, Field(lt=0)],
+    energy_kev: Annotated[StrictFloat | StrictInt, Field(gt=0)],
+    photon_absorption_factor_per_unit_length: Annotated[
+        StrictFloat | StrictInt, Field(gt=0)
+    ],
+    energy_dependence_exponent: Annotated[StrictFloat | StrictInt, Field(lt=0)],
 ) -> float:
     """Calculates mass attenuation per unit length.
 
+    See for example: https://en.wikipedia.org/wiki/Mass_attenuation_coefficient.
+
     Args:
         energy_kev (StrictFloat): X-ray energy in keV (positive values only).
-        photon_absorption_factor_per_unit_length (StrictFloat): Factors of e,
-            in X-ray flux reduction, per unit depth of absorber (positive values only).
+        photon_absorption_factor_per_unit_length (StrictFloat): Logarithmic constant, scaled in factors of e.
         energy_dependence_exponent (StrictFloat): Roll off of absorption factor (negative values only).
 
     Returns:
@@ -36,8 +129,8 @@ def photon_mass_attenuation_per_unit_length(
 
 @validate_call
 def attenuation_at_depth_cm(
-    depth_cm: Annotated[StrictFloat, Field(ge=0)],
-    absorption_coefficient_per_cm: Annotated[StrictFloat, Field(ge=0)],
+    depth_cm: Annotated[StrictFloat | StrictInt, Field(ge=0)],
+    absorption_coefficient_per_cm: Annotated[StrictFloat | StrictInt, Field(gt=0)],
 ) -> float:
     """Calculates attenuation in Barnett units, where 1000 Bn equivalent to 1/e,
     0Bn to 1 and 2000 Bn to 1/(e^2).
@@ -61,8 +154,8 @@ def attenuation_at_depth_cm(
 
 @validate_call
 def thickness_cm_required_to_attenuate(
-    target_attenuation_bn: Annotated[StrictFloat, Field(ge=0)],
-    absorption_coefficient_per_cm: Annotated[StrictFloat, Field(gt=1.0e-8)],
+    target_attenuation_bn: Annotated[StrictFloat | StrictInt, Field(ge=0)],
+    absorption_coefficient_per_cm: Annotated[StrictFloat | StrictInt, Field(gt=1.0e-8)],
 ) -> float:
     """Calculates material depth in cm.
 
@@ -72,7 +165,6 @@ def thickness_cm_required_to_attenuate(
         absorption_coefficient_per_cm (StrictFloat): Factors of e per cm, reduction in flux.
             N.B. This coefficient is positive (and greater than a lower bound for realism).
 
-
     Raises:
         ValueError: if attenuation is below zero,
             or absorption is less than the lower bound set at the round magnitude 1.0e-8,

src/dodal/common/general_maths/transmission_interconversion.py

@@ -1,15 +1,15 @@
 import math
+from typing import Annotated
 
-from pydantic import StrictFloat, validate_call
+from pydantic import Field, StrictFloat, validate_call
 
 _CANONICAL_BARNETT_CONVERSION = -1.0e3
 _REVERSE_BARNETT_CONVERSION = -1.0e-3
 
 
 @validate_call
 def attenuation_from_natural_log_of_transmission(ln_t: StrictFloat) -> float:
-    """Converts from natural log of transmission fraction into Barnett attenuation units
-    .
+    """Converts from natural log of transmission fraction into Barnett attenuation units.
 
     Args:
         ln_t (StrictFloat): natural log of transmission fraction
@@ -21,7 +21,9 @@ def attenuation_from_natural_log_of_transmission(ln_t: StrictFloat) -> float:
 
 
 @validate_call
-def attenuation_from_transmission(transmission_as_fraction: StrictFloat) -> float:
+def attenuation_from_transmission(
+    transmission_as_fraction: Annotated[StrictFloat, Field(ge=0, le=1)],
+) -> float:
     """Converts from transmission fraction into Barnett attenuation units.
 
     Args:
@@ -35,25 +37,28 @@ def attenuation_from_transmission(transmission_as_fraction: StrictFloat) -> floa
 
 
 @validate_call
-def natural_log_of_transmission_from_attenuation(attenuation_bn: StrictFloat) -> float:
-    """Converts from Barnett attenuation units into natural log of transmission fraction
-    .
+def natural_log_of_transmission_from_attenuation(
+    attenuation_bn: Annotated[StrictFloat, Field(ge=0)],
+) -> float:
+    """Converts from Barnett attenuation units into natural log of transmission fraction.
 
     Args:
-        attenuation_bn (StrictFloat): Barnett attenuation units
+        attenuation_bn (StrictFloat): Barnett attenuation units - zero or a positive number
 
     Returns:
-        (float): natural log of transmission fraction
+        (float): natural log of transmission fraction - zero or, more usually, a negative number
     """
     return _REVERSE_BARNETT_CONVERSION * attenuation_bn
 
 
 @validate_call
-def transmission_from_attenutation(attenuation_bn: StrictFloat) -> float:
+def transmission_from_attenutation(
+    attenuation_bn: Annotated[StrictFloat, Field(ge=0)],
+) -> float:
     """Converts from Barnett attenuation units into transmission fraction.
 
     Args:
-        attenuation_bn (StrictFloat): Barnett attenuation units
+        attenuation_bn (StrictFloat): Barnett attenuation units - zero or a positive number
 
     Returns:
         (float): transmission fraction