Carlson's Quadrature Set A

When the radiative transfer equation is solved in three dimensions, it is necessary to discretize the angular dependence of the specific intensity. In principle, the continuous integral over solid angle is replaced by a discrete sum over selected rays. This is the often used technique known as the discrete ordinates method (DOM). However, there is no unique solution for the discretization of the solid angle.

The state-of-art radiative magnetohydrodynamical codes (eg. MURaM and BIFROST) use a quadrature known as Set A by Carlson (1963). The advantage of this particular quadrature is that it is invariant with respect to 90^o axis rotation.

Terminology

The quadrature is defined by directional cosines $\mu_l$ and their weights $w_l$. They are combined into sets that define rays $\Omega_m = (\mu_x, \mu_y, \mu_z)$ and their weights $w_m$. The order of quadrature $n$ specifies the number of $\mu$ levels in the [-1, 1] interval. Since the quadrature is symmetric in respect to axis, it is sufficient to find values of the angles in half of that range. Thus the number of $\mu$ levels is $l_{max} = n/2$. The number of rays per octant in this scheme is $n(n+2)/8$. Each ray is labeled by three numbers indexing $\mu$ level along one of the axes. The sum of indices of a ray is the same of all rays, $l_{max}+2$. For example, when $n=8$ and $l_{max}=4$, possible indices are 1, 2, 3 and 4. Their possible variations are those that have sum of 6. The number of these variations (and the rays) is $n(n+2)/8 = 10$. The number of different classes (variation without repetition with the specified sum) is $l_{max}-1 = 3$. In our example: (4, 1, 1), (3, 2, 1) and (2, 2, 2).

Requirements

The quadrature is constructed to fulfill the following requirements: $$\sum_{l=1}^{n/2} w_l = 1\;\;\;\;\;\;\sum_{l=1}^{n/2} w_l\mu_l = 0\;\;\;\;\;\;\sum_{l=1}^{n/2} w_l\mu_l^2 = \frac{1}{3}.$$

The actual computation of the quadrature starts from the values $\mu_1$ and $w_1$. The full derivation is available in Carlson (1963). Here I just summarize the algorithm. The directional cosines $\mu_1$ are given by: $$\mu_l^2 = \mu_1^2 + 2\frac{l-1}{n-1},\;\;\;\;\;l=1,...,n/2-1$$ where $\mu_1^2 = 1/(3(n-1))$. There is a similar equation for the weights ($W_l = w_1 + \cdots +w_l$): $$W_l^2 = W_1^2 + 2\frac{l-1}{n-1},\;\;\;\;\;l=1,...,n/2-1.$$ However, the weights have to fulfill also: $$\sum_{l=1}^{n/2-1} W_l = \frac{n-2}{3}.$$ These two equations give us the equation for $w_1$: $$\sqrt{w_1^2} + \sqrt{w_1^2 + 2\frac{l-1}{n-1}} + \cdots + \sqrt{w_1^2 + \frac{n-4}{n-1}} = \frac{n-2}{3}.$$ This can be solved by using the Euler summation formula (see Carlson, 1963). Here I solve it numerically using the Newton-Raphson method. Once $w_1$ and $\mu_1$ are known, it is trivial to find other values of $\mu_l$ and $w_l$. It still remains to find the weights of the rays. The procedure is simple. One linear equation can be written for every $\mu$ level: the sum of the weights of the rays at the level $l$ is equal to the weight of that level. The number of equations is equal to the number of levels, but the last equation is redundant. For example, for $n=8$ we write the following equations: $$l=4: w_{411} = w_4,$$ $$l=3: w_{321} + w_{312} = w_3,$$ $$l=2: w_{231} + w_{222} + w_{213} = w_2,$$ where the weights of all rays that belong to one class is same ($w_{321} = w_{312} = w_{231} = w_{213}$). The system can be solved by any of the standard methods.

IDL code

I coded the described recipe for the construction of the Set A quadrature of Carlson in IDL. The main procedure carlson_quadrature_a.pro calls two other w1nr_a.pro and variations3.pro (click on the procedure names to download them). The result of the procedure is equivalent to the numbers in Table 2 of Carlson (1963). For example,

carlson_quadrature_a, 8, wl=wl, mu=mu, wm=wm

gives the following values (Set A, n=8):

IDL> print, wl
      0.43672982      0.25352174      0.18276705      0.12698138
IDL> print, wm
     0.070754692     0.091383524      0.12698138
IDL> print, mu
      0.21821789      0.57735027      0.78679579      0.95118973
IDL> print, variations
       1       1       1       1       2       2       2       3       3       4
       1       2       3       4       1       2       3       1       2       1
       4       3       2       1       3       2       1       2       1       1

Note: The procedure in Carlson's original paper is a bit cryptic. The paper of Bruls et al (1999) is much more clear, although they take just an approximative expression for $w_1$.