Upgraded subroutines.generate_s_poly_coeffs to be able to compute emissivities if omega_arr is supplied.

Slight improvements to docs and input checks.

Author: dhjx1996

docs/Pythonic-DISORT.ipynb

@@ -1035,7 +1035,7 @@
     "\n",
     "$$B(T(\\tau)) = a_{l0} + a_{l1} \\tau$$\n",
     "\n",
-    "for constants $a_{l0}, a_{l1}$ to be specified. The convenience function `PythonicDISORT.subroutines.generate_s_poly_coeffs` is provided to generate `s_poly_coeffs` such that the outputs of PythonicDISORT will match those of Stamnes' DISORT given the same inputs. PythonicDISORT allows much more flexibility in choosing blackbody emission profiles than DISORT."
+    "for constants $a_{l0}, a_{l1}$ to be specified. The convenience function `PythonicDISORT.subroutines.generate_s_poly_coeffs` is provided to generate `s_poly_coeffs` such that the outputs of PythonicDISORT will match those of Stamnes' DISORT given the same inputs. PythonicDISORT allows much more flexibility in choosing blackbody emission profiles than DISORT though."
    ]
   },
   {
@@ -1056,7 +1056,8 @@
     "\n",
     "##################################################################\n",
     "\n",
-    "emissivity = 1 # This should actually be determined using Kirchoff's law of thermal radiation"
+    "emissivity = 1 # In practice this should be determined using Kirchoff's law\n",
+    "               # of thermal radiation, and Kirchoff's law is enforced in Stamnes' DISORT"
    ]
   },
   {
@@ -1071,7 +1072,7 @@
     "# This is the coefficient matrix \\mathscr{a}\n",
     "s_poly_coeffs = PythonicDISORT.subroutines.generate_s_poly_coeffs(\n",
     "    tau_arr, TEMPER, WVNMLO, WVNMHI\n",
-    ") * emissivity\n",
+    ") * emissivity[:, None]\n",
     "\n",
     "#####################################################################################################\n",
     "\n",
@@ -5996,7 +5997,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.2"
+   "version": "3.11.8"
   }
  },
  "nbformat": 4,

pydisotest/8_test.ipynb

@@ -307,7 +307,7 @@
       "\n",
       "Upward (diffuse) fluxes\n",
       "Difference = 1.8563882056565895e-08\n",
-      "Difference ratio = 3.5919432662810775e-07\n",
+      "Difference ratio = 0.9999999953171316\n",
       "\n",
       "Downward (diffuse) fluxes\n",
       "Difference = 6.096862970039751e-08\n",

pydisotest/9_test.ipynb

@@ -656,12 +656,12 @@
       "Max pointwise differences\n",
       "\n",
       "Upward (diffuse) fluxes\n",
-      "Difference = 2.5765318011072846e-08\n",
-      "Difference ratio = 6.125620222015996e-07\n",
+      "Difference = 2.576531828168971e-08\n",
+      "Difference ratio = 6.125620247828119e-07\n",
       "\n",
       "Downward (diffuse) fluxes\n",
-      "Difference = 1.1920928955078125e-07\n",
-      "Difference ratio = 3.967778597939065e-07\n",
+      "Difference = 1.1920928977282585e-07\n",
+      "Difference ratio = 3.967778612553841e-07\n",
       "\n",
       "Direct (downward) fluxes\n",
       "Difference = 0.0\n",
@@ -672,7 +672,7 @@
       "At tau = 0.0\n",
       "Max pointwise difference = 5.53211412651633e-08\n",
       "At tau = 1.05\n",
-      "Max pointwise difference ratio = 9.347970166516183e-07\n",
+      "Max pointwise difference ratio = 9.347970250378194e-07\n",
       "\n"
      ]
     }
@@ -1036,9 +1036,9 @@
       "Intensities\n",
       "\n",
       "At tau = 21.0\n",
-      "Max pointwise difference = 0.0027373875326694996\n",
+      "Max pointwise difference = 0.0027373888669937063\n",
       "At tau = 21.0\n",
-      "Max pointwise difference ratio = 0.0017543192767219193\n",
+      "Max pointwise difference ratio = 0.0017543201318550158\n",
       "\n"
      ]
     }
@@ -1129,7 +1129,7 @@
     {
      "data": {
       "text/plain": [
-       "<matplotlib.legend.Legend at 0x2406e552810>"
+       "<matplotlib.legend.Legend at 0x2a1c21ed810>"
       ]
      },
      "execution_count": 31,
@@ -1452,9 +1452,9 @@
       "Intensities\n",
       "\n",
       "At tau = 0.0\n",
-      "Max pointwise difference = 0.06855659986274398\n",
+      "Max pointwise difference = 0.06911426569533097\n",
       "At tau = 0.0\n",
-      "Max pointwise difference ratio = 0.03692288128517093\n",
+      "Max pointwise difference ratio = 0.03919935155907838\n",
       "\n"
      ]
     }
@@ -1487,8 +1487,8 @@
       "Max pointwise differences\n",
       "\n",
       "Upward (diffuse) fluxes\n",
-      "Difference = 0.004794243247684626\n",
-      "Difference ratio = 0.0009782676032947248\n",
+      "Difference = 0.004794243247683738\n",
+      "Difference ratio = 0.0009782676032945436\n",
       "\n",
       "Downward (diffuse) fluxes\n",
       "Difference = 0.00870651337438666\n",
@@ -1501,9 +1501,9 @@
       "Intensities\n",
       "\n",
       "At tau = 2.1\n",
-      "Max pointwise difference = 0.012109081699047897\n",
+      "Max pointwise difference = 0.010067514209913586\n",
       "At tau = 2.1\n",
-      "Max pointwise difference ratio = 0.007287281829096514\n",
+      "Max pointwise difference ratio = 0.005664629397768609\n",
       "\n"
      ]
     }

src/PythonicDISORT/_assemble_intensity_and_fluxes.py

@@ -592,6 +592,6 @@ def flux_down(tau, is_antiderivative_wrt_tau=False, return_tau_arr=False):
         # --------------------------------------------------------------------------------------------------------------------------
 
     if only_flux:
-        return flux_up, flux_down, u0#, GC_collect_0, K_collect_0, B_collect_0
+        return flux_up, flux_down, u0#, GC_collect_0, K_collect_0#, B_collect_0
     else:
         return flux_up, flux_down, u0, u

src/PythonicDISORT/pydisort.py

@@ -230,10 +230,6 @@ def pydisort(
     # The user must supply at least as many phase function Legendre coefficients as intended for use
     if not NLeg > 0:
         raise ValueError("The number of phase function Legendre coefficients must be positive.") 
-    if not np.all(omega_arr * Leg_coeffs_all[:, 0] == omega_arr):
-        raise ValueError("The first phase function Legendre coefficient must equal 1.") 
-    if not (np.all(-1 <= Leg_coeffs_all) and np.all(Leg_coeffs_all <= 1)):
-        raise ValueError("The phase function Legendre coefficients must all be between -1 and 1.") 
     if not NLeg <= NLeg_all:
         raise ValueError("`NLeg` cannot be larger than the number of phase function Legendre coefficients provided.")
     # Ensure that the first dimension of the following inputs corresponds to the number of layers
@@ -245,6 +241,11 @@ def pydisort(
         raise ValueError("The length of `f_arr` does not match the number of layers as deduced from the length of `tau_arr`.")
     if there_is_iso_source and not np.shape(s_poly_coeffs)[0] == NLayers:
         raise ValueError("The zeroth dimension of the shape of `s_poly_coeffs` does not match the number of layers as deduced from the length of `tau_arr`.")
+    # Value checks on the phase function Legendre coefficients
+    if not np.all(omega_arr * Leg_coeffs_all[:, 0] == omega_arr):
+        raise ValueError("The first phase function Legendre coefficient must equal 1.") 
+    if not (np.all(-1 <= Leg_coeffs_all) and np.all(Leg_coeffs_all <= 1)):
+        raise ValueError("The phase function Legendre coefficients must all be between -1 and 1.") 
     # Conditions on the number of quadrature angles (NQuad), Legendre coefficients (NLeg) and loops (NFourier)
     if not NQuad >= 2:
         raise ValueError("There must be at least two streams.")
@@ -353,9 +354,10 @@ def pydisort(
         scaled_s_poly_coeffs = (scaled_s_poly_coeffs / rescale_factor).copy()
     else:
         rescale_factor = np.max((I0, np.max(b_pos), np.max(b_neg)))
-        I0 = (I0 / rescale_factor).copy()
-        b_pos = (b_pos / rescale_factor).copy()
-        b_neg = (b_neg / rescale_factor).copy()
+        if rescale_factor != 0:
+            I0 = (I0 / rescale_factor).copy()
+            b_pos = (b_pos / rescale_factor).copy()
+            b_neg = (b_neg / rescale_factor).copy()
     # --------------------------------------------------------------------------------------------------------------------------
 
     if NT_cor and not only_flux and there_is_beam_source and np.any(f_arr > 0) and NLeg < NLeg_all:

src/PythonicDISORT/subroutines.py

@@ -384,11 +384,12 @@ def linear_spline_coefficients(x, y, check_inputs=True):
 
 
 
-def generate_s_poly_coeffs(tau_arr, TEMPER, WVNMLO, WVNMHI, **kwargs):
+def generate_s_poly_coeffs(tau_arr, TEMPER, WVNMLO, WVNMHI, omega_arr=None, **kwargs):
     """Generate DISORT-equivalent ``s_poly_coeffs`` for input into ``pydisort``.
     This convenience function is provided to help match the inputs for Stamnes' DISORT to those for PythonicDISORT.
-    Users will have to manually adjust emissivities, however, and they should note that PythonicDISORT allows much
-    more flexibility in choosing blackbody emission profiles than DISORT.
+    If ``omega_arr`` is specified, the coefficients will be multiplied by emissivity factors equal to ``1 - omega_arr`` 
+    per Kirchoff's law of thermal radiation, and Kirchoff's law is enforced in Stamnes' DISORT. 
+    Otherwise, the emissivities will be set to 1 and users can multiply their own emissivity factors.
     
     Parameters
     ----------
@@ -401,6 +402,8 @@ def generate_s_poly_coeffs(tau_arr, TEMPER, WVNMLO, WVNMHI, **kwargs):
         Lower bound of wavenumber interval with units m^-1. This variable is identically named in Stamnes' DISORT.
     WVNMHI : scalar
         Upper bound of wavenumber interval with units m^-1. This variable is identically named in Stamnes' DISORT.
+    omega_arr : optional, array or scalar
+        Single-scattering albedo of each atmospheric layer.
     **kwargs
         Keyword arguments to pass to ``scipy.integrate.quad_vec``.
         
@@ -414,12 +417,35 @@ def generate_s_poly_coeffs(tau_arr, TEMPER, WVNMLO, WVNMHI, **kwargs):
     """
     tau_arr = np.atleast_1d(tau_arr)
     if not len(TEMPER) == len(tau_arr) + 1:
-        raise ValueError("Missing temperature specification at some boundaries / interfaces.")
-    
+        raise ValueError(
+            "Missing temperature specification at some boundaries / interfaces."
+        )
+
     tau_arr_with_0 = np.insert(tau_arr, 0, 0)
-    blackbody_emission_at_each_boundary = sc.integrate.quad_vec(lambda WVNM: Planck(TEMPER, WVNM), WVNMLO, WVNMHI, **kwargs)[0]
-    
-    return linear_spline_coefficients(tau_arr_with_0, blackbody_emission_at_each_boundary, check_inputs=False)
+    blackbody_emission_at_each_boundary = sc.integrate.quad_vec(
+        lambda WVNM: Planck(TEMPER, WVNM), WVNMLO, WVNMHI, **kwargs
+    )[0]
+
+    if omega_arr == None:
+        return (
+            linear_spline_coefficients(
+                tau_arr_with_0, blackbody_emission_at_each_boundary, check_inputs=False
+            )
+        )
+    elif np.isscalar(omega_arr):
+        return (
+            linear_spline_coefficients(
+                tau_arr_with_0, blackbody_emission_at_each_boundary, check_inputs=False
+            )
+            * (1 - omega_arr)
+        )
+    else:
+        return (
+            linear_spline_coefficients(
+                tau_arr_with_0, blackbody_emission_at_each_boundary, check_inputs=False
+            )
+            * (1 - omega_arr)[:, None]
+        )
 
 
 
@@ -448,7 +474,9 @@ def generate_emissivity_from_BDRF(N, zeroth_BDRF_Fourier_mode):
         return 1 - zeroth_BDRF_Fourier_mode
     else:
         mu_arr_pos, W = Gauss_Legendre_quad(N)
-        return 1 - 2 * zeroth_BDRF_Fourier_mode(mu_arr_pos, mu_arr_pos) * mu_arr_pos[None, :] @ W
+    return (
+        1 - 2 * zeroth_BDRF_Fourier_mode(mu_arr_pos, mu_arr_pos) * mu_arr_pos[None, :] @ W
+    )
 
 
 
@@ -778,7 +806,7 @@ def mathscr_b(i):
 
 
 
-def _compare(results, mu_to_compare, reorder_mu, flux_up, flux_down, u=None):
+def _compare(results, mu_to_compare, reorder_mu, flux_up, flux_down, u=None, div_threshold=1e-15):
     """Performs pointwise comparisons between results from Stamnes' DISORT,
     which are stored in ``.npz`` files, against results from PythonicDISORT. Used in our PyTests.
 
@@ -807,7 +835,7 @@ def _compare(results, mu_to_compare, reorder_mu, flux_up, flux_down, u=None):
         diff_flux_up,
         flup,
         out=np.zeros_like(diff_flux_up),
-        where=flup > 1e-8,
+        where=flup > div_threshold,
     )
     print("Difference =", np.max(diff_flux_up))
     print("Difference ratio =", np.max(ratio_flux_up))
@@ -820,7 +848,7 @@ def _compare(results, mu_to_compare, reorder_mu, flux_up, flux_down, u=None):
         diff_flux_down_diffuse,
         rfldn,
         out=np.zeros_like(diff_flux_down_diffuse),
-        where=rfldn > 1e-8,
+        where=rfldn > div_threshold,
     )
     print("Difference =", np.max(diff_flux_down_diffuse))
     print(
@@ -836,7 +864,7 @@ def _compare(results, mu_to_compare, reorder_mu, flux_up, flux_down, u=None):
         diff_flux_down_direct,
         rfldir,
         out=np.zeros_like(diff_flux_down_direct),
-        where=rfldir > 1e-8,
+        where=rfldir > div_threshold,
     )
     print("Difference =", np.max(diff_flux_down_direct))
     print(
@@ -852,7 +880,7 @@ def _compare(results, mu_to_compare, reorder_mu, flux_up, flux_down, u=None):
             diff,
             uu[mu_to_compare],
             out=np.zeros_like(diff),
-            where=uu[mu_to_compare] > 1e-8,
+            where=uu[mu_to_compare] > div_threshold,
         )
         max_diff_tau_index = np.argmax(np.max(np.max(diff, axis=0), axis=1))
         max_ratio_tau_index = np.argmax(np.max(np.max(diff_ratio, axis=0), axis=1))