Contents

• Radiative transfer in moving media
• 1 Introduction
• 2 The equation of the radiative transfer
• 3 Method of discontinuous finite elements
• 4 One-dimensional planar equation in moving media
• 5 Test calculations
• 6 Conclusions
• References
• Online Material

A&A 401, 419-428 (2003)
DOI: 10.1051/0004-6361:20030144

Radiative transfer in moving media

I. Discontinuous finite element method for one-dimensional atmospheres

D. Korcáková^1,2 - J. Kubát²

1 - Ústav teoretické fyziky a astrofyziky PrF MU, Kotlárská 2, CZ-611 37 Brno, Czech Republic
2 - Astronomický ústav, Akademie ved Ceské republiky, CZ-251 65 Ondrejov, Czech Republic

Received 21 May 2002 / Accepted 17 January 2003

Abstract
We describe the solution of the one-dimensional radiative transfer equation using the discontinuous finite element method. Tests of this method prove its reliability for static media. The static solution is then applied to the solution of the transfer equation in a moving medium using a method which we call here the Local Lorentz Transformation (LLT) method. This new method is applicable also to non-monotonic velocity fields. Direct solution of the radiative transfer equation in the moving medium using the discontinuous finite element method is also presented for cases both with and without inclusion of the effect of aberration. All three techniques are compared and they give comparable results.

Key words: radiative transfer - methods: numerical

1 Introduction

The solution of the equation of radiative transfer belongs to the most important problems in astrophysics. There are many methods and computer codes for static media, but there are also many problems that do not allow ignoring of the motion of a medium, like studying atmospheres of hot stars with fast stellar winds in O stars or Wolf-Rayet stars, expanding nebulae, novae, supernovae, extragalactic radio sources, and fast changes in solar prominences. Consequently, the problem of the transfer of radiation is more complicated. The reasons for this are the coefficients of emissivity and opacity, which in the case of differential motion are angle dependent.

For continuum radiation we can simply use the static equation of radiative transfer, since the opacity and emissivity coefficients change only a little within the region of the Doppler frequency shift. However, this is not true for lines. Changes of the opacity caused by the Doppler shift may be large, and we have to take it into consideration. Generally, there are two basic classes of methods for its solution, namely the solution in the observer frame (usually a frame connected with the center of a star) and the solution in the comoving frame, which is the frame connected with (usually outflowing) material. A basic review and discussion of the applicability of both classes of methods can be found in Mihalas & Kunasz (1986).

Solution of the radiative transfer equation in the comoving frame is a commonly used possibility. It takes advantage of the fact that the opacity and emissivity coefficients are not affected by motion and may be treated as in the static case (i.e., if they are isotropic for the static case then they are isotropic also in a local reference frame in a moving atmosphere). However, the transfer equation becomes more complicated due to an additional term containing the frequency derivative of the intensity. This term describes the Doppler shift of radiation. The comoving-frame transfer equation is usually solved using Feautrier variables, which introduce additional restrictions to monotonic velocity fields (see Mihalas et al. 1975,1976a,1976b,1976c; Mihalas & Kunasz 1978). Feautrier variables have also been used by Hillier (1987). In addition, Feautrier variables are not usable for the solution of the relativistic transfer equation, where characteristic rays are not symmetric (Mihalas 1980) and use of short-characteristics is more appropriate (Hauschildt 1992).

Another effect, which is usually neglected, is the effect of aberration. The importance of aberration was discussed by Mihalas et al. (1976b); they concluded that for stellar winds (where $v/c \la 0.01$) the aberration and advection effects may be safely ignored. This was in principle confirmed by Hauschildt et al. (1995) who found the aberration and advection effects for novae to be negligible. On the other hand, Hauschildt et al. (1991) showed that for $v\approx 0.1 c$, omission of aberration leads to errors larger than 50% and Baron et al. (1996) found non-negligible effects of advection for velocities larger than $\sim$ $2000~{\rm km}~{\rm s}^{-1}$. Another study of the effects of aberration and advection was performed by Peraiah (1987, 1991a), who found significant differences in the mean intensity and flux for maximum velocities between  $1000~{\rm km}~{\rm s}^{-1}$and  $5000~{\rm km}~{\rm s}^{-1}$ both for planar and spherical geometry. He also obtained changes in spectral lines (Peraiah 1991a), but it is not clear what the real effect of aberration is.

In the first part of this paper we discuss briefly the one-dimensional equation of radiative transfer, which we will use in the subsequent parts. Here we use the slightly overlooked method of discontinuous finite elements, which has been introduced into astronomy by Castor et al. (1992), and Dykema et al. (1996). The solution of the static radiative transfer equation using this method is described in Sect. 3. In Sect. 4 we describe two new methods of solution of the one-dimensional equation of radiative transfer in moving media in a comoving frame based on the discontinuous finite elements method. The solution of the radiative transfer equation where the effect of aberration is not neglected is presented as well. In the last part we discuss the advantages and disadvantages of the method of discontinuous finite elements in these cases.

2 The equation of the radiative transfer

The equation of the radiative transfer is one of the most important equations which can be found in the description of the stellar atmosphere. A general form of this equation can be written as

 

$$\left( \frac{1}{c} \frac{\partial}{\partial t} + \vec{n} \cdot \n... ...(\vec{r},\vec{n},\nu,t ) - \chi(\vec{r},\vec{n},\nu,t) I(\vec{r},\vec{n},\nu,t)$$ (1)

where $I(\vec{r},\vec{n},\nu,t)$ is the specific intensity of radiation, $\vec{n}$ is the propagation vector, $\nu$ is the radiation frequency, $\eta(\vec{r},\vec{n},\nu,t)$ is the emissivity, and $\chi(\vec{r},\vec{n},\nu,t)$ is the opacity. In some cases it is possible to use a much simpler form of Eq. (1) - the one-dimensional static plane-parallel equation,

$$\mu \frac{{\rm d}I(z,\mu,\nu)}{{\rm d} z}=\eta(z,\nu)-\chi(z,\nu)I(z,\mu,\nu).$$ (2)

The situation is a bit more complicated in moving media. The reason is that emissivity and opacity of the stellar material depend in this case on an angle. We can circumvent this annoyance by using the comoving frame, where these coefficients are angle independent. If we transform Eq. (2) to the comoving frame and take into account only terms up to first order in v/c, we obtain

 

$$\left[ \left(\mu_{0}+\frac{v}{c}\right) \frac{\partial}{\partial ... ...{c}\frac{\partial v}{\partial z_{0}} \frac{\partial}{\partial \mu_{0}}- \right.$$

$$\frac{\nu_{0} \mu_{0}^{2} }{c} \frac{\partial v}{\partial z_{0}} ... ...{2}}{c} \frac{\partial v}{\partial z_{0}} \right] I_{0}(z_{0},\mu_{0},\nu_{0})=$$

$$\eta_{0}(z_{0},\nu_{0})-\chi_{0}(z_{0},\nu_{0}) I_{0}(z_{0},\mu_{0},\nu_{0}).$$ (3)

The subscript ₀ indicates quantities expressed in the comoving frame. A detailed derivation of Eq. (3) can be found, for example, in the textbook by Mihalas (1978, Eq. (14.129)) and we will not repeat it here.

If the velocity of the medium is not very high with respect to the observer frame, we can neglect aberration, because the influence of the Doppler effect is much stronger. In this case, Eq. (3) reduces to (see Mihalas 1978, Eq. (14.99))

 

$$\mu_{0}\frac{\partial I^{0}(z_{0},\mu_{0},\nu_{0})}{\partial z_{0... ...= \eta^{0}(z_{0},\nu_{0})-\chi^{0}(z_{0},\nu_{0}) I^{0}(z_{0},\mu_{0},\nu_{0}).$$ (4)

We have to supply the Eqs. (3) and (4) with boundary conditions.

3 Method of discontinuous finite elements

 

In the past, several classes of methods for the solution of the radiative transfer equation have been used. Possibly the most successful class of methods that have been used in a number of NLTE model atmosphere codes (e.g. Mihalas et al. 1975; Werner 1986; Hubeny 1988; Kubát 1994, 2002) is the solution of the transfer equation as a differential equation using the Feautrier (1964) method. Another class of methods is based on an integral method (see, e.g., Kalkofen 1974; Avrett & Loeser 1984).

In this paper, we deal with the in astrophysical radiative transfer slightly overlooked method of discontinuous finite elements (DFE). This method was originally developed for a study of neutron transfer by Goldin (1964). It is very similar to the finite element method, which has been used very often in various calculations in physics in the last few decades (see, e.g., Rektorys 1980). It was at first used for a solution of the radiative transfer equation by Castor et al. (1992). Further improvements and generalization of this method to more spatial dimensions were published by Dykema et al. (1996). Recently, a finite element method was used by Richling et al. (2001) for the solution of the monochromatic static radiative transfer equation in three dimensions.

The discontinuous finite element method is a generalization of the finite element method (see, e.g., Rektorys 1980). This method is based on the Galerkin method. To understand the discontinuous finite element method, we first describe briefly the Galerkin method.

Let us take a differential equation

 

Au=f, (5)

where A is a differential operator, u is the unknown function, and f is the right hand side vector. The operator A is defined on D(A), which is dense in the real Hilbert space H. We choose the finite K-dimensional subspace V_K of the Hilbert space H. We seek the approximate solution u_K in V_K in the form

$$u_{K}=\sum_{k=1}^{K} a_{k} \phi_{k}.$$ (6)

For this solution the equation

$$(Au_{K}-f,\phi_{k})=0, \hspace{1cm} (k=1\dots K)$$ (7)

must hold. Here, the round brackets mean a scalar product, which is defined as

$$(u,v)=\int u(x) v(x) {\rm d}x.$$ (8)

The functions $\phi_{k}$ are the base functions in the subspace V_K. However, since we use a finite dimensional subspace V_K instead of the real space H, the difference Au_K-f is no longer zero, it has a nonzero value called the residue (cf. Rektorys 1980). If the dimension K of the subspace V_K is large enough, the residue Au_K-f becomes negligible.

Appropriate base functions for the finite element method may be found, e.g., according to Rektorys (1980). We section the region, where we seek the solution u, to cells. We consider the base function to be nonzero only in a single cell. We choose the base functions to be as simple as possible, e.g., polynomials or splines. Let us assume that the subspace in a particular cell is L-dimensional, i.e. we have L nonzero base functions w_l ( $l=1,\dots,L$) there ( L=K/number of cells). Then we solve the following system of algebraic equations for every cell,

 

$$\left(A(a_{1} w_{1}+\dots+a_{L} w_{L}\right)-f, w_{1})=0,$$

$$\dots,$$

$$\left(A(a_{1} w_{1}+\dots+a_{L} w_{L}\right)-f, w_{L})=0 ,$$ (9)

to obtain the coefficients $a_l, l=1, \dots, L$for each cell, which together with the Eq. (6) describe the solution.

Contrary to the above described Galerkin method, in the discontinuous finite element method we do not neglect the residue  $\Delta u_{K}$

$$\Delta u_{K}= Au_{K}-f .$$ (10)

Instead of (9) we obtain the system of equations

 

$$(A(a_{1} w_{1}+\dots+a_{L} w_{L})-f, w_{1})=(\Delta a_{1}, w_{1}),$$

$$\dots,$$

$$(A(a_{1} w_{1}+\dots+a_{L} w_{L})-f, w_{L})=(\Delta a_{L}, w_{L}) .$$ (11)

Here, the $\Delta a_{l}, l=1\dots L$ are the components of the vector of the residue  $\Delta u_{L}$ in the current cell. The vector of residue can be interpreted as the difference between solution from the previous cell and the current cell. We solve the system (11) of algebraic equations for every cell and obtain the coefficients a_L, which is the solution in the grid points. This technique allows us to choose the simplest base functions and reduce the number of cells compared to the original finite element method. Although the obtained solution from one cell at the given grid point is not the true solution, it is very close to it. To smooth the solution we perform the average of the solutions from all cells that meet in the given grid point. For a solution of radiative transfer problems in stellar atmospheres it is preferable to use the DFE method rather than the finite element method, since basic physical quantities like opacity, emissivity, and density vary by several orders of magnitude throughout the atmosphere.

   
3.1 One-dimensional static plane-parallel equation

In this section we describe the solution of Eq. (2). The method basically follows Castor et al. (1992). We present this simple case here, because we shall use this method later in a more complicated situation of moving atmosphere.

First, we introduce the discretization of the geometrical depth $z=\left\{z_{d}\right\}, (d=1,\dots,D)$. Note that the subscript d=1 corresponds to the deepest point. We assume a linear dependence of the intensity on the depth within each cell (see Fig. 1).

Figure 1: Schematic plot of intensity in one-dimensional method of discontinuous finite elements.

The intensity at a depth point $z \in \langle z_{d},z_{d+1} \rangle$ can be expressed as

$$I(z)=\frac{z_{d+1}-z}{z_{d+1}- z_{d}} I_{1,d} +\frac{z-z_{d}}{z_{d+1}- z_{d}} I_{2,d+1},$$ (12)

so the intensity in the interval $\langle I_{1,d}, I_{2,d+1} \rangle$is a linear combination of space functions with coefficients equal to the boundary intensities I_1,d and I_2,d+1. We denote these functions of z as

 

$$w_{1,d}(z)=\frac{z_{d+1}-z}{z_{d+1}-z_{d}}$$

$$w_{2,d+1}(z)=\frac{z-z_{d}}{z_{d+1}-z_{d}}\cdot$$ (13)

In addition, we assume that all terms in the Eq. (12) change as the intensity, i.e.

 

$$\eta(z) = w_{1,d}(z) \eta_{d} + w_{2,d+1}(z) \eta_{d+1}$$ (14a)

$$\chi(z) I(z) = w_{1,d}(z) \chi_{d} I_{1,d} + w_{2,d+1}(z) \chi_{d+1} I_{2,d+1}.$$ (14b)

Inserting Eqs. (12) and (14) into the Eq. (2) we obtain

 

$$\mu \left(I_{1,{\rm d}}\frac{{\rm d} w_{1,d}}{{\rm d} z_{d}} + I_... ...1,d}+ \eta_{d+1}w_{2,d+1}- \chi_{d}I_{1,d}w_{1,d}-\chi_{d+1}I_{2,d+1}w_{2,d+1}.$$ (15)

Now we have to perform a scalar product of the Eq. (15) with the test functions (13). We first multiply it by a corresponding w and we obtain a system of equations

 

$$\mu I_{1,d}w_{1,d}\frac{{\rm d} w_{1,d}}{{\rm d} z_{d}}+ \mu I_{2,d+1}w_{1,d}\frac{{\rm d} w_{2,d+1}}{{\rm d} z_{d}}$$

$$=\eta_{d}w_{1,d}w_{1,d}+\eta_{d+1}w_{1,d}w_{2,d+1}$$

$$-\chi_{d}I_{1,d}w_{1,d}w_{1,d}-\chi_{d+1}I_{2,d+1}w_{1,d}w_{2,d+1}$$ (16a)

$$\mu I_{1,d}w_{2,d+1}\frac{{\rm d} w_{1,d}}{{\rm d} z_{d}}+ \mu I_{2,d+1}w_{2,d+1}\frac{{\rm d} w_{2,d+1}}{{\rm d} z_{d}}$$

$$=\eta_{d}w_{1,d}w_{2,d+1}+ \eta_{d+1}w_{2,d+1}w_{2,d+1}$$

$$-\chi_{d}I_{1,d}w_{1,d}w_{2,d+1}-\chi_{d+1}I_{2,d+1}w_{2,d+1}w_{2,d+1}$$

$$-\left(I_{1,d}-I_{2,d}\right)w_{2,d+1}.$$ (16b)

In the last parenthesis of the Eq. (16b) we introduced the jump of intensity as discussed in the previous section (cf. Eq. (11)). Note that the vector of residue is nonzero only for the second of the Eqs. (16), because we choose the difference between the solution in the previous and the current cell to be the vector of the residue. The difference between the current and following cells is considered in the following cell. We integrate these equations over the cell and obtain a system of linear algebraic equations

$${\mathcal A}\vec{I}=\vec{b},$$ (17)

where

 

$${\mathcal A} \left(\begin{array}{cc} -\displaystyle\frac{\mu}{2}+\frac{\Delta ... ... \displaystyle\frac{\mu}{2}+\frac{\Delta z_{d}}{3}\chi_{d+1} \end{array}\right)$$ = (18)

$$\vec{b} \left(\begin{array}{c} \displaystyle\eta_{d}\frac{\Delta z_{d}}{3... ...a_{d}\frac{\Delta z_{d}}{6}+\eta_{d+1}\frac{\Delta z_{d}}{3} \end{array}\right)$$ =

and

$$\vec{I}= \left(\begin{array}{c} I_{1,d}\\ I_{2,d+1} \end{array}\right)\cdot$$ (19)

Note the presence of the jump in intensities in the expression for the matrix  $\mathcal A$ and vector $\vec{b}$ (Eq. (18)). The system (17) must be solved for each cell to obtain intensities I_1,d and I_2,d+1. The resulting intensity at a depth point d is then

$$I_{d}=\frac{I_{1,d}+I_{2,d}}{2}\cdot$$ (20)

To complete the solution of the transfer equation we must supply the boundary condition to specify the radiation that enters the medium. For a semi-infinite stellar atmosphere for upward direction we may use the diffusion approximation as the lower boundary condition,

$$I_{2,1}=B_{\nu}(z=R_{0})+\left.\mu\frac{\partial B_{\nu}}{\partial z} \right\vert _{z=R_{0}}$$ (21)

where $B_{\nu}$ is the Planck function. For the opposite downward direction, we may assume no incident radiation at the upper boundary,

I_1,D = 0. (22)

In order to check the DFE method, we compared results of the solution of the radiative transfer Eq. (2) using DFE method with the frequently used solution using Feautrier variables (Feautrier 1964, see also Mihalas 1985). Since the DFE method, unlike the Feautrier method, is a first order method, the easiest way of comparison is using the Feautrier variable u (the average of the specific intensities in opposite directions). Results of the comparison are in Fig. 2, and we see that solutions found by both methods match perfectly. The maximum relative difference between these two methods is less than 2% in continuum and less than 1% in the center of the line.

Figure 2: Comparison of the Feautrier method (dashed line) and the method of discontinuous finite elements (solid line) for a static planar atmosphere.

4 One-dimensional planar equation in moving media

 

In this section we describe the application of the method of discontinuous finite elements to the solution of the transfer equation with velocity fields, Eqs. (3) and (4).

   
4.1 Solution using the Local Lorentz transformation

This method is based on the fact that when we solve a radiative transfer equation in a medium, which moves with a constant uniform velocity with respect to a static observer, the static solution of the radiative transfer equation is fully sufficient. This is obviously not true if there is a velocity gradient present in the medium. However, if we split the differentially moving medium into cells with constant uniform velocities moving with respect to each other, the transfer equation inside these cells may still be considered a static one thus preserving the simplicity of a static radiative transfer Eq. (2) locally. We proceed as follows.

We split the region into cells as described in Sect. 3.1. In every cell we assume a constant velocity of the medium. The velocity changes only at cell boundaries. For quantities that do not change within the cell, but which change on cell boundaries (grid points) we use "half indexes'', e.g., $d+\frac{1}{2}$, which means that this cell is between grid points d and d+1. Since we consider individual cells as local inertial frames, we must transform both intensity and frequency at the boundaries between cells. Transformation laws for these quantities are

 

$$I_{1,d}=\left(\frac{\nu_{d+\frac{1}{2}}}{\nu_{d-\frac{1}{2}}}\right)^{3} I_{2,d}$$ (23)

$$\nu_{d+\frac{1}{2}}=\nu_{d-\frac{1}{2}}\left(1-\mu \frac{\Delta v_{d}}{c}\right) \left(1-\frac{\Delta v_{d}^{2}}{c^{2}}\right)^{-\frac{1}{2}},$$ (24)

$\Delta v_{d} =v_{d+\frac{1}{2}}-v_{d-\frac{1}{2}}$ is the relative velocity of moving neighboring cells. If the velocity of the medium is not relativistic, we can set  $(\nu_{d+\frac{1}{2}}/\nu_{d-\frac{1}{2}})^{3}\approx1$.

In practice we use the same frequency points as in the static case. We solve the Eq. (2) for each cell for these frequency points. At the cell boundary, we calculate the Doppler shift of frequency, Eq. (24), and interpolate the intensity I_2,d (Sect. 3.1 and Fig. 1) to these new "rest'' frequencies. We tested both linear and parabolic polynomial interpolation and found the difference at the sixth order only. Because the equation of radiative transfer is Lorentz invariant, we solve the static Eq. (2) in each cell using the DFE method described in Sect. 3.1, i.e. we solve the Eq. (17) with $\mathcal A$ and $\vec{b}$ given by (18). The resulting intensity is given by Eq. (20). We chose the geometrical depth as an independent variable. Due to the strong dependence of the opacity on frequency in a spectral line the geometrical depth is more suitable as an independent variable than the optical depth.

We tested this method for several particular cases. Here we show the result of calculations for a stellar wind with a velocity expressed by the beta law (Cassinelli & Lamers 1999, Eq. (2.3))

$$v_{d+1/2}=v_\infty \left\{1-\left[1- \left(\frac{v_{R}}{v_\in... ...{\frac{1}{\beta}}\right] \frac{R} {z_{d}}\right\}^{\beta}\cdot$$ (25)

In Fig. 3 we show the H$\alpha$ line profile compared to a line profile in a static case. This line was computed for $v_\infty = 2000~{\rm km}~{\rm s}^{-1}$, velocity in photosphere $v_{R}=200~{\rm km}~{\rm s}^{-1}$ and parameter $\beta=1$. The line looks as if it were only shifted in frequency. In fact, its shape has also changed, but this effect is almost negligible due to the wind parameters which we chose (which correspond to the thin wind of Be stars).

Figure 3: Comparison of H$\alpha$ line profiles calculated for the static case using the simple DFE method (dashed line) and for a stellar wind with $v_\infty = 2000~{\rm km}~{\rm s}^{-1}$ and $v_{0}=200 ~{\rm km}~{\rm s}^{-1}$ using the Local Lorentz Transformation together with DFE method (full line).

This method is most effective in a case where the velocity gradient is not "too large''. We simply have to guarantee that the sets of new and old frequency points at a cell boundary significantly overlap. More exactly, the frequency shift must be smaller than the Doppler halfwidth. With increasing velocity gradient the number of cells has also to increase and, consequently, one needs more computing time to solve the equation. However, for "very'' large velocity gradients we may use the Sobolev method.

The local Lorentz transformation method is extremely useful for media with small and medium velocity gradients where it enables easy and fast implementation of the velocity field while preserving the basic simplicity of a static formulation of the radiative transfer equation.

   
4.2 Solution of the transfer equation in a comoving frame without aberration

Figure 4: The diagram of the cell for computation of the intensity.

Now we come to the direct solution of the equation of radiative transfer (4)

 

$$\mu_{0}\frac{\partial I^{0}(z,\mu_{0},\nu_{0})}{\partial z}- \fra... ...tial \nu_{0}} = \eta^{0}(z,\nu_{0})-\chi^{0}(z,\nu_{0})I^{0}(z,\mu_{0},\nu_{0})$$ (26)

using the method of discontinuous finite elements. Since this equation is valid for non-relativistic velocities, $z \simeq z_{0}$, we may omit the subscript ₀ at z. We discretize the geometrical depth $z=\left\{z_d\right\}, d=1,\dots,D$ and the frequency $\nu_0=\left\{\nu_{0,n}\right\}, n=1,\dots,N$. The intensity inside every cell is expressed as a linear combination of intensities in grid points and space functions, which are determined as follows. In every cell we place a plane, which passes through the intensity grid points of the cell (see Fig. 4). We can choose the space functions as

$$\begin{array}{ll} T_{d}=\displaystyle\frac{z_{d+1}-z}{\Delta ... ...style\frac{\nu_{0}-\nu_{0,n}}{\Delta \nu_{0,n}} \end{array}$$

where $\Delta z_{d}=z_{d+1}-z_{d}$ and $\Delta \nu_{0,n}=\nu_{0,n+1}-\nu_{0,n}$. In the following we use a simplifying notation

$$T_{1}=T_{d} \quad\quad T_{2}=T_{d+1} \quad\quad N_{1}=N_{n} \quad\quad N_{2}=N_{n+1}.$$

We can express the intensity in a point in the given cell as follows (see Fig. 4). First, we write an expression for the intensity  $I_1^\prime$between I₁₁ and I₂₁,

$$I_1^\prime=I_{11}\frac{z_{d+1}-z}{\Delta z_{d}}+ I_{21}\frac{z-z_{d}}{\Delta z_{d}} =T_{1}I_{11}+T_{2}I_{21}.$$ (27)

Similarly, the expression for $I_2^\prime$ between I₁₂ and I₂₂ is

$$I_2^\prime=I_{12}\frac{z_{d+1}-z}{\Delta z_{d}}+ I_{22}\frac{z-z_{d}}{\Delta z_{d}} =T_{1}I_{12}+T_{2}I_{22}.$$ (28)

The intensity $I_{z,\nu_0}$ between $I_1^\prime$ and $I_2^\prime$ is

$$I_{z,\nu_0}=I_1^\prime\frac{\nu_{0,n+1}-\nu_{0}}{\Delta \nu_{... ...nu_{0,n}}{\Delta \nu_{0,n}} =N_{1}I_1^\prime+N_{2}I_2^\prime.$$ (29)

Combining Eqs. (27)-(29) we obtain

$$I_{z,\nu_{0}}=T_{1}N_{1}I_{11}+T_{2}N_{1}I_{21}+ T_{1}N_{2}I_{12}+T_{2}N_{2}I_{22}.$$ (30)

For easier description we identify

$$\begin{array}{ll} I_{11}\longrightarrow I_{1}&\hspace{12mm} ... ...}&\hspace{12mm} \quad I_{22}\longrightarrow I_{4} \end{array}$$

and

$$\begin{array}{ll} w_{1}=T_{1}N_{1}& \hspace{12mm}\quad w_{2... ...}=T_{1}N_{2}& \hspace{12mm}\quad w_{4}=T_{2}N_{2} \end{array}$$

so the expression for $I_{z,\nu_0}$ is simplified to

$$I_{z,\nu_0} =w_{1}I_{1}+w_{2}I_{2}+w_{3}I_{3}+w_{4}I_{4}= \sum_{k=1}^{4} w_{k} I_{k}.$$ (31)

We insert this expression for the intensity and similar expressions for $\eta$ and $\chi I$ (cf. Sect. 3.1) into the equation of radiative transfer (26) and we obtain for each zone

 

$$\mu_{0} \sum_{k=1}^{4} I_{k} \frac{\partial w_{k}}{\partial z}- \... ...ial \nu_{0}} = \sum_{k=1}^{4} w_{k} \eta_{k} -\sum_{k=1}^{4} w_k \chi_{k} I_{k}$$ (32)

where $\eta_{k}$ and $\chi_{k}$ are the emissivity and opacity in the grid points of the cell.

Figure 5: The scheme for the expression of the jump of the intensity at the boundary of the cell (see Eq. (34)). The thick square denotes the current cell.

We again perform a scalar product with the functions w_l, i.e. we multiply the Eq. (32) by functions w_l and integrate it over the volume of the cell,

 

$$\sum_{k=1}^{4}\left[ \mu_{0} I_{k} \iint w_{l} \frac{\partial w_{... ... w_{l} \frac{\partial w_{k}}{\partial \nu_{0}} {\rm d}z {\rm d} \nu_{0} \right]$$

$$=\sum_{k=1}^{4} \left[\eta_{k} \iint w_{l} w_{k} {\rm d}z {\rm d} \nu_{0} - \chi_{k}I_{k} \iint w_{l} w_{k} {\rm d}z {\rm d} \nu_{0} \right.$$

$$\left.+(I-I^{*}) \iint w_{l}~ {\rm d}z ~{\rm d} \nu_{0} \right].$$ (33)

The superscript ^* denotes the intensity from the previous cell. Thus, the last term in the square bracket on the right hand side represents the jump in intensity at the boundary of the cell. Its expression can be derived using Fig. 5.

$$(\vec{I}-\vec{I}^{*})= \left(\begin{array}{c} I_{1}-\left(I_... ...I_{d,n+1,2}+I_{d,n+1,4}\right)/2\\ [2mm] 0 \end{array}\right).$$ (34)

For each cell we have a system of four Eqs. (33) for a vector of unknown intensities $\vec{I}$

$$\vec{I}~=~ \left(\begin{array}{c} I_{1}\\ [2mm] I_{2}\\ [2mm]... ...,n+1})\\ [2mm] I(z_{d+1},\nu_{0,n+1}) \end{array} \right)\cdot$$ (35)

For each cell we solve the Eq. (17)

$${\mathcal A}\vec{I}=\vec{b}.$$ (36)

Detailed expressions for the matrix  $\mathcal A$ and for the vector $\vec{b}$ are in Appendix A.1 (only in the electronic version).

The solution is performed from deep layers to the stellar surface to obtain the emergent intensity and in the opposite direction for the incoming intensity. The solution starts at the first (either uppermost or lowest) depth point and the intensity for all frequencies is determined. After the intensity is known for all frequencies, we continue with the next depth point. Lower boundary conditions can be obtained, e.g., from the diffusion approximation (21) and for the upper boundary condition we can assume, for example, no incoming radiation. The boundary conditions for the frequency follow directly from this solution of transfer equation. In the first cell $\{d,1\}$ the change of opacity and emissivity caused by the Doppler effect is very small, so it is sufficient to solve the static equation there. The intensity in a grid point $\{d,n\}$ we obtain as in the previous case (20),

$$I_{d,n} = \frac{1}{4}\sum_{k=1}^4 I_{d,n,k}.$$ (37)

   
4.3 Solution of the transfer equation in a comoving frame including aberration

Figure 6: The scheme for the expression of the jump at a boundary of a cell.

In this section we describe the solution of the radiative transfer equation with the inclusion not only of the Doppler effect, but now also with aberration. We solve Eq. (3),

 

$$\left[ \left(\mu_{0}+\frac{v}{c}\right) \frac{\partial}{\partial ... ...2}-1)}{c}\frac{\partial v}{\partial z} \frac{\partial}{\partial \mu_{0}}\right.$$

$$- \left.\frac{\nu_{0} \mu_{0}^{2} }{c} \frac{\partial v}{\partial... ...{\mu_{0}^{2}}{c} \frac{\partial v}{\partial z} \right] I_{0}(z,\mu_{0},\nu_{0})$$

$$=\eta_{0}(z,\nu_{0})-\chi_{0}(z,\nu_{0})I_{0}(z,\mu_{0},\nu_{0}).$$ (38)

We omit the subscript ₀ of z as in the previous case.

We introduce grids $z=\left\{z_{d}\right\}$, $\mu_{0}=\left\{\mu_{i}\right\}$, $\nu_{0}=\left\{\nu_{0,n}\right\}$. For each cell we choose following space functions as base functions

$$\begin{array}{ll} T_{1}=\displaystyle\frac{z_{d+1}-z}{\Delta... ...playstyle\frac{\nu_{0}-\nu_{0,n}}{\Delta \nu_{0,n}} \end{array}$$

where $\Delta z_{d}=z_{d+1}-z_{d}$, $\Delta \mu_{i}=\mu_{i+1}-\mu_{i}$, and $\Delta \nu_{0,n}=\nu_{0,n+1}-\nu_{0,n}$. We express the intensity

$$I_{0}(z,\mu_{0},\nu_{0})$$ = R₁M₁N₁I₁₁₁+R₂M₁N₁I₂₁₁ + R₁M₂N₁I₁₂₁+R₂M₂N₁I₂₂₁ + R₁M₁N₂I₁₁₂+R₂M₁N₂I₂₁₂ + R₁M₂N₂I₁₂₂+R₂M₂N₂I₂₂₂.

As in the previous section we use following notation to enable easier description

$$\begin{array}{ll} I_{111}\longrightarrow I_{1}&\hspace{12mm... ...7}&\hspace{12mm}\quad I_{222}\longrightarrow I_{8} \end{array}$$

$$\begin{array}{ll} w_{1}=R_{1}M_{1}N_{1}&\hspace{12mm}\quad ... ...{2}N_{2}&\hspace{12mm}\quad w_{8}=R_{2}M_{2}N_{2} \end{array}$$

and the expression for $I_{0}(z,\mu_{0},\nu_{0})$ is similar to the case without aberration

$$I_{0}(z,\mu_{0},\nu_{0})=\sum_{k=1}^{8} w_{k}I_{k}.$$ (39)

After inserting this expression for the intensity into Eq. (38), multiplying it by functions w_l, and integrating it over the volume of the cell we obtain

 

$$\sum_{k=1}^{8} \left[ \mu_{0,k} I_{k} \int w_{l}\frac{\partial w_... ...c{\partial w_{k}}{\partial z} {\rm d}z ~{\rm d}\mu_{0} ~ {\rm d}\nu_{0} \right.$$

$$\left.+\frac{1}{c}\frac{\partial v}{\partial z}~ \mu_{0,k} \left(... ...tial w_{k}}{\partial \mu_{0}}~ {\rm d}z ~{\rm d}\mu_{0} ~ {\rm d}\nu_{0}\right.$$

$$\left.- \frac{1}{c}\frac{\partial v}{\partial z}~ \nu_{0,k}~ \mu_... ...tial w_{k}}{\partial \nu_{0}}~ {\rm d}z ~{\rm d}\mu_{0} ~ {\rm d}\nu_{0}\right.$$

$$+\left.\frac{3}{c}\frac{\partial v}{\partial z} ~ \mu_{0,k}^{2}~ I_{k} \int w_{l}~w_{k}{\rm d}z ~{\rm d}\mu_{0} ~ {\rm d}\nu_{0} \right]$$

$$= \sum_{k=1}^{8} \left[ \eta_{k} \int~ w_{l}~w_{k} {\rm d}z ~{\rm... ...i_{k} I_{k} \int ~w_{l}~w_{k} {\rm d}z ~{\rm d}\mu_{0} ~ {\rm d}\nu_{0} \right.$$

$$\left. + (I_{k}-I_{k}^{*}) \int w_{l} {\rm d}z ~{\rm d}\mu_{0} ~ {\rm d}\nu_{0} \right].$$ (40)

The last term in the Eq. (40) describes the jump at a boundary of a cell. Its expression can be derived using Fig. 6,

$$(\vec{I}-\vec{I}^{*})= \left(\begin{array}{c} I_{1}-(I_{d,i,... ..._{d,i+1,n+1,6}\\ + I_{d,i+1,n+1,8})/4\\ 0 \end{array}\right).$$ (41)

We obtained a system of eight linear algebraic equations for eight unknowns I_k that must be supplied with boundary conditions, e.g. the diffusion approximation (21) at the lower boundary and no incident radiation at the upper boundary. The solution is performed in the same manner as in the simpler case without aberration. First, the equation is solved for all frequencies, then for all angles, and finally, for all depth points. The resulting intensity in a grid point $\{d,i,n\}$ we obtain similarly as in the previous cases (20), (37),

$$I_{d,i,n} = \frac{1}{8}\sum_{k=1}^8 I_{d,i,n,k}.$$ (42)

Figure 7: Comparison of the DFE methods of solution of the radiative transfer equation (described in Sects. 3.1 and 4) for the case of a zero velocity gradient. Left panel: the dependence of the Feautrier variable u on the optical depth for a continuum frequency. The solid line (denoted as static) indicates the static solution (2), the dashed line ( LLT) denotes solution using the local Lorentz transformation, the fine dashed line ( with-aberration) means solution including aberration, Eq. (38), and the dotted line ( without-aberration) means the solution without aberration, Eq. (26). Note that all curves coincide. Right panel: the dependence of the intensity in the H$\alpha$ line on frequency. Types of lines have the same meaning as in the left panel.

5 Test calculations

5.1 The static case

Test calculations for the simplest case of the static radiative transfer equation have been described in the Sect. 3.1. Here we would like to discuss the choice of an independent variable for the solution of this static equation.

We found that it is better to solve this equation with geometrical depth as an independent variable. If we choose the optical depth as an independent variable, we encounter problems, since the models converge to a wrong solution. The reason lies in the numerical calculation of the system of algebraic equation, which we solve using the LU decomposition. Although the geometrical depth does not change too much throughout the atmosphere, the optical depth changes by several orders of magnitude and its use as an independent variable results in numerical inaccuracy.

We should stress that contrary to what we have said here, Castor et al. (1992) used the optical depth as an independent variable. However, they report a reduced stability for the Galerkin prescription of the mass matrix. They overcome the stability problems with the help of replacing the mass matrix  $\mathcal A$ of the system of algebraic Eqs. (17) by a unit matrix. By this step one sets the term Au_K-f in the Eq. (7) to be orthonormal with respect to the chosen set of the base functions $\phi_{k}$, which is not generally true. We would like to avoid such a step since in doing it we may risk that the method will converge to a wrong solution, because it will be stable and we will not be able to identify errors.

5.2 Zero velocity gradient

In order to perform a basic test of our methods in moving media, we compared the solution of the static radiative transfer equation with the solution of the transfer equation in a moving atmosphere for the special case of a zero velocity gradient. All calculations have been done using the DFE method described in Sects. 3.1 and 4. We performed a test solution for a set frequencies from the H$\alpha$ line assuming the Doppler profile. Opacities have been calculated for state parameters $n_{\rm e}$ (electron density) and T (temperature) from an LTE model atmosphere calculated with a code of Kubát (2002). This model atmosphere is presented in the Table 1 (only in electronic form). Results of these calculations are shown in Fig. 7. In the left panel there is a dependence of the Feautrier variable u on an optical depth in four cases - static equation of radiative transfer, equation written in local frequency, equation including aberration, and calculation using local Lorentz transformation. As we can see, all curves coincide, which means that the basic algorithm of the solution of the transfer equation is correct. This result is supported also by a profile of the H$\alpha$ line emerging from this model atmosphere.

5.3 Nonzero velocity gradient

In Fig. 8 we plotted H$\alpha$ line profiles for a nonzero velocity gradient. The velocity field is described by the beta law (25) with the terminal velocity $v_\infty = 2000~{\rm km}~{\rm s}^{-1}$ and the velocity in the photosphere $v_{0}=200 ~{\rm km}~{\rm s}^{-1}$. The H$\alpha$ line profile is computed for the same condition as in the previous case. Three situations described in preceding sections are shown in Fig. 8 - the local Lorentz transformation, discontinuous finite element method with aberration, and without it.

The profiles shown are in a comoving frame. For comparison, we plotted also the static solution. We can see small deviation between the LLT method and methods using the equation of the radiative transfer written in comoving frame. The reason for it is in the assumption of a constant velocity in the cell. We can influence this deviation by changing the density of the geometrical depth grid, which is discussed later. On the other hand, results obtained with the inclusion of aberration and without it are in a good agreement.

Figure 8: The dependence of the intensity on frequency in the case of a stellar wind of B stars. The solid line indicates the static solution, Eq. (2), the dashed line the solution using the local Lorentz transformation, the fine dashed line the solution including aberration, Eq. (38), and the dotted line the solution without aberration, Eq. (26).

5.4 Influence of aberration

In order to see the maximum possible effect of aberration, we solved the transfer equation for a small value of $\mu$. For this case, we obtained a small difference between the solution including aberration and without it (cf. Fig. 9) as expected. Bearing in mind that no effect of aberration was found in the previous section, we may state that in the region of validity of Eqs. (38) and (26), i.e., in non-relativistic velocity fields, we can safely neglect aberration and use the simpler Eq. (26) in accordance with previous results of Mihalas & Hummer (1976b).

5.5 The tests of the grid

The grid tests show a linear dependence of the computing time on the number of depth points (see Fig. 10, left panel). In the right panel we plotted the dependence of the accuracy on the number of depth grid points. As we can see, the optimum number of depth points is about one hundred. For this choice the solution is accurate enough and the computing time is not too high. For the velocity field, introduced before, the dependence of computing time on the number of grid points is linear, too (see Fig. 11, left panel). It is necessary to take enough frequency grid points in the dynamic case, because we solve the differential equation with respect to the variable $\nu$ now (see Fig. 11, right panel). We performed the calculations with the same number of grid points in the geometrical grid and in the frequency grid in the line (D=N=100).

Figure 9: The dependence of the intensity on frequency for a large angle $\theta$ ( $\cos \theta = \mu = 4\times 10^{-2}$). The solid line indicates the static solution, the dashed line the solution without aberration, and the dotted line the solution including aberration.

Figure 10: Tests of the depth grid. Left panel: the dependence of the computing time on the number of the depth grid points D. Linear fits of these dependences are also plotted. Right panel: the dependence of the accuracy of the solution on the number of the depth grid points D. Here, $\Delta I=(I_{D=1000}-I_{D=x})/I_{D=1000}$, where x= 50, 100, 200 and 600.

Figure 11: Tests of the frequency grid. Left panel: the dependence of the computing time on a number of frequency grid points within a line. Linear fits of these dependences are also plotted. Right panel: the comparison of the dependence of the line profile on the number of frequency grid points.

6 Conclusions

In this paper we studied an application of the discontinuous finite element method for the solution of the radiative transfer equation both for static and differentially moving media.

First, we solved the static equation. We confirmed the fact that the discontinuous finite element method is comparable to the standard Feautrier method in the static case (Fig. 2).

Using the DFE method for the static case we developed a new method for the solution of the radiative transfer equation in moving media, which uses the Lorentz transformation of intensity and frequency between subsequent cells. Since this method uses local inertial frames we call it the local Lorentz transformation method. This method takes advantage of the possibility of a discontinuity at the cell boundary, which in this case comes from the Lorentz transformation.

We solved the transfer equation using the DFE method in the comoving frame. First, we neglected aberration and then we solved it more exactly including both the Doppler shift and aberration. Results of both methods are in good agreement, see Figs. 7 and 8.

In the region of validity of these equations (in non-relativistic velocity fields) it is not necessary to account for aberration, as we expected. This enables faster calculations, since one has to solve the set of only four algebraic equations instead of eight. On the other hand, if we aim at a more precise study of the radiation field, the possibility of taking aberration into account may become very useful.

In addition, we found that it is better to solve these equations with geometrical depth as an independent variable, since choosing the optical depth as an independent variable causes stability problems.

Another important point is a sufficiently fine frequency grid to enable good line profile resolution. If we do not resolve the line sufficiently, we obtain wrong line profiles.

The method of discontinuous finite elements may be applied to the solution of rather complex radiative transfer problems both in static and moving media. We tested it on a problem of expanding atmospheres of hot B stars. However, using this method in other types of moving media is possible as well. We plan to employ our results in future studies of Be stars.

Acknowledgements

The authors would like to thank the referee for valuable comments on the manuscript. This research has made use of NASA's Astrophysics Data System. This work was supported by a grant GA CR 205/01/0656 and by projects K2043105 and AV 0Z1 003909.

References

 
• Auer, L. H., & Mihalas, D. 1970, MNRAS, 149, 65 NASA ADS  
• Avrett, E. H., & Loeser, R. 1984, in Methods in Radiative Transfer, ed. W. Kalkofen (Cambridge Univ. Press), 341  
• Baron, E., Hauschildt, P. H., & Mezzacappa, A. 1996, MNRAS, 278, 763 NASA ADS  
• Cassinelli, J. P., & Lamers, H. J. G. L. M. 1999, Introduction to Stellar Winds (Cambridge Univ. Press)  
• Castor, J. I. 1970, MNRAS, 149, 111 NASA ADS  
• Castor, J. I., Dykema, P. G., & Klein, R. I. 1992, ApJ, 387, 561 NASA ADS  
• Dykema, P. G., Klein, R. I., & Castor, J. I. 1996, ApJ, 457, 892 NASA ADS  
• Feautrier, P. 1964, C. R. Acad. Sci. Paris, 258, 3189  
• Goldin, V. Ya. 1964, USSR Comp. Math. and. Math. Phys., 4, 136  
• Hadrava, P. 1992, A&A, 256, 519 NASA ADS  
• Hauschildt, P. H. 1992, JQSRT, 47, 433 NASA ADS  
• Hauschildt, P. H., Best, M., & Wehrse, R. 1991, A&A, 247, L21 NASA ADS  
• Hauschildt, P. H., Starrfield, S., Shore, S., Allard, F., & Baron, E. 1995, ApJ, 447, 829 NASA ADS  
• Hillier, D. J. 1987, ApJS, 63, 947 NASA ADS  
• Hubeny, I. 1988, Comput. Phys. Commun., 52, 103  
• Kalkofen, W. 1974, ApJ, 188, 105 NASA ADS  
• Kubát, J. 1994, A&A, 287, 179 NASA ADS  
• Kubát, J. 2002, in Stellar Atmospheres Modelling, ed. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., submitted  
• Mihalas, D. 1978, Stellar Atmospheres, 2nd ed. (San Francisco: W. H. Freeman & Comp.)  
• Mihalas, D. 1980, ApJ, 237, 534  
• Mihalas, D. 1985, J. Comput. Phys., 57, 1 NASA ADS  
• Mihalas, D., Heasley, J. N., & Auer, L. H. 1975, NCAR-TN/STR-104, NCAR Boulder  
• Mihalas, D., & Kunasz, P. B. 1978, ApJ, 219, 635 NASA ADS  
• Mihalas, D., & Kunasz, P. B. 1986, J. Comput. Phys., 64, 1 NASA ADS  
• Mihalas, D., Kunasz, P. B., & Hummer, D. G. 1975, ApJ, 202, 465 NASA ADS  
• Mihalas, D., Kunasz, P. B., & Hummer, D. G. 1976a, ApJ, 203, 647 NASA ADS  
• Mihalas, D., Kunasz, P. B., & Hummer, D. G. 1976b, ApJ, 206, 515 NASA ADS  
• Mihalas, D., Kunasz, P. B., & Hummer, D. G. 1976c, ApJ, 210, 419 NASA ADS  
• Peraiah, A. 1987, ApJ, 317, 271 NASA ADS  
• Peraiah, A. 1991a, ApJ, 371, 683  
• Peraiah, A. 1991b, ApJ, 380, 212 NASA ADS  
• Rektorys, K. 1980, Variation methods in Mathematics, Boston, Massachussets  
• Richling, S., Meinköhn, E., Kryzhevoi, N., & Kanschat, G. 2001, A&A, 380, 776 NASA ADS  
• Werner, K. 1986, A&A, 161, 177 NASA ADS

 

Online Material

 

 

Table 1: The LTE model atmosphere parameters used in our test calculations.

d	z(d)	T(d)	$n_{\rm e}(d)$	d	z(d)	T(d)	$n_{\rm e}(d)$
	[cm]	[K]	[cm-3]		[cm]	[K]	[cm-3]
1	2.2862507E+11	10869.33	1.2087334E+08	51	2.2725908E+11	11798.41	3.5495587E+13
2	2.2859749E+11	10871.11	1.5565563E+08	52	2.2723026E+11	11914.26	4.5219122E+13
3	2.2856988E+11	10872.91	2.0052332E+08	53	2.2720105E+11	12068.55	5.7414772E+13
4	2.2854223E+11	10874.71	2.5839916E+08	54	2.2717131E+11	12270.27	7.2621213E+13
5	2.2851456E+11	10876.52	3.3305214E+08	55	2.2714089E+11	12528.18	9.1470967E+13
6	2.2848688E+11	10878.34	4.2934243E+08	56	2.2710964E+11	12849.62	1.1471892E+14
7	2.2845919E+11	10880.16	5.5353549E+08	57	2.2707735E+11	13240.32	1.4328346E+14
8	2.2843149E+11	10881.97	7.1370837E+08	58	2.2704383E+11	13703.47	1.7833501E+14
9	2.2840379E+11	10883.77	9.2026996E+08	59	2.2700886E+11	14240.48	2.2141124E+14
10	2.2837608E+11	10885.55	1.1866316E+09	60	2.2697219E+11	14851.56	2.7455037E+14
11	2.2834839E+11	10887.31	1.5300680E+09	61	2.2693356E+11	15535.93	3.4030908E+14
12	2.2832070E+11	10889.04	1.9728199E+09	62	2.2689270E+11	16290.47	4.2150812E+14
13	2.2829302E+11	10890.72	2.5435123E+09	63	2.2684940E+11	17113.95	5.2092714E+14
14	2.2826536E+11	10892.36	3.2789661E+09	64	2.2680352E+11	18003.13	6.4184323E+14
15	2.2823772E+11	10893.95	4.2265177E+09	65	2.2675499E+11	18952.82	7.8885477E+14
16	2.2821010E+11	10895.49	5.4469899E+09	66	2.2670373E+11	19961.21	9.6801721E+14
17	2.2818252E+11	10896.97	7.0184887E+09	67	2.2664966E+11	21036.43	1.1863688E+15
18	2.2815497E+11	10898.38	9.0412928E+09	68	2.2659261E+11	22192.59	1.4518842E+15
19	2.2812747E+11	10899.70	1.1644098E+10	69	2.2653238E+11	23450.71	1.7733060E+15
20	2.2810002E+11	10900.88	1.4992139E+10	70	2.2646867E+11	24833.30	2.1604747E+15
21	2.2807264E+11	10901.84	1.9297694E+10	71	2.2640111E+11	26364.98	2.6243590E+15
22	2.2804532E+11	10902.44	2.4833720E+10	72	2.2632926E+11	28071.95	3.1771988E+15
23	2.2801808E+11	10902.53	3.1951684E+10	73	2.2625257E+11	29981.45	3.8328005E+15
24	2.2799092E+11	10902.03	4.1104624E+10	74	2.2617042E+11	32123.75	4.6068802E+15
25	2.2796385E+11	10900.97	5.2877084E+10	75	2.2608201E+11	34535.64	5.5183242E+15
26	2.2793687E+11	10899.61	6.8023408E+10	76	2.2598631E+11	37271.26	6.5926417E+15
27	2.2790998E+11	10898.41	8.7516950E+10	77	2.2588178E+11	40401.01	7.8689375E+15
28	2.2788316E+11	10897.91	1.1261283E+11	78	2.2576656E+11	43881.84	9.4007882E+15
29	2.2785642E+11	10898.67	1.4492799E+11	79	2.2563988E+11	47382.47	1.1278833E+16
30	2.2782974E+11	10901.14	1.8654342E+11	80	2.2550264E+11	50646.24	1.3655364E+16
31	2.2780310E+11	10905.84	2.4013307E+11	81	2.2535622E+11	53737.15	1.6662339E+16
32	2.2777650E+11	10913.35	3.0912736E+11	82	2.2520149E+11	56787.39	2.0422591E+16
33	2.2774991E+11	10924.34	3.9791981E+11	83	2.2503875E+11	59889.69	2.5085568E+16
34	2.2772333E+11	10939.45	5.1213412E+11	84	2.2486796E+11	63106.12	3.0837464E+16
35	2.2769673E+11	10959.12	6.5897143E+11	85	2.2468887E+11	66483.65	3.7905999E+16
36	2.2767009E+11	10981.80	8.4778709E+11	86	2.2450106E+11	70063.15	4.6565851E+16
37	2.2764376E+11	10745.77	1.1152910E+12	87	2.2430399E+11	73884.18	5.7145280E+16
38	2.2761763E+11	10829.46	1.4273560E+12	88	2.2409699E+11	77987.71	7.0033853E+16
39	2.2759129E+11	10929.20	1.8242781E+12	89	2.2387926E+11	82417.70	8.5691233E+16
40	2.2756469E+11	11037.59	2.3300047E+12	90	2.2364989E+11	87221.97	1.0465732E+17
41	2.2753782E+11	11148.33	2.9755770E+12	91	2.2340781E+11	92452.26	1.2756469E+17
42	2.2751070E+11	11255.31	3.8015103E+12	92	2.2315185E+11	98162.81	1.5515582E+17
43	2.2748332E+11	11352.59	4.8608780E+12	93	2.2288064E+11	104406.52	1.8831074E+17
44	2.2745573E+11	11434.12	6.2236334E+12	94	2.2259275E+11	111227.09	2.2809544E+17
45	2.2742798E+11	11495.12	7.9815046E+12	95	2.2228665E+11	118644.74	2.7585136E+17
46	2.2740010E+11	11536.20	1.0251265E+13	96	2.2196086E+11	126638.05	3.3334260E+17
47	2.2737213E+11	11567.58	1.3173932E+13	97	2.2161412E+11	135133.44	4.0294841E+17
48	2.2734408E+11	11602.42	1.6919765E+13	98	2.2124550E+11	144013.67	4.8785993E+17
49	2.2731593E+11	11649.40	2.1701062E+13	99	2.2085454E+11	153151.40	5.9221375E+17
50	2.2728761E+11	11712.65	2.7785210E+13	100	2.2044116E+11	162451.27	7.2115492E+17

Appendix A: