create od2conc function and update beer_lambert

src/cedalion/nirs.py

@@ -51,34 +51,53 @@ def channel_distances(amplitudes: xr.DataArray, geo3d: xr.DataArray):
 
 
 def int2od(amplitudes: xr.DataArray):
-    """Calculate optical density from intensity data."""
+    """Calculate optical density from intensity amplitude  data."""
     od = - np.log( amplitudes / amplitudes.mean("time") )
     return od
 
 
-def beer_lambert(
-    amplitudes: xr.DataArray,
+def od2conc(
+    od: xr.DataArray,
     geo3d: xr.DataArray,
     dpf: xr.DataArray,
     spectrum: str = "prahl",
 ):
-    """Calculate concentration changes from amplitude data."""
-    validators.has_channel(amplitudes)
-    validators.has_wavelengths(amplitudes)
+    """Calculate concentration changes from optical density data."""
+    validators.has_channel(od)
+    validators.has_wavelengths(od)
     validators.has_wavelengths(dpf)
     validators.has_positions(geo3d, npos=3)
 
-    E = get_extinction_coefficients(spectrum, amplitudes.wavelength)
+    E = get_extinction_coefficients(spectrum, od.wavelength)
 
     Einv = xrutils.pinv(E)
 
-    dists = channel_distances(amplitudes, geo3d)
+    dists = channel_distances(od, geo3d)
     dists = dists.pint.to("mm")
 
-    optical_density = -np.log(amplitudes / amplitudes.mean("time"))
     # conc = Einv @ (optical_density / ( dists * dpf))
-    conc = xr.dot(Einv, optical_density / (dists * dpf), dims=["wavelength"])
+    conc = xr.dot(Einv, od / (dists * dpf), dims=["wavelength"])
     conc = conc.pint.to("micromolar")
     conc = conc.rename("concentrations")
 
     return conc
+
+
+def beer_lambert(
+    amplitudes: xr.DataArray,
+    geo3d: xr.DataArray,
+    dpf: xr.DataArray,
+    spectrum: str = "prahl",
+):
+    """Calculate concentration changes from amplitude data."""
+    validators.has_channel(amplitudes)
+    validators.has_wavelengths(amplitudes)
+    validators.has_wavelengths(dpf)
+    validators.has_positions(geo3d, npos=3)
+
+    # calculate optical densities
+    od = int2od(amplitudes)
+    # calculate concentrations
+    conc = od2conc(od, geo3d, dpf, spectrum)
+
+    return conc