\documentstyle[12pt,epsfig]{article}

\setlength{\textwidth}{6.5in}
\setlength{\textheight}{8.5in}
%\setlength{\footskip}{0.35in}
\setlength{\topmargin}{0.0in}
\setlength{\oddsidemargin}{0in}
\setlength{\parindent}{5ex}
\setlength{\parskip}{.12in}
\newcommand{\detw}{\mbox{$\frac{\Delta E}{2}$}}
\newcommand{\hs}{\hspace*{5ex}}


%\def\decayright#1{\kern#1em\raise1.1ex\hbox{$|$}\kern-.5em\rightarrow}

\begin{document}

\begin{center}
{\Large
Summary of a Simulation of the ANITA Detector
}
~\\
Amy Connolly\\
\today\\

~\\
~\\
{\bf Abstract}\\

\parbox[t]
{5.5in}
{
\hs We outline our code written to simulate the ANITA
detection system.
}
\end{center}

\tableofcontents


\section{Introduction}

The code outlined here is meant to be a simulation of the
ANITA and ANITA-lite detection systems.  It should provide
the collaboration with a tool to study the impact of various
effects on the experiment's sensitivity.  It will also 
provide the energy-dependent sensitivity that will be used
for setting limits on models for UHE neutrinos from
ANITA-lite data and ultimately from the full ANITA flights.
It also serves as a consistency check with the Monte Carlo
simulation being developed by our collaborators at the
University of Hawaii.

\section{Overall Strategy}
\label{sec:strategy}
\hs The simulation models interactions from neutrinos
of a specified energy.
The Askaryan pulse is parameterized by the current
theoretical model.
We require an interaction to occur within
the Antarctic ice volume within the balloon's horizon.  
We generate two different types of events.  The first is where the
ray seen by the balloon is direct; the ray is emitted from the
interaction upward.  The second is where the ray is emitted downward
and is reflected from the ice-rock interface before being reaching
the surface.  For half of the neutrinos, we force the ray to be the
first type and for the other half, the second type.
For a given type of ray, we find the unique path along which
an RF signal would travel from the interaction to the balloon, snelled
through ice layers 
at the ice-air surface.  Next, we pick a direction for the neutrino
path, only considering directions such that the Cerenkov cone is
close enough to the unique ray from interaction to balloon that
the signal is still detectable under a best-case scenario.
Depth-dependent attenuation lengths and indices of refraction
in the ice are based on recent South Pole measurements.  Surface
slopeyness is taken into account under a simple model.
Frequency-dependent 
antenna response is based on the manufacturer's specifications.
The trigger configuration is based on the most
current design.
The restrictions in neutrino phase space 
that come about in our selection of neutrino interaction positions
and neutrino directions 
are corrected for in the
final calculation (see Section~\ref{sec:volstr}).


\subsection{Coordinate System}


\hs Where Cartesian coordinates are used, we define $x$,$y$ and $z$ to 
be zero at the center of the earth.  The $+z$ axis 
runs from the center of the earth
to the south pole (so, our coordinate system is defined so that
the earth is ``upside-down'').  The $x$ and $y$ axes lie in the
0$^{\circ}$ latitude plane, with $+x$ pointing to
the 90$^{\circ}$ E longitudinal line, and $+y$ points to
0$^{\circ}$ longitude.  This coordinate system is right-handed.

The $\phi$ coordinate is defined as usual, zero along the $+x$
axis, and increasing moving counter-clockwise from the view 
looking down on the
$+z$ axis.  The $\theta$ coordinate is measured relative to the
$+z$ axis.

Bedrock and ice surface elevations are measured radially from
the center of the earth, and are quoted relative to sea level.
Sea level $r(\theta)$ is latitude-dependent and is quoted as a
distance from the center of the earth.  
We use a geoid shape given by~\cite{geoid}:
\begin{equation}
r(\theta)=\frac{r_{min}\cdot r_{max}}{\sqrt{r_{min}^2-(r_{min}^2-r_{max}^2) \cos^2{\theta}}}  
\end{equation}
where $r_{min}=6356.752$ km and $r_{max}=6378.137$ km.

To accommodate the global crust model 2.0 \cite{crust2}
for thicknesses of sediment layers through the crust as well as ice
thicknesses, which is on a $2^{\circ} \times 2^{\circ}$ grid, we bin
the earth's  surface into 180 bins in longitude and
90 bins in latitude.  Note that the resulting bins are not uniform
in area.
 The bins in latitude are
numbered such that the 0$^{th}$ bin is closest to the $+z$ axis (at the
south pole), and the 89$^{th}$ bin is closest to the north pole ($-z$).


\subsection{CSEDI Crust 2.0}
\label{sec:crust2}
\hs Crust 2.0 is the latest model of the Earth's interior near the 
surface published as the result of
an initiative called Cooperative Studies of the Earth's Deep
Interior (CSEDI).  It is based on seismological data, and
the model gives thicknesses and densities of each of seven
layers in $2^{\circ} \times 2^{\circ}$ bins:  ice, water, soft sediments, hard sediments, upper crust,
middle crust, and lower crust.  The ice thicknesses in Crust 2.0 are 
claimed to be within 250 m of the true ice thickness. Sediment thicknesses 
in each cell are to within 1.0 km of the true sediment thickness 
and crustal thickness are within 5 km of the true crustal thicknesses. 




%\hs  The user may input a list of
%bedrock and ice surface altitudes along with their coordinates
%in longitude and latitude, and the gaps are filled in with a 
%smoothing process.  In other words, for the bins where no data
%is given,  the altitudes for the bedrock and ice surface are
%interpolated.  The data that we are using now for bedrock and 
%ice surface elevations comes from: \\
%\centerline{http://classroomantarctica.aad.gov.au/materials/Elevations.pdf }
%which gives rockbed and ice surface elevations only along a single
%line across the continent. 

%The smoothing procedure is as follows.  The longitude and
%latitude coordinates for the inputted data points are translated
%to their appropriate $\phi$ and $\cos{\theta}$ bins and the
%bedrock and ice surface altitudes are set for these bins according
%to the data.  Also, the bins along the edge of the circle that
%defines our coordinate system (latitude=-60$^{\circ}$ and all bins
%in $\phi$) are set to be at sea level (altitudes set to zero
%relative to the earth radius).  All of the bins that have been
%set up to this point are called \emph{reference bins}.

%For a given bin that is not a reference bin, we need to
%set its bedrock and ice surface altitude.  First, we find
%the closest reference bin that is non-trivial.  A trivial reference
%bin is one set to sea level on the outer edge of the circle.
%Next, we find the closest reference bin (this time allowing trivial
%reference bins) that is on the other side of the bin of interest.
%So we want to be able to draw a straight line whose endpoints are
%reference points (at least one non-trivial) and that includes the
%bin of interest somewhere along the line.  Then, the bedrock altitude
%for the bin of interest is set to a value between those of the
%two reference points, determined by the bin's relative distance 
%from each of the reference points.  The same is done for the ice
%surface altitude.

Figures~\ref{fig:upperlayers} and~\ref{fig:lowerlayers} show the 
elevations of the seven crustal layers included in Crust 2.0.
\begin{figure}
\begin{center}
\epsfig{figure=upperlayers.eps,height=3in}
\caption{Altitude of the upper four layers given in Crust 2.0 along the
the $75^{\circ}$ S latitude line.  The horizontal axis is degrees in
longitude.
\label{fig:upperlayers}}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\epsfig{figure=lowerlayers.eps,height=3in}
\caption{Altitude of the lower three layers given in Crust 2.0 along the
the $75^{\circ}$ S latitude line.
\label{fig:lowerlayers}}
\end{center}
\end{figure}

We compute
the total Antarctic ice volume by summing the product of ice
thickness and surface area for each bin within the Antarctic continent.
For the area of each bin, we use:
\begin{equation}
\int_{\phi_1}^{\phi_2}\int_{\theta_1}^{\theta_2}\sin{\theta}~d\theta~d\phi=\left(\phi_1-\phi_2\right)\times\left(\cos{\theta_1}-\cos{\theta_2}\right)
\end{equation}
Where the limits of the integrals define the edges of the bin in
latitude and longitude.
We find $2.7 \times 10^{16}$ m$^2$ of Antarctic ice in this model.
Compare this to the $3.01098 \times 10^{16}$ m$^3$ volume of ice
in Antarctica reported by the US Geological Survey~\cite{usgs};
they are different by 10\%.
Unfortunately, Crust 2.0 calls any ice that sits above water (ice
shelves) just water.  Therefore, the Ross Ice Shelf needs to be
added into the model by hand.  Ice shelves make up 2.4\% of the
ice volume in Antarctica, so it won't account for the entire
$10\%$ difference. 

%\begin{center}
%\begin{tabular}{|l|l|l|l|}
%\hline\hline
%term      &   symbol &  with arguments & units for neutrinos\\
%\hline\hline
%neutrinos &
%  $N_{\nu}$    & 
%  $N_{\nu}$ &
%  (unitless) \\
%\hline
%neutrino intensity & 
%  $I$ &
%  $I(\theta,\phi)$  &
%  [length]$^{-2}$ [sr]$^{-1}$ [time]$^{-1}$\\
%\hline
%neutrino brightness, {\it or} & 
%  $I_{E}$, &
%  $I_{E}(E, \theta, \phi)$ &
%  [length]$^{-2}$ [sr]$^{-1}$ [time]$^{-1}$ [energy]$^{-1}$\\
%neutrino specific intensity  &
%  &
%  &
%  \\
%\hline
%net flux & 
%  $F_{E}$ &
%  $F_{E} (E)$ &
%  [length]$^{-2}$ [time]$^{-1}$ [energy]$^{-1}$\\
%\hline 
%total integrated flux &
%  $F$ &
%  $F$ &
%  [length]$^{-2}$ [time]$^{-1}$\\
%\hline\hline
%\end{tabular}
%\end{center}


\section{The Askaryan Signal}

\subsection{Magnitude of the Pulse}

We use the parameterization outlined in~\cite{jaime} for the peak of the
Askaryan signal:
\begin{equation}
\label{eq:vmmhz}
\frac{\cal E^{(\mathrm{@ 1m})} }{\mathrm{V/m/MHz}}=2.53\times 10^{-7}\cdot \frac{\sin{\theta_v}}{\sin{\theta_c}}\cdot  \frac{E_{em}}{\mathrm{TeV}} \cdot \frac{\nu}{\nu_0} \cdot \frac{1}{1.+\left( \frac{\nu}{\nu_0} \right)^{1.44}}
\end{equation}
where $\nu$ is the frequency, $\nu_0=1.15$ GHz, $E_{em}$ is the shower energy, $\theta_v$ is the viewing angle and $\theta_c$ is the Cerenkov angle. 

\subsection{Electromagnetic and Hadronic Components of the Shower}
 We assume flavor democracy; that by the time the neutrinos get
here, the flavors are fully mixed.  We also keep track of each
flavor individually so the sensitivity to each flavor may be
quoted separately.

We pick an inelasticity $y$ according to~\cite{ghandi},
which we have approximated by a double-exponential.
If we call the electromagnetic component $f_e$ and the 
hadronic component $f_h$, then the dependence of $f_e$ and
$f_h$ on $y$ depends on the whether it is a charged current
interaction that occurred or a neutral current interaction.
We choose
$70.64\%$ of the events at random to be charged-current
and the remainder to be neutral current.
If the incident neutrino was a $\nu_e$ and it is a charged current
event, then $f_e=1-y$ and $f_h=y$.
If the incident neutrino was $\nu_\tau$ or $\nu_\mu$, or
if the $\nu_e$ interacts through a neutral current, then
the electromagnetic component is for now treated as negligible, and
$f_h=y$.  A future version of the code will include $\mu/\tau$
bremsstrahlung and photo-nuclear interactions.


\subsection{Width of Cerenkov Cone}
\label{sec:cerenkov_width}

The user can select one of two different parameterizations of the
width of the Cerenkov cone in inputs.txt.  Here we refer to the
two parameterizations as the ``old'' parameterization and 
``new'' parameterization, for lack of a better terminology.
The user should select 
0 (1) for the old (new) parameterizations in the appropriate line in
the inputs.txt file.  In both parameterizations, the width
of the Cerenkov cone is 
different for the electromagnetic
and hadronic components of the shower.

\subsubsection{Old Parameterization}


For this parameterization, we model the cone widths as 
described in~\cite{jaime97} for electromagnetic showers and
and in~\cite{jaime98} for hadronic showers.
  
The width of the electromagnetic
component in degrees is characterized by:
\begin{equation}
\label{eq:deltheta_em}
\Delta\theta_{em}(\nu)=2.7\cdot \frac{\nu_0}{\nu}\cdot \left( \frac{E_{LPM}}{ 0.14 E_{\nu}+E_{LPM}} \right)^{0.3}
\end{equation}
where $E_{LPM}$ is energy above which the LPM effect becomes important.  
The LPM effect causes
the bremsstrahlung interaction to become suppressed because
the momentum transfer ($\propto k/E^2$) becomes so small 
that the Heisenberg uncertainty
causes the interaction to occur over many
scattering centers, resulting in destructive interference.  This
effect reduces the width of the Cerenkov cone, but not
the magnitude of the electric field at the Cerenkov angle.
For $E_{LPM}$, we use~\cite{pdg}:
\begin{equation}
%E_{LPM}=7.7 \mathrm{\hspace{0.05in}TeV/cm} \times \chi_{\mathrm{ice}}=0.3 \mathrm{~PeV}
  %double elpm=2.E15*(X0DEPTH/X0ICE);  // this is what Jaime uses.  see caption under figure 4 of 0003315.
E_{LPM}=2\times10^{15} \frac{\chi_{\mathrm{medium}}}{\chi_{\mathrm{ice}}}
\end{equation}
where $\chi_{\mathrm{ice}}=39.5$ cm is the radiation length of ice and
$\chi_{\mathrm{medium}}$ is the radiation length of the medium of
interest.


The radio pulse from the hadronic component of the shower is parameterized
as in ~\cite{jaime98}, defining $\epsilon=\log_{10}{E_{\nu}/1\mathrm{TeV}}$:
%    If $\epsilon>=0 && \epsilon<2$:	 

\begin{eqnarray}
&\Delta\theta_{had}=\frac{\nu_0}{\nu}\cdot\left(2.07-0.33\cdot\epsilon+0.075\cdot\epsilon^2 \right)& \epsilon\geq0 \mathrm{~and~} \epsilon<2\\
&\Delta\theta_{had}=\frac{\nu_0}{\nu}\cdot\left(1.744-0.0121\cdot\epsilon\right) & 2\leq\epsilon<5 \\
&\Delta\theta_{had}=\frac{\nu_0}{\nu}\cdot\left(4.23-0.785\cdot\epsilon+0.055\cdot\epsilon^2\right) & 5\leq\epsilon<7\\
\nonumber&\Delta\theta_{had}=\frac{\nu_0}{\nu}\cdot\left(4.23-0.785\cdot7.0+0.055\cdot7.0^2\right)\\
&\cdot\left[1.0+(\epsilon-7.0)\cdot0.075\right] &\epsilon\geq7
\end{eqnarray}
According to~\cite{jaime}, the angular distribution of the 
radio pulse is expected to exhibit a diffraction
pattern centered at the Cerenkov angle.  However, for now 
we approximate the
shape of the electromagnetic or hadronic 
pulse with a double Gaussian distribution with
standard deviations given by $\Delta\theta_{em}$ and 
$\Delta\theta_{had}$.
%  Due to the tails in the diffraction pattern at large angles,
%replacing the Gaussian distribution with a diffraction pattern
%would lead to a higher predicted sensitivity at the highest energies.

In this ``old'' parameterization, the signal strength at viewing
angle $\theta$ 
away from the Cerenkov angle $\theta_C$ is given as in~\cite{jaime97}
and~\cite{jaime98}:
\begin{equation}
E(\theta)=\frac{\sin{\theta}}{\sin{\theta_C}} \cdot
E(\theta_C)\cdot \exp\left[-\ln{2} \left( \frac{\theta-\theta_C}{\Delta\theta_{\mathrm{em,had}}} \right)^2 \right]
\end{equation}



\subsubsection{New Parameterization}

For this parameterization, we model the cone widths guided
by~\cite{jaime05}.

Electromagnetic showers in ice 
are modelled just as described in
the previous section, but here we allow for the possibility
of scaling width for other media, in case one would like to use
this same code for a detector in salt, for example.  We
use Equation~\ref{eq:deltheta_em} for the cone width 
for the electromagnetic component of the shower and
then scale for other media,  inspired by Equation 9 
in~\cite{jaime05},
according to:
\begin{equation}
\Delta \theta_{\mathrm{em, medium}}=\Delta \theta_{\mathrm{em, ice}} \cdot
	\frac{\rho_{\mathrm{medium}}}{\rho_{\mathrm{ice}} } \cdot
	\frac{\mathrm{K}_{\Delta,\mathrm{medium}}}{\mathrm{K}_{\Delta,\mathrm{ice}}} \cdot
	\frac{n_{\mathrm{medium}}^2-1}{n_{\mathrm{ice}}^2-1}
%=RHOMEDIUM/RHOICE
%	/KDELTA_MEDIUM*KDELTA_ICE
%	/X0MEDIUM*X0ICE
%	/sqrt(NMEDIUM*NMEDIUM-1)*sqrt(NICE*NICE-1);
\end{equation}

Hadronic showers are modelled as laid out in Equation 9
in~\cite{jaime05}:
\begin{equation}
\Delta \theta_{\mathrm{had, medium}}=\frac{c}{\nu} \cdot
	\frac{\rho_{\mathrm{medium}}}{\rho_{\mathrm{ice}} } \cdot
	\frac{\mathrm{K}_{\Delta,\mathrm{medium}}}{\mathrm{K}_{\Delta,\mathrm{ice}}} \cdot
	\frac{n_{\mathrm{medium}}^2-1}{n_{\mathrm{ice}}^2-1}
%=RHOMEDIUM/RHOICE
%	/KDELTA_MEDIUM*KDELTA_ICE
%	/X0MEDIUM*X0ICE
%	/sqrt(NMEDIUM*NMEDIUM-1)*sqrt(NICE*NICE-1);
\end{equation}

In this ``new'' parameterization, the signal strength at viewing
angle $\theta$ 
away from the Cerenkov angle $\theta_C$ is given as in~\cite{jaime05}:
\begin{equation}
E(\theta)=\frac{\sin{\theta}}{\sin{\theta_C}} \cdot
E(\theta_C)\cdot \exp\left[-\left( \frac{\theta-\theta_C}{\Delta\theta_{\mathrm{em,had}}} \right)^2 \right]
\end{equation}

\subsection{Additional factor of  root 2/2}






\section{Event Geometry}

\subsection{Picking Interaction Point and Direction}
\label{sec:pickneutrino}
We force the neutrino interaction to occur within the balloon's
horizon, given the balloon position for that event.  
Within that horizon, we stack the ice thicknesses of all 
of the bins in longitude and latitude.  The height of each
bin in the stack is the difference between the ice surface 
and soft sediment elevations for that bin, which is the ice thickness.
We pick a random point along the height of the stack to 
determine which bin the interaction occurs in, and at what
elevation within that bin.  The exact position of the interaction 
in the longitudinal and latitudinal directions 
is chosen at random within the bin.  The position of the
interaction is denoted here as $\vec{x}_{int}$. (wufan.From $\vec(x)_(int)$ we get 
the mirror position of the interaction below the local bottom of the ice where we assume 
the local bottom of the ice is an approximately flat mirror.)

The event is then given a weight $d_i^{pos}$ 
that is the volume of ice within
the horizon divided by the total volume of ice in Antarctica.
The volume of ice for each balloon position (in 100 discrete
bins) is calculated once at the beginning of the code so that
it does not need to be calculated for each event. 

The event is rejected if the
interaction-to-balloon
ray traverses rock.  
Since this is essentially the same as requiring that 
the balloon be within the 
horizon, it would never reject an event if it weren't for
surface tilt and the fact that we generate events to just beyond
the horizon to be safe. 

We pick a neutrino direction at random in $\cos{\theta}$ and $\phi$.  It
is denoted here as $\vec{n}_{\nu}$.  The neutrino momentum (non-normalized)
is denoted $\vec{p}_{\nu}$ and its magnitude simply $p_{\nu}$.
If we choose $\theta$ and $\phi$ among all possibilities,
the interaction-to-balloon ray may be far enough off of the Cerenkov cone
that the signal is undetectable, even if all other conditions are ideal.
Therefore, we only generate neutrino directions such that $|\theta_{\nu,ice}-\theta_{c}|<\theta_{th}$ where $\theta_{\nu,ice}$ is the angle 
between the neutrino 
direction($\vec{n}_{\nu}$) and (wufan the ray out from the interaction point in the ice($\vec{n}_{ray}$). 
For direct rays $\vec{n}_{ray}$ is the direction of interaction-to-balloon ray,but for reflected rays it is the
mirror vector of the 
mirror\_interaction\_point-to-balloon,
which means here $\vec{n}_{ray}$ is downward).
)  
The Cerenkov angle 
$\theta_{c}=\mathrm{Cos}^{-1}(1/n)$, see Section~\ref{sec:indexofrefraction} for how the index of refraction is chosen.  
The angle $\theta_{th}$ is the maximum
angle that the ray may diverge from the center of the Cerenkov cone
for the interaction to still be detectable, given the parameters
of the event.

To find $\theta_{th}$, 
first we find the maximum electric field (in V/m/MHz) at 1~m that is 
possible to be emitted from the interaction within the antenna bandwidth  
using Equation~\ref{eq:vmmhz}.  For
${\cal E^{(\mathrm{@1m})}}={\cal E^{(\mathrm{@1m})}}_{max}$,
we take $\nu=\nu_{max}=1200$ MHz.  Here, we use
$\sin{\theta_v}=1$ since it is the most conservative value 
and it allows us to solve for $\Delta\theta_{th}=\theta_v-\theta_c$.  
Then, we account for the
distance between the interaction and payload $r$ by dividing that
distance into ${\cal E^{(\mathrm{@1m})}}_{max}$ to get the maximum
electric field possibly seen at the payload, $\cal E$$_{max}$.  
From $\cal E$$_{max}$, we derive $\theta_{th}$.

When the shower is dominated by either electromagnetic or hadronic
energy deposition, the angular threshold (the angle beyond which
the signal could not possibly be detectable) is found by solving
for $\Delta\theta_{th}$ in the equation:
\begin{eqnarray}
\label{eq:definingthreshold}
\nonumber V^{max}_{noise} \cdot N_{\sigma}& = f_{em}\cdot {\cal E^{(\mathrm{@1m})}}_{max}(\sin{\theta_v}=1) \cdot h_{eff}^{max} \cdot B & \\ 
& \cdot \mathrm{exp}\left(-\frac{1}{2}\cdot\frac{\Delta\theta_{th}^2}{\Delta\theta\mathrm{(max)}_{em}^2}\right) & f_{em}\gg f_{had}\\
\nonumber &= f_{had}\cdot {\cal E^{(\mathrm{@1m})}}_{max}(\sin{\theta_v}=1) \cdot h_{eff}^{max} \cdot B \\ 
& \cdot \mathrm{exp}\left(-\ln{2}\cdot\frac{\Delta\theta_{th}^2}{\Delta\theta\mathrm{(max)}_{had}^2}\right) & f_{had}\gg f_{em}
\end{eqnarray}
Here, $f_{em,had}$ is the EM or hadronic fraction, if only one of them 
is non-negligible (otherwise, see the next paragraph). 
The maximum effective height $h_{eff}^{max}$ of the
antennas over the antenna bandwidth is 63 cm (from the antenna specs).
The width of one bandwidth slice within the larger antenna bandwidth
is denoted $B$.  The variables 
$\Delta \theta\mathrm{(max)}_{em,had}$ are the maximum widths 
of the Cerenkov cone due to the electromagnetic and hadronic
showers over the bandwidth available to the antennas.  
Section~\ref{sec:cerenkov_width} describes how these widths
are determined.  The largest widths are at the lowest frequency, 200 MHz;
we find $\theta\mathrm{(max)}_{had}~=~3.7^{\circ}$ and 
$\theta\mathrm{(max)}_{em}~=~1.1^{\circ}$ for deep ice, higher within the
firn (see Section~\ref{sec:cerenkov_width}).
$V^{max}_{noise}$ is the maximum expected noise level in a single bandwidth slice
among the 4 layers of antennas, given their cant angles (see Section~\ref{sec:noise}). 
$N_{\sigma}=2.3$ defines the
voltage threshold.
Explicitly, $\theta_{th}$ is (when electromagnetic or hadronic energy deposition 
dominates):
\begin{equation}
\Delta\theta_{th}=\sqrt{-2\cdot\Delta\theta\mathrm{(max)}^2_{em,had}\cdot \ln{\frac{V^{max}_{noise} \cdot N_{\sigma}}{f_{em,had}\cdot {\cal E ^{(\mathrm{@1m})}}_{max}(\sin{\theta_v}=1) \cdot h_{eff}^{max}  \cdot B}}}
\end{equation}

When the electromagnetic and hadronic components of the shower are both
non-negligible, then we find the angular threshold by stepping
by increments of 0.5~$\Delta\mathrm{(max)}\theta_{em}$ (we could
have chosen to step in either
$\Delta\mathrm{(max)}\theta_{em}$
or $\Delta\mathrm{(max)}\theta_{had}$) 
away from the center of the
Cerenkov cone.  When the ray is $\Delta\theta$ from the center
of the Cerenkov cone, the maximum voltage observed at the antenna is
derived as in Equation~\ref{eq:definingthreshold}, 
but adding the electromagnetic and
hadronic components (this time we use $\theta_v=\theta_c+\Delta\theta$ since
we are not inverting the equation in this case):
\begin{eqnarray}
\label{eq:definingthreshold_inquad}
\nonumber V_{rx} = 
	\left[ f_{em}\cdot \mathrm{exp}\left(-\frac{1}{2}\cdot\frac{\Delta\theta^2}{\Delta\theta\mathrm{(max)}_{em}^2}\right) + f_{had} \cdot exp\left(-\ln{2}\cdot\frac{\Delta\theta^2}{\Delta\theta\mathrm{(max)}_{had}^2}\right)\right] \\
\cdot {\cal E^{(\mathrm{@1m})}}_{max}(\theta_v=\theta_c+\Delta\theta) \cdot h_{eff}^{max} \cdot B 
\end{eqnarray}
Here, $V_{rx}$ is the maximum voltage seen at the antenna for a given
$\theta$, and 
all of the variables in Equation~\ref{eq:definingthreshold_inquad} 
are found as described for the previous case above.  The angular
threshold $\theta_{th}$ is the angle at which $V_{rx}$ falls below
$V^{max}_{noise} \times N_{\sigma}$.



\subsection{Finding Earth/Ice Entrance Points for the Neutrino}
\label{sec:entrance}
\hs Finding the entrance point for the neutrino $\vec{r}_{in}$,
given $\vec{p}_{\nu}$ and $\vec{x}_{int}$
is a simple geometry
problem.  Define:
%if (pow(R_EARTH,2)-pow(p,2)*(1-pow(costheta,2))>
\begin{equation}
    a\equiv p\times \cos{\theta}+\sqrt{R_{e}^2-p^2\cdot \sin^2{\theta}}
\end{equation}
where $\theta$ is the angle between $\vec{p}_{\nu}$ and $\vec{x}_{int}$.
Then, the entrance point for the neutrino is:
\begin{equation}
    \vec{r}_{in}=\vec{x}_{int}-a\times\vec{n}_{\nu}.
\end{equation}
The quantity inside the square root is negative for events where
the interaction occurred above sea level and the neutrino is
down-going.  These events are treated separately and the
entrance point is found iteratively.  First, the elevation $h^{(0)}_{int}$
is found at the latitude and longitude where the interaction
occurred, and then $R_e$ in the above equation is replaced with
$R_e+h^{(0)}_{int}$ and a first guess at an entrance point is found.  Then
the elevation at that latitude and longitude, $h^{(1)}_{int}$, is found
and $R_e$ is replaced with $R_e+h^{(1)}_{int}$.  After four iterations,
all but a few percent of the events have converged with a difference
between entrance points by successive iterations being less than
the interaction length for that energy.



\subsection{Balloon and Antenna Positions and Orientations}

\hs The balloon sits at an altitude of 37 km at 80$^\circ$ S latitude.
The position of the balloon in longitude is selected at random for
\emph{each} neutrino generated.

The payload design of interest may be selected in the input file.
For each design, the antennas are given the correct orientation in $\phi$
and inclination in $\theta$
at the right radius from the payload axis, spaced appropriately 
along the vertical axis.  The input file also asks for the 
number of layers to be included.  If the SMEX payload design
has been chosen, and the user sets the number of layers
to be 3, then only the top 3 layers of antennas out 
of the 4 layers in that design are included in the simulation.
This is so that the nadirs can easily be turned on and off.


\section{Ice}
\subsection{Attenuation in the Ice}


The user may set the attenuation length for radio in ice to be 
a constant value (currently 700~m, on the conservative side), 
or a depth-dependent attenuation
length based on (wufan.measurements performed at Ross Ice Shelf and the South Pole.~\cite{smex} )
The latter is the default.  
Fenfang Wu, working with David Goldstein,
 incorporated this feature into the code.

(wufan.
i. For Ice Shelf, we just treat Ronne Iceshelf the same way as Ross Ice Shelf and we find the 
average attenuation length for upward direct rays:
\begin{equation}
\begin{large}
\langle L \rangle = \frac{\int_{-d^{eff}_{int}}^0{1250*0.38*0.08886*\mathrm{exp}(-0.048827*(225.67460-86.517596*\mathrm{log_{10}}(x+848.870)))}dx}{d^{eff}_{int}} 
\end{large}
\end{equation}

For reflected rays which are first downward and then upward, the average attenuation length:

\begin{equation}
\begin{large}
\langle L \rangle = \frac{\int_{-d^{eff}_{int}}^{-d^{eff}_{max}}{1250*0.38*0.08886*\mathrm{exp}(-0.048827*(225.67460-86.517596*\mathrm{log}(x+848.870)))}dx+
\int_{-d^{eff}_{max}}^0{1250*0.38*0.08886*\mathrm{exp}(-0.048827*(225.67460-86.517596*\mathrm{log}(x+848.870)))}dx}{2*d^{eff}_{max}-d^{eff}_{int}} 
\end{large}
\end{equation}

where the interaction performed from the height of the interaction
relative to the surface
$x=-d^{eff}_{int}$ (a negative number) to the surface ($x=0$).
Since the measurements were only carried out to a depth of 420 m at Ross Ice Shelf, an effective maximum depth and effective interaction depth
are utilized given by:
\begin{equation}
d^{eff}_{int} = d_{int}  \hspace{1.0in}(d_{max}<420 \mathrm{m}) 
\end{equation}  
\begin{equation}
d^{eff}_{max}=d_{max}    \hspace{1.0in}(d_{max}<420 \mathrm{m}) 
\end{equation} 

\begin{equation}
d^{eff}_{int} = d_{int}\times\frac{420 \mathrm{m}}{d_{max}} \hspace{0.3in} (d_{max}>420 \mathrm{m})
\end{equation}	

\begin{equation}
  d^{eff}_{max} = 420        \hspace{0.3in} (d_{max}>420 \mathrm{m})(wufan)
\end{equation}	

ii. For other places on Antarctica, 

For direct rays which is always upward, the average attenuation length:
\begin{equation}    
\begin{large}
\langle L \rangle = \frac{\int_{-d^{eff}_{int}}^0{1250*0.08886*\mathrm{exp}(-0.048827*(\mathrm{exp}(\frac{1214.15-x}{1004.40})-55.5889))}dx}{d^{eff}_{int}}
\end{large}
\end{equation}

For reflected rays which are first downward and then upward, the average attenuation length:
\begin{equation}
\begin{large}
\langle L \rangle = \frac{\int_{-d^{eff}_{int}}^{-d^{eff}_{max}}{1250*0.08886*\mathrm{exp}(-0.048827*(\mathrm{exp}(\frac{1214.15-x}{1004.40})-55.5889))}dx+
\int_{-d^{eff}_{max}}^0{1250*0.08886*\mathrm{exp}(-0.048827*(\mathrm{exp}(\frac{1214.15-x}{1004.40})-55.5889))}dx}{2*d^{eff}_{max}-d^{eff}_{int}} 
\end{large}
\end{equation}


Here the measurements were only carried out to a depth of
2810 m, so:
  
\begin{equation}
d^{eff}_{int} = d_{int}  \hspace{1.0in}(d_{max}<2810 \mathrm{m}) 
\end{equation}    
\begin{equation}
d^{eff}_{max} = d_{max}  \hspace{1.0in}(d_{max}<2810 \mathrm{m}) 
\end{equation}  
 
\begin{equation}
d^{eff}_{int} = d_{int}\times\frac{2810 \mathrm{m}}{d_{max}} \hspace{0.3in} (d_{max}>2810 \mathrm{m})
\end{equation}	
\begin{equation}
d^{eff}_{max} = 2810 (d_{max}>2810 \mathrm{m})
\end{equation}


In the above formulas,  $d_{max}$ is the real maximum depth of the ice at the longitude and
latitude position where the interaction occurred.
$d_{int}$ is depth of the interaction (a positive number). wufan)
Then, this attenuation is imposed on the signal:
\begin{equation}
E\rightarrow E \cdot e^{-D/\langle L\rangle}.
\end{equation}
where $E$ is the electric field of the signal and
$D$ is the distance that the RF signal travels 
through the ice.  $D$ is the magnitude of the 
difference between $\vec{x}_{int}$ and
the point where the RF signal exits the ice, the latter found by the
method described in Section~\ref{sec:raytrace}.
(wufan.For reflect rays, we use 10 percent of power reflection rate, so:
\begin{equation}
E\rightarrow 0.316227766E \cdot e^{-D/\langle L\rangle}.
\end{equation}
)
 
\subsection{Index of Refraction}
\label{sec:indexofrefraction}
The index of refraction of deep ice $n_{ice}$ is taken as 1.79, and at the surface the
index of fraction $n_{surf}$ is 1.325.  
At depths shallower than 150 m is the firn, which is 
layers of packed snow.  We use the following equation for
the index of fraction $n_{firn}$ as a function of altitude ($h_{int}$, a negative number, here in m) in the firn~\cite{peter}. 
\begin{equation}
n(h_{int})=1.325+0.463251\times e^{-0.140157\times h_{int}}
\end{equation}
This equation is a fit to the data measured by members of the RICE
Collaboration.  See figure~\ref{fig:nrice}.
\begin{figure}
\begin{center}
\epsfig{figure=nrice.eps,height=4.5in}
\caption{Fit to RICE data for index of refraction as a function of depth below the surface.  Plot is from Peter Gorham.
\label{fig:nrice}}
\end{center}
\end{figure}



\subsection{Ray Tracing}
\label{sec:raytrace}  

\hs   We need to determine the path of the ray that emerges from
the interaction point and arrives at the balloon's position (the
distances between antennas on the payload are small enough to ignore
compared to the distance between any interaction point and the balloon).

Our first guess for the RF exit point is the point on the ice surface that
is directly radially outward from the interaction point.  We call the
vector from the earth's center (the origin of the coordinate system)
to the RF exit point
$\vec{r}_{exit,0}$ and we call the normalized vector from the interaction point
to the RF exit point $\vec{n}_{refr,0}$.
The surface normal at this exit point would also
point radially away from the center of the earth were it not for the
``tilt'' at the surface.  We take this into account by rotating 
the radial vector in the $\phi$ and $\theta$ directions to get the
surface normal $\vec{n}_{surf,0}$.  The magnitude
of the rotation in each direction is the inverse tangent of the slope
along that dimension determined from the difference in the altitudes
of its two neighboring bins in latitude and longitude.

Once $\vec{n}_{surf,0}$ is found, its direction is randomized
to account for surface ``slopeyness.''  To each component of 
$\vec{n}_{surf,0}$ we add a number picked from a Gaussian
distribution of width 0.012.  Then, $\vec{n}_{surf,0}$ is
normalized again.

Then, $\vec{n}_{e-b,0}$ is the vector running from
the first guess RF exit point to the balloon's position.

Using $\vec{n}_{e-b,0}$, we arrive at our second guess for
the direction of the ray in the ice $\vec{n}_{refr,1}$ 
using Snell's law, in its
unfamiliar and not-so-appealing form
\begin{equation}
\label{eq:snell}
\vec{n}_{refr,1}=-\left( \frac{n_{air}}{n_{surface}}\cos{\theta_i}-\sqrt{1-\left( \frac{n_{air}}{n_{surface}}\sin{\theta_i} \right)^2} \right) \vec{n}_{surf,0}+\frac{n_{air}}{n_{surface}}\vec{n}_{e-b,0}
\end{equation}
to get the refracted ray at the surface (the top of the firn).  Then,
we use Snell's law again, replacing $n_{air}$ with $n_{surface}$ and
$n_{surface}$ with $n_{depth}$ to get the refracted ray at the
depth where the interaction occurred.  
This ray is called $\vec{n}_{refr,1}$.  The value 
$n_{depth}$ is found by the method described in Section~\ref{sec:indexofrefraction}.

Using this new direction for the RF in ice, we find a new exit position
$\vec{r}_{exit,1}$ and a new
exit-to-balloon direction, $\vec{n}_{e-b,1}$. (wufan.Here for reflect rays, we use the mirror interaction point to get the new exit position and a new exit-to-balloon directio because the exit rays looks like coming from the mirror interaction point.)  
Instead of finding a new $\vec{n}_{surf,1}$, however, we
still use $\vec{n}_{surf,0}$ due to the slopeyness that was imposed.
Finally, we do one last
iteration and use Equation~\ref{eq:snell} twice again, changing subscripts on
the RHS from $0 \rightarrow 1$, to find $\vec{n}_{refr,2}$, and
from that deduce $\vec{r}^{exit,2}$, $\vec{n}_{surf,2}$, $\vec{n}_{e-b,2}$.

Figures~\ref{fig:converging1} and~\ref{fig:converging2} 
shows the difference between the
exit points from successive iterations, showing that the method
does converge.
\begin{figure}
\begin{center}
\epsfig{figure=diff2.eps,height=3in}
\caption{Difference in meters between RF exit points along the ice-air boundary
from the first and second ray-tracing iterations.
\label{fig:converging1}}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\epsfig{figure=diff3.eps,height=3in}
\caption{Difference in meters between RF exit points along the ice-air boundary
from the second and third ray-tracing iterations.
\label{fig:converging2}}
\end{center}
\end{figure}

 


\subsection{Slopeyness}
We take the surface normal at the longitude and latitude where
the interaction occurs and randomize each component, allowing
it to change by as much as 1.2\%.  This gives 0.01 rad as the
mean magnitude of the slope.  

\subsection{Fresnel Factor}

When the rays from the Cerenkov cone reach the ice-air boundary
they just obey boundary conditions and thus the 
magnitude of the electric field
is altered and the signal traverses the boundary~\cite{dawn}.
This Fresnel factor, as it is called, for the ``pokey'' case is:
\begin{equation}
f_{\parallel,\perp}=\sqrt{\frac{\tan{\theta_I}}{\tan{\theta_T}}(1-r_{\parallel,\perp}^2)}
\end{equation}
where $r_{\parallel}$ is:
\begin{equation}
r_{\parallel}=\frac{\tan(\theta_I-\theta_T)}{\tan(\theta_I+\theta_T)}.
\end{equation}
and $r_{\perp}$ is the Fresnel coefficient for the 
``slappy'' case:
\begin{equation}
r_{\perp} =\frac{\sin(\theta_I-\theta_T)}{\sin(\theta_I+\theta_T)}.
\end{equation}
For the pokey case, the polarization vector is perpendicular
to the plane of incidence, and for the slappy case the polarization
vector is in the plane of incidence.  The general case is
composed of both components.

After dividing up the electric field emerging from the ice 
into its ``pokey'' and
``slappy'' components, $E^{ice}_{\parallel}$ and 
$E^{ice}_{\perp}$, the
electric field that emerges from the ice $E_{air}$ is
\begin{equation}
E_{air}=\sqrt{\left( E^{ice}_{\parallel}\cdot r_{\parallel} \right)^2 + \left( E^{ice}_{\perp} \cdot r_{\perp} \right)^2 }
\end{equation}
We are in the process of verifying that the magnification (separate
from the Fresnel coefficient) is also treated properly in
our simulation.





\section{Antenna Response}

\subsection{Noise}
\label{sec:noise}
An antenna that sees only the sky will see noise due to the temperature of the sky (we take $T_{sky}=15$ K) and the temperature of the front-end electronics ($T_{syst}=200$ K).  An antenna that points downward will see less sky and more ice ($T_{ice}=240$ K).  For a given antenna declination, we find the effective temperature by performing a weighted average of the power seen by the antenna at each declination angle, taking into account the antenna beam widths.  The effective temperature $T_{eff}$ for a cant angle $\theta_{ant}$ is:
\begin{equation}
%T_{eff}=\frac{\int_{-\pi}^{\theta_{horiz}-\theta_{ant}}{(V_{sky}+V_{system})^2 \cdot e^{-ln{2}\cdot ($\theta/\theta_0$)^2}}}{1}
T_{eff}=\frac{\int_{-\pi}^{\theta_{horiz}-\theta_{ant}}{(T_{sky}+T_{syst}) \cdot } e^{-2\cdot \ln{2} \cdot (\theta/\theta_0)^2} d\theta 
+ \int_{\theta_{horiz}-\theta_{ant}}^{+\pi}{(T_{ice}+T_{syst}) \cdot } e^{-2\cdot \ln{2} \cdot (\theta/\theta_0)^2} d\theta
}
{\int_{-\pi}^{+\pi}{e^{-2\cdot \ln{2} \cdot (\theta/\theta_0)^2} d\theta} }
\end{equation}
where $\theta_{horiz}$ is the declination at which the balloon at altitude 37 km sees the horizon 700~km away.  Taking into account the earth's curvature, we find $\theta_{horiz}=\sqrt(2*h_{b}/g_{b})=0.11\mathrm{rad}\approx 6.2^\circ$ where $h_{b}=37.0$ km is the altitude of the balloon and $g_{b}$ is the
radius of the geoid at the balloon's latitude.  At $80^{\circ}$S, 
$g_{b}=6357.7$ km.



Based on this calculation, an antenna canted by $10^{\circ}$ sees an effective temperature of 338.8 K and an antenna canted by $55^{\circ}$ an effective temperature of 428.3 K.


\subsection{Trigger}

We divide up the total bandwidth of the system (200 MHz to 1200 MHz)
into four sub-bands whose central values are:  265, 435, 650 and 
980 MHz and whose widths are:
130, 160, 250 and 370 MHz.  
A trigger channel is one polarization (left- or right-
circularly polarized) of one antenna for in one sub-band in frequency.

The trigger in this simulation is the one described in
Gary Varner's note~\cite{gary}.  The top two layers of eight antennas
on the payload are combined into one layer of 16 antennas for
the purpose of triggering.  So, the trigger sees three layers
of antennas:  two layers of 16 antenns, and the third layer is 
the eight nadir antennas.
At Level 1, an antenna passes the trigger if three out of
8 channels (2 polarizations, 4 bandwidth slices) surpass
the threshold, which is $N_{sigma}=2.3$ times the expected noise
voltage for that channel.
At Level 2, 
trigger channels are considered
in groups of contiguous antennas, and the trigger is fired if
enough antennas in that group pass the Level 1 trigger.
For the top two trigger layers, there are 5 antennas in 
a group and at least 2 must pass.  For the nadir antennas,
the requirement is 2 out of 3.  There is a separate trigger
that considers nadir antennas only, and that requirement
is 3 out of 3 contiguous nadir antennas is a group. 
At Level 3, either two different layers are required to 
have passed at Level 2, or the nadir-only Level 2
trigger must have fired. 
    The noise voltage $V_{noise}$ is 
chosen at random from a Gaussian distribution with
a standard deviation equal to the expected noise for that
channel's bandwidth and RF temperature.

The antenna parameters are inferred from the manufacturer's
specifications for the antennas used in the Anita-lite flight.
The manufacturer is Seavey Engineering Associates, Inc. and
antenna used is the 0.3-1.5 GHz Dual Linear Polarized
Quad Ridge Antenna Model 0312-810.  The specs give the
antenna gain for 5 values of frequency from 300 MHz to 1200 MHz
and the beam width for the E-plane and H-plane, and
for the vertical and horizontal polarization of each.
We interpolate between these points to find the gains
the beam widths at the frequencies between those values given.
For the region below 200 MHz, we use the values given in the
specs for 300 MHz.
Gains are all between around 8 and 12 dBi and are denoted $G$.
The beam widths
vary from around 35$^\circ$ to 65$^\circ$ are are denoted $b^h_E(\nu)$,
$b^v_E(\nu)$, $b^h_H(\nu)$
and $b^v_H(\nu)$, where the subscripts and superscripts denote the
horizontal and vertical polarizations in the E-plane and H-plane.  
At a given
frequency $\nu$, the effective height of the antenna is related
to the gain by:
\begin{equation}
h_{eff}(\nu)= 2\times \sqrt{\frac{G \cdot c^2}{4\pi \cdot \nu^2} \cdot \frac{Z_{rx}}{Z_{air}}}
\end{equation}
where $Z_{air}=377~\Omega$ is the impedance 
of free space and $Z_{rx}=50~\Omega$ is the
impedance of the antenna.

For a signal $F(\nu)$ incident on the $i^{th}$ antenna at an angle
$\theta_H$ off of the H-plane and an angle $\theta_E$ off
of the E-plane, the 0$^{th}$ polarization ($j=0$)
and 1$^{st}$ polarization ($j=1$) signal in the 
$k^{th}$ sub-band from the antenna is:
\begin{equation}
\label{eq:antennavoltage}
V_{ijk}=\sum_{p=p_{min}(k)}^{p_{max}(k)}{\frac{1}{\sqrt{2}}\cdot F(\nu_p) \cdot 0.5 \cdot h_{eff}(\nu_p) \cdot \Delta \nu \cdot
		   c_j\cdot
			    \mathrm{exp}(-\ln{2}\cdot \left[ \left( \frac{\theta_E}{b^j_E(\nu_p)} \right)^2 + \left(\frac{\theta_H}{b^j_H(\nu_p)} \right)^2 \right])}
\end{equation}
Frequency bins run across all four sub-bands between
$\nu_1=200$ MHz and $\nu_2=1200$ MHz,
so that $\Delta \nu=10$ MHz, but the $k^{th}$ sum is over 
frequency bins within the $k^{th}$ sub-band, $p_{min}$ 
to $p_{max}$.    We take $c_0=c_E$ and $c_1=c_H$, where 
$c_E$ ($c_H$) is the component of the incident ray parallel to the E-plane 
(H-plane) and
perpendicular to the normal plane.  If $c_N$ is the component parallel to the
normal plane and perpendicular to both $c_E$ and $c_H$, then the vector
$\left(c_E,c_H,c_N\right)$ is a unit vector.  The magnitudes of the 
left and right
circularly polarized signals are each:
\begin{equation}
V_{i(L,R)k}=\sqrt{V_{i0k}^2+V_{i1k}^2}.
\end{equation}






\section{Volume $\times$ Steradians}
\label{sec:volstr}

\hs The water-equivalent $[V\Delta \Omega]$ is computed using the following formula:

\begin{equation}
[V\Delta \Omega]=V\times \frac{\rho_{ice}}{\rho_{H_2O}} \times   4\pi  \times \frac{\Sigma{w_i}}{N} \times \frac{\sigma_{\rm{SM}}}{\sigma}
\end{equation}

The variables are:

\begin{itemize}
\item{$V$ is the entire volume of Antarctic ice, computed by summing over all bins the product of the ice thickness and the cross-sectional area of the bin.}
\item{$\rho_{ice}$ and $\rho_{H_2O}$ are the densities of ice and water, respectively.}
\item{$N$ is the number of neutrinos that the user asks to be generated.}
\item{$w_i$ is the weight given to each event, described in Section~\ref{sec:weights}.}
\item{$\sigma_{\rm{SM}}$ is the Standard Model $\nu$N cross section}
\item{$\sigma$ is the $\nu$N cross section used for this run}
\end{itemize}
\begin{figure}
\begin{center}
\epsfig{figure=volstr_gzk.eps,height=4.in}
\caption{Effective Volume times steradians from this simulation.  We
also show the predicted GZK flux with a reasonable set of parameters
from~\cite{ess}.
\label{fig:volstr}}
\end{center}
\end{figure}
Figure~\ref{fig:volstr} shows the volume $\times$ steradians
attained from this simulation.  


\subsection{Weighting the Events}
 
In addition to accounting for phase space reduction described
previously,
we weight each event according to the probability that it
interacted in the earth before it reached the ice.


As a neutrino moves through the earth, it encounters varying
densities as it passes through layers of the earth's interior,
and thus differing interaction lengths.

The probability that a neutrino interacts in the rock is:
\begin{equation}
\label{eq:weight}
w=\prod_{i=0}^{n}{e^{-x_i/L_i}}=\prod_{i=0}^{n}{e^{-x_i \rho_i/\ell }}=e^{\frac{1}{\ell}\sum_{i=0}^{n}{x_i \rho_i}}
\end{equation}
where $n$ is the number of layers the neutrino traverses.  A ``layer''
can be the crust, the mantle, or one seven layers defined in Crust 2.0.

Then $x_i$ is the distance the neutrino travels through the 
$i^{th}$ layer in meters, $\rho_i$ is the density of the $i^{th}$ layer,
$L_i$ is the interaction length 
of that layer in meters,
and $\ell$ is the interaction length in kg/m$^2$.

To find the column density of the chord in $kg/m^2$ (the sum 
$\sum_{i=0}^{n}{x_i\cdot \rho_i}$) traversed by the neutrino,
we have two different options in the code, but describe the default
method here.

%\subsection{Method 1}

%\hs Here we use the simple picture of the earth as being divided
%into three regions, the core, the mantle and the crust.
%The outer radii of the core, mantle and crust are denoted $r_1$, $r_2$
%and $r_3$, respectively.  We take $r_1^2=1.2 \times 10^{13}$ m$^2$,
%$r_2=R_{earth}-40$ km, and $r_3=R_{earth}$.
%Similarly, their densities are denoted
%$\rho_1$, $\rho_2$ and $\rho_3$.  In this simple method, we use
%the a spherical model of the earth instead of a geoid.

%We set $\cos{\theta_3}=(\vec{x}_{int}-\vec{r}_{in})\cdot \vec{r}_{in}$.  
%If $\sin{\theta_3}>r_2/r_3$,
%then the neutrino only
%skims the earth's crust before reaching the ice and never reaches
%the earth's mantle.  Then the total chord (in $kg/m^2$) traversed by the
%neutrino is just $|\vec{r}_{in}-\vec{x}_{int}|\cdot \rho_3$.

%However, if $\sin{\theta_3}<r_2/r_3$, then the neutrino reaches
%the mantle.  The distance it travels through the crust before reaching
%the mantle is:
%\begin{equation}
%\label{eq:distance}
%d_3=r_3\cdot \cos{\theta_3}-\sqrt{r_2^2-r_3^2\times \sin^2{\theta_3}}
%\end{equation}
%Some of these neutrinos reach the core and some do not.
%If $\vec{r}_{in}^{mantle}$ is the point at which the neutrino enters
%the mantle, then just as before, we define 
%$\cos{\theta_2}=(\vec{x}_{int}-\vec{r}_{in}^{mantle})\cdot \vec{r}_{in}^{mantle}$ 
%and then if $\sin{\theta_2}>r_1/r_2$, then the neutrino does not
%reach the core and only traverses the mantle before exiting the earth
%through the crust.  Then the chord that the neutrino travels through the
%mantle (in $km/m^2$) is given by:  $d_2=2  r_2 \cos{\theta_2} \cdot \rho_2$.

%If $\sin{\theta_2}<r_1/r_2$, however, then the neutrino enters the core.
%Then, the distance that the neutrino travels in the mantle is given by
%Equation~\ref{eq:distance} with $2 \rightarrow 1$ and $3 \rightarrow 2$.
%For these neutrinos, the distance traveled through the core can
%be found when we define $r_{in}^{core}$ to be the point where the
%neutrino enters the core, $\cos{\theta_1}=(\vec{x}_{int}-\vec{r}_{in}^{core})
%\cdot \vec{r}_{in}^{core}$, and then the distance (in $km/m^2$) is 
%$d_1=2 r_1 \cos{\theta_1} \cdot \rho_1$.
%Then, the appropriate chord lengths are input into Equation~\ref{eq:weight}.

%\subsection{Method 2}

The total chord length in kg/m$^2$ is calculated as follows.  
We first step in 50 m steps
along the neutrino's path as it enters the earth through the 
crust,
summing over $x_i \rho_i$ where $x_i$ is the length of the $i^{th}$ 
step and $\rho_i$ is the density the layer that contains that step.
We continue stepping until we 
are too deep to be in the earth's crust based on Crust 2.0.
At each point, the depth along the neutrino's path is calculated
relative to the height of sea level according to the geoid model.
Next, the distance that the neutrino travels through the
mantle is found through a simple geometrical calculation.
The density of mantle is taken to be 3400 kg/m$^2$.  The
product of the path in the mantle and the mantle density
is added to the sum.
Then, we step through the remaining neutrino path to the
interaction position, finishing the summation.




\subsection{Phase Space Factor}
\label{sec:results.phasespace}


The phase space factor $d_i$ for the $i^{th}$ event is the
product of the phase space factor that comes about when
we select the neutrino interaction position ($d^{pos}_i$)
and the one that comes about when we select the neutrino
direction ($d^{dir}_i$).  Therefore,  $d_i=d^{pos}_i \times d^{dir}_i$.

The position of the neutrino interaction is forced to 
be in the volume of ice that is 
within the horizon of the balloon.  The balloon's longitudinal 
position along
the 80$^\circ$ S latitude line is chosen at random for each event.
The phase space factor $d^{pos}_i$ is then just the
total volume of ice in Antarctica divided by the
volume of ice within the horizon for the $i^{th}$ event.
The volume of ice within the horizon is calculated once at the
start of the program for 100 equally
spaced balloon positions along its circular path.  For each event,
we take the pre-calculated volume for the 
balloon position that is closest to the balloon position for that event.

The neutrino direction is chosen such that the its Cherenkov cone
lies close enough in solid angle 
to the interaction-to-balloon line of sight that
it is possible that the event is observable 
under the most optimistic of
circumstances, as described in Section~\ref{sec:pickneutrino}.
Since the direction for the $i^{th}$ neutrino is chosen from the intersection
of a cone of width $2\sin{\theta_{th}}_i$ and a unit sphere, the corresponding phase space factor is:
\begin{equation}
d^{dir}_i=\frac{4\pi}{ \sin{\theta_{c}} \cdot 2  \sin{\theta_{th}}_i \cdot 2\pi}
\end{equation}

\section{Cross Section Enhancements}

\subsection{Error on the Effective Volume $\times$ Steradian}

Since $[V\Delta\Omega]$ is derived from a sum of weights, the
statistical error is non-trivial.   Imagine we histogram all of the
weights that are summed into $M$ equally spaced bins between 0 and 1.
Then, the total number of events that pass our Monte-Carlo simulation (which
only differs from $[V\Delta\Omega]$ by a constant) is given by:
\begin{equation}
\sum_{i=0}^{M}{f_i \cdot w_i}
\end{equation}
where $f_i$ is the number of events fall in the $i^{th}$ bin centered
on weight $w_i$.  We define $\epsilon$ to be 
the error on $f_i$.  When $f_i<20$, we use the appropriate
Poisson errors from~\cite{cousins}, keeping track of asymmetric errors.   
When $f_i>20$, $\epsilon=\sqrt{f_i}$.
Then, the total error on the sum of weights is:
\begin{equation}
\sqrt{\sum_{i=0}^{M}{(\epsilon_i \cdot w_i)^2}}
\end{equation}



\section{Coding}

\subsection{Random Number Generation}


The {\tt gRandom} variable is a global variable set up by Root that points to a random number generator instance of the type TRandom3.  At the beginning of the program, we 
set gRandom to an instance of TRandom3 with a seed that we set to a constant.  One can change the value of that seed to a different constant through an option in the input file.  You may want to use this option if, for example, you want to simultaneously run multiple jobs with different seeds and be able to combine the results for increased statistics.

For the sake of repeatability, we set the random number seed to a constant instead of, for example, setting it to the current time as is often done.  We would like to be able to get identical outputs for equivalent code with the same inputs.


% N=\int_{E_{1}}^{E_{2}} \alpha(E) I(E) dE \Delta t
% N=\int_{E_{1}}^{E_{2}} \alpha(E) \cdot I(E)\cdot E \cdot \frac{dE}{E} \cdot \Delta t
% N=\int_{E_{1}}^{E_{2}} \alpha(E) \cdot F(E) \cdot d(\log_{10}{E}) \cdot (\log_{10}e) \cdot \Delta t
% N= (\log_{10}e) \Delta t \sum_{j} \alpha_{j} \cdot F_{j} \cdot d(\log_{10}E_{j}) n_{j}= (\log_{10}e) \Delta t \cdot \alpha_{j} \cdot F_{j} \cdot d(\log_{10}E_{j}) 
% F_{j}= \frac{1}{(\log_{10}e) \Delta t \cdot \alpha_{j} \cdot d(\log_{10}E_{j})} 
% F_{j}= \frac{2.3}{\Delta t \cdot \alpha_{j} \cdot d(\log_{10}E_{j})} 





\section{Future Improvements}

The authors are aware that there are many effects that still
need to be included and improvements that need to be made for 
this simulation.  These include:
\begin{itemize}
\item{Surface roughness needs to be treated properly.}
\item{Double-bang taus need to be included.}
\end{itemize}


%\newpage

%\begin{figure}
%\begin{center}
%\epsfig{figure=areasr1.eps,height=6in}
%\caption{Comparison of $[A\Delta \Omega]_{eff} /N_{90}$ for
%different neutrino telescopes for electron-type neutrinos.
%\label{areasr1}}
%\end{center}
%\end{figure}


%\newpage 

%\begin{figure}
%\begin{center}
%\epsfig{figure=areasr2.eps,height=6in}
%\caption{Comparison of $[A\Delta \Omega]_{eff} /N_{90}$ for
%different neutrino telescopes for electron-type neutrinos.
%\label{areasr2}}
%\end{center}
%\end{figure}

%\newpage

%\begin{figure}
%\begin{center}
%\epsfig{figure=areasr3.eps,height=6in}
%\caption{Comparison of $[A\Delta \Omega]_{eff} /N_{90}$ for
%different neutrino telescopes for electron-type neutrinos.
%\label{areasr3}}
%\end{center}
%\end{figure}



%\newpage

%\begin{figure}
%\begin{center}
%\epsfig{figure=volsr1.eps,height=6in}
%\caption{Comparison of $[V\Delta \Omega]_{eff} /N_{90}$ for
%different neutrino telescopes for electron-type neutrinos.
%\label{volsr1}}
%\end{center}
%\end{figure}

%\newpage

%\begin{figure}
%\begin{center}
%\epsfig{figure=volsr2.eps,height=6in}
%\caption{Comparison of $[V\Delta \Omega]_{eff} /N_{90}$ for
%different neutrino telescopes for muon-type neutrinos.
%\label{volsr2}}
%\end{center}
%\end{figure}


%\newpage

%\begin{figure}
%\begin{center}
%\epsfig{figure=volsr3.eps,height=6in}
%\caption{Comparison of $[V\Delta \Omega]_{eff} /N_{90}$ for
%different neutrino telescopes for tau-type neutrinos.
%\label{volsr3}}
%\end{center}
%\end{figure}


\begin{thebibliography}{99}

\bibitem{peter} Peter Gorham, personal communication.

\bibitem{crust2} Bassin, C., Laske, G. and Masters, G., The Current Limits of Resolution for Surface Wave Tomography in North America, EOS Trans AGU, 81, F897, 2000, http://mahi.ucsd.edu/Gabi/rem.html.

\bibitem{smex} Anita SMEX proposal, www.phys.hawaii.edu/$\sim$gorham/ANITA/AnitaCSR.pdf, p. 15.

\bibitem{pdg} Eidelman et al., Physics Letters B592, 1 (2004)

\bibitem{ess} R. Engel, D. Seckel, T. Stanev, ``Neutrinos from Propagation of Ultrahigh Energy Protons,’’  Phys.Rev. D64:093010

\bibitem{jaime} J. Alvarez-Muniz, R.A. Vazquez and E. Zas, ``Calculation
Methods for Radio Pulses from High Energy Showers,'' astro-ph/0003315.

\bibitem{ghandi}  Ghandhi,Reno,Quigg,Sarcevic, hep-ph/9512364.


\bibitem{cousins}Gary J. Feldman, Robert D. Cousins, Phys.Rev.D57:3873-3889,1998, physics/9711021.

\bibitem{geoid}http://www.pol.ac.uk/psmsl/puscience/

\bibitem{gary}http://www.phys.hawaii.edu/$\sim$idlab/publications/ANITA\_GTM\_Note.pdf

\bibitem{usgs} Williams, Richard S. Jr. and Jane G. Ferrigno. Estimated present-day area and volume of glaciers and maximum sea level rise potential. Satellite Image Atlas of Glaciers of the World. US Geological Survey (USGS).  http://pubs.usgs.gov/fs/fs133-99/gl\_vol.html


\bibitem{jaime97}  J. Alvarez-Muniz, Phys.Lett.B411:218-224 (1997) astro-ph/9706064. 

\bibitem{jaime98} J. Alvarez-Muniz, Phys.Lett.B434:396 (1998) astro-ph/9806098.
 
\bibitem{jaime05} J. Alvarez-Muniz, {\em et al.}, 
Phys.\ Rev.\ D {\bf 74}, 023007 (2006) astro-ph/0512337.



\bibitem{dawn} Dawn Williams, Ph.D. Thesis, http://www.physics.ucla.edu/$\sim$dawnwill/

\end{thebibliography} 







\end{document}










