From 541d3c8cd654e9077d8b7e54e585d5ac04f4051f Mon Sep 17 00:00:00 2001
From: Ashesh <ashesh276@gmail.com>
Date: Sun, 30 Jun 2024 11:44:39 +0200
Subject: [PATCH 1/3] doc

---
 docs/emission.md                 | 47 ++++++++++++++++++++++++++++++++
 src/microsim/schema/_emission.py |  4 +++
 2 files changed, 51 insertions(+)
 create mode 100644 docs/emission.md

diff --git a/docs/emission.md b/docs/emission.md
new file mode 100644
index 0000000..a25ab4f
--- /dev/null
+++ b/docs/emission.md
@@ -0,0 +1,47 @@
+# Overview
+This document outlines the mechanism of how the excitation and emission of the fluoropore is simulated in this package. The code is present at [emission.py](../src/microsim/schema/_emission.py). 
+
+    Ref: https://en.wikipedia.org/wiki/Absorption_cross_section,
+        "https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Time_Dependent_Quantum_Mechanics_and_Spectroscopy_%28Tokmakoff%29/07:_Interaction_of_Light_and_Matter/7.05:_Absorption_Cross-Sections"
+## Excitation
+  The task is to estimate the number of photons absorbed by the fluorophores for each `1nm` wavelength interval. For this, following steps are followed:
+
+1. Compute the irradiance for every `1nm` wavelength interval of the light source. `irradiance = ex_filter_spectrum * light_power`.  
+2. Calculate the absorption cross section given the extinction coefficient. See [section below](#absorption-cross-section) for more details. This is inturn used to calculate the absorbtion cross section for every `1nm` wavelength interval of the fluorophore excitation spectrum.
+3. Now that we have the abrobption cross section and irradiance, we calculate the power absorbed by the fluorophores for every `1nm` wavelength interval.
+4. Finally, we compute the absorbed photon count by simply dividing the power absorbed by the energy of a single photon.
+
+## Emission 
+1. We scale up the emission spectrum by a factor to account for light power, extinction coefficient and overlap between fluorophore excitation spectrum and input light spectrum. `scaling_factor = ex_rate.integral() / fluor_em_spectrum.integral()`
+2. Multiply the scaled emission spectrum by the quantum yield to get the number of emission events. `fluor_em_rate = fluor_em_rate * fluor.quantum_yield`
+3. Apply the emission filter(s) to the emission spectrum.
+
+ 
+### Absorption Cross Section
+Given the extinction coefficient, we calculate the absorption cross section using the following formula:
+
+> cross section = 1000 * ec * np.log(10) / $N_{A}$,
+
+where $N_{A}$ is the Avogadro number. The factor `1000` appears because we convert the concentration from `mol/L` to `mol/cm3`. The natural logarithm is converted to base 10 logarithm by multiplying with `np.log(10)`.
+
+
+$$A = ec \times C $$
+where $A$ is abosorbance cross section and $C$ is sample concentration in mol/L.
+
+$$C~(mol/L) = \frac{N}{N_A} \times 1000~({cm}^{-3})$$
+where $N$ is number of molecules in sample per cm3 and $N_A$ is the Avogadro number.
+Finally, base 10 is used commonly instead of the natural base and assuming that, we have a $ln(10)$ factor. 
+So, expression for A becomes 
+$$A = ec \times ln(10) \times 1000 \times \frac{N}{N_A}$$
+
+$\sigma =1-10^{-\frac{\frac{\text{absorption}}{N_{\text{A}}}}{}}$
+
+Simplifying this expression, we get:
+<!-- $\sigma =\frac{\ln(10)\times 10^3}{N_{\text{A}}}\times \epsilon$ -->
+This is the expression for the absorption cross section \sigma in terms of the extinction coefficient \epsilon .
+
+Power absorbed by the fluorophores is given by $A \times E $, where $E$ is the irradiance of the light source.
+At this point we get the number of photons absorbed by the fluorophores by dividing the power absorbed by the energy of a single photon.
+$$N_{abs} = \frac{A \times E \times \lambda}{h \times c}$$
+where $h$ is the Planck constant, $c$ is the speed of light and $\lambda$ is the wavelength.
+
diff --git a/src/microsim/schema/_emission.py b/src/microsim/schema/_emission.py
index 89c88be..b5164b0 100644
--- a/src/microsim/schema/_emission.py
+++ b/src/microsim/schema/_emission.py
@@ -76,7 +76,11 @@ def get_excitation_rate(
         )
 
     # TODO: derive light power from model
+    # TODO: Does it make sense to normalize the ex_filter_spectrum? In that case, we
+    # can say that light power is distributed over the spectrum.
     irradiance = ex_filter_spectrum * _ensure_quantity(light_power, "W/cm^2")
+    # cross section has units of cm2
+    # ec_to_cross_section is a number with units of cm2
     cross_section = fluor_ex_spectrum * ec_to_cross_section(ext_coeff)
     power_absorbed = cross_section * irradiance
     excitation_rate = power_absorbed / energy_per_photon(power_absorbed.wavelength)

From c8630cc3102fc92c0e1fb3ede522363fe066cb6e Mon Sep 17 00:00:00 2001
From: Ashesh <ashesh276@gmail.com>
Date: Sun, 30 Jun 2024 12:28:08 +0200
Subject: [PATCH 2/3] bug in cross section fixed

---
 docs/{emission.md => ex_em_methodology.md} | 27 +++-------------------
 src/microsim/schema/_emission.py           | 18 ++++-----------
 2 files changed, 8 insertions(+), 37 deletions(-)
 rename docs/{emission.md => ex_em_methodology.md} (54%)

diff --git a/docs/emission.md b/docs/ex_em_methodology.md
similarity index 54%
rename from docs/emission.md
rename to docs/ex_em_methodology.md
index a25ab4f..ad5c5fc 100644
--- a/docs/emission.md
+++ b/docs/ex_em_methodology.md
@@ -1,8 +1,6 @@
-# Overview
+# Excitation and Emission Methodology
 This document outlines the mechanism of how the excitation and emission of the fluoropore is simulated in this package. The code is present at [emission.py](../src/microsim/schema/_emission.py). 
 
-    Ref: https://en.wikipedia.org/wiki/Absorption_cross_section,
-        "https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Time_Dependent_Quantum_Mechanics_and_Spectroscopy_%28Tokmakoff%29/07:_Interaction_of_Light_and_Matter/7.05:_Absorption_Cross-Sections"
 ## Excitation
   The task is to estimate the number of photons absorbed by the fluorophores for each `1nm` wavelength interval. For this, following steps are followed:
 
@@ -12,6 +10,7 @@ This document outlines the mechanism of how the excitation and emission of the f
 4. Finally, we compute the absorbed photon count by simply dividing the power absorbed by the energy of a single photon.
 
 ## Emission 
+  For emission, we need to estimate the emission flux for every `1nm` wavelength interval which is done as follows:
 1. We scale up the emission spectrum by a factor to account for light power, extinction coefficient and overlap between fluorophore excitation spectrum and input light spectrum. `scaling_factor = ex_rate.integral() / fluor_em_spectrum.integral()`
 2. Multiply the scaled emission spectrum by the quantum yield to get the number of emission events. `fluor_em_rate = fluor_em_rate * fluor.quantum_yield`
 3. Apply the emission filter(s) to the emission spectrum.
@@ -22,26 +21,6 @@ Given the extinction coefficient, we calculate the absorption cross section usin
 
 > cross section = 1000 * ec * np.log(10) / $N_{A}$,
 
-where $N_{A}$ is the Avogadro number. The factor `1000` appears because we convert the concentration from `mol/L` to `mol/cm3`. The natural logarithm is converted to base 10 logarithm by multiplying with `np.log(10)`.
+where $N_{A}$ is the Avogadro number. The factor `1000` appears because we convert the concentration from `mol/L` to `mol/cm3`. Note that this conversion is taken care by [pint](https://pint.readthedocs.io/en/stable/) and therefore is not explicit in the code. The natural logarithm is converted to base 10 logarithm by multiplying with `np.log(10)`. [Reference 1](https://en.wikipedia.org/wiki/Absorption_cross_section), [Reference 2](https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Time_Dependent_Quantum_Mechanics_and_Spectroscopy_%28Tokmakoff%29/07:_Interaction_of_Light_and_Matter/7.05:_Absorption_Cross-Sections)
 
 
-$$A = ec \times C $$
-where $A$ is abosorbance cross section and $C$ is sample concentration in mol/L.
-
-$$C~(mol/L) = \frac{N}{N_A} \times 1000~({cm}^{-3})$$
-where $N$ is number of molecules in sample per cm3 and $N_A$ is the Avogadro number.
-Finally, base 10 is used commonly instead of the natural base and assuming that, we have a $ln(10)$ factor. 
-So, expression for A becomes 
-$$A = ec \times ln(10) \times 1000 \times \frac{N}{N_A}$$
-
-$\sigma =1-10^{-\frac{\frac{\text{absorption}}{N_{\text{A}}}}{}}$
-
-Simplifying this expression, we get:
-<!-- $\sigma =\frac{\ln(10)\times 10^3}{N_{\text{A}}}\times \epsilon$ -->
-This is the expression for the absorption cross section \sigma in terms of the extinction coefficient \epsilon .
-
-Power absorbed by the fluorophores is given by $A \times E $, where $E$ is the irradiance of the light source.
-At this point we get the number of photons absorbed by the fluorophores by dividing the power absorbed by the energy of a single photon.
-$$N_{abs} = \frac{A \times E \times \lambda}{h \times c}$$
-where $h$ is the Planck constant, $c$ is the speed of light and $\lambda$ is the wavelength.
-
diff --git a/src/microsim/schema/_emission.py b/src/microsim/schema/_emission.py
index b5164b0..0c38f04 100644
--- a/src/microsim/schema/_emission.py
+++ b/src/microsim/schema/_emission.py
@@ -20,11 +20,12 @@
 C = c * ureg.meter / ureg.second
 
 
-def _ensure_quantity(value: Any, units: str) -> pint.Quantity:
+def _ensure_quantity(value: Any, units: str, strict:bool=False) -> pint.Quantity:
     """Helper function to ensure that a value is a pint Quantity with `units`."""
     if isinstance(value, pint.Quantity):
         quant = value
     else:
+        assert not strict, f"Expected a pint.Quantity with units {units}, got {value}"
         quant = ureg.Quantity(value)
     if quant.dimensionless:
         quant *= ureg.Quantity(units)
@@ -35,13 +36,8 @@ def _ensure_quantity(value: Any, units: str) -> pint.Quantity:
 
 
 def ec_to_cross_section(ec: Any) -> pint.Quantity:
-    """Gives cross section in cm^2 from extinction coefficient in M^-1 * cm^-1."""
-    ec = _ensure_quantity(ec, "cm^2/mol")
-    # x1000?
-    # this came from calculations elsewhere, and looking at wikipedia
-    # and looking at Nathan Shaner's code
-    # need to double check whether it's still correct with our units
-    ec = ec * 1000
+    """Gives cross section in cm^2 from extinction coefficient in cm^2 * mol^-1."""
+    ec = _ensure_quantity(ec, "cm^2/mol", strict=True)
     return (ec * np.log(10) / AVOGADRO).to("cm^2")  # type: ignore [no-any-return]
 
 
@@ -68,7 +64,7 @@ def get_excitation_rate(
         raise NotImplementedError("Fluorophore has no excitation spectrum.")
 
     if (ext_coeff := fluor.extinction_coefficient) is None:
-        ext_coeff = _ensure_quantity(55000, "cm^2/mol")
+        ext_coeff = _ensure_quantity(55000, "cm^-1/M")
         warnings.warn(
             "No extinction coefficient provided for fluorophore, "
             "using 55,000 M^-1 * cm^-1.",
@@ -76,11 +72,7 @@ def get_excitation_rate(
         )
 
     # TODO: derive light power from model
-    # TODO: Does it make sense to normalize the ex_filter_spectrum? In that case, we
-    # can say that light power is distributed over the spectrum.
     irradiance = ex_filter_spectrum * _ensure_quantity(light_power, "W/cm^2")
-    # cross section has units of cm2
-    # ec_to_cross_section is a number with units of cm2
     cross_section = fluor_ex_spectrum * ec_to_cross_section(ext_coeff)
     power_absorbed = cross_section * irradiance
     excitation_rate = power_absorbed / energy_per_photon(power_absorbed.wavelength)

From 65b21c790644b5fc4726bb34ff8ccd7c14b7989a Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Sun, 30 Jun 2024 10:54:06 +0000
Subject: [PATCH 3/3] style(pre-commit.ci): auto fixes [...]

---
 src/microsim/schema/_emission.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/microsim/schema/_emission.py b/src/microsim/schema/_emission.py
index 0c38f04..164542b 100644
--- a/src/microsim/schema/_emission.py
+++ b/src/microsim/schema/_emission.py
@@ -20,7 +20,7 @@
 C = c * ureg.meter / ureg.second
 
 
-def _ensure_quantity(value: Any, units: str, strict:bool=False) -> pint.Quantity:
+def _ensure_quantity(value: Any, units: str, strict: bool = False) -> pint.Quantity:
     """Helper function to ensure that a value is a pint Quantity with `units`."""
     if isinstance(value, pint.Quantity):
         quant = value