add FluorescenceData object

lib/imaging.py

@@ -11,6 +11,9 @@ import toml
 import PIL.Image
 import lib.plotdefs as pd
 
+H = 6.626070040e-34
+C = 2.99792458e+8
+
 DEF_DX = 3.45 # pixel size for Flir Grasshopper [um]
 
 def _fresh_filename(path: pathlib.Path, overwrite: bool=False) -> pathlib.Path:
@@ -23,7 +26,7 @@ def _fresh_filename(path: pathlib.Path, overwrite: bool=False) -> pathlib.Path:
         print(f"Write instead to {_path}")
     return _path
 
-class AbsorptionData:
+class ImagingData:
     def __init__(self, name: str, arrays: dict[str, np.ndarray],
             config: dict[str, ...], comments: str=None,
             results: dict[str, ...]=None):
@@ -124,6 +127,10 @@ class AbsorptionData:
         if printflag: print("  Done.")
         return self
 
+    def compute_results(self, dA=(2*3.45)**2, printflag: bool=True):
+        raise NotImplementedError
+
+class AbsorptionData(ImagingData):
     def compute_results(self, dA=(2*3.45)**2, printflag: bool=True):
         _dA = self.config.get("bin_size", 1)**2 \
                 * (DEF_DX**2 if dA is None else dA)
@@ -144,6 +151,40 @@ def compute_od(shadow: np.ndarray, bright: np.ndarray, dark: np.ndarray):
     OD[~np.isfinite(OD)] = 0.0
     return OD
 
+class FluorescenceData(ImagingData):
+    def compute_results(self, dA=(2*3.45)**2, printflag: bool=True):
+        _dA = self.config.get("bin_size", 1)**2 \
+                * (DEF_DX**2 if dA is None else dA)
+        if printflag: print("Compute results")
+        N = compute_mot_number(
+            self.arrays["image"],
+            0.76,
+            self.config["gain"],
+            self.config["exposure_time"],
+            (0.0127**2 * np.pi) / (4 * np.pi * 0.1**2)
+            62e6 * 2 * np.pi,
+            0.1
+        )
+        sx, sy = fit_gaussian(self.arrays["image"], _dA)
+        self.results = dict() if self.results is None else self.results
+        self.results.update({"N": N, "sx": sx, "sy": sy})
+        return self
+
+def compute_mot_number(image, QE, gain, exposure_time, solid_angle,
+        detuning, intensity_parameter):
+    if image.max() >= 2**16:
+        print("WARNING: image may contain clipping")
+    electron_rate = im.sum() * 16
+    photon_rate = electron_rate / gain / QE / solid_angle / exposure_time
+
+    Y = 29.1e6 * 2 * np.pi # transition linewidth; s^-1
+    l = 399e-9 # transition wavelength; m
+    I_sat = (Y * H * C * np.pi) / (3 * l**3)
+    scatter_rate = lambda D, I: \
+            (Y / 2) * I / (1 + 4 * (D / Y)**2 + I)
+    N = photon_rate / scatter_rate(detuning, intensity_parameter)
+    return N
+
 def _fit_gaussian(A, dA):
     """
     Assumes uniform, square dA.