swift/README.first

******************************************************************************
                PRELIMINARY USER'S MANUAL FOR SWIFT
******************************************************************************
Authors : Martin Duncan and Hal Levison
Date:     June 15/93
Last Revisions: 11/29/00 HFL


                                INTRODUCTION

     The SWIFT subroutine package provided here is designed to
integrate a set of mutually gravitationally interacting bodies
together with a group of test particles which feel the gravitational
influence of the massive bodies but do not affect each other or the
massive bodies.  Four integration techniques are included:
            1) Wisdom-Holmam Mapping (WHM).  This is the n-body mapping 
               method of Wisdom & Holman~(1991; AJ, 102, 1528).
            2) Regularized Mixed Variable Symplectic (RMVS) method. This
               is an extension of the WHM that handles close approachs 
               between test particles and planets. See Levison & Duncan 
               (1994; Icarus, 108, 18).
            3) A fourth order T+U Symplectic (TU4) method.  This is the 
               T+U method of Candy & Rozmus (1991; J. Comput. Phy., 92, 230).
               Also see Gladman, Duncan, & Candy(1991; CeMDA, 52, 221.)
            4) A Bulirsch-Stoer method.  This is a Bulirsch-Stoer from Press,
               Teukolsky, Vetterling, & Flannery (1992. `Numerical Recipes in 
               FORTRAN'). 
            5) SyMBA.  This is a symplectic N-body that handles close
               encounters between massive bodies from Duncan et al
               (1998, aj, 116, 2067). SyMBA is run a little
               differently from the rest of Swift.  See README.symba
               for more details. 
 The package is designed so that the calls to each of these look
identical so that it is trivial to replace one with another.

     In order to get to this file you have presumably obtained,
uncompressed and untar-ed the swift.tar.Z file found on the Web at
www.boulder.swri.edu or in the anonymous ftp directory on
gort.boulder.swri.edu.  You will now find various files and
subdirectories.  Most of the directories contain the subroutines that
make up SWIFT.  However, fully functional examples of working drivers
are found in the subdirectory "main". They all have the same i/o
except for minor differences noted below. The subroutines used are
grouped by function in subdirectories with what we hope to be good
internal documentation.  There is a Tex file describing the WHM force
calculations in the subdir. "mvs/getacch". The subdir. called
"example" contains a working example for integrating the Sun plus the
four giant planets as well as Pluto (the latter treated as a test
particle) using the integrator "swift_whm". Descriptions of the
drivers and the input and output files follow.


                             COMPILING SWIFT   

    There are two files that must be edited in order to complile
SWIFT. In the top SWIFT directory, you will find a file called
"@make".  Edit this file.  You will have to change the values of:
     SWIFT_DIR which is the top directory in which you put SWIFT.
     FORTRAN   which is the command on your machine to run the Fortran 
                 compliler
     FFLAGS    which are the flags you use for your Fortran. You MUST
                 include a flag that will tell Fortran to make an object
                 file, but does not try to link it.  This is a "-c" on
                 most machines.  (In addition for SyMBA users, the
                 flags must be set so that recursion will work on your
                 machine.)  
      PRECOMP  is the full path of the C pre-compiler on your machine.
                 This program is usually called cpp.
      CPPFLAGS are the flags for the cpp pre-compiler.  See the next section
                 for a discussion.

   Then type "@makeall".  This will compile all the subroutines in the
package.  Swift will create a library called "libswift.a".  This file
must be copied into the "/usr/local/lib" directory on your machine.
You must be root in order to do this.  Once this is completed you can
link a program to this package by adding the flag "-L/usr/local/lib -lswift" 
to your Fortran statements.

   Now you must complile the drivers.  Change to the directory main
(enter: " cd main").  You will have to edit "@make" and change the
same variables as you did above.  In this case, you do not want the -c
option of the Fortran compiler.  Now enter "@make".  The executable
files in the directory are now ready for use.

   You can now remove swift.tar, if you want.


                               THE CPP FLAGS
 
   We have made a significant effort making sure that SWIFT is ANSI
standard Fortran 77.  Therefore, in principal, the Fortran source code
should compile on any UNIX platform.  Unfortunately, there are two
places where this ideal breaks down: in the OPEN function and in
recursive routines.  In these cases, we make use of the C
pre-compiler, cpp, which is universally available on UNIX machines (at
least we hope).  Cpp translates files ending with .F into Fortran .f
files.  In SWIFT, cpp needs two flags that will tell it what kind of
fortran files to create.  These are:

1) OPEN function: In several places, SWIFT needs to open files and
   append to the end of them.  So, far we have found the there are two
   ways in which Fortran compilers do that.  The first is with OPEN
   specifier "position='append'" and the second with "access='append'".
   If your machine is of the first type then use the CPPFLAG of
   "-D_OPEN_POSITION".  If it is of the second use "-U_OPEN_POSITION".
   We have found that most machines are of the second type, however IBM
   RS6000's and HP's Fortran 90 are of the first type.  

2) RECURSIVE ROUTINES: (Only SyMBA users need to worry about this!)
   Fortran 77 and Fortran 90 compilers handle recursion differently. As
   far as our experience goes, F77 allows any subroutine to be recursive,
   although you most likely will have to set compiler flags to accomplish
   this.  F90 compilers require recursive subroutines to declared in the
   source code as `recursive subroutine X(....)'.  If your compiler
   requires this than use the CPPFLAG of "-D_RECUR_SUB".  If not use
   "-U_RECUR_SUB"

3) FXDR package: XDR is a platform independent binary format that
   Swift can use to write bin.dat files that can move between
   systems. If requires FXDR which can be found at:
       http://meteora.ucsd.edu/~pierce/fxdr_home_page.html
   However, Swift is written so that it can be run without it.  If you
   have the fxdr library then set the -D_FXDR_AVAIL.  If you do not
   want to use it then set -U_FXDR_AVAIL.  If you try to write an XDR
   without the FXDR package Swift will complain and stop.

If you don't know which is correct for you, we have supplied a test. 
Change to the directory cpp_test (enter: " cd cpp_test").  Follow the
instructions in the README file.

                                THE DRIVERS

        All of the drivers are designed to take essentially the same
input files and to produce the same type of output.  The basic step
for all of them begins with heliocentric positions and velocities and
advances them a timestep dt.  Some methods (such as the Bulirsch-Stoer
integrator) have internal substeps or iterations etc. but this is all
hidden from the user.  A description of the basic i/o files is given
in the next section.  The differences between drivers is given next.

SWIFT_WHM : The basic WHM integrator when particles are to be
            removed at the time of close planetary encounters.
            Arbitrarily close solar encounters can occur.
            (This was called swift_mvs in previously releases of Swift.) 

SWIFT_RMVS3: The Regularized Mixed Variable Symplectic (RMVS) integrator. 
            This handles close approachs between test particles and planets.

SWIFT_BS  :  A Bulirsch-Stoer integrator with the same I/O as
            WHM except that in addition to the 3 input files
            the user is prompted for a tolerance EPS (1.e-8 is often
            used). The tolerance could be added to the end of
            the file mvs.in if redirection is used at input (see I/O below).
            The BS routine CAN get through close encounters, although
            it is usually about 10 times slower then SWIFT_RMVS.  
            One can imagine writing a hybrid driver to switch from the
            WHM stuff to BS during occasional close encounters.

SWIFT_LYAP: Computes Lyap. exponents for a set of test particles
           by integrating a set of shadow particles (one per test
           particle) with initial offsets in pos. and vel. read in
           from a file called 'shadow.in' with same format as
           'tp.in' (see I/O description below).  At time intervals
           spaced by dtout, the shadows are renormalized to their
           initial phase space distances from the real particles.
           An output file called 'lyap.out'  has for each
           output time the log (base 10) of the Lyap. exponent
           gamma for each particle.  NOTE:  We don't recommend
           that you use this without talking to one of us first!

SWIFT_LYAP2: Computes Lyap. exponents for a set of test particles
           by integrating the difference equations.  Besides the usual
           par.in, pl.in, and tp.in, the program asked you how often
           to output the phase space distances.  It also asks for a
           phase space vector representing the initial separation
           between a particle and a fictitious particle.  Note that
           the orbit of the fictitious particle is not integrated, 
           the difference equations are integrated.  An output file called
           'lyap2.out'  has for each output time the phase space distance
           for each particle.  NOTE: WE ARE STILL TESTING THIS, SO BE CAREFULL.

SWIFT_TU4 : A fourth order symplectic integrator based not on the WHM
          decomposition of the Hamiltonian but rather on the traditional
          Kinetic (T) plus Potential (U) form. Used mostly as a check
          on the other routines.


SWIFT_SYMBA5: SyMBA.

                                INPUT/OUTPUT

        The simplest way to use e.g. 'swift_whm' is to use the file called
        'mvs.in' in the 'examples' directory, which has just 3 lines
        corresponding to the names of 3 input files : param.in pl.in  tp.in.
        Thus one types

                        swift_whm < mvs.in

        The input parameters are as follows:

PARAM.IN : This file has all the run-time parameters of the form :

        t0, tstop, dt
        dtout, dtdump
        L1 L2 L3 L4 L5 L6
        rmin, rmax, rmaxu, qmin, lclose
        binary_outputfile
        status_flag_for_open_statements

   where
        t0 is the initial time
        tstop is the time to stop the integration
        dt is the timestep

        dtout is the time interval between outputs (see below)
        dtdump is time interval between dumps of the 3 basic files (see below)

        L1-L6 are logicals and are set to T or F
          L1:  .true.==> Include J2 and J4 for central body.  These values
                         must be included in the pl.in file.
          L2:  .true.==> Does the various checks for coming too close to Sun
                         or planet, or too far from Sun.  If particles are 
                         removed a file called 'discard.out' will contain a 
                         set of entries for each removed particle that gives 
                         the time, the particle's posn and vel and istat flags
                         and rstat flags (see a file called README.stat for a 
                         discription of the status flags), as well as all the 
                         planetary posns and vels. at the time of removal.
                         So you will have a complete discription of the 
                         planetary system at removal time if you need it.
                 NOTE: IF THIS BIT IS NOT SET YOU SHOULD OMIT THE LINE
                 IN PARAM.IN THAT CONTAINS RMIN, RMAX ETC. AND HAVE ONLY
                 THREE LINES IN ALL :  THE LAST HAVING THE NAME OF
                 THE BINARY OUTPUT FILE.
          L3:  .true.==>  Computes Jacobi integral for the tps (only useful
                          if have only one planet on circ. orbit).  Writes to 
                          a file called jacobi.out
          L4:  .true.==> Computes energy and ang. momentum and writes out 
                         to a file called energy.out. (of course this only 
                         checks the accuracy of the planetary integrations).
          L5:  .true.==> Write to a binary output file the heliocentric
                         orbital elements of all bodies and test particles.
                         The elements are written as REAL*4. The elements are 
                         a,e,i,OMEGA,omega and M, where OMEGA is the long. of 
                         ascending node, omega is arg. of peri. and M is mean 
                         anomaly. ROUTINES TO READ FROM THE BINARY FILE ARE 
                         GIVEN IN FOR EXAMPLE THE ROUTINE "FOLLOW" IN THE 
                         "TOOLS" DIRECTORY. TO USE FOLLOW SIMPLY RESPOND TO 
                         THE PROMPTS (e.g. the last query is for which 
                         particle # to follow - use a positive integer for a 
                         tp or a negative number if you want to follow a 
                         planet (e.g. -3 for the third massive body ie the 
                         second planet).
          L6:  .true.==> Same as L5 except that Swift write an XDR file.
                 NOTE: L5 and L6 cannot both be set to .true.

        rmin is the heliocentric distance at which a tp is stopped because
             it is deemed to be too close to the central body.
             SET TO -1. if you want to ignore this check
        rmax is the heliocentric distance at which a tp is stopped because
             it is deemed to be too far from the central body
             SET TO -1. if you want to ignore this check
        rmaxu is the heliocentric distance at which a tp is stopped because
             it is deemed to be too far from the central body AND it
             is also unbound with respect to the central body
             SET TO -1. if you want to ignore this check
        qmin is the perihelion distance at which a tp is stopped because
             it is deemed to be too close to the central body.
             SET TO -1. if you want to ignore this check
	     Note:  The orbital elements computed are the osculating 
                    (instantaneous) HELIOCENTRIC values. Because of the 
                    motion of the Sun induced by the planets, these values
                    are unrepresentative of the barycentric values for 
                    particles which have semi-major axes beyond a few 
                    hundred AU.  In particular, values for the pericentric 
                    distance and the inclination will vary considerably. As 
                    a result we STRONGLY RECOMMEND that a value of RMAX of 
                    a few hundred AU be used if a nonzero value of QMIN is 
                    used.
        lclose is a logical.  If .true. then the code will check for close 
             approaches between tp and planets.  If tp gets within a distance 
             rpl of the planet then it is stopped.  rpl for each planet must 
             be included in the pl.in if this flag is set.  NOTE: we recommend 
             that rpl be greater than a Hill sphere if you are using WHM.
        binary_outputfile is the name of a binary file to store the
            orbital elements at time intervals dtout (see below).
        status_flag_for_open_statements is the status flag for the open
        statements of all the output files but lyap.dat.  Must be one
        of the following:
                         new      (program dies if the file exists)
                         append   (data is appended to the file)
                         unknown  (data is written over existing data).

NOTES : 1) Every dtout time increments, the code outputs various quantities,
         depending on the value of the logical flags L1-L6 (read in from
         param.in). 

        2) At time increments set by dtdump, the code dumps all of
        the information needed to resume the integration at that time
        in case of power failures or in case one wishes to resume an
        integration from its endpoint. The info. is in 3 files called
        dump_pl.dat,  dump_tp.dat, and dump_param.dat.  The format is
        identical to pl.in, tp.in and param.in respectively.  Note
        that t0 in dump_param.dat records the time of the dump and
        that the status_flag_for_open_statements is changed to
        `append', so that these files can be used to restart a stopped
        integration.  Depending on the situation, one may wish to 
        increase tstop in order to extend an integration.   The files 
        are overwritten with each dump so only the most recent one is 
        preserved.

        NOTE: We strongly suggest that you set dtdump=dtout.  This will 
              make it easier to restart the program if your computer
              goes down.

PL.IN : The code requires units in which the grav. constant G is unity.
        Any combo of lengths, masses and times that keeps that true is O.K.
        For example, the pl.in file given in the examples dir. uses lengths
        in AU and time in days (where one year is exactly 365.25 days). This
        forces the Solar mass to be about 2.96e-4.  In the example, we
        give initial conditions from the classic paper of CHO for the
        integration of the outer planets.  One could instead use units
        in which lengths are in AU, and the Solar Mass is unity but
        then the orbital period of a test particle at one AU would
        be 2*pi. The format is simple : first the # of bodies on the first
        line  (INCLUDING the Sun) and then 3 lines for each body giving
        mass on the first line, heliocentric x,y,z on the next and
        heliocentric vx,vy,vz on the third.  NOTE:  The x,y,z and
        vx,vy,vz for the Sun MUST be 0!!

        In addition:  if L1 = .true. then next to the Sun must be the
        values of J_2 * R^2 and J_4 R^4, where R is the radius of the
        central body.  

        If L2=.true. and lclose=.true. than the lines that contain the masses 
        of the planets must also include the stopping radius.

TP.IN : In the same units as PL.IN, the first line is # of test particles.
        The test particles are assumed massless so for each particle there
        are 6 lines giving heliocentric x,y,z on first, helioc. vx,vy,vz on
        second and then 13 istat integers and 13 rstat integers (See 
        README.stat for the definitions of the status flags).  For most 
        runs one would begin with ststud flags equal to zero
        unless it is a resumption of an earlier run. The code may
        alter the arrays as it runs.

********NOTE TO USERS WHO MAY MODIFY THE CODE : ***************************
        For reasons of efficiency and simplicity, in some of the subroutines
        which use the array istat, istat is dimensioned and used as if it were
        a vector rather than a 2D array. These subroutines need only know
        whether or not istat(i,1) is zero and we don't pass in or out the
        other elements of the matrix. This is standard practice in Fortran
        and we note clearly in the description in each subroutine if this
        is the case.
**************************************************************************