Fixed major bugs with respect to the isotropic source (Thank you Rich…

…ard Larsson for catching these bugs!).

Added an interpolation subroutine with discussion and accuracy tests in the Jupyter Notebook.

README.md

@@ -7,8 +7,10 @@ On the other hand, PythonicDISORT should be easier to install, use, and modify t
 
 PythonicDISORT is based on Stamnes' FORTRAN DISORT (see References, in particular [2, 3, 8]) and has its main features: multi-layer solver, 
 delta-M scaling, Nakajima-Tanaka (NT) corrections, only flux option, direct beam source, isotropic internal source (blackbody emission), 
-Dirichlet boundary conditions (diffuse flux boundary sources), Bi-Directional Reflectance Function (BDRF) for surface reflection, and more.
-In addition, subroutine to calculate actinic fluxes was added to satisfy a user request and further feature requests as well as feedback are welcome.
+Dirichlet boundary conditions (diffuse flux boundary sources), Bi-Directional Reflectance Function (BDRF) for surface reflection, 
+interpolation with respect to polar angle and more.
+In addition, we added a subroutine to calculate actinic fluxes to satisfy a user request, and integration with respect to tau was also added.
+Further feature requests as well as feedback are welcome.
 
 You may contact me, Dion, through dh3065@columbia.edu.
 
@@ -29,8 +31,8 @@ Not only are there verification tests in `Pythonic-DISORT.ipynb`,
 most of the test problems in Stamnes' `disotest.f90` (download DISORT 4.0.99 from http://www.rtatmocn.com/disort/) have also been recreated.
 In these tests, the solutions from PythonicDISORT are compared against solutions 
 from a F2PY-wrapped Stamnes' DISORT (version 4.0.99; wrapper inspired by https://github.com/kconnour/pyRT_DISORT). With PyTest installed, execute the console command `pytest` 
-in the `pydisotest` directory to run these tests. The `pydisotest` directory also contains Jupyter Notebooks that correspond to each test
-to show how it was implemented. These notebooks double up as examples of how to use PythonicDISORT.
+in the `pydisotest` directory to run these tests. The `pydisotest` directory also contains Jupyter Notebooks to show the implementation of each test.
+These notebooks double up as examples of how to use PythonicDISORT.
 
 # Installation
 

docs/conf.py

@@ -9,7 +9,7 @@
 project = 'Pythonic DISORT'
 copyright = '2023, HO Jia Xu Dion'
 author = 'Dion HO Jia Xu'
-release = '0.8.0'
+release = '0.9.0'
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

pyproject.toml

@@ -8,7 +8,7 @@ dependencies = [
   "scipy>=1.8.0",
 ]
 name = "PythonicDISORT"
-version = "0.8.0"
+version = "0.9.0"
 authors = [
   { name="Dion HO Jia Xu", email="dh3065@columbia.edu" },
 ]

src/PythonicDISORT/_assemble_intensity_and_fluxes.py

@@ -233,8 +233,8 @@ def u(tau, phi, return_Fourier_error=False, is_antiderivative_wrt_tau=False):
                     K_collect_0,
                     G_inv_collect_0,
                     mu_arr,
-                    _autograd_compatible,
-                    is_antiderivative_wrt_tau
+                    is_antiderivative_wrt_tau,
+                    _autograd_compatible
                 )
 
                 if _autograd_compatible:
@@ -386,8 +386,8 @@ def u0(tau, is_antiderivative_wrt_tau=False, _return_act_dscale_for_reclass=Fals
                 K_collect_0,
                 G_inv_collect_0,
                 mu_arr,
-                _autograd_compatible,
-                is_antiderivative_wrt_tau
+                is_antiderivative_wrt_tau,
+                _autograd_compatible
             )
             u0 = u0 + _mathscr_v_contribution
 
@@ -439,8 +439,8 @@ def flux_up(tau, is_antiderivative_wrt_tau=False):
                 K_collect_0,
                 G_inv_collect_0,
                 mu_arr,
-                _autograd_compatible,
-                is_antiderivative_wrt_tau
+                is_antiderivative_wrt_tau,
+                _autograd_compatible
             )[:N, :]
         else:
             _mathscr_v_contribution = 0
@@ -518,8 +518,8 @@ def flux_down(tau, is_antiderivative_wrt_tau=False):
                 K_collect_0,
                 G_inv_collect_0,
                 mu_arr,
-                _autograd_compatible,
-                is_antiderivative_wrt_tau
+                is_antiderivative_wrt_tau,
+                _autograd_compatible
             )[N:, :]
         else:
             _mathscr_v_contribution = 0

src/PythonicDISORT/_solve_for_coeffs.py

@@ -124,7 +124,8 @@ def _solve_for_coeffs(
         # _mathscr_v_contribution
         if there_is_iso_source and m_equals_0:
             _mathscr_v_contribution_top = -_mathscr_v(
-                        np.array([0]), np.array([0]),
+                        np.array([0]), 
+                        np.array([0]),
                         s_poly_coeffs[0:1, :],
                         Nscoeffs,
                         G_collect_m[0:1, N:, :],
@@ -139,38 +140,51 @@ def _solve_for_coeffs(
                 _mathscr_v_contribution_middle = (
                     _mathscr_v(
                         tau_arr[:-1],
-                        indices,
-                        s_poly_coeffs[:-1, :],
+                        indices + 1,
+                        s_poly_coeffs,
                         Nscoeffs,
-                        G_collect_m[:-1, :, :],
-                        K_collect_m[:-1, :],
-                        G_inv_collect_0[:-1, :, :],
+                        G_collect_m,
+                        K_collect_m,
+                        G_inv_collect_0,
                         mu_arr,
                     )
                     - _mathscr_v(
                         tau_arr[:-1],
                         indices,
-                        s_poly_coeffs[1:, :],
+                        s_poly_coeffs,
                         Nscoeffs,
-                        G_collect_m[1:, :, :],
-                        K_collect_m[1:, :],
-                        G_inv_collect_0[1:, :, :],
+                        G_collect_m,
+                        K_collect_m,
+                        G_inv_collect_0,
                         mu_arr,
                     )
-                ).ravel()
+                ).ravel(order='F')
 
-            _mathscr_v_contribution_bottom = _mathscr_v(
-                    tau_arr[-1:], np.array([0]),
+            _mathscr_v_contribution_bottom = -_mathscr_v(
+                    tau_arr[-1:], 
+                    np.array([0]),
                     s_poly_coeffs[-1:, :],
                     Nscoeffs,
-                    G_collect_m[-1:, N:, :],
+                    G_collect_m[-1:, :N, :],
                     K_collect_m[-1:, :],
                     G_inv_collect_0[-1:, :, :],
-                    mu_arr).ravel()
+                    mu_arr
+                ).ravel()
             if NBDRF > 0:
                 _mathscr_v_contribution_bottom = (
-                    R + np.eye(N)
-                ) @ _mathscr_v_contribution_bottom
+                    _mathscr_v_contribution_bottom
+                    + R
+                    @ _mathscr_v(
+                        tau_arr[-1:],
+                        np.array([0]),
+                        s_poly_coeffs[-1:, :],
+                        Nscoeffs,
+                        G_collect_m[-1:, N:, :],
+                        K_collect_m[-1:, :],
+                        G_inv_collect_0[-1:, :, :],
+                        mu_arr
+                    ).ravel()
+                )
 
             _mathscr_v_contribution = np.concatenate(
                 [

src/PythonicDISORT/pydisort.py

@@ -82,7 +82,7 @@ def pydisort(
     -------
     mu_arr : array or float
         All `mu` (cosine of polar angle) quadrature nodes.
-    Fp(tau) : function, function
+    Fp(tau) : function
         (Energetic) Flux function with argument `tau` (type: array or float) for positive (upward) `mu` values.
         Returns the diffuse flux magnitudes (same type and size as `tau`).
         Pass `is_antiderivative_wrt_tau = True` (defaults to `False`)
@@ -104,7 +104,7 @@ def pydisort(
         Pass `is_antiderivative_wrt_tau = True` (defaults to `False`)
         to switch to an antiderivative of the function with respect to `tau`.
     u(tau, phi) : function, optional
-        Intensity function with arguments `(tau, phi)` of types `(array or float, array or float)`.
+        Intensity function with arguments `(tau, phi)` each of type array or float.
         Returns an ndarray with axes corresponding to variation with `mu, tau, phi` respectively.
         Pass `return_Fourier_error = True` (defaults to `False`) to return the 
         Cauchy / Fourier convergence evaluation (type: float) for the last Fourier term.

src/PythonicDISORT/subroutines.py

@@ -1,7 +1,6 @@
 import numpy as np
 import scipy as sc
 from math import pi
-from numpy.polynomial.legendre import leggauss
 
 
 def transform_interval(arr, c, d, a, b):
@@ -71,7 +70,7 @@ def calculate_nu(mu, phi, mu_p, phi_p):
 
     Returns
     -------
-    ndarray
+    nu : ndarray
         Cosine of scattering angles which axes capture variation with `mu, phi, mu_p, phi_p` respectively.
 
     """
@@ -104,7 +103,7 @@ def Gauss_Legendre_quad(N, c=0, d=1):
         Quadrature weights.
 
     """
-    mu_arr_pos, W = leggauss(int(N))
+    mu_arr_pos, W = np.polynomial.legendre.leggauss(int(N))
 
     return transform_interval(mu_arr_pos, c, d, -1, 1), transform_weights(W, c, d, -1, 1)
 
@@ -224,8 +223,10 @@ def atleast_2d_append(*arys):
 
 def generate_diff_act_flux_funcs(u0):
     """Generates the up and down diffuse actinic flux functions respectively.
-    This is a use case of the ``u0`` function that is an output of ``pydisort``.
+    This is a use case of the ``u0`` function which is an output of ``pydisort``.
     The reclassification of delta-scaled actinic flux is automatically performed.
+    The computation of actinic fluxes is similar to that of energetic fluxes
+    and the latter is discussed in section 3.8 of the Comprehensive Documentation.
 
     Parameters
     ----------
@@ -235,12 +236,12 @@ def generate_diff_act_flux_funcs(u0):
 
     Returns
     -------
-    function
+    Fp_act(tau) : function
         Actinic flux function with argument `tau` (type: array) for positive (upward) `mu` values.
         Returns the diffuse flux magnitudes (type: array).
         Pass `is_antiderivative_wrt_tau = True` (defaults to `False`)
         to switch to an antiderivative of the function with respect to `tau`.
-    function
+    Fm_act(tau) : function
         Actinic flux function with argument `tau` (type: array) for negative (downward) `mu` values.
         Returns the diffuse flux magnitudes (type: array).
         Pass `is_antiderivative_wrt_tau = True` (defaults to `False`)
@@ -262,6 +263,65 @@ def flux_act_down_diffuse(tau, is_antiderivative_wrt_tau=False):
     return flux_act_up, flux_act_down_diffuse
 
 
+def interpolate(mu_arr_pos, u):
+    """Polynomial (Barycentric) interpolation with respect to mu. The output 
+    is a function that is continuous and variable in all three arguments: mu, tau and phi.
+    Discussed in sections 3.7 and 6.3 in the Comprehensive Documentation.
+
+    Parameters
+    ----------
+    mu_arr_pos : array
+        Positive `mu` (cosine of polar angle) quadrature nodes.
+        Equivalent to `mu_arr[:N]`.
+    u : function
+        Non-interpolated intensity function as outputted by `pydisort`.
+
+    Returns
+    -------
+    u_interpol : function
+        Intensity function with arguments `(mu, tau, phi)` each of type array or float.
+        Returns an ndarray with axes corresponding to variation with `mu, tau, phi` respectively.
+        Pass `return_Fourier_error = True` (defaults to `False`) to return the 
+        Cauchy / Fourier convergence evaluation (type: float) for the last Fourier term.
+        Pass `is_antiderivative_wrt_tau = True` (defaults to `False`)
+        to switch to an antiderivative of the function with respect to `tau`.
+
+    """
+    N = len(mu_arr_pos)
+    u_pos_interpol = sc.interpolate.BarycentricInterpolator(mu_arr_pos)
+    u_neg_interpol = sc.interpolate.BarycentricInterpolator(-mu_arr_pos)
+
+    def u_interpol(mu, tau, phi, return_Fourier_error=False, is_antiderivative_wrt_tau=False):
+        assert np.all(np.abs(mu) <= 1)
+        mu = np.atleast_1d(mu)
+
+        if return_Fourier_error:
+            u_cache, Fourier_error = u(tau, phi, True, is_antiderivative_wrt_tau)
+        else:
+            u_cache = u(tau, phi, False, is_antiderivative_wrt_tau)
+
+        u_shape =  np.shape(u_cache)
+        assert N == u_shape[0] // 2 # Check that mu_arr_pos (of length N) matches the function u
+
+        results = np.empty((len(mu),) + u_shape[1:])
+        mask_pos = mu > 0
+        mask_else = ~mask_pos
+
+        if np.any(mask_pos):
+            u_pos_interpol.set_yi(u_cache[:N])
+            results[mask_pos] = u_pos_interpol(mu[mask_pos])
+        if np.any(mask_else):
+            u_neg_interpol.set_yi(u_cache[N:])
+            results[mask_else] = u_neg_interpol(mu[mask_else])
+
+        if return_Fourier_error:
+            return np.squeeze(results)[()], Fourier_error
+        else:
+            return np.squeeze(results)[()]
+
+    return u_interpol
+
+
 def _mathscr_v(tau,                             # Input optical depths
                 l,                              # Layer index of each input optical depth
                 s_poly_coeffs,                  # Polynomial coefficients of isotropic source
@@ -270,8 +330,8 @@ def _mathscr_v(tau,                             # Input optical depths
                 K,                              # Eigenvalues
                 G_inv,                          # Inverse of eigenvector matrix
                 mu_arr,                         # Quadrature nodes for both hemispheres
+                is_antiderivative_wrt_tau=False, # Switch to an antiderivative of the function?
                 _autograd_compatible=False,      # Should the output functions be compatible with autograd?
-                is_antiderivative_wrt_tau=False # Switch to an antiderivative of the function?
                 ):
     """Particular solution for isotropic internal sources.
     Refer to Section 3.6.1 of the Comprehensive Documentation.
@@ -289,8 +349,8 @@ def _mathscr_v(tau,                             # Input optical depths
     | `K`                         | `NLayers<= x NQuad`                   |
     | `G_inv`                     | `NLayers<= x NQuad x NQuad` or `None` |
     | `mu_arr`                    | `NQuad`                               |
-    | `_autograd_compatible`      | boolean                               |
     | `is_antiderivative_wrt_tau` | boolean                               |
+    | `_autograd_compatible`      | boolean                               |
     
     Notable internal variables of _mathscr_v
     |     Variable     |                Shape                |