This is the README file for the code Geokerr, which calculates all coordinates of null geodesics in 
a Kerr spacetime, as presented in Dexter & Agol (2009). Please cite this paper if you make use of 
this code in your research.

The package contains the FORTRAN source file geokerr.f, as well as example files input_example.f and 
abgrid.in. To get started, compile geokerr.f, for example with gfortran:

gfortran -ffixed-line-length-132 -O3 geokerr.f -o geokerr

The -ffixed-line-length-132 option relaxes the FORTRAN line length limit, while -O[1-3] are 
optimization options. The code has also been tested with ifort and g77.

QUICK START AND DESCRIPTION OF PROGRAM TESTGEOPHIT

For a quick start, the included program TESTGEOPHIT can be run using one of the example input files
(.in), for instance abgrid.in. By default the code is set to read from standard input and write to
standard output. So to write the geodesic data for this example to a file 'abgrid.out', use

./geokerr < abgrid.in > abgrid.out

To run the other example, first compile inputex.f. Then run

./geokerr < example.in > example.out

To reproduce Fig. 1 of Dexter & Agol (2009), plot dt vs. uf for each geodesic.

  OUTPUT FORMAT

The first line in the output file contains a summary of the rest of the file, including the total 
number of geodesics, mu0, a, and u0. Then each geodesic is written with the first line specifying 
alpha, beta (or q2, l if provided instead), the number of points on the geodesic, and the case 
number corresponding to Tables 1,2 of Dexter & Agol (2009). Finally, the geodesic coordinates are 
output with the columns uf, muf, dt, dphi, lambda,  mu turning points, u turning points. A template 
is given here:

NGEO MU0 A U0
ALPHA BETA NUP NCASE
UFI(1) MUFI(1) DTI(1)        DPHI(1) LAMBDAI(1)            TPMI(1) TPRI(1)
UFI(2) MUFI(2) DTI(2)-DTI(1) DPHI(2) LAMBDAI(2)-LAMBDAI(1) TPMI(2) TPRI(2)
UFI(3) MUFI(3) DTI(3)-DTI(1) DPHI(3) LAMBDAI(3)-LAMBDAI(1) TPMI(3) TPRI(3)
...
ALPHA2 BETA2 NUP2 NCASE2
...

With slight modification of the program TESTGEOPHIT, the output variables or format can be changed. 
By changing the variables INUNIT (OUTUNIT), code input (output) can use the file INFILE (OUTFILE) 
instead of being redirected via standard input (output) as is the default.

  INPUT FILES

The input format is given in abgrid.in, and can be easily read from the source of PROGRAM 
TESTGEOPHIT. It's as follows:

STANDARD
MU0
A
RCUT
NROTYPE
A1 A2 B1 B2
NRO NPHI NUP

where RCUT is the outer boundary of \sqrt{ALPHA^2+BETA^2} for a circular grid (NROTYPE=1), and 
[A1,A2], [B1,B2] are the boundaries in alpha and beta respectively for a rectangular grid 
(NROTYPE=2). RCUT is not used when NROTYPE > 1, but must be present anyway. NRO is the number of 
geodesics along the alpha axis for a rectangular grid, or along the \SQRT{ALPHA^2+BETA^2} axis for a 
circular grid. NPHI is the number of geodesics along the beta axis for a rectangular grid, or along the
PHI=ATAN(ALPHA/BETA) axis for a circular grid. Finally, NUP is the number of points along each 
geodesic.

The program can also accept non-standard inputs with STANDARD=0. The format for such a file can be 
seen from looking at the source of input_example.in. A template is as follows:

STANDARD
QL
NGEO
A
NUP
ALPHA BETA U0 MU0 UF SU TPR
ALPHA2 BETA2 U02 MU02 UF2 SU2 TPR2
...

In this case, STANDARD=0. QL is a variable stating whether (=1) or not (=0) the inputs will be in 
terms of Q2 and L instead of ALPHA and BETA. If QL=1, Q2 replaces ALPHA and L replaces BETA in the 
list above. All other variables are the same as above.

SUBROUTINE GEOKERR

The main program works well for writing geodesic information to a file. However, it may be desirable 
to calculate null geodesics and store them in memory for use with other routines. This can be 
accomplished by calling the main driver subroutine geokerr, which inputs constants of the motion q2, 
l, a as well as starting positions u0, mu0 and either uf, tpr (r turning points) or muf, tpm (mu 
turning points) and outputs arrays of coordinates along the desired geodesic.

There are some subtleties with calling geokerr:

1) USEGEOR=.FALSE. (u as independent variable)
UPDATE 5/18/2009: A new boolean variable, MUFILL, has been added to GEOKERR which can provide extra resolution 
near u turning points, where small steps in u lead to large steps in mu and phi. It works by switching the 
independent variable to MU a user defined KEXT steps before the turning point, computing NPTS points evenly spaced 
in MU between KEXT steps before the turning point to KEXT steps after it. Both KEXT and NPTS can be changed in the
code itself, according to the number of extra points required and their desired locations. Setting MUFILL=.FALSE., 
as is done in the main program by default, reverts to the behavior of the last version.
--The variable UOUT sets a radial boundary on the portion of geodesic to trace. For the whole thing, 
set UOUT=U0.
--For uf as the independent variable, the value of TPR=0 means not to trace past a turning point if 
one exists. Likewise, TPR=1 doesn't mean a turning point must exist.
--Although it may seem redundant, U0, UF, TPR and SU must all be specified separately for using uf. 
This is to avoid ambiguity -- if su=1 with uf > u0 could mean tpr=0 or su=-1 with tpr=1.
--To trace a full geodesic, for instance from infinity to the horizon or to the turning point and 
back, set u0 small (for instance according to the formula used to TESTGEOPHIT), uf=uplus and tpr=1.

2) USEGEOR=.TRUE. (mu as independent variable)
UPDATE 5/18/2009: The caveat below is no longer present. GEOKERR can handle an arbitrary number of mu turning points in the region to be traced, as long as the total number is known in advance and supplied to the routine as TPM.
--The main caveat with this case is that at the moment, geokerr is not set up to handle tracing segments
  of geodesics which cross turning points. This involves a procedure similar to what's done in the 
USEGEOR=.FALSE. section and isn't difficult to implement, but it's not done yet. It's not an issue for 
any of the applications shown in Dexter & Agol (2009).
--The difficulty with this case in general is that the number of mu turning points must be input ahead 
of time, which requires prior knowledge of the geodesic. In the case of a thin disk or even one with angular spread
as shown in Dexter & Agol (2009), the number of turning points can be calculated given MU0 and BETA.

GEOMU, GEOR, GEOPHITIME

These routines and supporting functions contain the formulas from Dexter & Agol (2009). Unless the 
subroutine geokerr is unsatisfactory for the desired application, these don't need to be called 
directly.

OTHER

The supporting routines consist of Carlson elliptic integral routines based on his papers 
(ellquarticreal, ellcubicreal, ellquarticcomplex, elldoublecomplex, ellcubiccomplex) as well as 
routines to compute mu component integrals (ellphitmu, calcphitmuasym, calcphitmusym, etc.). All of 
these routines use Carlson elliptic integral functions (rj, rf, rc, rd).

WHEN YOU FIND BUGS OR HAVE QUESTIONS

Please e-mail jdexter@u.washington.edu