Radiative intensity is defined for any direction Ω, because the angular space is continuous. In order to handle the radiative intensity equation numerically, the angular space is discretized.

The discrete ordinates method (or SN approximation) provides a discretization of angular space into n = N(N + 2) in 3D (or n = N(N + 2)/2 in 2D) discrete directions. It consists in a set of directions and quadrature weights. Several sets are available in the literature. A set should satisfy first, second, and third moments (see Ref. 20); it is also recommended that the quadrature fulfills the half moment for vectors of Cartesian basis. Since it is not possible to fulfill exactly all these conditions, accuracy should be improved when N increases.

Following the conclusion of Ref. 21, the implementation uses the LSE symmetric quadrature for S2, S4, S6, and S8. The LSE symmetric quadrature fulfills the half, first, second, and third moments.

Depending on the value of N, a set of n dependent variables has to be defined and solved for I1,  I2, …, In.

where Si is the i-th discrete ordinate, with the following boundary condition