\documentclass{article}
\usepackage{amsmath}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% ***** macro definitions etc *****
%
\bibliographystyle{alpha}

% misc text stuff
\def\code#1{{\tt #1}}			% for formatting code
\def\defn#1{``#1''}			% definition of a term
\def\eg{\hbox{eg.\hbox{}}}
\def\ie{\hbox{i.e.\hbox{}}}
\def\ahf{\code{AHFinderDirect}}		% our own name
\def\Strahlkoerper{Strahl\-k\"{o}rper}           % Minkowski's term

% get size/spacing of "++" right, cf online C++ FAQ question 35.1
\def\Cplusplus{\hbox{C\raise.25ex\hbox{\footnotesize ++}}}

% misc math mode stuff
\def\ltsim{\lesssim}
\def\gtsim{\gtrsim}
\def\tfrac#1#2{{\textstyle\frac{#1}{#2}}}
\def\dfrac#1#2{{\displaystyle\frac{#1}{#2}}}
\def\thalf{\tfrac{1}{2}}
\def\const{{\rm const}}
\def\ij{{ij}}
\def\uv{{uv}}
\def\del{\nabla}

\def\A{{\cal A}}                % 2-D (continuous) domain of angular coords

\def\I{{\text{\scriptsize I}}}                  % grid-point index
\def\J{{\text{\scriptsize J}}}                  % grid-point index
\def\K{{\text{\scriptsize K}}}                  % grid-point index
\def\M{{\text{\scriptsize M}}}                  % molecule index

\def\Jac[#1]{{\bf J} \Bigl[ #1 \Bigr]}          % discrete Jacobian for fn #1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% ***** title/author/abstract *****
%

\begin{document}
\title{The \ahf{} Thorn}
\author{Jonathan Thornburg\quad{}\code{<jthorn@aei.mpg.de>}}
%
% We want CVS to expand the Id keyword on the next line, but we don't
% want TeX to go into math mode to typeset the expansion (because that
% wouldn't look as nice in the output), so we use the "$ $" construct
% to get TeX out of math mode again when typesetting the expansion.
%
\date{$ $Id$ $}
\maketitle

\abstract{
This document describes the \ahf{} thorn.  This thorn locates
apparent horizons in a numerically computed slice using a direct
method, posing the apparent horizon equation as an elliptic PDE
on angular-coordinate space.  This is very fast, but requires a
``reasonable'' initial guess.
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

The basic algorithm of this thorn is that of
\cite{Thornburg-1996-horizon-finding}:
We pose the apparent horizon equation
\begin{equation}
H \equiv \del_i n^i + K_\ij n^i n^j - K = 0
							\label{eqn-horizon}
\end{equation}
(where $n^i$ is the outward-pointing unit normal to the horizon)
as an elliptic PDE on angular-coordinate space.

We assume that a local coordinate origin has been chosen such that
the horizon is a \defn{\Strahlkoerper}, defined by Minkowski as
   ``a region in $n$-dimensional Euclidean space containing
     the origin and whose surface, as seen from the origin,
     exhibits only one point in any direction''
(Ref.~\cite[p.~108]{Schroeder-1986-number-theory}).
Introducing generic angular coordinates $(\rho,\sigma)$, we can thus
parameterize the horizon's shape by $r = h(\rho,\sigma)$ for some
single-valued \defn{horizon shape function} $h$ defined on the
2-dimensional domain $\A$ of the angular coordinates $(\rho,\sigma)$.

Finite differencing~\eqref{eqn-horizon}, we obtain a system of
simultaneous nonlinear algebraic equations for $h$ at the angular
grid points.  We solve this system of equations by Newton's method
(or a variant with improved convergence).

This algorithm has 3 main parts:
\begin{itemize}
\item	Evaluation of the ``horizon function'' $H(h)$.
\item	Evaluation of the Jacobian $\Jac[H(h)]$ of $H(h)$.
\item	Solving the nonlinear equations $H(h) = 0$ by Newton's method
	or a variant.
\end{itemize}
(These all actually apply to the finite differenced equations.)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Notation}

We use the following notation:
\begin{center}
\begin{tabular}{ll}
$x^i \equiv (x,y,z)$
		& 3-dimensional Cartesian coordinates			\\
$ijkl$		& indices ranging over $x^i$ coordinates		\\
$r$		& the Euclidean radius $(x^2 + y^2 + z^2)^{1/2}$	\\
$y^u \equiv (\rho,\sigma)$
		& the 2-dimensional angular coordinates on $S^2$	\\
$uvw$		& indices ranging over $y^u$ coordinates		%%%\\
\end{tabular}
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Evaluating the Horizon Function}

The key problem here is that in Cactus we only know the 3-metric $g_\ij$
and extrinsic curvature $K_\ij$ on a Cartesian ($xyz$) grid.  Following
\cite{Thornburg-1996-horizon-finding}, we define the function
$F(x^i) = r - h(\rho,\sigma)$, then define an outward-pointing
normal field (in general {\em not\/} of unit norm) on the horizon,
$s_i \equiv \del_i F$.  Then by equations~\hbox{(14)} and~\hbox{(15)}
of \cite{Thornburg-1996-horizon-finding}, we have
\begin{subequations}
							\label{eqn-ABCD(s-d)}
\begin{equation}
H = \frac{A}{D^{3/2}} + \frac{B}{D^{1/2}} + \frac{C}{D} - K
					\,\text{,}	% text punctuation
\end{equation}
where the subexpressions $A$, $B$, $C$, and $D$ are given by
\begin{eqnarray}
A       & = &   {}
		- (g^{ik} s_k) (g^{jl} s_l) \partial_i s_j
		- \thalf (g^{ij} s_j) \Bigl[ (\partial_i g^{kl}) s_k s_l \Bigr]
									\\
B       & = &   (\partial_i g^{ij}) s_j
		+ g^{ij} \partial_i s_j
		+ (\partial_i \ln \sqrt{g}) (g^{ij} s_j)
									\\
C       & = &   K^{ij} s_i s_j
									\\
D       & = &   g^{ij} s_i s_j
					\,\text{.}	% text punctuation
									%%%\\
\end{eqnarray}
\end{subequations}

The problem of evaluating $H(h)$ thus reduces to that of evaluating
$s_i$ and $\partial_i s_j$, given $h$ on a $(\rho,\sigma)$ grid
and given $g_\ij$ and $K_\ij$ on an $(x,y,z)$ grid.  To do this, we
observe that given 3-D coordinates $x^i$,
\begin{eqnarray}
s_i	& \equiv &
		\del_i F
									\\
	& = &	\del_i r - \del_i h(\rho,\sigma)
									\\
	& = &	\frac{\partial r}{\partial x^i}
		-
		\frac{\partial h}{\partial y^u}
		\left. \frac{\partial y^u}{\partial x^i} \right|_{x^j}
									\\
	& = &	\frac{x^i}{r}
		-
		X^u{}_i \partial_u h
					\,\text{,}	% text punctuation
							\label{eqn-s-i(h)}
									%%%\\
\end{eqnarray}
where we define the coefficients
\begin{equation}
X^u{}_i \equiv \frac{\partial y^u}{\partial x^i}
					\,\text{.}	% text punctuation
\end{equation}
We now have
\begin{eqnarray}
\partial_i s_j
	& \equiv &
		\frac{\partial s_i}{\partial x^j}
									\\
	& = &	\frac{\partial}{\partial x^j}
		\left(
		\frac{x^i}{r}
		-
		X^u{}_i \frac{\partial h}{\partial y^u}
		\right)
									\\
	& = &	\frac{\partial (x^i/r)}{\partial x^j}
		-
		(\partial_j X^u{}_i ) \frac{\partial h}{\partial y^u}
		-
		X^u{}_i \left(
			\frac{\partial}{\partial y^w}
			\frac{\partial h}{\partial y^u}
			\frac{\partial y^w}{\partial x^j}
			\right)
					\,\text{.}	% text punctuation
									%%%\\
\end{eqnarray}
A straightforward calculation gives the first term as
\begin{equation}
\frac{\partial (x^i/r)}{\partial x^j}
	= \begin{cases}
	  \dfrac{\sum_{k \neq i} (x^k)^2}{r^3}	& \text{if $i = j$}
								\\[2ex]
	  - \dfrac{x^i x^j}{r^3}	& \text{if $i \neq j$}	%%%\\
	  \end{cases}
\end{equation}
For the second term, we define the coefficients
\begin{equation}
X^u{}_\ij \equiv \frac{\partial y^u}{\partial x^i \partial x^j}
\end{equation}
The final result is that
\begin{equation}
\partial_i s_j
	= \left.
	  \begin{cases}
	  \dfrac{\sum_{k \neq i} (x^k)^2}{r^3}	& \text{if $i = j$}
								\\[2ex]
	  - \dfrac{x^i x^j}{r^3}	& \text{if $i \neq j$}	%%%\\
	  \end{cases}
	  \right\}
	  - X^u{}_\ij \partial_u y
	  - X^u{}_i X^v{}_j \partial_\uv h
						\label{eqn-partial-i-s-j(h)}
\end{equation}

Using~\eqref{eqn-ABCD(s-d)}, \eqref{eqn-s-i(h)},
 and~\eqref{eqn-partial-i-s-j(h)}, we can evaluate $H(h)$ using
only angular-coordinate derivatives of $h$ and $xyz$-coordinate
derivatives of $g^\ij$ and $\ln \sqrt{g}$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Implementation Notes}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Overview}

The \ahf{} thorn is written primarily in \Cplusplus{}, calling
C and Fortran~77 numerical libraries.  \Cplusplus{} has proven to be
a powerful and flexible language for this type of programming, and
I think this thorn (particularly the interpatch interpolation) would
have been significantly harder to write in a lower-level language
like C or Fortran~77.

The implementation is only loosely coupled to Cactus -- in general
the code uses its own types and classes, which are implemented on
top of the Cactus ones.  Notably, all floating point arithmetic is
done using the type \code{fp}, a typedef for \code{CCTK\_REAL}.

The code should be fairly portable to modern \Cplusplus{} compilers.
Templates are used, but only template classes templated on floating-point
or integer types, and these templates are always instantiated explicitly.
\code{bool}, \code{mutable}, and \code{typename} are used, as are the
new \code{for}-loop scope rules.  The code has been compiled and run
successfully using gcc~2.95.2 on x86 GNU/Linux systems and using
Digital C~V5.6 on Digital Unix~V4.0.  The code won't compile using
badly broken compilers like the current (mid-2001) version of Microsoft
Visual \Cplusplus.)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Maple Code}

The relativity code is all machine-generated from Maple code, which
also uses a Maple preprocessor I wrote in Perl.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliography{jt}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}