Replicate Matlab functions

.github/workflows/pylint.yml

@@ -11,7 +11,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.8", "3.9", "3.10"]
+        python-version: ["3.11"]
     steps:
     - uses: actions/checkout@v3
     - name: Set up Python ${{ matrix.python-version }}

pyproject.toml

@@ -8,17 +8,19 @@ description = """\
     Zenodo, Jul. 01, 2022. doi: 10.5281/zenodo.6025377. \
     """
 requires-python = ">=3.8"
-keywords = ["python", "nasa", "desci"]
+keywords = ["python", "nasa", "desci", "simulation", "lunar"]
 authors = [
     {name = "Philip Linden", email = "phil@moondao.com"},
     {name = "Daniel R. Cremons"},
 ]
-license = "CC-BY 4.0"
+license = {text = "CC-BY 4.0"}
 dependencies = [
     "matplotlib",
+    "numpy",
+    "pint",
+    "scipy",
 ]
 
 [build-system]
 requires = ["setuptools"]
 build-backend = "setuptools.build_meta"
-

repro/data_io.py

@@ -0,0 +1,2 @@
+import os
+import numpy as np

repro/main.py

@@ -1,9 +1,17 @@
-import .hapke
+from matplotlib import pyplot as plt
+import numpy as np
+from pprint import pprint as print
 
-
-class Hydrated_Regolith_Spectra_Generator_and_Retrieval:
-    pass
+from simulator import Hapke, get_sample_grid
 
 
 if __name__ == '__main__':
-    pass
+    WLS = get_sample_grid()
+    Refl = np.array([np.random.randn()*x for x in WLS])
+    R = Hapke(Refl, WLS)
+    print(WLS)
+    print(Refl)
+    print(R)
+    plt.figure()
+    plt.plot(WLS,Refl,WLS,R)
+    plt.show()

repro/simulator.py

@@ -0,0 +1,161 @@
+from functools import partial
+from typing import Any, Callable
+import numpy as np
+from scipy.optimize import fmin as fminsearch
+
+
+def get_sample_grid(
+    min_wavelength: float = 1, max_wavelength: float = 6, spacing: float = 0.005
+):
+    """Create a linear spectral simulation grid.
+
+    Defines the discrete-space spectrum to be simulated. All units are in microns.
+
+    Args:
+        min_wavelength (float, optional): Lowest simulated wavelength. Defaults to 1 um.
+        max_wavelength (float, optional): Highest simulated wavelength. Defaults to 6 um.
+        spacing (float, optional): Spectral distance between grid samples. Defaults to 0.005 um.
+
+    Returns:
+        array: Simulation sample grid
+
+    """
+    start = min_wavelength
+    stop = max_wavelength
+    step = spacing
+    grid = np.linspace(start, stop, int((stop - start) / step + 1))
+    return grid
+
+
+def ordinary_least_squares(x: Any, y: Callable, yx: Any):
+    '''Ordinary Least Squares Function.
+
+    y (Callable): The estimator function.
+    x (Any): The argument to y.
+    yx (Any): The observation. Must be the same type as the result of y.
+    '''
+    return sum((y(x) - yx) ** 2)
+
+
+def angular_width(filling_factor=0.41):
+    '''Angular width parameter, see Equation 3.
+
+    Args:
+        filling_factor (float, optional): Filling factor. Defaults to 0.41 for the
+            lunar regolith (Bowell et al., 1989).
+    '''
+    return (-3 / 8) * np.log(1 - filling_factor)
+
+
+def backscattering(h, g):
+    '''Backscattering function, see Equation 2.
+
+    This describes the opposition effect.
+
+    Args:
+        h (float): angular width parameter.
+        g (float): phase angle in radians.
+    '''
+    return 1 / (1 + (1 / h) * np.tan(g / 2))
+
+
+def R(SSA, P, mu, mu0, B):
+    """Bidirectional reflectance, see Equation 1.
+
+        R = (ω/4) (μ₀ / (μ + μ₀)) {(1 + B)P + H(ω)H₀(ω) - 1}
+
+    This equation is HUGE so I broke it down into smaller terms for code readability.    
+        let:
+            R = r1 * r2 * (r3 + r4 - 1)
+
+        where:
+            r1 = (ω/4)
+            r2 = (μ₀ / (μ + μ₀))
+            r3 = (1 + B)P
+            r4 = H(ω) * H₀(ω)
+
+    Equivalent to src/Hapke_Lidar_R_function.m
+
+    Args:
+        SSA (_type_): single-scattering albedo, aka ω
+        P (float): scattering phase function
+        mu (float): cosine of emission angle
+        mu0 (float): cosine of incident angle
+        B (float): backscattering function
+
+    Returns:
+        _type_: _description_
+    """
+
+    def Hfunc(SSA, mu, mu0):
+        '''Ambartsumian-Chandrasekhar H functions.
+
+        Computed using the approximation from equation 8.57 from Hapke (2012).
+
+        Args:
+            SSA (float): single-scattering albedo, aka ω
+            mu (float): cosine of emission angle
+            mu0 (float): coside of incident angle
+        '''
+        gamma = np.sqrt(1 - SSA)
+        r0 = (1 - gamma) / (1 + gamma)
+
+        h1 = np.log((1 + mu0) / mu0)
+        h2 = (r0 + ((1 - 2 * r0) * mu0) / 2) * h1
+
+        H = (1 - SSA * mu0 * h2) ** -1
+
+        h3 = np.log((1 + mu) / mu)
+        h4 = (1 - 2 * r0 * mu)
+        h5 = r0 + h3 * (h4 / 2)
+
+        H0 = (1 - SSA * mu * h5) ** -1
+        return H, H0 
+
+    H, H0 = Hfunc(SSA, mu, mu0)
+
+    r1 = SSA / 4
+    r2 = mu0 / (mu0 + mu)
+    r3 = (1 + B) * P
+    r4 = H * H0
+    return r1 * r2 * (r3 + r4 - 1)
+
+
+def Hapke(Refl, WLS, P=0.15, emission_angle=0, incident_angle=30, phase_angle=30, filling_factor=0.41):
+    """Convert reflectance spectrum to single-scattering albedo (SSA)
+
+    Uses scipy.optimize.fmin (equivalent to MATLAB fminsearch) to minimize 
+    ordinary least squares distance between SSA obtained from the supplied
+    reflectance, R, and the SSA from the estimated reflectance at each
+    sample point in WLS.
+
+    Hapke, B. (2012). Theory of reflectance and emittance spectroscopy (2nd ed.).
+        Cambridge University Press.
+
+    Equivalent to src/Hapke_Lidar_SSA_function.m
+    Default parameter values replicate src/Hapke_Inverse_Function_Passive.m
+
+    Args:
+        R (array[float]): Bidirectional reflectance, see Equation 1.
+        WLS: simulation sample grid?
+        p (float, optional): Scattering phase function. Defaults to 0.15 for ansiotropic
+            scattering on the modeled mean particle phase function for lunar soil
+            (Goguen et al., 2010).
+        emission_angle (float, optional): Emission angle in degrees. Defaults to 0.
+        incident_angle (float, optional): Incident angle in degrees. Defaults to 30.
+        phase_angle (float, optional): Phase angle in degrees. Defaults to 30.
+        filling_factor (float, optional: Particle filling factor. Defaults to 0.41.
+    """
+    mu = np.cos(np.deg2rad(emission_angle))
+    mu0 = np.cos(np.deg2rad(incident_angle))
+    g = np.deg2rad(phase_angle)
+    h = angular_width(filling_factor=filling_factor)
+    B = backscattering(h, g)
+
+    w = []
+    for m, x in zip(Refl, WLS):
+        w0 = 0.5  # initial guess, ω₀
+        y = partial(R, P=P, mu=mu, mu0=mu0, B=B) # turn R() into the form y=f(x)
+        OLS = partial(ordinary_least_squares, y=y, yx=m) # turn least squares into the form y=f(x)
+        w.append(fminsearch(OLS, w0, args=x, disp=False, maxiter=10_000, maxfun=10_000, ftol=1e-7, xtol=1e-7))
+    return np.concatenate(w)

repro/spectra.py

@@ -0,0 +1,24 @@
+from dataclasses import dataclass
+from typing import Tuple
+
+
+@dataclass
+class Endmember:
+    name: str
+    abundance_range: Tuple[float, float]
+    max_abundance: float
+    density: float
+    grain_size_range: Tuple[float, float]
+
+
+class HydratedMorbGlass(Endmember):
+    ''' Hydrated mid-ocean-ridge basalt (MORB) glass '''
+    sample_label = ''
+    spectrum_label = 'MORB D38A'
+    density = 2.8
+    grain_size = 69e-6
+
+
+class Regolith(Endmember):
+    density = 1.8
+    grain_size = 32e-6