The HR engine solves the radiative transfer equation in a region of interest using the successive orders technique on a true spherical Earth. The theoretical basis for the technique can be found in [Bourassa2008] and upgrades specific to the HR engine are described in [Zawada2015]. The engine is well-suited for calculating limb radiances for EUV to NIR (nominally 200 nm to 5 micron) wavelengths measured from high altitudes balloons, aircraft and spacecraft as well as ground-based slant angle observations. It is particularly suitable for conditions where the curvature of the Earth is a significant factor and solar illumination is an important term. The model features the following highlights:

  • Support to calculate the weighting functions of target species.
  • Adapative integration within optically thick regions.
  • Polarized or scalar radiance calculations.
  • Support for gradients across a region of interest using high resolution diffusive fields.
  • Support for terminator conditions.
  • Support for thermal emissions.
  • Support for elevated surface topology.



Here is a simple example that shows basic usage of the radiative transfer engine:

engine             = skif.ISKEngine('HR')                                                 # Create the engine
albedo             = 0.3                                                                  # Set the albedo to 0
usermjd            = 54832.000000000000;                                                  # Set the time of the measurements
satpos             = [366082.10014671006,  1005802.3038196121,  -6874430.9971755166]      # Set the satellite position
lookv              = [0.28845686317656621, 0.79252871806432690, 0.53729960834682389]      # Set the look unit-vector
wavelen            = [294.0,       295.0,       296.0]                                    # Set the wavelengths in nanometers
rayleigh           = skif.ISKOpticalProperty('RAYLEIGH')                                  # Get Rayleigh scattering optical properties of dry air
atmosphere         = skif.ISKClimatology    ('MSIS90')                                    # Get a climatology of air number density
o3_opticalprops    = skif.ISKOpticalProperty('O3_DBM')                                    # Get ozone cross-sections
o3numberdensity    = skif.ISKClimatology    ('O3LABOW')                                   # Get a climatology of ozone number density

engine.AddSpecies( 'SKCLIMATOLOGY_AIRNUMBERDENSITY_CM3', atmosphere, rayleigh );          # Add dry air to the radiative transfer engine
engine.AddSpecies( 'SKCLIMATOLOGY_O3_CM3', o3numberdensity, o3_opticalprops);             # add ozone to the radiative transfer engine
engine.SetAlbedo ( albedo );                                                              # Set the Lambertian albedo of the surface
engine.SetWavelengths( wavelen );                                                         # Add the three wavelengths

engine.AddLineOfSight( usermjd,  satpos, lookv);                                          # Add one line of sight, add more if necessary
ok, radiance = engine.CalculateRadiance();                                                # Calculate the radiance for all lines of sight

It is possible to customize and adjust the performance of the radiative transfer model but you must specify all of the settings,other than wavelength and species, before calling CalculateRadiance() or CalculateStokesVector() as all settings for each instance, other than wavelength and species, are locked in at this time. You must create a new instance of the engine if you wish to change any parameters after this time. The ability to change change wavelengths and species at any time allows the engine to be efficiently adjusted by atmopsheric retrieval codes.


The HR model should be used with caution for:

  • Extremely optical thick conditions like clouds. The code does not handle the large optical depths very well.
  • Nadir geometries. They are supported but do not adapt well to the strengths of the HR model.
  • Scattering and surface BRDF properties with sharp spikes at angles other than forward-scatter are not well supported.

Full Property List

All properties are set through:

ok = engine.SetProperty('PropertyName', value)

The property name is not case sensitive, and value may be a scalar, list, or array.

A complete list of properties is provided for reference.

Geophysical Properties

A list of geophysical properties


SetSun (array v) [Default: Determined from MJD]

Specifies the sun direction in the model. v is a unit vector from the Earth to the Sun. By default the sun location is determined from the average MJD of the input lines of sight.


SurfaceHeight (double d) [Default: d=0]

The height of the Earth’s surface can be changed to account for topography. This value is the altitude of the surface above mean sea level expressed in meters. The model uses this height for the whole of the Earth’s surface across the ray tracing region used in the model. It is useful for consistently high elevation areas such as as Greenland, Antarctica and the Himalayan plateau.

We currently have not implemented any internal topological maps such as Etopos in to Sasktran and the user must choose the surface height. As a guide we note that rays used for limb geometry in the Earth’s stratospheric region span several degrees of angle around the Earth and this scale size is appropriate for selecting a surface height from a topography model.

This method must be called during the initialization phase of the HR RT model. It cannot be set after an instance has called CalculateRadiance().


TOAHeight (double d) [Default: d=100000.0]

The height of the top of the atmosphere in meters. It is assumed there is no atmosphere above this altitude. All climatologies used in the radiative transfer model should support valid (not NaN) values up to this height. We have experienced stange warning messages during the calculations when this requirement was not met.

This method must be called during initialization of the HR RT model. It cannot be set after an instance has called CalculateRadiance().


EarthRadius (double d) [Default: Osculating Sphere Calculation]

Manually sets the radius of the Earth in meters. Normally not required but useful if you want to try and use the same setting as another radiative transfer model. This value cannot be set after an instance has called CalculateRadiance().

Weighting Functions

The HR engine supports analytical computation of weighting functions for absorbing and scattering species. Prior to calculating weighting functions it is necessary to tell the engine that is should calculate weighting functions, CalcWF, and also specify which species need weighting functions, WFSpecies. More details can be found in this description of the HR Weighting Functions.

The default configuration does not calculate weighting functions.


CalcWF (int n) [Default: n=0]

Sets the weighting function calculation mode. For more details see Weighting Functions.

n Setting
0 No weighting function calculation
1 Line of sight weighting function calculation (1D)
2 User specified (1D)
3 Two Dimensional weighting function calculation


WFPrecision (int n) [Default: n=0]

Weighting functions are calculated considering line of sight effects, direct solar effects, and higher order terms. Setting this option to 1 disables the direct solar and higher order terms, greatly improving weighting function calculation speed at the cost of accuracy.

n Setting
0 Weighting functions are calculated with full precision/
1 Weighting functions are calculated neglecting direct solar and higher order effects.


WFHeights (array v) [Default: linspace(500, 99500, 100)]

See Weighting Functions.


WFWidths (array v) [Default: ones(100) * 1000]

See Weighting Functions.


WFSpecies (ISKOpticalProperty* o) [Default: None]

Sets the species to calculate weighting functions for. Species must have scattering or absorbing cross sections. See Weighting Functions for detailed information on weighting functions.

Diffuse Field

The primary option that controls the speed and accuracy of the HR radiative transfer calculation is the number of diffuse profiles. Each diffuse profile represents a single solar zenith angle where the multiple scattering source function is calculated. By default, only one diffuse profile is used in the calculation and is placed at the solar zenith angle of the average tangent point. Whether or not this is appropriate is heavily dependant on a variety of factors. More details at HR :ref: hrdiffuse,


DiffusePlacementType (int n) [Default: See Text]

Changes the placement of diffuse profiles inside the engine. The default value changes depending on how the model is initialized. It is recommended not to change this setting unless there is good reason to.

n Setting
0 Diffuse profiles are placed linearly in angle along the line of sight plane
1 Diffuse profiles are placed linearly in solar zenith angle, forcing one profile at the tangent point
2 Diffuse profiles are placed linearly in solar zenith angle


ForceV21DiffuseIncoming (int n) [Default: n=0]

Changes the spacing of incoming rays on diffuse spheres. SASKTRANV21 by default had a set of hard coded incoming rays, but this could be changed to have higher resolution spacing immediately before the horizon. HR by default has higher resolution incoming rays both above and below the horizon.

n Setting
0 Incoming rays have high resolution both above and below the horizon
1 Incoming rays are hard coded to the SASKTRANV21 values
2 Incoming rays have high resolution only below the horizon


NumDiffuseProfilesInPlane (int n) [Default: n=1]

Sets the number of diffuse profiles in either the solar zenith angle or line of sight plane, depending on the setting of DiffusePlacementType. See Diffuse Profiles for more information.


NumDiffuseOutgoing (int n) [Default: 169]

Sets the number of outgoing rays on diffuse spheres. n must be a perfect square and is limited to up to approximately 1000.


DiffuseIncomingResolution (array v) [Default: See Text]

Sets the incoming ray resolution for diffuse spheres

Array Index Setting
0 Number of zenith discretizations before (above) the horizon
1 Number of zenith discretizations in the horizon region
2 Number of zenith discretizations below the horizon (towards the ground)
3 Number of azimuth discretizations in all regions

Default value is v=[6, 8, 10, 12].


ManualDiffuseHeights (array v) [Default: See Text]

Sets the altitude for diffuse points. This setting applies to every profile. It is recommended to place a diffuse point just above the ground. The default setting for this is:

v = np.concatenate(([10^-6], np.linspace(500, 99500, 100))]


NumDiffusePlanes (int n) [Default: n=1] [EXPERIMENTAL]

Diffuse profiles can also be placed outside the line of sight plane. This option sets the number of planes to create, they are spacing based off of DiffuseMaxAngleOffPlane.


DiffuseMaxAngleOffPlane (double d) [Default: d=5.0] [EXPERIMENTAL]

Sets the maximum angle in degrees that diffuse planes are placed off of the line of sight plane.


HorizonSize (double n) [Default: n=20]

Changes the horizon size in degrees for the incoming resolution on diffuse points. HR places more incoming rays above and below the horizon.

Adaptive Integration

Properties that set the adaptive integration parameters.


IntegrationTechnique (int n) [Default: n=1]

HR has added a new line integration technique where the source function is quadratically splined over each ray. Cells with areas of large extinction are also recursively split through an adaptive integration procedure. This option allows the integration improvements to be turned off, reverting to the old technique of assuming the source function is constant within each cell. Turning this improvement off offers a small performance increase, however unless necessary it is not recommended to do so.

n Setting
0 Uses the legacy SASKTRANV21 integration technique where the source function is assumed to be constant inside each cell. Adaptive integration is disabled.
1 The source function is quadratically splined along the line of sight. Adaptive integration is enabled.


MaxOpticalDepthOfCell (double d) [Default: n=1000000]

HR supports adaptive splitting of cells in the model when certain conditions are met. For a cell to be split, the optical depth of the cell must be less than MaxOpticalDepthOfCell and the ratio of scattering extinction at the start of the cell to the end of the cell must be less than MinExtinctionRatioOfCell. Note that the high default value means the by default, adaptive integration is essentially disabled.


MinExtinctionRatioOfCell (double d) [Default: d=0.9]

HR supports adaptive splitting of cells in the model when certain conditions are met. For a cell to be split, the optical depth of the cell must be less than MaxOpticalDepthOfCell and the ratio of scattering extinction at the start of the cell to the end of the cell must be less than MinExtinctionRatioOfCell.

Ray Tracing Properties

Properties that affect how the engine traces rays through the atmosphere.


UseShellRayTracer (int n) [Default: n=0]

The HR engine uses a new ray tracer which calculates intersections with a set of arbitrary geometry primitives (for one dimensional atmospheres this is a set of spherical shells). The old ray tracer from the SO engine can be used by calling this option with n=1, however there is no reason to.

n Setting
0 Uses the standard generic ray tracer
1 Uses the old shell ray tracer


SolarRayTracingShells (double d) [Default: d=1000]

Sets the vertical shell spacing used for tracing solar rays in [m]. This setting should be changed if higher vertical resolution calculations are desired.


RayTracingShells (double d) [Default: d=1000]

Sets the vertical shell spacing used for tracing rays other than the solar rays in [m]. This setting should be changed if higher vertical resolution calculations are desired.

Optical Table Properties

Properties that affect the internal optical property table.


OpticalTableType (int n) [Default: n=0]

Sets the type of atmosphere used internally in the model. This does not change how the user specifies the atmosphere through climatologies, this setting only affects how the model handles the atmosphere it is given.

n Setting
0 One dimensional atmosphere in height
1 Three dimensional atmosphere using a Delaunay triangulation.
2 Two dimensional atmosphere in height and angle along the line of sight plane.
3 Two dimensional atmosphere in solar zenith angle and height [EXPERIMENTAL]


ManualOpticalHeights (array v) [Default: See Text]

Sets the heights in the model where cross sections and number densities are specified. Typically this is twice the resolution of the ray tracing grid. The default value is v = linspace(0, 100000, 201).


OpticalNormalAndReference (array v) [Default: See Text]

Only has an effect when OpticalTableType is 2. v is a 6 element array where v[0:3] is the normal and v[3:6] is the reference. The normal vector is a unit vector normal to the optical table plane. The reference is the x axis of the plane, and the normal cross the reference is the y axis.


OpticalAngleGrid (array v) [Default: See Text]

Sets the angular grid used when OpticalTableType is set to 2 in degrees. The default value is v = linspace(-10, 10, 30)


ThreeDOpticalTableParam (array v) [Default: See Text]

When OpticalTableType is set to 1, a Delaunay triangulation is constructed around the reference point creating the three dimensional optical properties table. The construction is done through placing v[1] concentric cones around the tangent point spaced v[0] degrees apart. Then, v[2] points are placed uniformly on each cone. The default value is v=[1, 10, 10].


ForceOpticalCacheUpdates (int n) [Default: n=0] [EXPERIMENTAL]

Debugging option. If set to 1 the engine forces all climatologies and optical properties to update their caches every time the model is evaluated. This is useful for debugging purposes but will slow the code down substantially.

Miscellaneous Properties

Miscellaneous properties


NumOrdersOfScatter (int n) [Default: n=50]

Sets the number of scatter orders to n. For more details see Number of Scatter Orders.


NumThreads (int n) [Default: n=0]

The number of threads the model will run on. A setting of n=0 indicates to run on as many threads as the machine has logical cores.

Deprecated Properties

A list of properties that have been deprecated. The properties are still supported but users are urged to use alternatives


GroundShiftAlt (double d) [Default: d=0] [DEPRECATED]

This name is deprecated. Please use SurfaceHeight instead.


A.E. Bourassa, D.A. Degenstein, E.J. Llewellyn, SASKTRAN: A spherical geometry radiative transfer code for efficient estimation of limb scattered sunlight, Journal of Quantitative Spectroscopy and Radiative Transfer, Volume 109, Issue 1, January 2008, Pages 52-73, ISSN 0022-4073,
Zawada, D. J., Dueck, S. R., Rieger, L. A., Bourassa, A. E., Lloyd, N. D., and Degenstein, D. A.: High resolution and Monte Carlo additions to the SASKTRAN radiative transfer model, Atmos. Meas. Tech. Discuss., 8, 3357-3397, doi:10.5194/amtd-8-3357-2015, 2015.