%\documentstyle[11pt]{article} 

%\input{rep12}
%\input{symbollga}

%\begin{document}
%\setcounter{page}{2}
\newpage
 
\begin{center}
\Large{\bf II. RADIATION HYDRODYNAMICS}
\end{center}   
\addcontentsline{toc}{section}{\protect\numberline{II.}{\bf ~Radiation Hydrodynamics}}
 
\vspace*{5mm}
{\flushleft\large{\bf A. The Equations}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{II.A.}{~The Equations}}
 
\begin{center}
\it Reference
\end{center}
 
{\flushleft
\hspace{10mm}\rm D. Mihalas and B. W. Mihalas, ``Foundations of Radiation\\ 
\hspace{15mm}Hydrodynamics'', Oxford University Press, New York, 1984}\\
\vspace{2mm}
 
\noindent\rm 
The equations of Newtonian gray radiation hydrodynamics for a single fluid in the 
comoving-frame of reference taking into account terms of $O\Bigl({v\over c}\Bigr)$ are
formulated according to Mihalas \& Mihalas. The two meaningful geometries for radiation
transfer are slab symmetry, designated by the parameter $\mu = 0$, and spherical symmetry,
$\mu = 2$. The hydrodynamical equations by themselves can also be solved in cylindrical 
symmetry, $\mu = 1$.\vspace*{2mm}

\noindent\rm 
The mass within a shell of radius $r$ (a slab of distance $r$) is given by
by the density:
$$
      m(r) = \int^r \rho(x) {d(x^{\mu+1})\over{\mu+1}}             \eqno{(1)}
$$
where $\rho$ is the density of the gas. The set of governing equations are, defining 
the comoving time derivative as $\Dt\equiv {\p\overp t} + u{\p\overp r}$:

\begin{description}
\item~{\sl continuity}
$$
      \Dt(\rho) + \rho {\p(r^\mu u)\over r^\mu\p r} = 0            \eqno{(2)}
$$

\item~{\sl momentum}
$$
     \rho\Dt(u) + {\p(\Pg)\overp r}
                + {\p\Bigl(r^{3\mu/2}\PQ\Bigr)\over{r^{\mu/2}\p r}}
                + \rho \Bigl[\Bigl(1-{\mu\over2}\Bigr) g
                            + 2\pi\mu {G m\over r^\mu}\Bigr]
                - {1\over c} \rho \kF \Fr = 0                      \eqno{(3)}
$$

\item~{\sl radiative momentum}
$$
     \rho\Dt\Bigl({\Fr\over{\rho c^2}}\Bigr)
                + \Bigl[{\p (\pr)\overp r} + {\mu\over2}{(3\pr - \er)\over r}\Bigr]
                + {\Fr\over c^2} {\p u\overp r}
                + {1\over c} \rho \kF \Fr = 0                      \eqno{(4)}
$$
 
\item~{\sl fluid plus radiative energy}
$$
     \rho\Dt\Bigl(\eg+{\er\over\rho}\Bigr) + {\p (r^\mu \Fr)\over{r^\mu\p r}}
                + (\Pg + \pr){\p(r^\mu u)\over r^\mu\p r}
                + \PQ \Bigl[{\p u \overp r}-{\mu\over2}{u\over r}\Bigr]
$$
$$
                + {\mu\over2}(\er - 3\pr){u\over r} = 0            \eqno{(5)}
$$

\item~{\sl radiative energy}
$$
   \rho\Dt\Bigl({\er\over\rho}\Bigr) + {\p (r^\mu \Fr)\over{r^\mu\p r}}
                + \pr{\p(r^\mu u)\over r^\mu\p r}
                + {\mu\over2}(\er - 3\pr){u\over r}
$$
$$
                + c \rho \Bigl[\kE\er-\kP\ar T^4\Bigr] = 0         \eqno{(6)}
$$
\end{description}

\noindent\rm 
In equation (5) $\eg$, the specific internal energy of the gas, and the radiative energy
$\er$, are combined to form a ``total internal'' energy. Correspondingly there are two 
different contributions to the work: that of the gas pressure~$\Pg$ and that of the 
radiative stress $\pr {\p(r^\mu u)\over r^\mu\p r} + {\mu\over2}(\er - 3\pr){u\over r}$. 
In addition we provide for the work done by an artificial stress $\PQ$. The radiant energy 
is transported via the radiative flux $\Fr$. We have neglected terms that combine radiation 
quantities with the gas acceleration.\vspace*{2mm}
 
\noindent\rm 
The gas momentum equation (3) allows for gravitational forces in two geometries, with 
$G$ standing for the Newtonian gravitational constant, and $g$ for a constant acceleration.
The radiation force, ${1\over c}\rho\kF\Fr$, the gradient of the gas pressure, and the 
artificial stress are exerted on the fluid elements, which move with velocity 
$u$.\vspace*{2mm}
 
\noindent\rm 
The equation for the radiative flux (4) can be interpreted as the photon momentum 
equation which couples to the gas momentum via the last term.\vspace*{2mm}

\noindent\rm 
The radiative energy equation (6) describes transport and work of the photons and 
constitutes, together with equation (4), the first two frequency and solid angle
integrated moments of the transport equation. The last term in (6) can be viewed
approximately as ${1\over {\tau_r}} \Bigl(\er - \ar T^4\Bigr)$, $\ar$ denoting the
radiation constant, and then be interpreted as a heat exchange on a time scale
$\tau_r \equiv {1\over {c\:\kappa\rho}}$, which is the lifetime of a photon traveling a
mean free path ${1\over{\kappa\rho}}$ between two scattering or absorption events. If
$\er > \ar T^4$ then radiative energy is converted into heat, and {\it vice versa}.
By assuming such a form for the radiative coupling term we restrict ourselves to a
radiating fluid in LTE.\vspace*{2mm}

\noindent\rm 
Equations (4) \& (6) are referred as ``full transport'' which is a correct description in 
all regimes for the one-dimensional case. In particular, the ``diffusion regime'' is recovered 
when (4) becomes equivalent to
$$
      \Fr = - {c \over 3\rho\chi}~{\p\er\overp r}.                     \eqno{(7)}
$$
Substituting formula (7) into equation (6) yields the non-equilibrium diffusion because
$\ar T_r^4 \equiv \er$ defines a radiative temperature field which is allowed to differ
from the fluid temperature $T$. Imposing the additional condition that 
$\kE\er = \kP\ar T^4$, or $T_r = T$ leads back to equilibrium diffusion.\vspace*{2mm}

\noindent\rm 
In spherical symmetry it is customary to define the luminosity as the surface 
brightness ${\cal L} = 4\pi r^2 \Fr$ of a shell with radius $r$.
 
\vspace*{5mm}
{\flushleft \large{\bf B. Adaptive Mesh Transformation}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{II.B.}{~Adaptive Mesh Transformation}}
 
\begin{center}
\it Reference
\end{center}
 
{\flushleft
\hspace{10mm}\rm Adaptive--Mesh Radiation Hydrodynamics. I. The Radiation\\
\hspace{15mm}Transport Equation in a Completely Adaptive Coordinate System,\\
\hspace{15mm}K.--H. Winkler, M.L. Norman, and D. Mihalas, \it J.Q.S.R.T.,\\
\hspace{15mm}\bf 31, \rm 473, 1984\\
\hspace{10mm}\rm W. M. Tscharnuter and K.--H. A. Winkler, \it Comput. Phys. \\
\hspace{15mm}Commun., \bf 18, \rm 171, 1979\\}
\vspace*{2mm}
 
\noindent\rm 
The mass definition (1) can be immediately rewritten in differential form as:
$$
      dm  = \rho~d\Bigl[{r^{\mu+1}\over{\mu+1}}\Bigr]     \eqno{(8)} $$\vspace*{2mm}

\noindent\rm 
The partial differential equations (2)--(6) should be represented in finite difference form 
on an arbitrary moving grid distribution. This distribution can be obtained by uniquely 
subdividing the whole computational domain into cells around each grid point $k$ with 
$k = 1, \dots N$. If (2)--(6) are integrated over those cells each encompassing a volume 
(length) $dV_k = {r^{\mu+1}\over{\mu+1}}\Big|_k$, they reduce to a set of $N$ ordinary
differential equations (in time). The proper transformation is achieved by the 
adaptive--mesh Reynolds transport theorem as described in Winkler~\etal Hereby the cells 
are assumed to move with a relative velocity $\ur|_k$ with respect to the fluid flow.\vspace*{2mm}

\noindent\rm 
For each cell $(k = 1, \dots N)$ the following equations hold:
\begin{description}

\item~{\sl continuity}
$$
   {d\over d t} \int_V \rho~dV + \oint_{\p V}~\rho\ur~dS = 0        \eqno{(9)}
$$

\item~{\sl fluid momentum}
$$
   {d\over d t} \int_V u \rho~dV - \oint_{\p V}~u \rho\ur~dS
       + \int_V r^\mu d(\Pg)
       + \int_V r^{-\mu/2} d[r^{3\mu/2}~\PQ] 
$$
$$
       + \int_V \Bigl[\Bigl(1-{\mu\over2}\Bigr) g + 2\pi\mu {G m\over r^\mu}
       - \kF {\Fr\over c}\Bigr]~\rho dV = 0                         \eqno{(10)}
$$

\item~{\sl radiative momentum}
$$
   {d\over d t} \int_V {\Fr\over c^2}~dV - \oint_{\p V}~{\Fr\over c^2}\ur~dS
       + \int_V r^\mu d(\pr)
       + \int_V~\Bigl[{\mu\over2}{(3\pr-\er)\over r}\Bigr]~dV      
$$
$$
       + \int_V {\Fr\over c^2} r^\mu d(u)
       + \int_V \Bigl[\kF {\Fr\over c}\Bigr]~\rho dV = 0            \eqno{(11)}
$$

\item~{\sl fluid plus radiative energy}
$$
   {d\over d t} \int_V \Bigl(\eg+{\er\over\rho}\Bigr) \rho~dV - \oint_{\p V}~\Bigl(\eg+{\er\over\rho}\Bigr)\rho\ur~dS
       + \int_V d[r^\mu~\Fr]
$$
$$
       + \int_V (\Pg + \pr)~d[r^\mu~u]
       + \int_V \PQ\Bigl[{\p u\overp r}-{\mu\over2}{u\over r}\Bigr]~dV
       + \int_V \Bigl[{\mu\over2}(\er-3\pr){u\over r}\Bigr]~dV = 0  \eqno{(12)}
$$

\item~{\sl radiative energy}
$$
   {d\over d t} \int_V \er ~dV - \oint_{\p V}~\er \ur~dS
       + \int_V d[r^\mu~\Fr]
$$
$$
       + \int_V \pr~d[r^\mu~u]
       + \int_V \Bigl[{\mu\over2}(\er-3\pr){u\over r}\Bigr]~dV
       + \int_V~c~\Bigl[\kE\er-\kP\ar T^4\Bigr]~\rho dV = 0         \eqno{(13)}
$$
\end{description}

\noindent\rm 
Here $\p V$ denotes the surface of the volume element $V$.
In the transformation onto the adaptive grid Neumann-Richtmyer's idea of introducing 
an artificial viscosity to treat shocks in the flow has been implemented according 
to Tscharnuter \& Winkler. A tensor formulation is used which introduces the 
artificial stress:
$$ \PQ = - {4\over3} \rho \mq \Bigl[{\p u\overp r} 
         - {\mu\over2}{u\over r}\Bigr]                              \eqno{(14)} 
$$
in conjunction with a linear and quadratic artificial viscosity, $\ellQ$ giving the 
shock width:
$$ \mq = C_{\scriptscriptstyle 1} \ellQ   \sqrt{\gamma\Pg/\rho}
       - C_{\scriptscriptstyle 2} \ellQ^2 
           \min\Bigl[~0, {\p(r^\mu u)\over{r^\mu\p r}}~\Bigr]       \eqno{(15)} $$

\vspace*{5mm}
{\flushleft \large{\bf C. Symbolic Notation}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{II.C.}{~Symbolic Notation}}
 
\begin{center}
\it Reference
\end{center}
 
{\flushleft
\hspace{10mm}\rm WH80s: Numerical Radiation Hydrodynamics, K.--H. Winkler and \\
\hspace{15mm}M.L. Norman,in \it Astrophysical Radiation Hydrodynamics,\\
\hspace{15mm}\rm (Dordrecht: Reidel), 71, 1986}\\
\vspace*{2mm}

\noindent\rm 
The volume (length) integrals in the equations above (9)--(13) can be approximated with 
finite differences in the simplest way if one assumes that the integrands remain constant. 
In general this is valid as long as the linear dimensions of the integral volumes are 
small compared to any physical length scale in the problem. In the case of smooth variations
of the physical variables this can be achieved by decomposing the physical domain into 
a sufficiently large number of numerical zones. In the case of ({\em in praxi}) 
discontinuities the demand is that the physical variables still are not allowed to vary 
significantly between neighboring zones. This can be satisfied only if the zones become 
very small; or in other words if the numerical grid {\em resolves} such nonlinear features. 
Both situations are taken care of by the implicit nature of the adaptive mesh.\vspace*{2mm}

\noindent\rm 
In a like spirit the spatial differential operators can interpreted in the simplest form 
by spatial differences denoted by $\Delta$ in symbolic form. Equation (8) valid for all 
zones $k = 1, \dots, N$ on the numerical domain then can be straightforwardly expressed as:
$$       \Delta m = \rho \Delta V                               \eqno{(16)} $$\vspace*{2mm}

\noindent\rm 
The temporal derivatives remain intact, the surface integrals are resolved as differences 
of the surface terms, and the spatial derivatives approximated as spatial differences the 
symbolic form of the equations (9)--(15) are again valid for all zones $k = 1, \dots, N$:
\begin{description}

\item~{\sl continuity}
$$
   \frac{d(\rho\Delta V)}{dt} + \Delta(\rho r^{\mu} \ur) 
          = \Delta\left(r^{\mu}\sr \frac{\Delta\rho}{\Delta r}\right)  \eqno{(17)} $$

\item~{\sl fluid momentum}
$$
   \frac{d}{dt}[u\Big<\rho\Delta V\Big>] - \Delta\left(\Big<\frac{dm}{dt}\Big> u\right) 
                   + r^{\mu}\Delta\Pg
                   + \left(\frac{2-\mu}{2}g + \frac{2\pi\mu G m}{r^{\mu}}\right) \Big<\rho\Delta V\Big> $$
$$                 = \phi_{Q}\Delta V 
                   + \frac{1}{c}<\!\!\kF\!\!>\Fr\Big<\rho\Delta V\Big> \eqno{(18)} $$

\item~{\sl radiative momentum}
$$
   \frac{d}{dt} \left[ \frac{\Fr}{c^2}\Big<\Delta V\Big> \right] 
                   - \Delta\left[\Big<\frac{dm}{dt}\Big>\left(\frac{\Fr}{c^2<\!\!\rho\!\!>}\right)\right] 
                   + r^{\mu}\left(\Delta\pr + \frac{\Fr\Delta u}{c^2}\right) 
                   + \frac{(3\pr - \er)}{r}\Big<\Delta V\Big> $$
$$                 = 
                   - \frac{1}{c}<\!\!\kF\!\!>\Fr\Big<\rho\Delta V\Big> \eqno{(19)} $$

\item~{\sl fluid plus radiative energy}
$$
   \frac{d}{dt}\left[\left(\eg +\frac{\er}{\rho}\right)\rho\Delta V\right] 
                   - \Delta\left[\frac{dm}{dt}\left(\eg +\frac{\er}{\rho}\right) 
                   - r^{\mu}\Fr\right] $$
$$                 + (\Pg+\pr)\Delta(r^{\mu}u) + (\er-3\pr)\frac{u}{r}\Delta V 
                   = \Delta\left(r^{\mu}\se\frac{\Delta\eg}{\Delta r}\right) 
                   + \epsilon_{Q}\Delta V                                      \eqno{(20)} $$

\item~{\sl radiative energy}
$$
   \frac{d}{dt} \left( \er \Delta V \right) 
                   - \Delta\left[\frac{dm}{dt}\left(\frac{\er}{\rho}\right) 
                   - r^{\mu}\Fr \right] 
                   + \pr\Delta(r^{\mu}u) + (\er - 3\pr)\frac{u}{r}\Delta V $$
$$                 = c \Bigl[\kP\ar T^4 - \kE\er\Bigr]\rho\Delta V             \eqno{(21)} $$
\end{description}

\noindent\rm 
The bracketed quantities are averaged in this staggered mesh representation, \cf~the
discussion in section III A. This set of equations corresponds to the equations $M1$, 
$C1$, $GM1$, $RM1$, $FE1$, $RE1$ in the TITAN Code Reference Manual.\vspace*{2mm}
\noindent\rm 
The original advection operator $\Delta(\rho r^{\mu} \ur)$ has been retained only in the 
continuity equation (17). In all other equations, for reasons discussed by Winkler and
Norman, the combination $\rho r^{\mu} \ur$ is replaced by $- \frac{dm}{dt}$. This is possible 
because the motion of the grid relative to the fluid results in an apparent mass transfer rate 
that can be computed according to:
$$             \frac{dm}{dt} = - \rho r^{\mu} \ur                              \eqno{(22)} 
$$
Since the correct amount of mass per zone is already determined by (16) the advection 
operator in all other ordinary differentail equations (17)--(21) can be expressed through 
$-\Delta\left(\frac{dm}{dt}\right)$.\vspace*{2mm}

\noindent\rm 
The artificial stress enters into the momentum equation (18) in the form of a force:
$$  
    \phi_Q \Delta V \equiv r^{-(\mu/2)} \Delta\left[\frac{4}{3} \rho \mq r^{(3\mu/2)} \left(\frac{\Delta u}{\Delta r} - \frac{\mu}{2}\Bigl<\frac{u}{r}\Bigr>\right)\right] \eqno{(23)} 
$$
and into the fluid plus radiative energy equation (20) in the form of a source:
$$  \epsilon_Q \equiv \frac{4}{3}\mq\rho\left(\frac{\Delta u}{\Delta r} 
                    - \frac{\mu}{2}\Bigl<\frac{u}{r}\Bigr>\right)^{2}          \eqno{(24)} 
$$
The artificial viscosity coefficient is evaluated according to:
$$  \mq \equiv C_{\scriptscriptstyle 1} \ellQ a_s
      - \min [ C_{\scriptscriptstyle 2} \ellQ^2 \Delta(r^{\mu}u) / \Delta V, 0]\eqno{(25)} 
$$
The shock width can be assumed as a constant $\ellQ = \ell_0$ or proportional to the
radius $\ellQ = \ell_1 <\!\!r\!\!>$.
Additional diffusion operators appear in the continuity and fluid plus radiative energy 
equations. Their purpose is to treat numerical artifacts that occur in contact continuities 
and when shock reflect from walls (wall heating). The constants $\sr$ and $\se$ 
control the length scale over which density or internal energy is allowed to diffuse.

\vspace*{5mm}
{\flushleft \large{\bf D. Boundary Conditions}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{II.D.}{~Boundary Conditions}}

\noindent\rm 
TITAN is equipped with a broad variety of boundary conditions, allowing for everything but 
very special-purpose boundaries. They naturally fall into two categories which describe the 
behavior of the radiation field and of the fluid ath the boundaries. The inner and outer 
(left and right) boundary conditions obey similar laws; hence is suffices to give just the 
generic formulation. More details are given in the TITAN Code Reference Manual.\\

\noindent (1) radiative boundary:\vspace{2mm}

\noindent
In general the radiative boundary can be posed by requiring the radiative energy to remain 
constant: $\Delta\er = 0$. We provide for the optically transmitting and reflecting cases. 
The first is defined through a flux which is computed from the free streaming formula. In 
addition the flux is allowed to be modified by a constant intensity field $I^\pm$. The second 
case is given through a zero flux condition. Furthermore we have the possibility of imposing 
a time dependent net flux.

\begin{description}
\item~{\sl transmitting}
$$
     \Fr = \pm c \gE \er \mp (2 \gE + 1) \pi I^\pm.         \eqno{(26)} $$
\item~{\sl reflecting~flux}
$$
     \Fr = 0.                                               \eqno{(27)} $$
\item~{\sl imposed~flux}
$$
     \Fr = \Phi(t).                                         \eqno{(28)} $$
\end{description}

\noindent
In equation (26) $\gE$ is a surface Eddington factor connecting $\er$ to $\Fr$.\\

\noindent (2) Eulerian boundary:\vspace{2mm}

\noindent
The Eulerian boundary is given by the fixed grid condition $\dot r = 0$. We allow for three
different cases. The transmitting condition demands that temperature, density, and
velocity all remain constant across the boundary. The zero flux condition necessitates
that the temperature and density stay constant while the velocity vanishes. For the nonzero
flux condition density and temperature can be functions of time.
 
\begin{description}
\item~{\sl transmitting}
$$
  \Delta T = 0~;~\Delta\rho = 0~;~\Delta u = 0.  \eqno{(29)} $$
\item~{\sl zero~flux}
$$
  \Delta T = 0~;~\Delta\rho = 0~;~u = 0.         \eqno{(30)} $$
\item~{\sl nonzero~flux}
$$
  T = \Theta(t)~;~\rho = P(t)~;~u = 0.           \eqno{(31)} $$
\end{description}

\noindent (3) Lagrangean boundary:\vspace{2mm}

\noindent
The Lagrangean boundary is given by the fixed mass condition $\dot m = 0$. We allow for two
different forcings. In the piston driven case a time dependent velocity is assumed to
move the boundary. Otherwise we allow for a forcing through an external pressure.
 
\begin{description}
\item~{\sl external velocity}
$$    
    u = U_{\scriptscriptstyle ext}(t).           \eqno{(32)} $$
\item~{\sl external pressure}
$$
  \Pg = P_{\scriptscriptstyle ext}(t).           \eqno{(33)} $$
\end{description}

\noindent\rm 
In addition to the boundaries we need to provide for phantom zones so that the difference
stencil stays well defined. These phantom zones have to be chosen to reflect the physical
properties of the boundary conditions.
 
\vspace*{5mm}
{\flushleft \large{\bf E. Constitutive Relations}}\\
\vspace*{2mm}
\addcontentsline{toc}{subsection}{\protect\numberline{II.E.}{~Constitutive Relations}}
 
\begin{center}
\it References
\end{center}
 
{\flushleft
\hspace{10mm}\rm D. Mihalas, W. D\"appen, D. Hummer, and B. W. Mihalas,\\
\hspace{15mm}\it Ap. J., \bf 331, \rm 794, 1988; \it Ap. J., \bf 331, \rm 815, 1988;\\
\hspace{15mm}\it Ap. J., \bf 332, \rm 261, 1988; \it Ap. J., \bf 350, \rm 350, 1990}\\
\vspace{2mm}

\noindent\rm (1) fluid equation of state:\vspace{2mm}

\noindent\rm 
Arbitrary equations of state as functions of the density and gas temperature are permitted:
\begin{description}
\item~{\sl caloric}
$$
   \eg = \eg(\rho,T),                 \eqno{(34)} $$
\item~{\sl thermal}
$$
   \Pg = \Pg(\rho,T).                 \eqno{(35)} $$
\end{description}

\noindent\rm 
They can be given in form of tables or as analytical formulae. In the first case 
monotonic bicubic hermite polynomials are employed to perform the interpolations. 
Important thermodynamical properties for the fluid are derived through the Gibbs form 
$d\eg(\rho,T) = T ds(\rho,T) + \Pg(\rho,T) d({1\over\rho})$, where $s$ denotes the specific
entropy. Among those of particular interest is the adiabatic index
$\gamma = {\partial \ln\Pg \overp {\ln\rho}}\Big|_s$.\vspace{2mm}

\noindent\rm 
For the purpose of the test problems a perfect gas is assumed so that:
$\Pg = {{\it R}\over\overline{\mu}}\rho T$, with ${\it R}$ standing for the universal gas
constant and $\overline{\mu}$ for the mean molecular weight. Prescribing a constant
adiabatic index $\gamma$ then relates the gas pressure to the gas energy:
$\eg = {1\over (\gamma-1)}{\Pg\over\rho}$. \\

\noindent\rm (2) radiative equation of state:\vspace{2mm}

\noindent
The two-equation formalism for the radiative transport, equations (6) \& (4), is closed
via the variable Eddington factor~$\fe$. The radiative pressure~$\pr$ in this description
is simply proportional to $\er$:
$$
   \pr = \fe \er.                      \eqno{(36)}
$$
In order to evaluate $\fe$ we integrate the static transport equation for the radiative
intensity as a function of angle $\phi$, $I(x=\cos\phi)$, using a ray tracing technique. 
The Eddington factor is then computed as the ratio of the second to the zeroth moment:
$$
    \fe = {{\int_{-1}^{+1} I(x) x^2 dx}\over
           {\int_{-1}^{+1} I(x) dx}}.  \eqno{(37)}
$$
At the boundaries of the radiation field a corresponding factor, called $\gE$, is evaluated.
$\fe$ describes the geometric properties of the radiation field. For an isotropic intensity
distribution, $I(\phi) = I_o$, we obtain $\fe = {1\over3}$. This corresponds to the
diffusion limit of the two-equation formalism. If the intensity distribution peaks along
the line of sight, \ex~$I(\phi) \propto \cos^2\phi$, we get $\fe = {3\over5} > {1\over3}$.
The stronger the intensity field is pointed along the line of sight the closer
$\fe \rightarrow 1$. This corresponds to the free streaming limit of the full transport
formalism. On the other hand, if the intensity distribution peaks sideways, \ex~
$I(\phi) \propto \sin^2\phi$, we find $\fe = {3\over15} < {1\over3}$. This situation is
often found in radiating shock fronts. \\

\noindent\rm (3) opacities:\vspace{2mm}

\noindent
In general the frequency dependent opacity $\chi(\nu)$ is given by contributions from
absorption, $\chi^a(\nu)$, and scattering, $\chi^s(\nu)$. The gray formulation of the
radiation hydrodynamics equations introduces three different opacity means: the flux mean,
the absorption (or energy) mean, and the Planck mean: \vspace*{2mm}
\begin{description}
\item {\sl flux~mean}
$$
    \kF = \int_{0}^{\infty} [\chi^a(\nu)+\chi^s(\nu)]~{\Fr(\nu)\over\Fr}~d\nu. \eqno{(38)} $$
\item {\sl absorption~mean}
$$
    \kE = \int_{0}^{\infty} \chi^a(\nu)~{\er(\nu)\over\er}~d\nu.               \eqno{(39)} $$
\item {\sl Planck~mean}
$$
    \kP = \int_{0}^{\infty} \chi^a(\nu)~{B(\nu,T)\over \sigma T^4/\pi}~d\nu.   \eqno{(40)} $$
\end{description}
The flux and the energy mean both depend on the spectral distribution of the radiation
field, whereas the Planck mean is given by the local gas temperature and density via the
Planck distribution $B$. The evaluation of the first two opacity means therefore can
be carried out correctly only if one also computes the frequency dependent, \ie~non-gray,
radiation field. This is often done with a multigroup approach. For the moment TITAN is 
restricted to the gray transport. \vspace{2mm}

\noindent\rm 
Defining the Rosseland mean opacity
$$
     \kR^{-1} = \int_{0}^{\infty} [\chi^a(\nu)+\chi^s(\nu)]^{-1}~
                              {{\p B(\nu,T)/\p T}\over 4\sigma T^3/\pi}~d\nu,  \eqno{(41)}
$$
one can show that in diffusion approximation it equals to the flux mean, \ie~$\kF = \kR$.
For the purpose of our test problems we assume a constant $\kR$. In order to simulate the
effects of a scattering or absorption dominated fluid we set both the energy mean and the
Planck mean equal to a fraction $0 \le \fr \le 1$ of the Rosseland mean,
\ie~$\kE = \kP = \fr \kR$.\\

\noindent\rm
Both the fluid equation of state and the opacities are available in table lookup form.
The source is described by Mihalas~\etal in a series of articles.
In this case monotonized bicubic hermite interpolation is used to extract the desired
quantity $x(T,\rho)$ along with its partial derivatives ${\p\ln x\overp \ln T}$ and 
${\p\ln x\overp \ln\rho}$.  
 
%\end{document}
.
