### Initialize multi-layer matrices

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Example setup for collecting random matrices for layers and quadrature points.

```python
G_collect_m = np.random.random((NLayers, NQuad, NQuad))
K_collect_m = np.random.random((NLayers, NQuad))
```

--------------------------------

### Setup in _solve_for_gen_and_part_sols.py

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Initializes variables for Legendre polynomial calculations and weighted coefficients.

```python
M_inv = 1 / mu_arr_pos

ells = np.arange(m, NLeg)
degree_tile = np.tile(ells, (N, 1)).T
poch = sc.special.poch(ells + m + 1, -2 * m)
signs = np.empty(NLeg - m)
signs[::2] = 1
signs[1::2] = -1

asso_leg_term_pos = sc.special.lpmv(m, degree_tile, mu_arr_pos)
asso_leg_term_neg = asso_leg_term_pos * signs[:, None]
asso_leg_term_mu0 = sc.special.lpmv(m, ells, -mu0)


weighted_Leg_coeffs_l = ((2 * np.arange(NLeg) + 1) * Leg_coeffs[l, m:])
omega_l = omega_arr[l]
omega_times_Leg_coeffs = (omega_l / 2) * weighted_Leg_coeffs_l
```

--------------------------------

### Setup evaluation timing

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Prepares various solver outputs and test grids to measure evaluation performance.

```python
# Setup
flux_up, flux_down, u0, u = PythonicDISORT.pydisort(
    tau_arr, omega_arr,
    NQuad,
    Leg_coeffs_all,
    mu0, I0, phi0,
)[1:]

u_1 = PythonicDISORT.pydisort(
    tau_arr[-1], omega_arr[0],
    NQuad,
    Leg_coeffs_all[0, :],
    mu0, I0, phi0,
)[-1]

u_NT = PythonicDISORT.pydisort(
    tau_arr, omega_arr,
    NQuad,
    Leg_coeffs_all,
    mu0, I0, phi0,
    f_arr=f_arr,
    NT_cor=True,
)[-1]

u_NT_1 = PythonicDISORT.pydisort(
    tau_arr[-1], omega_arr[0],
    NQuad,
    Leg_coeffs_all[0, :],
    mu0, I0, phi0,
    f_arr=f_arr[0],
    NT_cor=True,
)[-1]

u_iso = PythonicDISORT.pydisort(
    tau_arr, omega_arr,
    NQuad,
    Leg_coeffs_all,
    mu0, 0, phi0,
    s_poly_coeffs=s_poly_coeffs,
)[-1]

Ntau_time = 1000
Nphi_time = 100
tau_test_arr_time = np.sort(np.random.random(Ntau_time) * tau_arr[-1])
phi_arr_time = np.random.random(Nphi_time) * 2 * pi
```

--------------------------------

### Setup test points and phi quadrature

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Initializes random optical depth test points and calculates the phi quadrature grid points for flux verification.

```python
Ntau = 100 # Number of tau test points
tau_test_arr = np.random.random(Ntau) * tau_arr[-1]

# Number of phi grid points
# This selection should ensure that the phi quadrature is at least as accurate as the mu quadrature
Nphi = int((NQuad * pi) // 2) * 2 + 1
phi_arr, full_weights_phi = PythonicDISORT.subroutines.Clenshaw_Curtis_quad(Nphi)
```

--------------------------------

### Configure direct beam source parameters

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Example configuration for the direct beam source parameters including mu0, I0, and phi0.

```python
#################### SUPPLIED TO PYDISORT ########################

# Cosine of polar angle (UMU0) -- MOST IMPORTANT
mu0 = 0.6
# Intensity (FBEAM) -- Simply a scale factor
I0 = 100
# Azimuthal angle (PHI0) -- No effect on fluxes
phi0 = 0.9 * pi

##################################################################
```

--------------------------------

### Calculate reflectance, transmittance, and absorptance

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Example calculation using PythonicDISORT to determine flux and derive optical properties.

```python
# Example calculation
flux_up, flux_down = PythonicDISORT.pydisort(
    tau_arr[-1], omega_arr[0],
    NQuad,
    Leg_coeffs_all[-1, :],
    mu0, I0, phi0,
    only_flux=True
)[1:3]

Ref = flux_up(0) / (I0 * mu0)
Trans = np.sum(flux_down(tau_arr[-1]), axis=0) / (I0 * mu0)

print("Reflectance =", Ref)
print("Transmittance =", Trans)
print("Absorptance =", 1 - Ref - Trans)
```

--------------------------------

### Example output for reflectance calculation

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

The expected output values resulting from the reflectance, transmittance, and absorptance calculation.

```text
Reflectance = 0.6295114216333902
Transmittance = 0.3704723849448305
Absorptance = 1.619342177927585e-05
```

--------------------------------

### Execute DISORT or load results

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Runs the DISORT solver if installed, otherwise loads pre-saved results from a file.

```python
if disort_is_installed:
    # Run disort, putting DFDT, UAVG, and UU in a, b, and c, respectively
    rfldir, rfldn, flup, dfdt, uavg, uu_interpol, albmed, trnmed = disort.disort(usrang, usrtau, ibcnd, onlyfl, prnt, plank, lamber, deltamplus, do_pseudo_sphere, dtauc, ssalb,
                            pmom, temper, wvnmlo, wvnmhi, utau, umu0, phi0 * 180/pi, umu, phi * 180/pi, fbeam, fisot, albedo, btemp, ttemp,
                            temis, earth_radius, h_lyr, rhoq, rhou, rho_accurate, bemst, emust, accur, header, rfldir,
                            rfldn, flup, dfdt, uavg, uu_interpol, albmed, trnmed)
else:
    results = np.load("section6_testresults2.npz")
    # Load saved results from Stamnes' DISORT
    uu_interpol = results["uu_interpol"]
    # Load comparison points
    mu_test_arr_pos = results["mu_test_arr_pos"]
    test_tau = results["test_tau"]
    test_phi = results["test_phi"]

    Nmu = len(mu_test_arr_pos)
    mu_test_arr_neg = -np.flip(mu_test_arr_pos)
```

--------------------------------

### Installation commands

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Commands to install the PythonicDISORT package using pip, conda, or mamba.

```bash
pip install PythonicDISORT
conda install -c conda-forge pythonicdisort
mamba install -c conda-forge pythonicdisort
```

--------------------------------

### Execute or load Stamnes DISORT results

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Runs the DISORT solver if installed, otherwise loads pre-saved results from a file.

```python
if disort_is_installed:
    # Run disort, putting DFDT, UAVG, and UU in a, b, and c, respectively
    rfldir, rfldn, flup, dfdt, uavg, uu, albmed, trnmed = disort.disort(usrang, usrtau, ibcnd, onlyfl, prnt, plank, lamber, deltamplus, do_pseudo_sphere, dtauc, ssalb,
                            pmom, temper, wvnmlo, wvnmhi, utau, umu0, phi0 * 180/pi, umu, phi * 180/pi, fbeam, fisot, albedo, btemp, ttemp,
                            temis, earth_radius, h_lyr, rhoq, rhou, rho_accurate, bemst, emust, accur, header, rfldir,
                            rfldn, flup, dfdt, uavg, uu, albmed, trnmed)
else:
    results = np.load("section6_testresults1.npz")
    # Load saved results from Stamnes' DISORT
    uu = results["uu"]
    flup = results["flup"]
    rfldn = results["rfldn"]
    rfldir = results["rfldir"]
    # Load comparison points
    tau_test_arr = results["tau_test_arr"]
    phi_arr = results["phi_arr"]
```

--------------------------------

### Check for Stamnes DISORT installation

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Attempts to import the disort module to determine if reference solutions can be generated locally.

```python
disort_is_installed = True
try:
    import disort
except ImportError:
    disort_is_installed = False
```

--------------------------------

### Delta-M scaling of tau_arr

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Example demonstrating how to perform Delta-M scaling on an optical depth array (tau_arr) by calculating layer thicknesses and applying the scaling factor.

```python
# Example: Proper delta-M scaling of `tau_arr`
thickness_arr = np.append(tau_arr[0], np.diff(tau_arr))
scaled_thickness_arr = (1 - omega_arr * f_arr) * thickness_arr
# Include 0 in the array for computational reasons
scaled_tau_arr_with_0 = np.append(0, np.cumsum(scaled_thickness_arr))
print("Delta-M scaled `tau_arr` with 0 =", scaled_tau_arr_with_0)
```

```text
Delta-M scaled `tau_arr` with 0 = [0.         0.49498871 0.98997741 1.48496612 1.97995483 2.47494354
 2.96993224 3.46492095 3.95990966 4.45489836 4.94988707 5.44487578
 5.93986449 6.43485319 6.9298419  7.42483061 7.91981931]
```

--------------------------------

### Initialize quadrature and test points

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Sets up the phi quadrature grid and random optical depth test points.

```python
# Number of phi grid points
# This selection should ensure that the phi quadrature is at least as accurate as the mu quadrature
Nphi = int((NQuad * pi) // 2) * 2 + 1
phi_arr, full_weights_phi = PythonicDISORT.subroutines.Clenshaw_Curtis_quad(Nphi)

Ntau = 100  # Number of tau test points
tau_test_arr = np.random.random(Ntau) * tau_arr[-1]

#tau_test_arr_small = np.random.choice(tau_test_arr, size=10, replace=False)
#phi_arr_small = np.random.choice(phi_arr, size=10, replace=False)

MU_ARR, PHI_ARR = np.meshgrid(phi_arr, mu_arr_RO)
```

--------------------------------

### Solving the System Naively

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Demonstrates the naive approach to solving the system, which results in high condition numbers for large optical thicknesses.

```python
if m < NBDRF:
    BDRF_LHS_contribution = R @ G_neg
else:
    BDRF_LHS_contribution = 0

LHS = np.vstack((G_neg, (G_pos - BDRF_LHS_contribution) * El[None, :]))

print("Condition number =", np.linalg.cond(LHS))
# The system is too ill-conditioned to be solved if the optical thickness is large
# xi = np.linalg.solve(LHS, RHS)
```

```text
Condition number = 2.0702755673804135e+163
```

--------------------------------

### Assemble coefficient matrix and source vector in _solve_for_gen_and_part_sols.py

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Constructs the final coefficient matrix A and the direct beam source vector X_tilde.

```python
alpha = M_inv[:, None] * (D_pos * W[None, :] - np.eye(N))
beta = M_inv[:, None] * D_neg * W[None, :]
A = np.vstack([np.hstack([-alpha, -beta]), np.hstack([beta, alpha])])
if I0 > 0:
    X_tilde = np.concatenate([-M_inv * X_pos, M_inv * X_neg])
```

--------------------------------

### Verify particular solution for direct beam source

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Sets up random optical depth test points for verifying the direct beam source solution.

```python
Ntau = 100

tau_test_arr = np.random.random(Ntau) * tau_arr[-1]
```

--------------------------------

### Print solver configuration parameters

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Displays the current configuration values for NQuad, NLeg, NFourier, and NLayers.

```python
print("NQuad, NLeg, NFourier, NLayers =", NQuad, NLeg, NFourier, NLayers)
```

--------------------------------

### Configure Stamnes DISORT parameters

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Initializes the input arrays and configuration parameters required for the Stamnes DISORT solver.

```python
nlyr = NLayers
nmom = NLeg
nstr = NQuad
numu = NQuad
nphi = Nphi
ntau = Ntau
usrang = False
usrtau = True
ibcnd = 0
onlyfl = False
prnt = np.array([False, False, False, False, False])  # Prints to cmd instead of this notebook
plank = False
lamber = False
deltamplus = False
do_pseudo_sphere = False
dtauc = np.append(tau_arr[0], np.diff(tau_arr))
ssalb = omega_arr
pmom = Leg_coeffs_all.T
temper = np.zeros(nlyr + 1)
wvnmlo = 0
wvnmhi = 0
utau = tau_test_arr
umu0 = mu0
phi0 = phi0
umu = mu_arr_RO
phi = phi_arr
fbeam = I0
fisot = 0
albedo = 0
btemp = 0
ttemp = 0
temis = 0
earth_radius = 6371
h_lyr = np.zeros(nlyr + 1)
rhoq = np.zeros((nstr // 2, nstr + 1, nstr))
rhou = np.zeros((numu, nstr // 2 + 1, nstr))
rho_accurate = np.zeros((numu, nphi))
bemst = np.zeros(nstr // 2)
emust = np.zeros(numu)
accur = 0
header = "Test Problem 3:  Henyey-Greenstein Scattering, g = 0.75 (Compare To Ref. VH2, Table 37)"
rfldir = np.zeros(ntau)
rfldn = np.zeros(ntau)
flup = np.zeros(ntau)
dfdt = np.zeros(ntau)
uavg = np.zeros(ntau)
uu = np.zeros((numu, ntau, nphi))
albmed = np.zeros(numu)
trnmed = np.zeros(numu)
```

--------------------------------

### Import libraries for testing and exposition

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Additional libraries used specifically for testing and demonstration purposes within notebooks.

```python
# Required only in this notebook for tests and exposition
import autograd as ag
import autograd.numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy import constants
```

--------------------------------

### Timing Stamnes' DISORT with %timeit

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Uses %timeit to measure execution time of disort.disort for intensity and flux-only modes. The code checks if disort is installed; if not, it prints the recorded timing results. Intensity computation is much slower (1.85 s) than flux-only (13.9 ms) because NT corrections are applied to every point.

```python
if disort_is_installed:
    print("Intensity")
    %timeit disort.disort(usrang, usrtau, ibcnd, False, prnt, plank, lamber, deltamplus, do_pseudo_sphere, dtauc, ssalb, pmom, temper, wvnmlo, wvnmhi, utau, umu0, phi0 * 180/pi, umu, phi * 180/pi, fbeam, fisot, albedo, btemp, ttemp, temis, earth_radius, h_lyr, rhoq, rhou, rho_accurate, bemst, emust, accur, header, rfldir, rfldn, flup, dfdt, uavg, uu, albmed, trnmed)
    print()

    print("Only fluxes")
    %timeit disort.disort(usrang, usrtau, ibcnd, True, prnt, plank, lamber, deltamplus, do_pseudo_sphere, dtauc, ssalb, pmom, temper, wvnmlo, wvnmhi, utau, umu0, phi0 * 180/pi, umu, phi * 180/pi, fbeam, fisot, albedo, btemp, ttemp, temis, earth_radius, h_lyr, rhoq, rhou, rho_accurate, bemst, emust, accur, header, rfldir, rfldn, flup, dfdt, uavg, uu, albmed, trnmed)
else:
    print("RECENT TIMING TESTS")
    print()
    print("Intensity")
    print("1.85 s ± 121 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)")
    print()
    print("Only fluxes")
    print("13.9 ms ± 275 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)")
```

```text
Intensity
1.85 s ± 121 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Only fluxes
13.9 ms ± 275 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```

--------------------------------

### Configure boundary condition inputs

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Calculates boundary conditions using blackbody contributions and validates the input shapes for the solver.

```python
#################### SUPPLIED TO PYDISORT ########################

# At bottom of atmosphere
b_pos = (
    PythonicDISORT.subroutines.blackbody_contrib_to_BCs(T_Surf, WVNMLO, WVNMHI)
    * emissivity_surface
)

# At top of atmosphere
b_neg = (
    PythonicDISORT.subroutines.blackbody_contrib_to_BCs(T_CMB, WVNMLO, WVNMHI)
    * emissivity_top_boundary
)

##################################################################

# Code within `pydisort` to ensure that the BC inputs are of the correct shape
b_pos_is_scalar = False
b_neg_is_scalar = False
b_pos_is_vector = False
b_neg_is_vector = False

if len(np.atleast_1d(b_pos)) == 1:
    b_pos_is_scalar = True
elif len(b_pos) == N:
    b_pos_is_vector = True
elif not np.shape(b_pos) == (N, NFourier):
    raise ValueError("The shape of the bottom boundary condition is incorrect.")
if len(np.atleast_1d(b_neg)) == 1:
    b_neg_is_scalar = True
elif len(b_neg) == N:
    b_neg_is_vector = True
elif not np.shape(b_neg) == (N, NFourier):
    raise ValueError("The shape of the top boundary condition is incorrect.")
```

--------------------------------

### Configure Parameters for Large-Scale Timing

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Sets up the simulation parameters to compute solutions across a large number of tau and phi points.

```python
# Changed parameters
ntau = Ntau_time
utau = tau_test_arr_time
nphi = Nphi_time
phi = phi_arr_time
print("Number of tau points =", ntau, ", Number of phi points =", nphi)


rho_accurate = np.zeros((numu, nphi))
rfldir = np.zeros(ntau)
rfldn = np.zeros(ntau)
flup = np.zeros(ntau)
dfdt = np.zeros(ntau)
uavg = np.zeros(ntau)
uu = np.zeros((numu, ntau, nphi))
```

--------------------------------

### Define test parameters

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Sets up the number of streams and test arrays for mu, optical depth, and azimuth angle.

```python
Nmu = 50
mu_test_arr_pos = np.sort(np.random.random(size=Nmu))
mu_test_arr_neg = -np.flip(mu_test_arr_pos)
test_tau = 0.2
test_phi = pi
```

--------------------------------

### Configure DISORT input variables

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Initializes the required input arrays and parameters for the DISORT solver, specifically for Henyey-Greenstein scattering.

```python
# Test Problem 3a:  Henyey-Greenstein Scattering (Compare To Ref. VH2, Table 37)
nlyr = NLayers
nmom = NLeg
nstr = NQuad
numu = 2 * Nmu
nphi = 1
ntau = 1
usrang = True
usrtau = True
ibcnd = 0
onlyfl = False
prnt = np.array([False, False, False, False, False])  # Prints to cmd instead of this notebook
plank = False
lamber = True
deltamplus = False
do_pseudo_sphere = False
dtauc = np.append(tau_arr[0], np.diff(tau_arr))
ssalb = omega_arr
pmom = Leg_coeffs_all.T
temper = np.zeros(nlyr + 1)
wvnmlo = 0
wvnmhi = 0
utau = test_tau
umu0 = mu0
phi0 = phi0
umu = np.concatenate((mu_test_arr_neg, mu_test_arr_pos)) # Needs to be in increasing order
phi = test_phi
fbeam = I0
fisot = 0
albedo = 0
btemp = 0
ttemp = 0
temis = 0
earth_radius = 6371
h_lyr = np.zeros(nlyr + 1)
rhoq = np.zeros((nstr // 2, nstr + 1, nstr))
rhou = np.zeros((numu, nstr // 2 + 1, nstr))
rho_accurate = np.zeros((numu, nphi))
bemst = np.zeros(nstr // 2)
emust = np.zeros(numu)
accur = 0
header = "Test Problem 3:  Henyey-Greenstein Scattering, g = 0.75 (Compare To Ref. VH2, Table 37)"
rfldir = np.zeros(ntau)
rfldn = np.zeros(ntau)
flup = np.zeros(ntau)
dfdt = np.zeros(ntau)
uavg = np.zeros(ntau)
uu_interpol = np.zeros((numu, ntau, nphi))
albmed = np.zeros(numu)
trnmed = np.zeros(numu)
```

--------------------------------

### Prepare test grids and comparison filters

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Reorders mu values and defines a filter to exclude intensities near the direct beam.

```python
# Reorder mu_arr from smallest to largest
reorder_mu = np.argsort(mu_arr)
mu_arr_RO = mu_arr[reorder_mu]

# We may not want to compare intensities around
# the direct beam in the forward direction (the solar aureole)
deg_around_beam_to_not_compare = 0
mu_to_compare = (
    np.abs(np.arccos(np.abs(mu_arr_RO)) - np.arccos(mu0)) * 180 / pi
    > deg_around_beam_to_not_compare
)
```

--------------------------------

### Prepare test points for isotropic internal source verification

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Generates random optical depths and identifies the corresponding layers for testing the isotropic internal source solution.

```python
Ntau = 100

tau_test_arr = np.random.random(Ntau) * tau_arr[-1]
l = np.argmax(tau_test_arr[:, None] <= tau_arr[None, :], axis=1)
l_uniq, l_inv = np.unique(l, return_inverse=True)
```

--------------------------------

### Verify general solution against boundary conditions

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Checks if the general solution satisfies boundary conditions at the top and bottom of the atmosphere using np.allclose.

```python
# At top of atmosphere
assert np.allclose(um(0)[N:], b_neg_m)

# At bottom of atmosphere
if m < NBDRF:
    if I0 > 0:
        surface_reflection = R @ um(tau_arr[-1])[N:] + mathscr_X_pos * np.exp(
            -tau_arr[-1] / mu0
        )
    else:
        surface_reflection = R @ um(tau_arr[-1])[N:]
else:
    surface_reflection = 0

assert np.allclose(um(tau_arr[-1])[:N], b_pos_m + surface_reflection)


print("Passed all tests")
```

--------------------------------

### Configure Parameters for Single-Point Timing

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Adjusts the simulation parameters to compute the solution at a single point for benchmarking purposes.

```python
# Changed parameters
ntau = 1
utau = tau_test_arr_time[Ntau_time//2]
nphi = 1
phi = phi_arr_time[Nphi_time//2]


rho_accurate = np.zeros((numu, nphi))
rfldir = np.zeros(ntau)
rfldn = np.zeros(ntau)
flup = np.zeros(ntau)
dfdt = np.zeros(ntau)
uavg = np.zeros(ntau)
uu = np.zeros((numu, ntau, nphi))
```

--------------------------------

### Prepare mu grid for verification

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Reorders mu_arr from smallest to largest for presentation and efficiency.

```python
# mu_arr is arranged as it is for code efficiency and readability
# For presentation purposes we re-arrange mu_arr from smallest to largest
reorder_mu = np.argsort(mu_arr)

full_weights_mu_RO = full_weights_mu[reorder_mu]
mu_arr_RO = mu_arr[reorder_mu]

MU_ARR, PHI_ARR = np.meshgrid(phi_arr, mu_arr_RO)
```

--------------------------------

### Generate phase function terms in _solve_for_gen_and_part_sols.py

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Calculates phase function terms D and X, including shortcuts for insignificant multiple scattering.

```python
if np.any(np.abs(omega_times_Leg_coeffs) > 1e-8): # There are shortcuts if multiple scattering is insignificant
    D_temp = omega_times_Leg_coeffs[None, :] * (poch[None, :] * asso_leg_term_pos.T)
    D_pos = D_temp @ asso_leg_term_pos
    D_neg = D_temp @ asso_leg_term_neg

    if I0 > 0:
        # Generate X
        X_temp = (
            (omega_l * I0 * (2 - (m == 0)) / (4 * pi))
            * weighted_Leg_coeffs_l
            * (poch * asso_leg_term_mu0)
        )
        X_pos = X_temp @ asso_leg_term_pos
        X_neg = X_temp @ asso_leg_term_neg
        X_tilde = np.concatenate([-M_inv * X_pos, M_inv * X_neg])
```

--------------------------------

### Benchmark solver performance

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Uses %timeit to measure the execution time of various pydisort configurations.

```python
print("Intensity")
%timeit PythonicDISORT.pydisort(tau_arr, omega_arr, NQuad, Leg_coeffs_all, mu0, I0, phi0)
print()

print("Intensity with only one layer")
%timeit PythonicDISORT.pydisort(tau_arr[-1], omega_arr[0], NQuad, Leg_coeffs_all[0, :], mu0, I0, phi0)
print()

print("Intensity with blackbody emission and Lambertian BDRF")
%timeit PythonicDISORT.pydisort(tau_arr, omega_arr, NQuad, Leg_coeffs_all, mu0, I0, phi0, s_poly_coeffs=s_poly_coeffs, BDRF_Fourier_modes=BDRF_Fourier_modes)
print()

print("Only fluxes")
%timeit PythonicDISORT.pydisort(tau_arr, omega_arr, NQuad, Leg_coeffs_all, mu0, I0, phi0, only_flux=True)
print()

print("Only fluxes with delta-M scaling")
%timeit PythonicDISORT.pydisort(tau_arr, omega_arr, NQuad, Leg_coeffs_all, mu0, I0, phi0, only_flux=True, f_arr=f_arr)
```

--------------------------------

### Validate direct beam source solution with autograd

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Uses the autograd package to differentiate the particular solution and verify it against the differential equation, printing the maximum pointwise error ratio.

```python
if I0 > 0:
    mathscr_v_direct = lambda tau: B[:, None] * np.exp(-tau[None, :] / mu0)

    LHS = np.sum(ag.jacobian(mathscr_v_direct)(tau_test_arr), axis=-1)
    RHS = A @ mathscr_v_direct(tau_test_arr) + X_tilde[:, None] * np.exp(
        -tau_test_arr[None, :] / mu0
    )

    print("Max pointwise error ratio =", np.max(np.abs((RHS - LHS) / LHS)))
else:
    print("There is no direct beam source")
```

```text
Max pointwise error ratio = 1.65962699110508e-06
```

--------------------------------

### Integrate intensity to test flux

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Uses einsum to perform numerical integration of the intensity cache to derive upward and downward fluxes.

```python
flux_up_test = np.einsum(
    "itp, i, p -> t",
    mu_arr_pos[:, None, None] * u_cache[:N, :],
    W,
    full_weights_phi,
    optimize=True,
)
flux_down_test = np.einsum(
    "itp, i, p -> t",
    mu_arr_pos[:, None, None] * u_cache[N:, :],
    W,
    full_weights_phi,
    optimize=True,
)
```

--------------------------------

### Gauss-Legendre Quadrature for Mu

Source: https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html

Calculates quadrature weights and points for the mu domain using Double Gauss-Legendre quadrature.

```python
# For positive mu values (the weights are identical for both domains)
mu_arr_pos, W = PythonicDISORT.subroutines.Gauss_Legendre_quad(N) # mu_arr_neg = -mu_arr_pos
mu_arr = np.concatenate([mu_arr_pos, -mu_arr_pos])
full_weights_mu = np.concatenate([W, W])
```