%\documentstyle[11pt]{article} 
%\begin{document}
%\setcounter{page}{65}
 
\newpage

\begin{center}
\Large{\bf X. EDDINGTON FACTORS}
\addcontentsline{toc}{section}{\protect\numberline{X.}{\bf ~Eddington Factors}}
\end{center}   
 
\vspace*{5mm}
{\flushleft \large{\bf A. Plane Geometry}}\\
\addcontentsline{toc}{subsection}{\protect\numberline{X.A.}{ ~Plane Geometry}}

\vspace*{2mm}
\begin{center}
\it(1) Solution of the Transfer Equation
\addcontentsline{toc}{subsubsection}{\protect\numberline{X.A.1.}{\sl ~Solution of the Transfer Equation}}
\end{center}   

\begin{flushleft}
\rm The planar transfer equation is
$$\mu\frac{\partial I}{\partial r} = \rho(\eta - \chi I) \hspace{85mm}\rm (EF1)$$\\

\rm or \\
$$\mu\frac{\partial I}{\partial \tau} = I - S \hspace{94mm}\rm (EF2)$$\\

\vspace*{2mm}
\rm where \\

\vspace*{2mm}
$S \equiv \eta/\chi = [(\rho \chi - n_{e}\sigma_{e})B + n_{e}\sigma_{e}(cE/4 \pi)]/\rho\chi \hfill \rm (EF3)$\\

\vspace*{2mm}
\rm and \\

\vspace*{2mm}
$d\tau = - \rho \chi dr  \hfill \rm (EF4)$\\

\vspace*{2mm}
\rm Discretizing, we have\\

\vspace*{2mm}
$\Delta\tau_{k} \equiv - \rho_{k}\chi_{k}(r_{k+1} - r_{k}) = \tau_{k} - \tau_{k+1} \hfill \rm (EF5)$\\

\vspace*{2mm}
\rm and\\

\vspace*{2mm}
$S_{k}(\tau) = [S_{k}- \frac{1}{2}(dS/d\tau)_{k}\Delta\tau_{k}] + (dS/d\tau)_{k}(\tau - \tau_{k+1}) \hfill $\\

\vspace*{2mm}
$\hspace{10mm}\equiv a_{k} + b_{k}(\tau - \tau_{k+1}) \hfill (\tau_{k+1}\leq\tau\leq\tau_{k}) \hspace{5mm} \rm (EF6)$\\
\end{flushleft}
 
\vspace*{2mm}
\noindent\rm where (dS/d$\tau$) is monotonized as described in section X.C below. Then for an outgoing ray,

\begin{flushleft}
$ I^{+}_{j,k+1} = a_{k}(1 - e^{-\Delta \tau_{k} / \mu_{j}}) + b_{k}\mu_{j}\{1-[1+(\Delta \tau_{k}/\mu_{j})]e^{-\Delta \tau_{k}/\mu_{j}}\} +I^{+}_{jk}e^{-\Delta\tau_{k}/\mu_{j}}\hfill$\\
$\hfill \rm (EF7) $\\

\rm For incoming rays,\\

\vspace*{2mm}
$I^{-}_{j,k+1} = I^{-}_{j,k+1}e^{-\Delta\tau_{k}/|\mu_{j}|} \hfill $\\
$\hfill + a_{k}(1 - e^{-\Delta \tau_{k} /|\mu_{j}|}) + b_{k}|\mu_{j}|[ e^{-\Delta \tau_{k}/|\mu_{j}|}+(\Delta \tau_{k}/|\mu_{j}|)-1] \hspace{2mm} \rm (EF8) $\\
\end{flushleft}

\noindent\rm In this formulation, the intensities are naturally centered at the interfaces $k$ and $k+1$. However, in order to compute Eddington factors, we actually need intensities at cell centers. Thus in practice we perform the integration in two steps, from the first interface to the cell center, then from the cell center to the next interface. This procedure greatly increases accuracy, at negligible cost. Therefore we augment the interfacial radial shells with a second set of shells through cell centers: $r_{k+\frac{1}{2}}$$\equiv$$\frac{1}{2} (r_{k} + r_{k+1})$. We index the combined sets of shells with a new index $\ell$$ \equiv$$ 2k - 1$, with $(k = 1,$$\ \frac{3}{2},$$\ 2,$$\ \frac{5}{2},$$ \ldots,\ N+1).$

\vspace*{2mm}
\begin{center}
\it(2) Boundary Conditions
\addcontentsline{toc}{subsubsection}{\protect\numberline{X.A.2.}{\sl ~Boundary Conditions}}
\end{center}   

\begin{flushleft}
\rm At the inner boundary ($k = 1$) take\\

\vspace*{2mm}
$I^{+}_{j,1} \equiv I^{+}_{L} \hfill \rm (EF9) $\\

\vspace*{2mm}
\rm for {\tt lribc} = 1, and

\vspace*{2mm}
$I^{+}_{j,1} \equiv \frac{1}{4\pi}(cE_{1} + \frac{3}{2}F_{1}\Delta\tau_{1}) + \frac{3}{4\pi} F_{1} \mu_{j} \hfill \rm (EF10) $\\

\vspace*{2mm}
\rm for {\tt lribc} = 3.\\

\vspace*{2mm}
\rm At the outer boundary ($k = N+1$) take\\

\vspace*{2mm}
$I^{-}_{j,N+1} \equiv I^{-}_{R} \hfill \rm (EF11) $\\

\vspace*{2mm}
\rm for {\tt lrobc} = 1, and

\vspace*{2mm}
$I^{-}_{j,N+1} \equiv I^{+}_{j,N+1} \hfill \rm (EF12) $\\
\end{flushleft}

\noindent \rm for {\tt lrobc} = 2. For \it both \tt lribc \rm = 2 \it and \tt lrobc \rm = 2, one must use a global- ly consistent solution, described in section X.B.3 below.

\vspace*{2mm}
\begin{center}
\it(3) Two Reflecting Boundaries
\addcontentsline{toc}{subsubsection}{\protect\numberline{X.A.3.}{\sl ~Two Reflecting Boundaries}}
\end{center}   

\vspace*{2mm}
\noindent \rm The physical problem is to find the ougoing intensity $I^{+}_{jk}$ $(k = 1,$ $\ldots,$ $ N+1)$ and the incoming intensity $I^{-}_{jk}$ $ (k = N+1,$ $ \ldots,$ $ 1)$ subject to reflecting boundary conditions at both the inner and outer boundaries. Because of the reflections, the intensities traveling in both directions are coupled.\\

\noindent \rm From equations (EF7), (EF10), and reflecting boundary conditions we get a linear system of the following form: Starting with the $jth$ outward traveling ray, the system is bidiagonal with elements on the principal diagonal and the first lower diagonal, $(k = 1,$$ \ldots,$$ N+1)$, representing the coupling of $I^{+}_{j,k+1}$ to $I^{+}_{jk}\,,$ and a column vector on the right hand side containing the contributions of the source terms in cell $(k, k+1)$. We then extend the system from $k = N+2$ to $k = 2N + 2$ to determine the inward depth variation of the incoming intensity $I^{-}_{jk'}$ where $k'\equiv  2N + 3 - k$. This part of the system is again bidiagonal, with the principal and the first lower diagonal. At the outer boundary $(k=N+2)$ we equate $I^{-}_{j,N+1}$ to $I^{+}_{j,N+1}$. Likewise, at the inner boundary $(k=1)$ we equate $I^{+}_{j1}$ to $I^{-}_{j1}$, which introduces one additional element in the last column of the first row.\\
 
\noindent \rm This system is easy to generate and solve. Let $d_{jk}$ denote an element on the main diagonal, $c_{jk}$ an element on the lower diagonal, $e_{jk}$ an element in the rightmost column, and $s_{jk}$ an element of the source vector on the right hand side. Then

\begin{flushleft}
$k = 1$ \\
$\hspace{28mm} d_{j1} =1, e_{j1}= -1, s_{j1}= 0 \hfill \rm (EF13)$\\

\vspace*{2mm}
$k = 2, \ldots, N+1 $ \\

\vspace*{2mm}
$c_{jk} = - e^{-\Delta\tau_{k}/\mu_{j}}, d_{jk} =1, e_{jk}=\hspace{2mm} 0, s_{jk}=$ \rm sources on rhs of (EF7). \hfill (EF14)\\

\vspace*{2mm}
$k = N+2$ \\

\vspace*{2mm}
$c_{j,N+2} =\ -1, d_{j,N+2} = 1, e_{j,N+2} = 0, s_{j,N+2k} = 0$ \hfill \rm (EF15)\\

\vspace*{2mm}
$k = 2N+3-k',$ where $k' = N, \ldots, 1 $ \\

\vspace*{2mm}
$c_{jk} = - e^{-\Delta\tau_{k}/|\mu_{j}|}, d_{jk} = 1, e_{jk}=\hspace{2mm} 0, s_{jk}=$ \rm sources on rhs of (EF10). \hfill (EF16)\\

\vspace*{2mm}
\rm To solve the system, carry out the forward elimination:\\

\vspace*{2mm}
$k = 2, \ldots, 2N+1 $ \\

\vspace*{2mm}
$e_{jk} = \hspace{7mm}- c_{jk}e_{j,k-1}$ \hfill \rm (EF17)\\
$s_{jk} = s_{jk} - c_{jk}s_{j,k-1}$ \hfill \rm (EF18)\\

\vspace*{2mm}
$k = 2N+2 $ \\

\vspace*{2mm}
$d_{jk} = \hspace{3mm} 1 - c_{jk}e_{j,k-1}$ \hfill \rm (EF19)\\
$s_{jk} = s_{jk} - c_{jk}e_{j,k-1}$ \hfill \rm (EF20)\\
$I^{-}_{j1} = s_{j1}/d_{j1} = I^{+}_{j1} \hfill \rm (EF21) $\\

\vspace*{2mm}
\rm Given this globally consistent value for $I^{+}_{j1}$, we can now go back to (EF5) - (EF12) for boundary and cell center intensities.
\end{flushleft}

\vspace*{5mm}
{\flushleft \large{\bf B. Spherical Geometry}}\\
\addcontentsline{toc}{subsection}{\protect\numberline{X.B.}{ ~Spherical Geometry}}
  
\vspace*{2mm}
\noindent \rm For each radial shell $(k, k+1)$ we know the cell-center values $S_{k+\frac{1}{2}}$ and can determine the monotonized derivative $(dS/d\tau)_ {k+\frac{1}{2}}$, see section X.C. From this information we can calculate values for $S_{L}$ and $S_{R}$ at the left and right interfaces. Thus intensities are naturally centered at the interfaces $k$ and $k+1$. However, in order to compute Eddington factors, we actually need the intensities at cell centers. Thus as in the planar case we perform the integration in two steps, from the first interface to the cell center, then from the cell center to the next interface.  Therefore we augment the interfacial radial shells with a second set of shells through cell centers: $r_{k+\frac{1}{2}}$$ \equiv $$ \frac{1}{2} (r_{k} + r_{k+1})$.  We index the combined sets of shells with a new index $\ell = 2k - 1$, with $(k = 1,$ $\ \frac{3}{2},$ $\ 2,$ $\ \frac{5}{2},$ $\ldots,$ $\ N+1).$

\begin{flushleft}
\vspace*{2mm}
\rm First compute all interface values of the source function:

\vspace*{2mm}
$S_{L\ell}= S_{k+\frac{1}{2}} + \frac{1}{2} \Delta \tau_{k}\left(\frac{dS}{d\tau}\right)_{k+\frac{1}{2}} \hspace{13mm}(=  S_{Lk}) \hfill \rm (EF22)$ \\

\vspace*{2mm}
$S_{R\ell}= S_{k+\frac{1}{2}} \hfill \rm (EF23)$\\

\vspace*{2mm}
$S_{L,\ell+1}= S_{k+\frac{1}{2}} \hfill \rm (EF24)$\\

\vspace*{2mm}
$S_{R,\ell+1}= S_{k+\frac{1}{2}} - \frac{1}{2} \Delta \tau_{k}\left(\frac{dS}{d\tau}\right)_{k+\frac{1}{2}} \hspace{10mm}(= S_{Rk}) \hfill \rm (EF25)$ \\

\vspace*{2mm}
$\Delta \tau_{k} = \tau_{k} - \tau_{k+1} \hfill \rm (EF26)$\\
\end{flushleft}

\vspace*{2mm}
\noindent\rm Then using the augmented set of shells, define the impact parameters of a parallel set of rays tangent to the spherical shells as

\vspace*{2mm}
\begin{flushleft}
$ p_{j} = (j - 1) r_{1}/{\tt ncore} \hspace{7mm} (1 \leq j \leq {\tt ncore}) \hfill \rm (EF27)$\\

\vspace*{2mm}
\rm and\\
$ p_{j} = r_{j - {\tt ncore}} \hspace{20mm} ({\tt ncore} + 1 \leq j \leq {\tt nray}) \hfill \rm (EF28)$\\

\vspace*{2mm}
\rm where\\
${\tt nray} = {\tt ncore} + 2N + 1 \hfill \rm (EF29)$\\
\end{flushleft}

\vspace*{2mm}
\noindent \rm Here {\tt nray} is the total number of rays, and {\tt ncore} is the number of rays penetrating {\it inside} the first (innermost) radial shell. In general, on shell $\ell$ we have {\it core rays} for $1 \leq j$$ \leq \tt ncore $, and {\it envelope rays} for ${\tt ncore} + 1$$ \leq j $ $\leq {\tt ncore} + \ell.$

\begin{flushleft}
\vspace*{2mm}
\rm For incoming rays we have\\

\vspace*{2mm}
$ x_{j\ell\hspace{5mm}} \equiv \left(r_{\ell\hspace{3mm}}^{2} - p_{j}^{2}\right)^{\frac{1}{2}} \hfill \rm (EF30) $\\

\vspace*{2mm}
$ x_{j,\ell+1} \equiv \left(r_{\ell+1}^{2} - p_{j}^{2}\right)^{\frac{1}{2}} \hfill \rm (EF31) $\\

\rm and\\

\vspace*{2mm}
$\Delta \tau_{j\ell} = \rho_{k}\chi_{k}(x_{j\ell} - x_{j,\ell+1}) \hfill \rm (EF32)$\\
\end{flushleft}

\vspace*{2mm}
\noindent\rm where $k = (\ell+1)/2$ using {\tt FORTRAN} integer arithmetic rules. For the interval $(\ell, \ell+1)$ take

$$S_{j\ell}(t) = S_{L\ell} + \left(\frac{S_{R\ell} - S_{L\ell}}{\Delta\tau_{j\ell}}\right)t = a_{\ell}+ b_{j\ell}t \hspace{44mm} \rm (EF33) $$

\begin{flushleft}
\rm where $t \equiv (\tau_{\ell} - \tau),\ \  0 \leq t \leq \Delta \tau_{j\ell}.$ Then\\

\vspace*{2mm}
$I^{-}_{j,\ell} = I^{-}_{j,\ell+1}e^{-\Delta\tau_{j\ell}} + a_{\ell}(1 - e^{-\Delta \tau_{j\ell}}) + b_{j\ell}[1 - (1 + \Delta \tau_{j\ell})e^{-\Delta \tau_{j\ell}}] \hfill \rm (EF34) $\\

\vspace*{2mm}
\rm and\\

\vspace*{2mm}
$\mu_{j\ell} = \left(1 - p_{j}^{2} / r_{\ell}^{2}\right)^{\frac{1}{2}}, \hspace{10mm} 0 \leq \mu_{jl} \leq 1 \hfill \rm (EF35)$\\

\vspace*{2mm}
\rm The outer boundary condition at $\ell = \ell_{max}$ is\\

\vspace*{2mm}
$ I_{j}^{-} = 0, \hspace{10mm}(j = 1, ..., {\tt nray})$ \hfill \rm (EF36) \\ 

\vspace*{2mm}
\rm For outgoing rays in the interval $(\ell-1, \ell)$ take\\

$$S_{j,\ell-1}(t) = S_{R,\ell-1} + \left(\frac{S_{L,\ell-1} - S_{R,\ell-1}}{\Delta\tau_{j,\ell-1}}\right)t = a_{\ell-1}+ b_{j,\ell-1}t \hspace{17mm} \rm (EF37) $$\\

\vspace*{2mm}
\rm where $t \equiv (\tau - \tau_{\ell}),\ \  0 \leq t \leq \Delta \tau_{j,\ell-1}.$ Then\\

\vspace*{2mm}
$I^{+}_{j\ell} = e^{-\Delta\tau_{j,\ell-1}}I^{+}_{j,\ell-1} \hfill $\\
$ \hfill + a_{\ell-1}(1 - e^{-\Delta \tau_{j,\ell-1}}) + b_{j,\ell-1}[1-(1+\Delta \tau_{j,\ell-1})e^{-\Delta \tau_{j,\ell-1}}] \hspace{5mm} \rm (EF38) $\\

\vspace*{2mm}
\rm At the inner boundary ($k = 1$) take\\

\vspace*{2mm}
$I^{+}_{1} \equiv I^{+}_{L} \hfill \rm (EF39) $\\

\vspace*{2mm}
\rm for {\tt lribc} = 1, and for {\tt lribc} = 3 take\\

\vspace*{2mm}
$I^{+}_{1} \equiv \frac{1}{4\pi}(cE_{1} + \frac{3}{2}F_{1}\Delta\tau_{1}) + \frac{3}{4\pi} F_{1} \mu \hfill \rm (EF40) $\\
\end{flushleft}

\noindent \rm where $\Delta \tau_{1}$ is the optical thickness of the innermost shell $(k,$ $k+1)$ $ = (1,$ $2)$. The $\Delta\tau$ term in equation (EF40) should be very small for a star: O($T_{\scriptscriptstyle eff}^{4}/T_{1}^{4}$) where $T_{\scriptscriptstyle eff}$ is the effective temperature of the star, and $T_{1}$ is the temperature at the inner boundary.

\vspace*{5mm}
{\flushleft \large{\bf C. Monotonic Interpolation of Slopes}}\\
\addcontentsline{toc}{subsection}{\protect\numberline{X.C.}{ ~Monotonic Interpolation of Slopes}}

\vspace*{2mm}
\noindent\rm Because the variation in cell size on an adaptive grid can be extreme, the physical accuracy of the integration along rays can be improved by using monotonized slopes of the source function on the ray. Thus define

\begin{flushleft}
\vspace*{2mm}
$\Delta x_{k} \equiv \frac{1}{2}(x_{k+2} - x_{k}) \hfill \rm (EF41)$\\

\vspace*{2mm}
$\Delta q_{k} \equiv \ \ \ q_{k+1} - q_{k} \hfill \rm (EF42)$\\

\vspace*{2mm}
\rm Then take\\
$$\left(\frac{dx}{dq}\right)_{k} = \frac{1}{2}\left(\frac{\Delta x_{k-1}}{\Delta q_{k-1}} + \frac{\Delta x_{k}}{\Delta q_{k}}\right) \hspace{63mm} \rm (EF43) $$\\

\rm or\\

$$ \overline{\left( \frac{dq}{dx} \right)}_{k} = \frac{2 \Delta q_{k-1} \Delta q_{k}}{ \Delta x_{k-1} \Delta q_{k} + \Delta x_{k} \Delta q_{k-1}} \hspace {55mm} \rm (EF44) $$\\
\end{flushleft}

\vspace*{2mm}
\noindent \rm which reduces to van Leer's formula for equal step sizes. If we use (EF44) when $\Delta q_{k-1}\Delta q_{k} > 0,$ and $\overline{(dq/dx)} = 0$ otherwise, we get monotonic interpolation. For example, at interface $k+1$

\begin{flushleft}
$$q_{I, k+1} = q_{k} + \frac{1}{2} \Delta x_{k} \left(\frac{dq}{dx}\right)_{k} = q_{k} + \frac{\Delta q_{k}}{\left(\frac{\Delta q_{k}}{\Delta q_{k-1}}\frac{\Delta x_{k-1}}{\Delta x_{k}} + 1 \right)} \hspace{27mm} \rm (EF45)$$\\
\end{flushleft}

\noindent \rm The term in the denominator lies on the range $(1, \infty)$. Negative values of the factor $(\Delta q_{k} / \Delta q_{k-1})$ are excluded because the filter sets $(dq/dx)_{k}$ $=$ $0$ when $(\Delta q_{k-1} \Delta q_{k})$ $\leq$ $0$. Therefore the interface value can never exceed $q_{k} + \Delta q_{k}$.
 
\vspace*{5mm}
{\flushleft \large{\bf D. Angle Quadrature}}\\
\addcontentsline{toc}{subsection}{\protect\numberline{X.C.}{ ~Angle Quadrature}}
 
\begin{flushleft}
\vspace*{2mm}
\rm On the interval $\mu_{j} \leq \mu \leq \mu_{j+1}$ represent the angular variation of the specific intensity as\\
 
$$I(\mu) = I_{j} + \frac{(I_{j+1} - I_{j})}{(\mu_{j+1} - \mu_{j})}\,(\mu - \mu_{j}) \equiv \alpha_{j} + \beta_{j} (\mu - \mu_{j}) \hspace{25mm} \rm (EF46) $$\\
 
\rm Then\\

$$\int_{\mu_{j}}^{\mu_{j+1}}[(\alpha_{j} - \beta_{j} \mu_{j}) + \beta_{j} \mu] \mu^{n} d\mu \hspace{75mm}$$\\

$$= (\alpha_{j}-\beta_{j}\mu_{j})\,\frac{(\mu_{j+1}^{n+1}-\mu_{j}^{n+1})}{(n+1)} + \beta_{j}\,\frac{(\mu_{j+1}^{n+2} - \mu_{j}^{n+2})}{(n+2)}\hspace{42mm}$$\\

$$= \frac{(I_{j}\mu_{j+1}-I_{j+1}\mu_{j})}{(n+1)}\frac{(\mu_{j+1}^{n+1}- \mu_{j}^{n+1})}{(\mu_{j+1}-\mu_{j})} + \frac{(I_{j}\mu_{j+1}-I_{j+1}\mu_{j})}{(n+2)}\frac{(\mu_{j+1}^{n+2}- \mu_{j}^{n+2})}{(\mu_{j+1}-\mu_{j})}$$\\
\hfill\rm (EF47)

\rm Hence\\

$${\cal I}_{j0} \equiv \int_{\mu_{j}}^{\mu_{j+1}} I(\mu)\hspace{4mm}d\mu = \frac{1}{2}(I_{j} + I_{j+1})(\mu_{j+1} - \mu_{j}) \hspace{30mm} \rm(EF48)$$\\

$${\cal I}_{j1} \equiv \int_{\mu_{j}}^{\mu_{j+1}} I(\mu)\mu\hspace{2mm}d\mu = \frac{1}{2}(I_{j}\mu_{j+1} - I_{j+1}\mu_{j})(\mu_{j} + \mu_{j+1}) \hspace{31mm}$$\\

$$\hspace{39mm}+ \frac{1}{3}(I_{j+1} - I_{j})(\mu_{j}^{2} + \mu_{j}\mu_{j+1} + \mu_{j+1}^{2}) \hspace{12mm} \rm(EF49)$$\\

$${\cal I}_{j2} \equiv \int_{\mu_{j}}^{\mu_{j+1}} I(\mu) \mu^{2} d\mu = \frac{1}{3}(I_{j}\mu_{j+1} - I_{j+1}\mu_{j})(\mu_{j}^{2} + \mu_{j}\mu_{j+1} + \mu_{j+1}^{2}) \hspace{15mm}$$\\

$$\hspace{39mm}+ \frac{1}{4}(I_{j+1} - I_{j})[(\mu_{j} + \mu_{j+1})(\mu_{j}^{2} + \mu_{j+1}^{2})] \hspace{11mm} \rm(EF50)$$\\

\rm Thus the Eddington factor at the center of cell $(k, k+1)$ is\\

$${\tt fedd(k)} = f_{k} = \sum_{j = 1}^{J} {\cal I}_{\ell j0} / \sum_{\ell j = 1}^{J} {\cal I}_{\ell j2} \hspace{59mm} \rm (EF51)$$\\

\rm where $\ell = 2k$.\\
\end{flushleft}

\noindent \rm For planar geometry the code uses a fixed set of angles $\{\mu_{j}\},$ $(1$ $\leq j$ $\leq$ $J+1)$ which may be changed by the user. For spherical geometry the code uses the set of $\{\mu_{j}\},$ $(1$ $\leq j$ $\leq$ $\ell + {\tt ncore})$, induced on each shell by the tangent rays that penetrate that shell.

%\end{document}   
.
