Producing the Coulsen Tables with the Low-Level Interface

In this example we will reproduce part of the Coulsen tables using the low level interface to SASKTRAN DO. These tables involve a single layer conservative scattering atmosphere with a Lambertian surface underneath.

%matplotlib inline
import sasktran.disco.lowlevel as lowlevel
import numpy as np

We can begin by setting up our configuration, we are performing a single calculation for 3 stokes parameters, a single layer, and we will set the number of streams to be a high number (40). We will explicitly disable the pseudo spherical correction since these tables correspond to a plane parallel atmosphere. For brevity we will only calculate 6 of the entries in the table.

We start by creating all of the low level input objects.

nstr = 40
nlyr = 1
nwavel = 1
nderiv = 0
nstokes = 3
nlos = 6

config = lowlevel.Config(nstr, nwavel, nlyr, nstokes, 0, use_pseudo_spherical=False)
weightingfunctions = lowlevel.WeightingFunctions(nstr, nlos, nstokes, nwavel, nderiv)
viewing_geometry = lowlevel.ViewingGeometry(nlos)
atmosphere = lowlevel.Atmosphere(nstr, nlyr, nwavel)

Next let’s setup our viewing geometry. We have to be slightly careful, DO uses the same definition of cosine solar zenith angle \(\mu_0\) (cos_sza), and viewing cosine zenith \(\mu\) (cos_vza), however we have a minus sign convention difference on relative azimuth angle \(\phi\) (saa).

viewing_geometry.cos_sza = 0.2

viewing_geometry.cos_vza[0] = 0.02
viewing_geometry.cos_vza[1] = 0.2
viewing_geometry.cos_vza[2] = 1.00
viewing_geometry.cos_vza[3] = 0.02
viewing_geometry.cos_vza[4] = 0.2
viewing_geometry.cos_vza[5] = 1.00

viewing_geometry.saa[0] = 0
viewing_geometry.saa[1] = 0
viewing_geometry.saa[2] = 0
viewing_geometry.saa[3] = -np.pi/3
viewing_geometry.saa[4] = -np.pi/3
viewing_geometry.saa[5] = -np.pi/3

And finally our atmosphere optical quantities. The table is for a single layer of optical depth 0.5 and single scatter albedo of 1. The table contains only Rayleigh scattering with no depolarization factor, of which the Legendre coefficients are easy to calculate.

atmosphere.od[0, :] = 0.5

atmosphere.ssa[0, :] = 1.0

atmosphere.a1[0, 0, :] = 1
atmosphere.a1[2, 0, :] = 0.5

atmosphere.a2[2, 0, :] = 3

atmosphere.a4[1, 0, :] = 3/2

atmosphere.b1[2, 0, :] = np.sqrt(6.0) * 0.5

Next we set the albedo to 0

atmosphere.albedo[0] = 0.0

And finally, even though the layer boundary altitudes are not used in this problem, we set it to an arbitrary number for good measure.

atmosphere.layer_boundary_altitudes[0] = 1000000

Next we perform the calculation and multiply by \(\pi\) since the tables assume a solar flux of \(\pi\) while DO assumes a solar flux of 1.

output = lowlevel.calculate(atmosphere, config, weightingfunctions, viewing_geometry).xarray().isel(wavelength=0)['radiance'] * np.pi

We can print the first element of the stokes vector’s and see that they are the same as those in the table to ~4 decimal places.

print(output.isel(stokes=0))

<xarray.DataArray 'radiance' (los: 6)>
array([0.44130148, 0.26939546, 0.05300458, 0.30091557, 0.19368246,
       0.05300458])
Dimensions without coordinates: los