Title\thanksref{label1}
Title\thanksref{label1}
This paper describes the Fortran 77 code {SIMU}, version 1.1, designed for numerical simulations of observational relations along the past null geodesic in the Lema\^{\i}tre-Tolman-Bondi (LTB) spacetime. {SIMU} aims at finding scale invariant solutions of the average density, but due to its full modularity it can be easily adapted to any application which requires LTB's null geodesic solutions. In version 1.1 the numerical output can be read by the {GNUPLOT} plotting package to produce a fully graphical output, although other plotting routines can be easily adapted. Details of the code's subroutines are discussed, and an example of its output is shown.
ms2.tex
% Template article for preprint document class `elsart'
% SP 2001/01/05

\documentclass{elsart}

% Use the option doublespacing or reviewcopy to obtain double line spacing
%\documentclass[doublespacing]{elsart}

% if you use PostScript figures in your article
% use the graphics package for simple commands
% \usepackage{graphics}
% or use the graphicx package for more complicated commands
% \usepackage{graphicx}
% or use the epsfig package if you prefer to use the old commands
% \usepackage{epsfig}

% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}

\topmargin 1cm

\begin{document}

\begin{frontmatter}

% Title, authors and addresses

% use the thanksref command within \title, \author or \address for footnotes;
% use the corauthref command within \author for corresponding author footnotes;
% use the ead command for the email address,
% and the form \ead[url] for the home page:
% \title{Title\thanksref{label1}}
% \thanks[label1]{}
% \author{Name\corauthref{cor1}\thanksref{label2}}
% \ead{email address}
% \ead[url]{home page}
% \thanks[label2]{}
% \corauth[cor1]{}
% \address{Address\thanksref{label3}}
% \thanks[label3]{}

\title{A Fortran Code for Null Geodesic Solutions in the
       Lema\^{\i}tre-Tolman-Bondi Spacetime}

% use optional labels to link authors explicitly to addresses:
% \author[label1,label2]{}
% \address[label1]{}
% \address[label2]{}

\author{Marcelo B.\ Ribeiro\thanksref{if}\thanksref{capes}}

\address{Vatican Observatory Group, Steward Observatory, University of
         Arizona, Tucson, AZ 85721, USA; e-mail: mribeiro@as.arizona.edu} 

\thanks[if]{On leave from Physics Institute, University of Brazil -- UFRJ,
            Rio de Janeiro.}
\thanks[capes]{Supported by Brazil's CAPES Foundation.}

\begin{abstract}
This paper describes the Fortran 77 code {\tt SIMU}, version 1.1, designed
for numerical simulations of observational relations along the past null
geodesic in the Lema\^{\i}tre-Tolman-Bondi (LTB) spacetime. {\tt SIMU}
aims at finding scale invariant solutions of the average density, but due
to its full modularity it can be easily adapted to any application which
requires LTB's null geodesic solutions. In version 1.1 the numerical
output can be read by the {\tt GNUPLOT} plotting package to produce a
fully graphical output, although other plotting routines can be easily
adapted. Details of the code's subroutines are discussed, and an example
of its output is shown.
\end{abstract}

\begin{keyword}
% keywords here, in the form: keyword \sep keyword
cosmology \sep numerical simulation; solution of equations \sep numerical
relativity
% PACS codes here, in the form: \PACS code \sep code
\PACS 98.80.-k \sep 02.60.Cb \sep 04.25.Dm
\end{keyword}
\end{frontmatter}

% main text
\section{Introduction}\label{int}

It is a basic result of relativity theory that light rays follow null
geodesics in curved 4-dimensional spacetimes, and this means that if one
is pursuing cosmological applications of solutions of Einstein's field
equations with the aim of comparing the model's theoretical predictions
with the actual data produced by astronomical observations, one requires
null geodesic solutions for the chosen metric. In fact, for cosmological
applications one only needs the {\it past} null geodesics, as what we
observe today are events which occurred in the past. Nevertheless, a
very basic problem arising in cosmology, and which hinders comparison
between theory and observations in other models than the standard
Friedmann-Lema\^{\i}tre-Robertson-Walker (FLRW), stems from the fact
that as soon as we depart from the simple FLRW models largely used by
observational cosmologists, the task of finding null geodesic solutions
quickly becomes an intractable analytical problem. In fact, this
situation is not fully appreciated by many of those dealing with
cosmological modelling of astronomical data, inasmuch as the simple
FLRW models largely applied in observational cosmology is the exception,
since it has a very simple analytical solution for the null geodesic
equation.

The Lema\^{\i}tre-Tolman-Bondi (LTB) spacetime is the most general
spherically symmetric dust solution of Einstein's field equations,
having the FLRW models as special sub-cases \cite{rib92b,rib94,kr}. By being
a spatially inhomogeneous model, that is, where $\rho=\rho(r,t)$, it allows
us to study scenarios which do not have FLRW's spatial homogeneity assumed
{\it a priori}. Due to this, it is the most widely used cosmological model
after the standard FLRW, having been applied in a wide range of cosmological
problems, from microwave background radiation studies to hierarchical
(fractal) modelling, just to quote a few (see \cite{kr}, and references
therein, for a large number of applications of LTB models in cosmology).
However, as is the case with almost all non-standard cosmologies, while
it is possible to solve analytically its Einstein's field equations,
LTB's null geodesic equation remains an intractable analytical problem,
unless we rewrite its geometry in terms of the so-called ``observational
coordinates''. That, however, has the handicap of adding a great deal of
mathematical complexity to the problem, since such an approach requires
complex mathematical calculations in curved spacetimes \cite{mhms},
something which may not be required in all LTB applications. Therefore,
due to such a wide range of application, a code which produces numerical
solutions of LTB null geodesics is desirable.

In this paper I will describe such a code. It deals with LTB geometry in
its full generality, being applicable to each of its special sub-cases,
parabolic, hyperbolic and elliptic \cite{bonnor74}, either separately or
together in a single problem encompassing all sub-cases, if desired. It
was originally designed to find scale invariant solutions of the average
density by means of numerical simulations \cite{rib92,rib93}, but as we
shall see below, without changing its core null geodesic calculations the
code can be easily adapted to do much more than this, due to its
modularity. The original results obtained by this code were recently
analytically confirmed, and extended, by other authors \cite{matra} by
means of the observational coordinates approach mentioned above. The
confirmation of the results of \cite{rib93} by \cite{matra} adds then
further reliability to the numerical results obtained with the code that
will be described below. Besides, those two approaches, numerical and
analytical, when applied to a difficult problem such as null geodesic
solutions in LTB spacetime, will in fact complement each other, rather
than exclude one another. In \S \ref{des} I will describe version 1.1 of
the code, and in \S \ref{io} the input and output will be discussed. The
results are summarized in \S \ref{conc}, where I also indicate where the
code can be obtained.

\section{Description of {\tt SIMU} 1.1}\label{des}

{\tt SIMU} solves simultaneously two ordinary differential equations
(ODE's) by means of the fourth order Runge-Kutta method with adaptive
stepsize control. One differential equation is required for the null
geodesic, while the other is needed for obtaining the redshift in this
cosmology \cite{rib92}. The subroutines that carry out the numerical
integration are from \cite{nr}, and some minor changes necessary for
taking the results back to the main program were made, but without
changing the actual numerical implementation of the Runge-Kutta method.
A full description of the LTB model and notation as used in {\tt SIMU}
is given in \cite{rib92}. The numerical simulations take advantage of
the fact that LTB spacetime geometry has three unknown
functions,\footnote{ \ Actually two, as a third  can be eliminated by
a coordinate transformation.}$F(r)$, $f(r)$, $\beta(r)$, respectively
representing the amount of gravitational mass within the radius coordinate
$r$, the overall curvature and dynamics of the model, and the time elapsed
since the big bang for each observer located at particular values of $r$
\cite{rib94,kr,bonnor74}. The simulations then take advantage of this
freedom, as one starts by choosing these functions, run the program
and analyse the results, concluding then whether or not the simulation
was successful, and if not choose another set of three functions (see
\cite{rib93} for a very detailed explanation of this procedure).

The code is implemented in double precision, and it has a built-in
methodology for checking the possible catastrophic accumulation of
round-off errors. That is done by monitoring the energy equation, as
derived from Einstein's field equations, and its derivative, since, by
theory, both must remain unchanged throughout the integration
\cite{rib93}. In addition, to check for accuracy and stability, after
the integration the program runs in reverse, that is, it takes the
final result and uses it as initial conditions to get another result
that can be compared with the original initial condition. The results
described in \cite{rib93} showed that double precision is enough for
maintaining the desired accuracy required by the problem under study.
{\tt SIMU} also calculates the errors associated with each observational
quantity evaluated, as well as the propagated errors. This is done by
taking the accuracy given by the adaptive method and using this value
as the input in the standard propagation error equations. All those
details are fully explained in the initial documentation of the code.

As mentioned above, the code uses some subroutines from \cite{nr},
namely the fourth order Runge-Kutta integration with adaptive stepsize
control, the root finding algorithm for the transcendental equations
appearing in LTB's elliptical and hyperbolic sub-cases \cite{rib92},
the extended trapezoidal rule quadrature algorithm, as well as Simpson's
rule to the desired accuracy. However, apart from those, the remaining
subroutines are new, as well as the actual way in which the whole set of
subroutines were linked to each other and organized. In addition, the
code is extensively commented and widely documented, also having a
structure chart showing the dependence among the subroutines and
functions. This detailed documentation is added in order to help
readers who might be interested in implementing {\tt SIMU} in their
own computer environments, or changing it to suit their specific
applications.

It must be mentioned that there are more ``state of the art'' ways
for finding numerical solutions of ODE's than using the adaptive
4th order Runge-Kutta method. The ODEPACK subroutine package is an
example \cite{r2,r3,r4}. Nevertheless, although the chosen ODE integrator
for {\tt SIMU} may have already been superseded by better methods, it
{\it does} solve the proposed problem, and at the desired accuracy.
Inasmuch as {\tt SIMU} is a tested code, with its results having already
been analytically confirmed by other authors \cite{matra}, there is no
need to change to another ODE integrator at {\tt SIMU}'s current version.
The reader must, however, be aware that such a change might be required
for other applications of {\tt SIMU} to LTB spacetime.

The code uses only standard Fortran 77 commands to avoid possible
compilation problems with machine dependent commands that are not
available universally. If a potential user is only interested in the
core of the program, that is, the part which integrates the null
geodesic and calculates the observational relations, he or she can
simply remove its ``tail'', where the observational relations are
manipulated according to the aims of the theory studied in
\cite{rib92,rib93}, and replace it with something else. This ``tail''
consists of the two last branches of the structure chart headed by
the subroutines {\tt FIT} and {\tt OUTPUT}, and are easily spotted
in the {\tt MAIN}. Finally, the functions $f(r)$, $F(r)$, $\beta(r)$
and their derivatives appear in six subroutines at the end of the
list so that they can be compiled separately.

The previous version 1.0 (formerly {\tt SIMU 5d}; see \cite{thesis}) had
the {\tt PGPLOT} plotting routines merged into the second half of the
{\tt OUTPUT} subroutine, while version 1.1 removed them and piped the
numerical results into eleven independent files, {\tt v1} to {\tt v11},
that are then {\it externally} read by the {\tt GNUPLOT} plotting
package in order to produce a graphical output of the results. This is
the only difference between versions 1.0 and 1.1, meaning that in
version 1.1 the graphical output is completely independent from the
code itself. Making this change was justified on two grounds: {\it (i)}
the {\tt GNUPLOT} package has greater plotting capabilities, even
allowing \LaTeX \ output formats that can be included directly in the
{\tt figure} environment of any \LaTeX \ macro. So, producing plots
with the results is done by feeding the numerical tables into {\tt
GNUPLOT} itself by means of a script file written for, say, a linux
environment; {\it (ii)} this feature of having the results outputted
into numerical files brings flexibility to the code, as its results
can be independently used by another program, if so desired or demanded
by another application of the LTB spacetime.

\section{Input and Output}\label{io}

The actual implementation for the original problem advanced in
\cite{rib92} was described at length in \cite{rib93}, where one can find
detailed plots with the results and analysis of the various
simulations,\footnote{ \ These results were summarized in \cite{rib94}.}
as well as detailed explanations of initial conditions and functions used
in each simulation. So, in here I shall limit the discussion to the code
itself, without dealing with any specific application. 

As mentioned above, we require a script file for running a
simulation and immediately producing a graphical output. An example of
such a script in linux is given as follows:

{\footnotesize \tt g77 simu1.1.for\\a.out < in.simu\\gnuplot < in.gnu\\latex
plots.tex\\xdvi plots.dvi}

Here the file {\tt in.simu} contains the maximum value for the radius
coordinate $r$ for the integration, while {\tt in.gnu} is a short
{\tt GNUPLOT} script file for inputting the {\tt v1} to {\tt v11} files
and producing \LaTeX \ files (see appendix \ref{app1}). Finally
{\tt plots.tex} is a file for running the {\tt .tex} \LaTeX \ files produced
by {\tt GNUPLOT} in order for the results be available for analysis
directly on the screen.

The details of each simulations, as well as accuracy monitoring are given
in an outputted file called {\tt s1}. Appendix \ref{app2} shows the
results of a simulation using the value 1.5 for {\tt in.simu}, and the
three functions given as follows: $f(r)=\cosh(r)$, $F(r)=\sinh^3(r)$,
$\beta(r)=0.7$. 

\section{Conclusion}\label{conc}

In this paper I have described the Fortran 77 code {\tt SIMU} version
1.1 for calculating solutions of the null geodesics equations in the
Lema\^{\i}tre-Tolman-Bondi spacetime geometry. I have discussed the
details of the code, input and output, as well as the ways in which it
can be possibly modified for other cosmological applications of this
geometry. The code is available for download at \cite{distro}.

\ack

I am grateful to William R.\ Stoeger for reading the original manuscript,
Newton Jos\'e de Moura Jr.\ for his assistance with the {\tt GNUPLOT}
package, and the referee for helpful comments. I also thank the Vatican
Observatory Research Group for their kind hospitality while writing this
paper.

\begin{thebibliography}{00}

% \bibitem{label}
% Text of bibliographic item

% notes:
% \bibitem{label} \note

% subbibitems:
% \begin{subbibitems}{label}
% \bibitem{label1}
% \bibitem{label2}
% If there is a note, it should come last:
% \bibitem{label3} \note
% \end{subbibitems}

\bibitem{rib92b} M.\ B.\ Ribeiro, {\it Astrophys.\ J.\ }{\bf 395} (1992) 29
\bibitem{rib94} M.\ B.\ Ribeiro, {\it in:} D.\ Hobill, A.\ Burd and A.\
        Coley, eds., {\it Deterministic Chaos in General Relativity}
        (New York, Plenum, 1994) 269-296
\bibitem{kr} A.\ Krasi\'nski, {\it Inhomogeneous Cosmological Models}
        (Cambridge University Press, 1997)
\bibitem{mhms} R.\ Marteens, N.\ P.\ Humphreys, D.\ R.\ Matravers and
        W.\ R.\ Stoeger, {\it Class.\ Quantum Grav.\ }{\bf 13} (1996)
	253; corrigenda {\bf 13} (1996) 1689; M.\ E.\ Ara\'ujo
	and W.\ R.\ Stoeger, {\it Phys.\ Rev.\ D} {\bf 60} (1999) 104020
\bibitem{bonnor74} W.\ B.\ Bonnor, {\it Monthly Not.\ Royal Astron.\
        Soc.\ }{\bf 167} (1974) 55 
\bibitem{rib92} M.\ B.\ Ribeiro, {\it Astrophys.\ J.\ }{\bf 388} (1992) 1 
\bibitem{rib93} M.\ B.\ Ribeiro, {\it Astrophys.\ J.\ }{\bf 415} (1993) 469
\bibitem{matra} N.\ P.\ Humphreys, D.\ R.\ Matravers and R.\ Marteens,
        {\it Class.\ Quantum Grav.\ }{\bf 15} (1998) 3041, gr-qc/9804025
\bibitem{nr} W.\ H.\ Press, B.\ P.\ Flannery, S.\ A.\ Teukolsky and W.\
        T.\ Vetterling, {\it Numerical Recipes}, (Cambridge University
	Press, 1986)
\bibitem{r2} A.\ C.\ Hindmarsh, {\it `ODEPACK, a Systematized Collection
        of ODE Solvers'}, {\it in:} R.\ S.\ Stepleman et al., eds., {\it
	Scientific Computing} (North-Holland, 1983) 55-64
\bibitem{r3} L.\ R.\ Petzold, {\it SIAM J.\ Sci.\ Sta.\ Comput.\ }{\bf 4}
        (1983) 136
\bibitem{r4} L.\ R.\ Petzold and A.\ C.\ Hindmarsh, ODEPACK subroutine
        package available from netlib at {\tt http://www.netlib.org/odepack/}
\bibitem{thesis}  M.\ B.\ Ribeiro, {\it On Modelling a Relativistic
        Hierarchical (Fractal) Cosmology by Tolman's Spacetime}, PhD
	thesis, (Queen Mary and Westfield College, University of London,
	1992); {\tt http://www.if.ufrj.br/\~{}mbr/publi.html}
\bibitem{distro} {\tt SIMU 1.1} can be downloaded from
        {\tt http://www.if.ufrj.br/\~{}mbr/codes} or
	{\tt http://arXiv.org/abs/gr-qc/0205}

\end{thebibliography}

\appendix
\section{{\tt in.gnu} Script}\label{app1}

{\footnotesize \tt
set terminal latex\\
set output 'v1.tex'\\
set size 1,1\\  
set xlabel "\$ $\backslash$log d\_l  (Gpc)\$"\\
set ylabel "\$ $\backslash$log  $\backslash$rho\$"\\
set nokey\\
plot 'v1' with points 1 8, 'v11' with lines 1\\ 
set output 'v2.tex'\\
set xlabel "\$d\_l  (Gpc)\$"\\
set ylabel "\$ $\backslash$rho\$"\\
plot 'v2' with points 1 8\\
set output 'v3.tex'\\
set xlabel "\$z\$"\\
set ylabel "\$d\_l\$ $\backslash$$\backslash$ \$(Gpc)\$"\\
plot 'v3' with points 1 8\\
set output 'v4.tex'\\
set xlabel "\$z\$"\\
set ylabel "\$N\_c\$"\\
plot 'v4' with points 1 8\\
set output 'v5.tex'\\
set xlabel "\$r\$"\\
set ylabel "\$Field\$ $\backslash$$\backslash$ \$EE\$"\\
set yrange [-10**-5:10**-5]\\
plot 'v5' with points 1 8\\
set output 'v6.tex'\\
set xlabel "\$r\$"\\
set ylabel "\$Deriv.\$ $\backslash$$\backslash$ \$Field\$
$\backslash$$\backslash$ \$EEDSH\$"\\
set yrange [-10**-5:10**-5]\\
plot 'v6' with points 1 8\\
set output 'v7.tex'\\
set title "\$Center  $\backslash$rightarrow Border\$"\\
set xlabel "\$r\$"\\
set ylabel "\$t\$"\\
set autoscale y\\
plot 'v7' with points 1 8\\
set output 'v8.tex'\\
set title "\$Border  $\backslash$rightarrow Center\$"\\
set xlabel "\$r\$"\\
set ylabel "\$t\$"\\
plot 'v8' with points 1 8\\
set output 'v9.tex'\\
set title "\$Center  $\backslash$rightarrow Border\$"\\
set xlabel "\$r\$"\\
set ylabel "\$I\$"\\
plot 'v9' with points 1 8\\
set output 'v10.tex'\\
set title "\$Border  $\backslash$rightarrow Center\$"\\
set xlabel "\$r\$"\\
set ylabel "\$I\$"\\
plot 'v10' with points 1 8\\
quit
}

\newpage
\section{Plots of a Single Simulation}\label{app2}

\begin{figure}[h]
    \input{v1}
\end{figure}
\begin{figure}[b]
    \input{v2}
\end{figure}
\begin{figure}[htb]
    \input{v3}
\end{figure}
\begin{figure}[htb]
    \input{v4}
\end{figure}
\begin{figure}[htb]
    \input{v5}
\end{figure}
\begin{figure}[htb]
    \input{v6}
\end{figure}
\begin{figure}[htb]
    \input{v7}
\end{figure}
\begin{figure}[htb]
    \input{v8}
\end{figure}
\begin{figure}[htb]
    \input{v9}
\end{figure}
\begin{figure}[htb]
    \input{v10}
\end{figure}

\end{document}


gr-qc0205095.extracted_bib
\begin{thebibliography}{00}

% \bibitem{label}
% Text of bibliographic item

% notes:
% \bibitem{label} \note

% subbibitems:
% \begin{subbibitems}{label}
% \bibitem{label1}
% \bibitem{label2}
% If there is a note, it should come last:
% \bibitem{label3} \note
% \end{subbibitems}

\bibitem{rib92b} M.\ B.\ Ribeiro, {\it Astrophys.\ J.\ }{\bf 395} (1992) 29
\bibitem{rib94} M.\ B.\ Ribeiro, {\it in:} D.\ Hobill, A.\ Burd and A.\
        Coley, eds., {\it Deterministic Chaos in General Relativity}
        (New York, Plenum, 1994) 269-296
\bibitem{kr} A.\ Krasi\'nski, {\it Inhomogeneous Cosmological Models}
        (Cambridge University Press, 1997)
\bibitem{mhms} R.\ Marteens, N.\ P.\ Humphreys, D.\ R.\ Matravers and
        W.\ R.\ Stoeger, {\it Class.\ Quantum Grav.\ }{\bf 13} (1996)
	253; corrigenda {\bf 13} (1996) 1689; M.\ E.\ Ara\'ujo
	and W.\ R.\ Stoeger, {\it Phys.\ Rev.\ D} {\bf 60} (1999) 104020
\bibitem{bonnor74} W.\ B.\ Bonnor, {\it Monthly Not.\ Royal Astron.\
        Soc.\ }{\bf 167} (1974) 55 
\bibitem{rib92} M.\ B.\ Ribeiro, {\it Astrophys.\ J.\ }{\bf 388} (1992) 1 
\bibitem{rib93} M.\ B.\ Ribeiro, {\it Astrophys.\ J.\ }{\bf 415} (1993) 469
\bibitem{matra} N.\ P.\ Humphreys, D.\ R.\ Matravers and R.\ Marteens,
        {\it Class.\ Quantum Grav.\ }{\bf 15} (1998) 3041, gr-qc/9804025
\bibitem{nr} W.\ H.\ Press, B.\ P.\ Flannery, S.\ A.\ Teukolsky and W.\
        T.\ Vetterling, {\it Numerical Recipes}, (Cambridge University
	Press, 1986)
\bibitem{r2} A.\ C.\ Hindmarsh, {\it `ODEPACK, a Systematized Collection
        of ODE Solvers'}, {\it in:} R.\ S.\ Stepleman et al., eds., {\it
	Scientific Computing} (North-Holland, 1983) 55-64
\bibitem{r3} L.\ R.\ Petzold, {\it SIAM J.\ Sci.\ Sta.\ Comput.\ }{\bf 4}
        (1983) 136
\bibitem{r4} L.\ R.\ Petzold and A.\ C.\ Hindmarsh, ODEPACK subroutine
        package available from netlib at {\tt http://www.netlib.org/odepack/}
\bibitem{thesis}  M.\ B.\ Ribeiro, {\it On Modelling a Relativistic
        Hierarchical (Fractal) Cosmology by Tolman's Spacetime}, PhD
	thesis, (Queen Mary and Westfield College, University of London,
	1992); {\tt http://www.if.ufrj.br/\~{}mbr/publi.html}
\bibitem{distro} {\tt SIMU 1.1} can be downloaded from
        {\tt http://www.if.ufrj.br/\~{}mbr/codes} or
	{\tt http://arXiv.org/abs/gr-qc/0205}

\end{thebibliography}