2023 11 04

repro/Hapke_Inverse_Function_Passive.m

@@ -5,12 +5,12 @@
 h=(-3/8)*log(1-0.41);
 B=1/(1+(1/h)*tand(g/2));
 for m=1:numel(WLS)
-Refl=R(m);
-y=@(SSA,x) (SSA./(4*(mu0+mu))).*(p.*(1+B)+(((1-SSA.*mu0.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu0)/2)*log((1+mu0)./mu0))).^-1).*((1-SSA.*mu.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu)/2)*log((1+mu)./mu))).^-1))-1);
-x=WLS(m);
-yx=Refl;
-w0=0.5; %Initial Guesses here
-OLS=@(SSA) sum((y(SSA,x)-yx).^2); %Ordinary Least Squares Function
-opts=optimset('MaxIter',10000,'MaxFunEvals',10000,'TolFun',1E-7,'TolX',1E-7);
-SSA(m,1) = fminsearch(OLS, w0, opts);      % Use ‘fminsearch’ to minimise the ‘OLS’ function
+    Refl=R(m);
+    y=@(SSA,x) (SSA./(4*(mu0+mu))).*(p.*(1+B)+(((1-SSA.*mu0.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu0)/2)*log((1+mu0)./mu0))).^-1).*((1-SSA.*mu.*(((1-sqrt(1-SSA))./(1+sqrt(1-SSA)))+((1-2.*((1-sqrt(1-SSA))./(1+sqrt(1-SSA))).*mu)/2)*log((1+mu)./mu))).^-1))-1);
+    x=WLS(m);
+    yx=Refl;
+    w0=0.5; %Initial Guesses here
+    OLS=@(SSA) sum((y(SSA,x)-yx).^2); %Ordinary Least Squares Function
+    opts=optimset('MaxIter',10000,'MaxFunEvals',10000,'TolFun',1E-7,'TolX',1E-7);
+    SSA(m,1) = fminsearch(OLS, w0, opts);      % Use ‘fminsearch’ to minimise the ‘OLS’ function
 end

repro/Hydrated_Regolith_Spectra_Generator_and_Retrieval.m

@@ -52,43 +52,43 @@
 clear MORB_D38A2
 %Perform Normalization of MORB spectra at 650, 700, 750, 800 C
 for i=1:4
-MORB_D38A2=fullfile(MORB_D38A(i).folder, MORB_D38A(i).name); %start from 650 C spectrum
-MORB_D38A2=readtable(MORB_D38A2);
-MORB_D38A2=table2array(MORB_D38A2);
-MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns
-MORB_D38A2(:,2)=MORB_D38A2(:,2)/100; %convert from percent to decimal
-MORB_D38A2=flipud(MORB_D38A2);
-%Normalize
-NormFactor=MORB_D38A2(12752,2)/NormalizationValue; %normalize at 2.6 micron
-MORB_D38A2(:,2)=MORB_D38A2(:,2)./NormFactor;
-%Stitch with low wavelength spectrum
-Morb_D38A_LowLam2=Morb_D38A_LowLam;
-Morb_D38A_LowLam2(:,2)=Morb_D38A_LowLam(:,2)*(MORB_D38A2(6901,2)/0.15);
-Morb_D38A_LowLam2=Morb_D38A_LowLam2(1:15,:);
-MORB_D38A2=cat(1,Morb_D38A_LowLam2,MORB_D38A2(6902:end,:));
-% Interpolate data to 5 nm spacing between 1 and 4 microns
-InterpMorb=interp1(MORB_D38A2(:,1),MORB_D38A2(:,2),WLS);
-if i>1 % Use 650 C spectrum below 2.6 microns to isolate 3 micron feature changes
-    InterpMorb(1:321)=RMORB(1:321,1);
-else
-end
-RMORB(:,i)=InterpMorb;
-%Use Hapke model to convert from laboratory reflectance to single
-%scattering albedo using reflectance and scattering asymmetry
-% factor (p) of 0.81 (see manuscript for details on Hapke model)
-SSAMORB(:,i)=Hapke_Inverse_Function_Passive(RMORB(:,i),0.81,WLS);
+    MORB_D38A2=fullfile(MORB_D38A(i).folder, MORB_D38A(i).name); %start from 650 C spectrum
+    MORB_D38A2=readtable(MORB_D38A2);
+    MORB_D38A2=table2array(MORB_D38A2);
+    MORB_D38A2(:,1)=1E4./MORB_D38A2(:,1); %convert from cm^-1 to microns
+    MORB_D38A2(:,2)=MORB_D38A2(:,2)/100; %convert from percent to decimal
+    MORB_D38A2=flipud(MORB_D38A2);
+    %Normalize
+    NormFactor=MORB_D38A2(12752,2)/NormalizationValue; %normalize at 2.6 micron
+    MORB_D38A2(:,2)=MORB_D38A2(:,2)./NormFactor;
+    %Stitch with low wavelength spectrum
+    Morb_D38A_LowLam2=Morb_D38A_LowLam;
+    Morb_D38A_LowLam2(:,2)=Morb_D38A_LowLam(:,2)*(MORB_D38A2(6901,2)/0.15);
+    Morb_D38A_LowLam2=Morb_D38A_LowLam2(1:15,:);
+    MORB_D38A2=cat(1,Morb_D38A_LowLam2,MORB_D38A2(6902:end,:));
+    % Interpolate data to 5 nm spacing between 1 and 4 microns
+    InterpMorb=interp1(MORB_D38A2(:,1),MORB_D38A2(:,2),WLS);
+    if i>1 % Use 650 C spectrum below 2.6 microns to isolate 3 micron feature changes
+        InterpMorb(1:321)=RMORB(1:321,1);
+    else
+    end
+    RMORB(:,i)=InterpMorb;
+    %Use Hapke model to convert from laboratory reflectance to single
+    %scattering albedo using reflectance and scattering asymmetry
+    % factor (p) of 0.81 (see manuscript for details on Hapke model)
+    SSAMORB(:,i)=Hapke_Inverse_Function_Passive(RMORB(:,i),0.81,WLS);
 end
 %% Interpolation of MORB Spectra
-%Set range and resolution for interpolation 
+%Set range and resolution for interpolation
 %(standard is 1 ppm spacing from 0 to 1666 ppm) With 30% maximum MORB
 %abundance this gives a maximum total water abundance of 500 ppm.
-WatInterp=linspace(0,1666,1667); 
+WatInterp=linspace(0,1666,1667);
 
-for m=1:numel(WLS) 
+for m=1:numel(WLS)
     %For each wavelength, fit a line to the ESPAT values for each heating step spectrum
     ESPAT(m,:)=polyfit(WaterPPM(1:4),(1-SSAMORB(m,1:4))./SSAMORB(m,1:4),1); %linear function
     %For each wavelength, create new synthetic spectra from the linear ESPAT relationship
-    TotalWatInterpESPAT(:,m)=1./(ESPAT(m,1).*WatInterp+ESPAT(m,2)+1); 
+    TotalWatInterpESPAT(:,m)=1./(ESPAT(m,1).*WatInterp+ESPAT(m,2)+1);
 end
 
 % Check match of interpolated spectra with real MORB spectra
@@ -99,7 +99,7 @@
     hold on
 end
 for k=1:4
-plot(WLS,TotalWatInterpESPAT(RealMorbIndices(k),:),'Color','black','linestyle','-','linewidth',2);
+    plot(WLS,TotalWatInterpESPAT(RealMorbIndices(k),:),'Color','black','linestyle','-','linewidth',2);
 end
 xlabel('Wavelength (\mum)','FontSize',16,'FontWeight','bold');
 ylabel('SSA','FontSize',16,'FontWeight','bold');
@@ -203,43 +203,43 @@
         PCTsMareImmature(j)=((HighAbunReg-LowAbunReg).*rand(1)+LowAbunReg);
         PCTsHighlandsImmature(j)=0;
         PCTsHighlandsMature(j)=0;
-    %Second half of spectra will be highlands spectra
+        %Second half of spectra will be highlands spectra
     else
         PCTsMareMature(j)=0;
         PCTsMareImmature(j)=0;
         PCTsHighlandsImmature(j)=((HighAbunReg-LowAbunReg).*rand(1)+LowAbunReg);
         PCTsHighlandsMature(j)=((HighAbunReg-LowAbunReg).*rand(1)+LowAbunReg);
     end
     PCTsPyroxene(j)=((HighAbunPyr-LowAbunPyr).*rand(1)+LowAbunPyr);
-%Select a random interpolated MORB spectrum to simulate  hydration
-RandWatAmt(j)=round(rand(1)*max(WatInterp));
-%Round to the nearest 1 ppm
-[M I]=min(abs(RandWatAmt(j)-WatInterp));
-%Actual water amount in mixture is morb abundance*interpolated MORB
-%spectrum used
-ActWatAmt(j)=WatInterp(I)*PCTsMORB(j);
-WtPercent(j,:)=[PCTsMORB(j),PCTsHighlandsMature(j),PCTsHighlandsImmature(j),PCTsMareMature(j),PCTsMareImmature(j),PCTsPyroxene(j)]; %Wt%
-if j<NumSpectra/2
-    WtPercent(j,4)=1-sum(WtPercent(j,1:3))-sum(WtPercent(j,5:6)); % Fill in with mature mare so sum=1
-else
-    WtPercent(j,2)=1-WtPercent(j,1)-sum(WtPercent(j,3:6)); % Fill in with mature highlands so sum=1
-end
-%Placing of SSAs, grain sizes (d) and densities (rho) into arrays
-SSAMORBChosen=TotalWatInterpESPAT(I,:).';
-AM=[SSAMORBChosen,SSAHighlandsMature,SSAHighlandsImmature,SSAMareMature,SSAMareImmature,SSAPyrox].';
-rho=1000*[densMORB,densReg,densReg,densReg,densReg,densPyrox].';
-d=[dMORB,dReg,dReg,dReg,dReg,dPyrox].';
-%Weighting of SSAs by cross-sectional area
-AV=AM(1:6,:)./(rho(1:6).*d(1:6))/(1/(mean(rho(1:6))*mean(d(1:6))));
-Numerator=WtPercent(j,:).'.*AM./(rho.*d);
-Denominator=WtPercent(j,:).'./(rho.*d);
-% Nonlinear mixing of SSAs to create intimate mixture via Hapke Radiative Transfer modeling
-SSAInt=sum(Numerator,1)/sum(Denominator); %SSA of intimate mixture
-% Need to convert to lidar geometry using SSA and scattering asymmetry
-% factor (p) of 1.5 (see manuscript for details on Hapke model)
-RInt=Hapke_Lidar_R_Function(SSAInt,1.5,WLS);
-SimulationLidarR(j,:)=RInt;
-SimulationSSA(j,:)=SSAInt;
+    %Select a random interpolated MORB spectrum to simulate  hydration
+    RandWatAmt(j)=round(rand(1)*max(WatInterp));
+    %Round to the nearest 1 ppm
+    [M I]=min(abs(RandWatAmt(j)-WatInterp));
+    %Actual water amount in mixture is morb abundance*interpolated MORB
+    %spectrum used
+    ActWatAmt(j)=WatInterp(I)*PCTsMORB(j);
+    WtPercent(j,:)=[PCTsMORB(j),PCTsHighlandsMature(j),PCTsHighlandsImmature(j),PCTsMareMature(j),PCTsMareImmature(j),PCTsPyroxene(j)]; %Wt%
+    if j<NumSpectra/2
+        WtPercent(j,4)=1-sum(WtPercent(j,1:3))-sum(WtPercent(j,5:6)); % Fill in with mature mare so sum=1
+    else
+        WtPercent(j,2)=1-WtPercent(j,1)-sum(WtPercent(j,3:6)); % Fill in with mature highlands so sum=1
+    end
+    %Placing of SSAs, grain sizes (d) and densities (rho) into arrays
+    SSAMORBChosen=TotalWatInterpESPAT(I,:).';
+    AM=[SSAMORBChosen,SSAHighlandsMature,SSAHighlandsImmature,SSAMareMature,SSAMareImmature,SSAPyrox].';
+    rho=1000*[densMORB,densReg,densReg,densReg,densReg,densPyrox].';
+    d=[dMORB,dReg,dReg,dReg,dReg,dPyrox].';
+    %Weighting of SSAs by cross-sectional area
+    AV=AM(1:6,:)./(rho(1:6).*d(1:6))/(1/(mean(rho(1:6))*mean(d(1:6))));
+    Numerator=WtPercent(j,:).'.*AM./(rho.*d);
+    Denominator=WtPercent(j,:).'./(rho.*d);
+    % Nonlinear mixing of SSAs to create intimate mixture via Hapke Radiative Transfer modeling
+    SSAInt=sum(Numerator,1)/sum(Denominator); %SSA of intimate mixture
+    % Need to convert to lidar geometry using SSA and scattering asymmetry
+    % factor (p) of 1.5 (see manuscript for details on Hapke model)
+    RInt=Hapke_Lidar_R_Function(SSAInt,1.5,WLS);
+    SimulationLidarR(j,:)=RInt;
+    SimulationSSA(j,:)=SSAInt;
 end
 %%
 %Add random Gaussian noise to reflectance based on the instrument SNR
@@ -250,9 +250,9 @@
 %% Plot a subset of the mixtures in Reflectance
 figure(2);
 if NumSpectra<200
-plot(WLS,SimulationLidarRNoisy(:,:),'color',[0,0,0,0.5]);
+    plot(WLS,SimulationLidarRNoisy(:,:),'color',[0,0,0,0.5]);
 else
-plot(WLS,SimulationLidarRNoisy(1:200,:),'color',[0,0,0,0.5]);    
+    plot(WLS,SimulationLidarRNoisy(1:200,:),'color',[0,0,0,0.5]);
 end
 xlabel('Wavelength (\mum)','FontSize',16,'FontWeight','bold');
 ylabel('Lidar Reflectance','FontSize',16,'FontWeight','bold');
@@ -286,29 +286,29 @@
 AV=AVolSolve./(rhoSolve.*dSolve)/(1/(mean(rhoSolve)*mean(dSolve)));
 %% Perform retrieval for each synthetic spectrum
 for j=1:NumSpectra
-%Downselect to laser wavelengths
-WLSCut=LaserLams.';
-ObsSpectrum=SimulationLidarRNoisy(j,:);
-%Convert from noisy lidar reflectance to noisy SSA
-ObsSSA=Hapke_Lidar_SSA_function(ObsSpectrum,1.5,WLS);
-for i=1:numel(LLindices)
-    %Select only measured and endmember SSAs at laser wavelengths
-    AVol(:,i)=AV(:,LLindices(i));
-    ObsIntimateSSACut(i)=ObsSSA(LLindices(i));
-end
-%Perform non-negative least squares algorithm using four endmember spectra
-%and observed SSA values of mixture
-MeasuredAbundances=lsqnonneg(AVol.',ObsIntimateSSACut.');
-MeasuredVolAmplitudes(j,:) = MeasuredAbundances;
-
-%Calculate total water amount using solved-for MORB endmember abundances at total
-%water amounts.
-AreaFactor=0.8428; %This factor is the ratio between the mean grain size and density difference
-%for the mixtures as generated (4 regolith, 1 Morb, 1 Pyroxene) and the retrieved
-%mixtures (2 regolith, 2 MORB.
-MeasuredWater(j,:)=MeasuredAbundances(1)/0.8428*WaterPPM(MorbSelection1)+MeasuredAbundances(2)/0.8428*WaterPPM(MorbSelection2); %retrieve two MORBs
-%Calculate total water error
-WatError(j,:) = MeasuredWater(j)-ActWatAmt(j);
+    %Downselect to laser wavelengths
+    WLSCut=LaserLams.';
+    ObsSpectrum=SimulationLidarRNoisy(j,:);
+    %Convert from noisy lidar reflectance to noisy SSA
+    ObsSSA=Hapke_Lidar_SSA_function(ObsSpectrum,1.5,WLS);
+    for i=1:numel(LLindices)
+        %Select only measured and endmember SSAs at laser wavelengths
+        AVol(:,i)=AV(:,LLindices(i));
+        ObsIntimateSSACut(i)=ObsSSA(LLindices(i));
+    end
+    %Perform non-negative least squares algorithm using four endmember spectra
+    %and observed SSA values of mixture
+    MeasuredAbundances=lsqnonneg(AVol.',ObsIntimateSSACut.');
+    MeasuredVolAmplitudes(j,:) = MeasuredAbundances;
+    
+    %Calculate total water amount using solved-for MORB endmember abundances at total
+    %water amounts.
+    AreaFactor=0.8428; %This factor is the ratio between the mean grain size and density difference
+    %for the mixtures as generated (4 regolith, 1 Morb, 1 Pyroxene) and the retrieved
+    %mixtures (2 regolith, 2 MORB.
+    MeasuredWater(j,:)=MeasuredAbundances(1)/0.8428*WaterPPM(MorbSelection1)+MeasuredAbundances(2)/0.8428*WaterPPM(MorbSelection2); %retrieve two MORBs
+    %Calculate total water error
+    WatError(j,:) = MeasuredWater(j)-ActWatAmt(j);
 end
 %% Plot the results showing scatter around 1:1 retrieval and histogram of total water error
 MareError=WatError(1:end/2);

repro/simulator.py

@@ -8,7 +8,7 @@
 """
 from collections import namedtuple
 from functools import partial
-from typing import Any, Callable
+from typing import Callable
 
 import numpy as np
 from numpy.typing import ArrayLike
@@ -27,12 +27,12 @@ def interpolate_series(data: Series, new_index: ArrayLike) -> Series:
     return Series(data=new_y, index=new_index)
 
 
-def ordinary_least_squares(x: Any, y: Callable, yx: Any):
+def ordinary_least_squares(x: float, y: Callable, yx: float) -> float:
     """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.
+    x (float): The argument to y.
+    yx (float): The observation. Must be the same type as the result of y.
     """
     return sum((y(x) - yx) ** 2)
 
@@ -61,7 +61,9 @@ def backscattering(h: float, g: float) -> float:
 def bidirectional_reflectance(
     SSA: float, P: float, mu: float, mu0: float, B: float
 ) -> float:
-    """Bidirectional reflectance, see Equation 1.
+    """Estimate bidirectional reflectance from single-scattering albedo.
+
+    see Equation 1 in the manuscript:
 
         R = (ω/4) (μ₀ / (μ + μ₀)) {(1 + B)P + H(ω)H₀(ω) - 1}
 
@@ -85,7 +87,7 @@ def bidirectional_reflectance(
         B (float): backscattering function
 
     Returns:
-        float: _description_
+        float: bidrectional reflectance
     """
 
     def h_func(SSA: float, mu: float, mu0: float) -> HFunction:
@@ -123,14 +125,14 @@ def h_func(SSA: float, mu: float, mu0: float) -> HFunction:
 
 
 def reflectance_to_ssa(
-    Refl: Series,
-    P: float = 0.15,
+    reflectance: Series,
+    asymmetry_factor: float = 0.81,
     emission_angle: float = 0,
     incident_angle: float = 30,
     phase_angle: float = 30,
     filling_factor: float = 0.41,
 ):
-    """Convert reflectance spectrum to single-scattering albedo (SSA)
+    """Estimate single-scattering albedo (SSA) from reflectance
 
     Uses scipy.optimize.fmin (equivalent to MATLAB fminsearch) to minimize
     ordinary least squares distance between SSA obtained from the supplied
@@ -144,10 +146,8 @@ def reflectance_to_ssa(
     Default parameter values replicate src/Hapke_Inverse_Function_Passive.m
 
     Args:
-        R (Series): Bidirectional reflectance, see Equation 1.
-        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).
+        reflectance (Series): Bidirectional reflectance, R. see Equation 1.
+        asymmetry_factor (float, optional): Scattering asymmetry factor, p.
         emission_angle (float, optional): Emission angle in degrees.
             Defaults to 0.
         incident_angle (float, optional): Incident angle in degrees.
@@ -163,13 +163,18 @@ def reflectance_to_ssa(
     B = backscattering(h, g)
 
     w = []
-    for index, value in Refl.items():
+    for index, value in reflectance.items():
         w0 = 0.5  # initial guess, ω₀
         # turn bidrectional_reflectance() into the form y=f(x)
-        y = partial(bidirectional_reflectance, P=P, mu=mu, mu0=mu0, B=B)
+        y = partial(
+            bidirectional_reflectance, P=asymmetry_factor, mu=mu, mu0=mu0, B=B
+        )
+
+        # formulate ordinary least squares between estimated reflectance, y
+        # and observed reflectance, yx, and arrange into the form y=f(x)
         OLS = partial(
             ordinary_least_squares, y=y, yx=value
-        )  # turn least squares into the form y=f(x)
+        )
         w.append(
             fmin(
                 OLS,