Photobleaching calculations

photobleach.py

import math
import pint

AVOGADRO = pint.Quantity("6.02214076e23 /mol")
PLANCK = pint.Quantity("6.62607015e-34 J*s")
C = pint.Quantity("299792458 m/s")


def _ensure_quantity(value: str | float | pint.Quantity, units: str) -> pint.Quantity:
    quant = pint.Quantity(value)
    if quant.dimensionless:
        quant *= pint.Quantity(units)
    _u = pint.Quantity(units).units
    if not quant.is_compatible_with(_u):
        raise pint.DimensionalityError(quant.units, _u)
    return quant


def ec_to_cross_section(ec: str | float | pint.Quantity) -> pint.Quantity:
    """Gives cross section in cm^2 from extinction coefficient in M^-1 * cm^-1."""
    ec = _ensure_quantity(ec, "cm^2/mol")  # CHECK ME ... x1000?
    return ec * math.log(10) * 10**3 / AVOGADRO


def energy_per_photon(wavelength: str | float | pint.Quantity) -> pint.Quantity:
    """Converts a wavelength to energy per photon."""
    wavelength = _ensure_quantity(wavelength, "1nm")
    return PLANCK * C / wavelength


def watts_to_photons_per_second(
    power: str | float | pint.Quantity, wavelength: str | float | pint.Quantity
) -> pint.Quantity:
    """Converts a power to photons per second at a given wavelength.

    If not provided as pint Quantities, power is assumed to be in Watts and
    wavelength is assumed to be in nanometers.
    """
    power = _ensure_quantity(power, "1W")
    return (power / energy_per_photon(wavelength)).to("1/s")


def photons_per_fluorophore_per_second(
    irradiance: str | float | pint.Quantity,
    wavelength: str | float | pint.Quantity,
    extinction_coefficient: str | float | pint.Quantity,
) -> pint.Quantity:
    irradiance = _ensure_quantity(irradiance, "1W/cm^2")
    E_photon = energy_per_photon(wavelength)
    cross_section = ec_to_cross_section(extinction_coefficient)
    # Calculate the number of photons hitting the fluorophore per second
    return (cross_section * irradiance / E_photon).to("1/s")


def saturating_irradiance(
    extinction_coefficient: str | float | pint.Quantity,
    wavelength: str | float | pint.Quantity,
    lifetime: str | float | pint.Quantity,
):
    """Determine the irradiance required to saturate a fluorophore."""
    E_photon = energy_per_photon(wavelength)
    cross_section = ec_to_cross_section(extinction_coefficient)
    lifetime = _ensure_quantity(lifetime, "1ns")
    return E_photon / (cross_section * lifetime)


if __name__ == "__main__":
    EC = 56000
    WAVE = 488
    LIFETIME = 2.6
    print(ec_to_cross_section(EC))
    # print(watts_to_photons_per_second(10, WAVE))
    print(photons_per_fluorophore_per_second(177.43564667323508, WAVE, EC).to("1/ns"))
    print(saturating_irradiance(EC, WAVE, LIFETIME).to("kW/cm^2"))