Getting started¶

This tutorial shows to new users how to use nanodisort. The nanodisort workflow is similar to what the CDISORT workflow would be using the C API.

The `DisortState` class¶

The main entry point to nanodisort is the DisortState class. It exposes the eponym C struct to Python through a set of bindings written with the nanobind library. To start running nanodisort, we instantiate this class:

import numpy as np

import nanodisort as nd
from nanodisort.utils import phase_functions as pf

ds = nd.DisortState()

At this point, our DisortState instance is empty. If we print it, we will see the sizes of the memory buffers that will be allocated to its input values. An allocated boolean variable indicates that memory has not been allocated yet:

ds

<DisortState nstr=0 nlyr=0 nmom=0 ntau=0 numu=0 nphi=0 allocated=False>

Memory allocation¶

Since nanodisort reflects the CDISORT interface with barely any changes, it is up to the user to trigger memory allocation for the input fields. Before calling allocate(), we must configure the problem dimensions:

nstr: number of computational streams (must be even; higher values yield better accuracy)
nlyr: number of atmospheric layers
nmom: number of Legendre moments for the phase function; should satisfy nmom >= nstr
ntau: number of optical thickness levels at which outputs are computed
numu: number of user-defined polar viewing directions (0 for flux-only computation)
nphi: number of user-defined azimuthal angles (0 for flux-only computation)

For this tutorial, we model a single Rayleigh-scattering atmospheric layer illuminated by a solar beam, and compute both fluxes and radiances at the layer top and bottom:

ds.nstr = 16  # 16 computational streams
ds.nlyr = 1  # single atmospheric layer
ds.nmom = 16  # nmom = nstr (minimum recommended)
ds.ntau = 2  # outputs at τ = 0 (TOA) and τ = 0.1 (surface)
ds.numu = 6  # 6 polar viewing directions
ds.nphi = 1  # 1 azimuthal angle
ds

<DisortState nstr=16 nlyr=1 nmom=16 ntau=2 numu=6 nphi=1 allocated=False>

With dimensions set, we can now allocate memory:

ds.allocate()
ds

<DisortState nstr=16 nlyr=1 nmom=16 ntau=2 numu=6 nphi=1 allocated=True>

Setting flags¶

The solver’s behavior is controlled by boolean flags. The most commonly used ones are:

usrtau: if True, outputs are computed at user-specified optical thickness (utau); otherwise at layer boundaries
usrang: if True, radiances are computed at user-specified polar and azimuthal directions
lamber: if True, the lower boundary is a Lambertian (perfectly diffuse) reflector
planck: if True, thermal emission is included in the calculation
onlyfl: if True, only fluxes are computed (radiances are skipped)
quiet: if True, diagnostic output from the C library is suppressed
intensity_correction: if True, the Nakajima-Tanaka correction is applied for improved radiance accuracy near the direct-beam direction

For our example, we enable user-specified output levels and angles, use a Lambertian surface, and request both fluxes and radiances:

ds.usrtau = True  # use user-specified output optical thicknesses
ds.usrang = True  # compute radiances at user-specified angles
ds.lamber = True  # Lambertian lower boundary
ds.onlyfl = False  # compute both fluxes and radiances
ds.quiet = True  # suppress C-level diagnostic output
ds.intensity_correction = True  # apply Nakajima-Tanaka correction

Specifying input¶

All array-valued inputs must be assigned after allocate() has been called.

Optical properties¶

Three arrays describe the atmospheric column:

dtauc: extinction optical thickness of each layer, shape (nlyr,)
ssalb: single-scattering albedo of each layer, shape (nlyr,); 0 = purely absorbing, 1 = purely scattering
pmom: Legendre expansion coefficients of the scattering phase function, shape (nmom+1, nlyr); the zeroth moment must be 1 by convention

The utils.phase_functions module provides helpers to generate moment arrays for common phase functions. Here we use Rayleigh scattering, whose phase function has only a non-trivial second moment (\(\beta_2 = 0.1\)):

ds.dtauc = np.array([0.1])  # extinction optical thickness of the layer
ds.ssalb = np.array([1.0])  # pure scattering, no absorption
ds.pmom = pf.rayleigh(ds.nmom).reshape(-1, 1)  # Rayleigh moments, shape (nmom+1, nlyr)

Output levels and viewing angles¶

We request output at the top of the layer (\(\tau = 0\)) and at the surface (\(\tau = 0.1\)). Six polar viewing directions cover the full hemisphere in one azimuthal plane.

In DISORT’s convention, \(\mu < 0\) denotes downward-propagating radiation and \(\mu > 0\) upward-propagating radiation.

ds.utau = np.array([0.0, 0.1])  # output optical thicknesses
ds.umu = np.array([-1.0, -0.5, -0.1, 0.1, 0.5, 1.0])  # cosines of polar angles
ds.phi = np.array([0.0])  # azimuthal angle [degrees]

Boundary conditions¶

We illuminate the layer with a collimated solar beam at \(\theta_0 = 60°\) from the zenith (\(\mu_0 = 0.5\)), using the standard radiometric convention \(F_\mathrm{beam} = \pi\). The surface is a black Lambertian reflector (albedo = 0). There is no diffuse illumination from above.

ds.fbeam = np.pi  # incident solar irradiance [W/m²]
ds.umu0 = 0.5  # cosine of solar zenith angle (θ₀ = 60°)
ds.phi0 = 0.0  # solar azimuth [degrees]
ds.albedo = 0.0  # Lambertian surface albedo
ds.fisot = 0.0  # isotropic top-of-atmosphere illumination

Running the solver¶

With all inputs set, we call solve() to run the DISORT algorithm:

ds.solve()

Collecting the output¶

After solve() returns, all output arrays are populated and ready to read.

Fluxes¶

The flux outputs are 1D arrays indexed by the ntau optical thickness levels:

rfldir: direct-beam downward irradiance
rfldn: diffuse downward irradiance
flup: diffuse upward irradiance

print(f"{'Level':>6}  {'τ':>6}  {'F_dir_dn':>12}  {'F_diff_dn':>12}  {'F_diff_up':>12}")
print("-" * 60)
for i, tau in enumerate(ds.utau):
    print(
        f"{i:>6}  {tau:>6.3f}  {ds.rfldir[i]:>12.6f}  {ds.rfldn[i]:>12.6f}  {ds.flup[i]:>12.6f}"
    )

 Level       τ      F_dir_dn     F_diff_dn     F_diff_up
------------------------------------------------------------
     0   0.000      1.570796     -0.000000      0.143026
     1   0.100      1.286059      0.141711     -0.000000

The direct-beam flux at TOA equals \(F_\mathrm{beam} \cdot \mu_0 = \pi \times 0.5 \approx 1.5708\). After propagation through the layer, Beer’s law gives \(F_\mathrm{dir} = F_\mathrm{beam} \cdot \mu_0 \cdot e^{-\tau/\mu_0} = \pi \times 0.5 \times e^{-0.2} \approx 1.287\).

Radiances¶

The uu array contains the diffuse radiance field with shape (numu, ntau, nphi). We inspect the upward radiances leaving the top of the atmosphere and the downward radiances reaching the surface:

print(
    f"Radiance array shape: {ds.uu.shape}  (numu={ds.numu}, ntau={ds.ntau}, nphi={ds.nphi})\n"
)

print("Upward radiances leaving the atmosphere at TOA (τ = 0.0, μ > 0):")
for i, mu in enumerate(ds.umu):
    if mu > 0:
        print(f"  μ = {mu:+.1f}:  {ds.uu[i, 0, 0]:.6f} W/m²/sr")

print()
print("Downward radiances reaching the surface (τ = 0.1, μ < 0):")
for i, mu in enumerate(ds.umu):
    if mu < 0:
        print(f"  μ = {mu:+.1f}:  {ds.uu[i, 1, 0]:.6f} W/m²/sr")

Radiance array shape: (6, 2, 1)  (numu=6, ntau=2, nphi=1)

Upward radiances leaving the atmosphere at TOA (τ = 0.0, μ > 0):
  μ = +0.1:  0.212784 W/m²/sr
  μ = +0.5:  0.047129 W/m²/sr
  μ = +1.0:  0.023936 W/m²/sr

Downward radiances reaching the surface (τ = 0.1, μ < 0):
  μ = -1.0:  0.023866 W/m²/sr
  μ = -0.5:  0.070090 W/m²/sr
  μ = -0.1:  0.225105 W/m²/sr