PythonicDISORT package

PythonicDISORT.pydisort module

PythonicDISORT.pydisort.pydisort(tau_arr, omega_arr, NQuad, Leg_coeffs_all, mu0, I0, phi0, NLeg=None, NFourier=None, b_pos=0, b_neg=0, only_flux=False, f_arr=0, NT_cor=False, BDRF_Fourier_modes=[], s_poly_coeffs=array([], shape=(1, 0), dtype=float64), use_banded_solver_NLayers=10, autograd_compatible=False)

Solves the 1D RTE for the fluxes, and optionally intensity, of a multi-layer atmosphere with the specified optical properties, boundary conditions and sources. Optionally performs delta-M scaling and NT corrections. Refer to the *_test.ipynb Jupyter Notebooks in the pydisotest directory for examples of use.

See https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html#1.-USER-INPUT-REQUIRED:-Choose-parameters for a more detailed explanation of each parameter. See https://pythonic-disort.readthedocs.io/en/latest/Pythonic-DISORT.html#2.-PythonicDISORT-modules-and-outputs for a more detailed explanation of each output. The notebook also has numerous examples of this function being called.

Parameters:

  • tau_arr (array or float) – Optical depth of the lower boundary of each atmospheric layer.

  • omega_arr (array or float) – Single-scattering albedo of each atmospheric layer.

  • NQuad (int) – Number of mu quadrature nodes, i.e. number of streams.

  • Leg_coeffs_all (2darray) – All available unweighted phase function Legendre coefficients. Each row pertains to an atmospheric layer (from top to bottom). Each coefficient should be between 0 and 1 inclusive.

  • mu0 (float) – Cosine of polar angle of the incident beam.

  • I0 (float) – Intensity of the incident beam.

  • phi0 (float) – Azimuthal angle of the incident beam.

  • NLeg (optional, int) – Number of phase function Legendre coefficients to use in the pre-correction solver.

  • NFourier (optional, int) – Number of Fourier modes to use to construct the intensity function.

  • b_pos (optional, 2darray or float) – Dirichlet condition at the bottom boundary for the upward direction. Each column pertains to a Fourier mode (ascending order).

  • b_neg (optional, 2darray or float) – Dirichlet condition at the top boundary for the downward direction. Each column pertains to a Fourier mode (ascending order).

  • only_flux (optional, bool) – Do NOT compute the intensity function?

  • f_arr (optional, array or float) – Fractional scattering into peak for each atmospheric layer. Each row pertains to an atmospheric layer (from top to bottom). We recommend setting f_arr to Leg_coeffs_all[NQuad], or Leg_coeffs_all[:, NQuad] for a multi-layer atmosphere.

  • NT_cor (optional, bool) – Perform Nakajima-Tanaka intensity corrections?

  • BDRF_Fourier_modes (optional, list of functions) – BDRF Fourier modes, each a float, or a function with arguments mu, -mu_p of type array which output has the same dimensions as the outer product of the two arrays.

  • s_poly_coeffs (optional, array) – Polynomial coefficients of isotropic internal sources. Each row pertains to an atmospheric layer (from top to bottom). Arrange coefficients from lowest order term to highest.

  • use_banded_solver_NLayers (optional, int) – At or above how many atmospheric layers should scipy.linalg.solve_banded be used?

  • autograd_compatible (optional, bool) – If True, the autograd package: https://github.com/HIPS/autograd can be used to compute the tau-derivatives of the output functions but pydisort will be less efficient.

Returns:

  • mu_arr (array) – All mu (cosine of polar angle) quadrature nodes.

  • 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) to switch to an antiderivative of this function with respect to tau. For example, Fp(tau_arr, True) - Fp(np.insert(tau_arr[:-1] + 1e-15, 0, 0), True) will produce an array of the tau-integral over each layer. Pass return_tau_arr to return tau_arr (defaults to False).

  • Fm(tau) (function) – (Energetic) Flux function with argument tau (type: array or float) for negative (downward) mu values. Returns a tuple of the diffuse and direct flux magnitudes respectively where each entry is of the same type and size as tau. Pass is_antiderivative_wrt_tau = True (defaults to False) to switch to an antiderivative of this function with respect to tau. Pass return_tau_arr to return tau_arr (defaults to False).

  • u0(tau) (function) – Zeroth Fourier mode of the intensity with argument tau (type: array or float). Returns an ndarray with axes corresponding to variation with mu and tau respectively. This function is useful for calculating actinic fluxes and other quantities of interest, but reclassification of delta-scaled flux and other corrections must be done manually (for actinic flux subroutines.generate_diff_act_flux_funcs will automatically perform the reclassification). Pass is_antiderivative_wrt_tau = True (defaults to False) to switch to an antiderivative of this function with respect to tau. Pass return_tau_arr to return tau_arr (defaults to False).

  • u(tau, phi) (function, optional) – 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 is_antiderivative_wrt_tau = True (defaults to False) to switch to an antiderivative of this function with respect to tau. Pass return_Fourier_error = True (defaults to False) to return the Cauchy / Fourier convergence evaluation (type: float) for the last Fourier term. Pass return_tau_arr to return tau_arr (defaults to False).

PythonicDISORT.subroutines module

PythonicDISORT.subroutines.Clenshaw_Curtis_quad(Nphi, c=0, d=6.283185307179586)

Generates Clenshaw-Curtis quadrature zero points and weights for integration from c to d.

Parameters:

  • Nphi (int) – Number of quadrature nodes.

  • c (float, optional) – Start of integration interval.

  • d (float, optional) – End of integration interval.

Returns:

  • array – Quadrature zero points.

  • array – Quadrature weights.

PythonicDISORT.subroutines.Gauss_Legendre_quad(N, c=0, d=1)

Generates Gauss-Legendre quadrature zero points and weights for integration from c to d.

Parameters:

  • N (int (will be converted to int)) – Number of quadrature nodes.

  • c (float, optional) – Lower bound of integration interval.

  • d (float, optional) – Upper bound of integration interval.

Returns:

  • array – Quadrature zero points.

  • array – Quadrature weights.

PythonicDISORT.subroutines.Planck(T, WVNM)

The Planck function for intensity leaving a blackbody surface, with units W / m^2 to match Stamnes’ DISORT.

Parameters:

  • T (float or array) – Temperatures in kelvin.

  • WVNM (float) – Wavenumber with units m^-1.

Returns:

Emitted power per unit area from a blackbody surface with units W / m^2.

Return type:

float or array

PythonicDISORT.subroutines.affine_transform_poly_coeffs(poly_coeffs, a_arr, b_arr)

Given a polynomial C_0 + C_1 x + … C_n x^n and the affine transformation y -> ax + b, determine the coefficients of the polynomial D_0 + D_1 y + … D_n y^n.

Parameters:

  • poly_coeffs (2darray) – Array with columns [C_0, C_1, …, C_n]. The rows correspond to different affine transformations.

  • a_arr (array) – Scale factors.

  • b_arr (array) – Translations.

Returns:

transformed_poly_coeffs – Array with columns [D_0, D_1, …, D_n]. The rows correspond to different affine transformations.

Return type:

2darray

PythonicDISORT.subroutines.atleast_2d_append(*arys)

View inputs as arrays with at least two dimensions. Dimensions are added as necessary to the back of the shape tuple rather than to the front.

This is exactly NumPy’s numpy.atleast_2d function but altered to add dimensions to the back of the shape tuple rather than to the front. View the documentation for numpy.atleast_2d at https://numpy.org/doc/stable/reference/generated/numpy.atleast_2d.html.

Parameters:

arys1, arys2, … (array_like) – One or more array-like sequences. Non-array inputs are converted to arrays. Arrays that already have two or more dimensions are preserved.

Returns:

res, res2, … – An array, or list of arrays, each with a.ndim >= 2. Copies are avoided where possible, and views with two or more dimensions are returned.

Return type:

ndarray

PythonicDISORT.subroutines.blackbody_contrib_to_BCs(T, WVNMLO, WVNMHI, **kwargs)

Compute blackbody contribution to each BC (b_pos or b_neg), i.e. the blackbody emission of that boundary, with units W / m^2. This convenience function is provided to help match the inputs for Stamnes’ DISORT to those for PythonicDISORT. Users will have to manually adjust emissivities but PythonicDISORT.subroutines.generate_emissivity_from_BDRF can help with that.

Parameters:

  • T (float or array) – Temperatures in kelvin.

  • WVNMLO (float) – Lower bound of wavenumber interval with units m^-1. This variable is identically named in Stamnes’ DISORT.

  • WVNMHI (float) – Upper bound of wavenumber interval with units m^-1. This variable is identically named in Stamnes’ DISORT.

  • **kwargs – Keyword arguments to pass to scipy.integrate.quad_vec.

Returns:

The blackbody emission of each boundary with units W / m^2.

Return type:

float or array

PythonicDISORT.subroutines.cache_BDRF_Fourier_modes(N, BDRF_Fourier_modes, mu0=0)

If the same BDRF and number of streams will be used repeatedly, consider using this function to cache BDRF_Fourier_modes and pass cached_BDRF_Fourier_modes to pydisort. It is also possible to cache BDRF_Fourier_modes with respect to a fixed mu0.

Parameters:

  • N (int) – Number of upper hemisphere quadrature nodes. Equal to NQuad // 2.

  • BDRF_Fourier_modes (list of functions and scalars) – BDRF Fourier modes, each a scalar representing a constant function, or a function with arguments mu, -mu_p of type array which output has the same dimensions as the outer product of the two arrays.

  • mu0 (float, optional) – Cosine of polar angle of the incident beam.

Returns:

cached_BDRF_Fourier_modes – Cached BDRF Fourier modes, each a function with arguments mu, -mu_p of type array and which output is a scalar (if the Fourier mode was a scalar) or has the same dimensions as the outer product of the two arrays.

Return type:

list of functions

PythonicDISORT.subroutines.calculate_nu(mu, phi, mu_p, phi_p)

Calculates the cosine of the scattering angle nu between incident angle (mu_p, phi_p) and scattering angle (mu, phi).

Parameters:

  • mu (array) – Cosine of outgoing polar angles.

  • phi (array) – Outgoing azimuthal angles.

  • mu_p (array) – Cosine of incident polar angles.

  • phi_p (array) – Incident azimuthal angles.

Returns:

nu – Cosine of scattering angles which axes capture variation with mu, phi, mu_p, phi_p respectively.

Return type:

ndarray

PythonicDISORT.subroutines.generate_FD_mat(Ntau, a, b)

Generates a sparse first derivative central difference (second-order accuracy) matrix on [a,b] in csr format with Ntau grid points. Second order forward or backward differences are used at the boundaries.

Parameters:

  • Nphi (int) – Number of grid nodes.

  • a (float, optional) – Start of diffentiation interval.

  • b (float, optional) – End of diffentiation interval.

Returns:

  • array – The differentiation grid.

  • sparse 2darray – The finite difference matrix.

PythonicDISORT.subroutines.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 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:

u0 (func) – Zeroth Fourier mode of the intensity. See the fourth “return” of the pydisort function.

Returns:

  • 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 this function with respect to tau. Pass return_tau_arr to return tau_arr (defaults to False).

  • 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) to switch to an antiderivative of this function with respect to tau. Pass return_tau_arr to return tau_arr (defaults to False).

PythonicDISORT.subroutines.generate_emissivity_from_BDRF(N, zeroth_BDRF_Fourier_mode)

Use Kirchoff’s law of thermal radiation to determine the (directional) emissivity of the surface given the zeroth Fourier mode of the Bi-Directional Reflectance Function (BDRF). This computation is automatic and internal to Stamnes’ DISORT. This function supplements PythonicDISORT.subroutines.blackbody_contrib_to_BCs.

Parameters:

  • N (int) – Number of upper hemisphere quadrature nodes. Equal to NQuad // 2.

  • zeroth_BDRF_Fourier_mode (function or float) – Zeroth BDRF Fourier mode with arguments mu, -mu_p of type array and which output has the same dimensions as the outer product of the two arrays. A scalar input represents a constant function and in that case the output is simply one minus the scalar input.

Returns:

Emissivity for the blackbody contribution to the lower boundary source b_neg.

Return type:

array or float

PythonicDISORT.subroutines.generate_s_poly_coeffs(tau_arr, TEMPER, WVNMLO, WVNMHI, **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. Note that the coefficients will be multiplied by emissivity factors equal to 1 - omega_arr inside pydisort, i.e. Kirchoff’s law of thermal radiation is enforced, just like in Stamnes’ DISORT.

Parameters:

  • tau_arr (array or float) – Optical depth of the lower boundary of each atmospheric layer.

  • TEMPER (array) – Temperature in kelvin at each boundary / interface from top to bottom. This variable is identically named in Stamnes’ DISORT.

  • WVNMLO (float) – Lower bound of wavenumber interval with units m^-1. This variable is identically named in Stamnes’ DISORT.

  • WVNMHI (float) – Upper bound of wavenumber interval with units m^-1. This variable is identically named in Stamnes’ DISORT.

  • **kwargs – Keyword arguments to pass to scipy.integrate.quad_vec.

Returns:

s_poly_coeffs – Polynomial coefficients of isotropic internal sources. Each row pertains to an atmospheric layer (from top to bottom). The coefficients are arranged from lowest order term to highest.

Return type:

2darray

PythonicDISORT.subroutines.interpolate(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 of the Comprehensive Documentation.

Parameters:

u or u0 (function) – Non-interpolated intensity function as given by pydisort.

Returns:

u_interpol – Intensity function with arguments (mu, tau, phi) (for u) or (mu, tau) (for u0) each an array or scalar. Returns an ndarray with axes corresponding to variation with each argument in the same order. Pass is_antiderivative_wrt_tau = True (defaults to False) to switch to an antiderivative of this function with respect to tau. Pass return_Fourier_error = True (defaults to False) to return the Cauchy / Fourier convergence evaluation (type: float) for the last Fourier term. Pass return_tau_arr to return tau_arr (defaults to False).

Return type:

function

PythonicDISORT.subroutines.linear_spline_coefficients(x, y, check_inputs=True)

Compute the coefficients of a linear spline interpolation.

Parameters:

  • x (array) – Array of x data points.

  • y (array) – Array of y data points.

Returns:

The coefficients with axes (segment, ascending polynomial order) which is required by pydisort.

Return type:

2darray

PythonicDISORT.subroutines.prepend(arr, arr_len, value)

Prepend value to array. Faster than np.insert(arr, 0, value).

Parameters:

  • arr (array) – The 1D array to which to prepend.

  • arr_len (int) – Length of original array (arr).

  • value (float) – Value to prepend.

Returns:

Array of length arr_len + 1 with value prepended.

Return type:

array

PythonicDISORT.subroutines.to_diag_ordered_form(A, Nsuperdiags, Nsubdiags)

Convert some matrix A to the diagonal ordered form required by scipy.linalg.solve_banded.

Parameters:

  • A (2darray) – The square matrix to be converted.

  • Nsuperdiags (int) – The number of super-diagonals.

  • Nsubdiags (int) – The number of sub-diagonals.

Returns:

The diagonal ordered form matrix as required by scipy.linalg.solve_banded.

Return type:

2darray

PythonicDISORT.subroutines.transform_interval(arr, c, d, a, b)

Affine transformation of an array from interval [a, b] to [c, d].

Parameters:

  • arr (array) – The 1D array to transform.

  • c (float) – Lower bound of interval [c, d].

  • d (float) – Upper bound of interval [c, d].

  • a (float, optional) – Lower bound of interval [a, b].

  • b (float, optional) – Upper bound of interval [a, b].

Returns:

The transformed 1D array.

Return type:

array

PythonicDISORT.subroutines.transform_weights(weights, c, d, a, b)

Transforms an array of quadrature weights from interval [a, b] to [c, d].

Parameters:

  • weights (array) – The weights to transform.

  • c (float) – Lower bound of interval [c, d].

  • d (float) – Upper bound of interval [c, d].

  • a (float, optional) – Lower bound of interval [a, b].

  • b (float, optional) – Upper bound of interval [a, b].

Returns:

The transformed weights.

Return type:

array