M A N U A L F O R T H E M E R C U R Y I N T E G R A T O R
=================================================================
P A C K A G E V E R S I O N 6
=================================
by John E. Chambers
(with some subroutines supplied by Hal Levison and Martin Duncan)
Dedicated to the memory of Fabio Migliorini
Many thanks to all of you for reporting bugs and suggesting
improvements. Special thanks to David Asher, Scott Manley and
Eugenio Rivera for your help.
< Last modified 1 March 2001 >
N.B. If you publish the results of calculations using MERCURY, please
=== reference the package using J.E.Chambers (1999) ``A Hybrid
Symplectic Integrator that Permits Close Encounters between
Massive Bodies''. Monthly Notices of the Royal Astronomical
Society, vol 304, pp793-799.
C O N T E N T S
===============
(1) Introduction
(2) Initial preparations
(3) How to do an integration
(4) Converting data to orbital elements
(5) Examining data on close encounters
(6) Continuing an integration from dump files
(7) Extending a previous integration
(8) Note for previous users: Changes from Mercury5.
------------------------------------------------------------------------------
(1) I N T R O D U C T I O N
=======================
MERCURY is a general-purpose software package for doing N-body
integrations. It is designed to calculate the orbital evolution of
objects moving in the gravitational field of a large central body.
For example MERCURY can be used to simulate the motion of the planets,
asteroids and comets orbiting the Sun; or a system of moons orbiting
a planet; or a planetary system orbiting another star.
MERCURY is written in Fortran 77. The code is slightly non standard
(e.g. it contains `end do' statements) but it should work using most
compilers.
MERCURY currently includes the following N-body algorithms:
o A second-order mixed-variable symplectic (MVS) algorithm
incorporating simple symplectic correctors (see J.Wisdom et al.
1996, Fields Instit. Commun. vol 10 pp217) - this is very fast
but it cannot compute close encounters between objects.
o A general Bulirsch-Stoer - slow but accurate in most situations.
You can use this when all else fails, or to test whether the other
algorithms are appropriate for the problem you want to study.
o Conservative Bulirsch-Stoer - twice as fast as the general BS
routine, but it will only work for conservative systems, in which
accelerations are a function of position only (e.g. Newtonian
gravity, but not General Relativity).
o Everhart's RA15 (RADAU) - about 2-3 times faster than the general
version of Bulirsch-Stoer. Usually reliable, except for very close
encounters or very eccentric (e.g. Sun grazing) orbits.
o Hybrid symplectic/Bulirsch-Stoer integrator - very fast but only
moderately accurate. This algorithm can compute close encounters.
N.B. The symplectic integrators may give spurious results if some
=== objects have highly eccentric orbits during an integration.
MERCURY includes the effects of Newtonian gravitational forces between
bodies that are assumed to be point masses. It can also calculate
non-gravitational forces for comets (using equations by B.Marsden et al.
1973. Astron. J. vol 78, pp211). You can include the effects of other
forces, using some algorithms, by modifying the subroutine mfo_user.for.
The MERCURY package consists of several drivers and subroutines,
together with a number of input files that you may want to alter to
suit your application.
The Drivers
-----------
1) mercury6_1.for
This is the basic integration programme. It contains all the
subroutines you need to carry out integrations using any of the
algorithms described above. mercury6_1.for produces some output files
that are in a machine-independent compressed format, and you will
need the following programmes to convert this output into a format
you can read.
2) element6.for
This programme converts an output file created by mercury6_1.for into
a set of files containing Keplerian orbital elements for each of the
objects in the integration. These files allow you to see how object's
orbits change with time, and can be used as the basis for making
graphs or movies using a graphics package.
3) close6.for
This programme converts an output file produced by mercury6_1.for
into a set of files containing details of close encounters between
objects during the integration.
Other files
-----------
1) mercury.inc
This contains constants and general parameters used by programmes in
the mercury package. You may want to alter some of the parameters before
you compile and run mercury6_1.for
2) swift.inc
This contains constants and parameters used in the subroutines written
by H.Levison and M.Duncan (1994, Icarus, vol 108, pp18). These
subroutines have names that begin with either `drift' or `orbel'.
N.B. If you change mercury.inc or swift.inc, you must recompile the
=== driver programmes (mercury6_1.for etc) for the changes to take
effect.
------------------------------------------------------------------------------
(2) I N I T I A L P R E P A R A T I O N S
=======================================
Before using the MERCURY package for the first time, you should do the
following:
a) Make sure you have copies of these files:
mercury6_1.for
mercury.inc
swift.inc
element6.for
close6.for
message.in
files.in
You will need some additional input files to run the programmes in the
MERCURY package, but you can create these yourself (see the following
sections for details).
b) Compile mercury6_1.for. The precise command you use will depend on
your Fortran compiler. I suggest you call the executable version
of the programme mercury6
e.g. On Linux systems, try g77 -o mercury6 mercury6_1.for
On DEC Unix systems, try f77 -o mercury6 mercury6_1.for
c) Compile element6.for. I suggest you call the executable element6
d) Compile close6.for. I suggest you call the executable close6
------------------------------------------------------------------------------
(3) H O W T O D O A N I N T E G R A T I O N
===============================================
a) Make sure the compiled version of mercury6_1.for exists.
b) Make sure each of the input files described below exists, and alter
them to suit your needs.
c) Run the compiled version of mercury6_1.for
e.g. On Unix or Linux systems, use the command: ./mercury6
For long integrations, you may want to run the programme in the
background or as a batch job.
d) Read the information summary file produced by mercury6_1.for to make
sure that no problems occurred during the integration. You can read
this file while mercury6_1.for is still running.
The Input Files
---------------
1) files.in
This should contain a list of 10 names for the input and output files
used by mercury6_1.for. List each file name on a separate line.
The first 3 names should be input files that already exist.
In order, these are:
i) A file containing initial data for the Big bodies in the integration
(e.g. planets in the Solar System).
ii) A file containing initial data for the Small bodies in the integration
(e.g. asteroids or comets that you want to include).
iii) A file containing parameters used by the integrator (e.g. start
and end times of the integration etc.).
The MERCURY package includes some sample versions of these files which
you can use as templates. I call these big.in, small.in and param.in
respectively.
The last 7 names in files.in are the names that mercury6_1.for will use
for its output files. In order, these will contain:
iv) Position and velocity information for the objects in the
integration, produced at periodic intervals.
v) Details of close encounters that occur during the integration.
vi) A summary of the integration parameters used in by the integrator,
and a list of any events that took place (e.g. collisions between
objects).
vii) A dump file containing data for the Big bodies. You can use this
to continue the integration if your computer crashes or the
programme is interupted.
viii) A dump file containing data for the Small bodies.
ix) A dump file containing the integration parameters.
x) An additional dump file containing other variables used by
mercury6_1.for
I usually call these 7 files xv.out, ce.out, info.out, big.dmp,
small.dmp, param.dmp and restart.dmp
2) big.in
This file contains the initial data for all the Big bodies in the
integration EXCEPT for the central body (i.e. the Sun if you are
integrating the Solar System). A Big body is defined as one that
perturbs and interacts with all other objects during the integration.
Any lines beginning with the ) character are assumed to be comments,
and will be ignored by mercury6_1.for, however, you should not delete
the first comment line beginning with )O+_06
o The first non-comment line should end with a word that tells the
programme what format you use for the initial data. You should specify
either
Cartesian = for xyz coordinates and velocities. Distances should be
in AU and velocities in AU per day (1 day = 86400 seconds).
Asteroidal = Keplerian orbital elements, in an `asteroidal' format.
i.e. a e I g n M, where
a = semi-major axis (in AU)
e = eccentricity
I = inclination (degrees)
g = argument of pericentre (degrees)
n = longitude of the ascending node (degrees)
M = mean anomaly (degrees)
Cometary = Keplerian orbital elements in a `cometary' format.
i.e. q e I g n T, where
q = pericentre distance (AU)
e,I,g,n = as above
T = epoch of pericentre (days)
o The next line should end with the epoch of osculation in days (i.e.
the time at which the initial coordinates/elements are valid).
E.g. the first few lines of big.in might look like this:
)O+_06 Big-body initial data (WARNING: Do not delete this line!!)
) Lines beginning with `)' are ignored.
)---------------------------------------------------------------------
style (Cartesian, Asteroidal, Cometary) = Asteroid
epoch (in days) = 2451544.5
o The remaining lines provide data for each Big body. The first line
for each body should begin with the body's name, having up to 8
characters.
After that you can include any of the following on the same line:
m = X where X is a real number, to indicate the body's mass in
Solar masses. If you don't specify a value the mass is
assumed to be 0.
r = X where X is a real number, to indicate the maximum
distance from the body (in Hill radii) that constitutes
a close encounter. If you don't include this the default
is r=1
d = X where X is a real number, to indicate the density of the
body in g/cm^3. If you don't include this the default is d=1
a1 = X where X is a real number, to indicate the A1 non-gravitational
force parameter for this body. Realistically this should be
zero for Big bodies (the default is 0).
a2 = X where X is a real number, to indicate the A2 non-gravitational
force parameter for this body (the default is 0).
a3 = X where X is a real number, to indicate the A1 non-gravitational
force parameter for this body (the default is 0).
E.g. the line might look something like this:
MARS m=3.22715144505386530E-07 d= 3.94
The next line(s) for a body should contain the 6 initial coordinates and
velocities or the 6 orbital elements, separated by one or more spaces or
carriage returns. After these numbers you should give the 3 components
of spin angular momentum for the body, in units of solar masses AU^2 per
day (if in doubt enter these as 0).
3) small.in
This file contains the initial data for all the Small bodies in the
integration. A Small body is defined as one that only perturbs and
interacts with Big bodies during the integration. Hence, Small bodies
ignore one another completely (i.e they do not perturb one another,
and they cannot collide with each other).
If you give these objects zero mass they will behave as test particles.
Any lines beginning with the ) character are assumed to be comments,
and will be ignored by mercury6_1.for, however, you should not delete
the first comment line beginning with )O+_06
o The first non-comment line should end with a word that tells the
programme what format you use for the initial data. The possible
formats are the same as those in big.in
o The remaining lines provide data for each Small body. These are
exactly analogous to the lines in big.in, except that you can also
specify an epoch of osculation for each Small body using
ep = X where X is a real number. If you don't include this the
default is X = the same as the epoch for the Big bodies.
Small bodies with differing epochs will be integrated
to the same epoch prior to the main integration.
E.g. the line might look something like this:
HALLEY Ep=2446480.5 a1=0.04d-8 A2 =0.0155d-8
Note that if any of the Small bodies have different epochs than
the Large bodies, the Small bodies must all have zero mass.
4) param.in
This file contains parameters that control how an integration
is carried out. Any lines beginning with the ) character are assumed
to be comments, and will be ignored by mercury6_1.for, however, you
should not delete the first comment line beginning with )O+_06
The file should contain the following items, one per line and in
this order (the programme actually searches for information after
an `=' sign, so you may change the text of the message before this).
o The integration algorithm. Choose one of the following:
mvs : second-order mixed-variable symplectic
bs : Bulirsch-Stoer (general)
bs2 : " " (conservative systems only)
radau : RA15 (RADAU)
hybrid : hybrid symplectic/Bulirsch-Stoer integrator
o The time that the integration should start (in days). This doesn't
have to be the same as the epoch of the Big bodies. Rather, the
integrator will start producing output at this date. If you are
integrating objects in the Solar System, you may want to measure
time in Julian Day numbers (e.g. 1st Jan 2000 = 2451544.5)
o The time at which the integration will finish (in days).
o The output interval (in days). This determines how often the
programme will store orbital elements for the bodies.
o The timestep used by the integrator (in days). For variable
timestep algorithms (e.g. Bulirsch-Stoer and Radau), this is the
stepsize used for the first timestep only. After that the programme
will choose its own timestep. Note that if you choose a large
initial timestep, the variable timestep algorithms may reduce it
in order to maintain the desired accuracy.
o A integration accuracy parameter. This is approximately how much
error per step the variable-timestep algorithms will tolerate.
It is also used by the hybrid algorithm during close approaches.
This number is ignored by the MVS algorithm but you should
provide a number anyway.
The next lines in param.in should contain options that you will only
want to change occasionally. If in doubt, you can use the same
options as the sample param.in file. The options are:
o Should the integrator stop if a close encounter occurs (yes or no)?
o Should the programme check for collisions and take appropriate
action if they occur (yes or no)? If you answer no, all the
bodies will behave as point masses.
o Should collisions result in fragmentation (yes or no)? This
version of MERCURY does not include fragmentation, so this is
ignored at present. You should still specify yes or no however.
o How should the time be expressed (days or years)? This option
controls how information messages produced by mercury6_1.for are
formatted.
o Should time be measured with respect to the integration start
time (yes or no)?
If you choose `years' and `no' for the previous option and this
one, the time will expressed as a Julian date before October
1582, and as a Gregorian date for later dates.
o What level of output precision (low, medium or high)? This
determines how many significant figures will be used to store
the orbital elements (roughly 4, 9 or 15).
o This line is no longer used. However, for backwards compatibility,
you should still include a line here in param.in, although it
doesn't matter what you put on this line.
o Include the effects of general relativity (yes or no)? This
version of MERCURY does not include relativity, so this is
ignored at present.
o Include the effects of the user-defined force routine (yes or
no). You can add additional forces to the integrator in the
subroutine mfo_user in mercury6_1.for
The remaining lines in param.in should contain some other parameters
that you will only need to change rarely. These are:
o The distance from the central body at which objects are
removed (in AU). These bodies are assumed to be so far from
the central body that they are no longer important.
Note that this number is used to scale the output (on a log
scale), so don't make it bigger than you need to.
o The radius of the central body (in AU). Objects coming closer
than this are assumed to collide with the central body.
This number is also used to scale the output (on a log scale).
o The mass of the central body (in solar masses).
o The J2 moment of the central body in units of its radius^2.
o The J4 moment of the central body in units of its radius^4.
o The J6 moment of the central body in units of its radius^6.
o A line which is not used at present. Write whatever you like
on this line.
o Another line which is not used at present.
o The changeover distance used by the hybrid integrator (in Hill
radii). This is the minimum separation between objects before
the hybrid (close encounter) part of the integrator takes effect.
o The maximum number of timesteps between data dumps. This also
controls how often mercury6_1.for notifies you of its progress.
o The number of timesteps between other periodic effects. At
present this controls how often mercury6_1.for checks for
ejections and recomputes objects' Hill radii.
5) message.in
N.B. Alter the contents of this file at your peril!!
===
This file contains the text of various messages output by MERCURY,
together with an index number and the number of characters in the string
(including spaces used for alignment).
------------------------------------------------------------------------------
(4) C O N V E R T I N G D A T A T O O R B I T A L E L E M E N T S
=====================================================================
After doing an integration you can see how the objects' orbits
varied over time. To do so, follow these steps:
a) Make sure the output files produced by the original integration still
exist.
b) Make sure the compiled version of element6.for exists.
c) Make sure each of the input files described below exists, and alter
them to suit your needs.
d) Run the compiled version of element6.for
e.g. On Unix/Linux systems, use the command: ./element6
The programme will produce a set of new files, one per object, containing
orbital elements. Each file has the name of the object with the extension
.aei
The Input Files
---------------
1) element.in
This file contains parameters and options used by element6.for
Any lines beginning with the ) character are assumed to be comments,
and will be ignored by element6.for
o In the first non-comment line, give the number N of compressed
data files you want to read from. Usually this will be 1.
(If you specify more than one file, element6.for will combine
the orbital elements from all the integrations into one set
of output files).
o The next N lines should contain the name(s) of the compressed
data file(s) produced by mercury6_1.for Put each file name on a
different line.
o At the end of the next line, indicate what origin you want to
use for your elements. Choose between Central (for elements
with respect to the central body), Barycentric (for barycentric
elements) or Jacobi (for Jacobi elements).
o On the next line, specify the minimum time interval between the
times at which you would like orbital elements. For example,
if the original integration stored data every 100 days, but you
are only interested in seeing orbital elements every 500 days,
put 500 on this line in elements.in
o On the next line say whether you would like time expressed in days
or years (write years or days).
o Then, on a new line, state whether you want to express the time
with respect to the integration start date (write yes or no).
o Next comes a line indicating which orbital elements you want,
and in what format. Each element is indicated by a code letter,
followed by a number indicating the desired number of digits and
decimal places. If the number of figures is followed by `e', then
the programme will use exponential notation for that element.
The code letters are:
a = semi-major axis (in AU)
b = apocentre distance (in AU, b is short for Big q)
d = density (g per cm^3)
e = eccentricity
f = true anomaly (degrees)
g = argument of perihelion (degrees)
i = inclination (degrees)
l = mean anomaly (degrees)
m = mass (solar masses)
n = longitude of ascending node
o = obliquity (degrees)
p = longitude of perihelion (degrees)
q = pericentre distance (AU)
r = radial distance (AU)
s = spin period (days)
x, y or z = Cartesian coordinates x, y or z
u, v or w = Cartesian velocities vx, vy or vz
E.g. a8.4 e8.6 i7.3 g7.3 n7.3 l7.3 m13e
indicates that you want the semi-major axis (8 digits including
4 decimal places), eccentricity (8 digits including 6 decimal
places) etc.... and mass (13 digits in exponential notation).
Note that if you choose to express an element using a large number
of significant figures, the last few digits might not be meaningful
if the output precision of the original integation was low or medium.
o The remaining lines in elements.in should contain the names of the
objects for which you want orbital elements. If you don't supply
any names, the programme assumes that you want elements for all the
objects.
2) message.in
This is the same file as used by mercury6_1.for
------------------------------------------------------------------------------
(5) E X A M I N I N G D A T A O N C L O S E E N C O U N T E R S
===================================================================
To examine details of close encounters that occurred during an
integration, follow these steps:
a) Make sure the output files produced by the original integration still
exist.
b) Make sure the compiled version of close6.for exists.
c) Make sure each of the input files described below exists, and alter
them to suit your needs.
d) Run the compiled version of close6.for
e.g. On Unix/Linux systems, use the command: ./close6
The programme will produce a set of new files, one per object, containing
details of close encounters with that object. Each file has the name of
the object with the extension .clo
The Input Files
---------------
1) close.in
This file contains parameters and options used by close6.for
Any lines beginning with the ) character are assumed to be comments,
and will be ignored by close6.for. Don't delete the first comment
line though.
o In the first non-comment line, give the number N of compressed
data files you want to read from. Usually this will be 1.
(If you specify more than one file, close6.for will combine
close-encounter details from all the integrations into one set
of output files).
o The next N lines should contain the name(s) of the compressed
data file(s) produced by mercury6_1.for Put each file name on a
different line.
o On the next line say whether you would like time expressed in days
or years (write years or days in close.in).
o Then, on a new line, state whether you want to express the time
with respect to the integration start date (write yes or no).
o The remaining lines in m_close.in should contain the names of the
objects for which you want close-encounter details. If you don't
supply any names, the programme assumes that you want details for
all the objects.
2) message.in
The same file as used by mercury6_1.for
N.B. When using the hybrid symplectic algorithm, only those close
=== encounters that are integrated using the Bulirsch-Stoer part of
the integrator will be saved. In practice, this means that some
distant ``close'' encounters will not be recorded.
------------------------------------------------------------------------------
(6) C O N T I N U I N G A N I N T E G R A T I O N F R O M
===========================================================
D U M P F I L E S
===================
If your computer crashes while MERCURY is doing an integration, all is
not lost. You can continue the integration from the point at which
mercury6_1.for last saved data in dump files, rather than having to
redo the whole calculation. Just follow these steps:
a) Make sure all of the input, output and dump files used by the original
integration are still present.
b) Make sure the filenames listed in files.in correspond to these files.
c) Run the compiled version of mercury6_1.for
If for some reason one of the dump files has become corrupted, look
to see if you still have a set of files with the extension .tmp
produced during the original integration (if you have subsequently
used mercury6_1.for to do another integration in the same directory,
you will have lost these unfortunately). These .tmp files are duplicate
copies of the dump files. Copy each one so that they form a set of
uncorrupted dump files (e.g. copy big.tmp to big.dmp etc.), and then
run the compiled version of mercury6_1.for
N.B. It is important that you replace all the dump files with the .tmp
=== files in this way, rather than just the file that is corrupted.
------------------------------------------------------------------------------
(7) E X T E N D I N G A P R E V I O U S I N T E G R A T I O N
===============================================================
If you want to extend an old integration that finished successfully
(i.e. not one that crashed), follow these steps:
a) Make sure all of the input, output and dump files used by the original
integration are still present.
b) Make sure the filenames listed in files.in correspond to these files.
c) Change the finish time in the parameter dump file (see section 3)
to the end point of the extended integration.
d) Run the compiled version of mercury6_1.for
------------------------------------------------------------------------------
(8) C H A N G E S F R O M M E R C U R Y 5
===========================================
The code has been restructured, primarily to make it easier to
incorporate new coordinate systems (binary-star coordinates, and
some others that I'm tinkering with, that are not included in this
version). All the coordinate change routines now accept the same
input and output (except going to/from Keplerian elements).
A major result of this is that the arrays for objects are now
indexed starting at 2 instead of 1 (the central body uses index 1).
The compressed output has changed again (!! sorry about that).
However, the new version is an improvement for a few reasons.
Firstly, it makes it possible to choose your coordinate origin
(central/barycentric/Jacobi) when you run the decompression
programme rather than when you do the original integration.
Secondly, it handles hyperbolic elements better - no more problems
with the hyperbolic mean anomaly. In addition, you no longer have
to specify emax or qmin in the param.in file - the new compressed
format can handle all osculating values of q and e. Finally, it is
now possible to produce .aei files when the number of objects in the
integration exceeds the number of files that can be open at the
same time for your operating system.
Rather than remodifying m_elem.for to cope with yet another type
of compressed output, I have written a new programme called elements
(which is also easier to pronounce). Use this for compressed output
generated by Mercury6, and use m_elem5 for the earlier versions of
Mercury.
A bug fix: Mercury5_2, and earlier, had a bug which occurred when
integrating test particles using the MVS algorithm. This incorrectly
modelled the perturbations due to the innermost massive body.
The central-body oblateness terms J2, J4 and J6 now function. These
were also included in a few older versions of Mercury, but were not
properly tested until now.
The MVS algorithm now incorporates simple symplectic correctors
which should improve the accuracy using this algorithm. The
correctors remove error terms proportional to eh^3 thru eh^6, where
h is the timestep and e is the ratio of the object masses to the
central mass. Hence, the leading error terms are proportional to
e^2h^3 and eh^7.
Close encounter details are now output in batches, rather than after
every encounter. Furthermore, the programme does not always do regular
data dumps in addition to data dumps at every output interval. These
changes are made in the hope that they will reduce delays caused by
accessing the hard disk.