Table of Contents

In addition to integrating the equations
of motion, the particles are allowed to loose some (**eta**) of their random
energy (smoothed over a certain ‘‘box’’ size (**cell**)), and turn it into mean
orbital motion.

Orbits can also be diffused: each timestep the velocity
vector can be rotated over the randomly gaussian distributed angle (**sigma**).

A special integration mode (mode=-1) will advance particles assuming their orbits are epicycles, and launched with velocities deviating from the angular velocity consistent with rotation. The supplied potential should be axisymmetric and numerical derivatives are used to calculate the angular, radial and vertical frequencies such that the orbital motions can be solved for analytically. from their guidance center.

**in=***in-file*- Initial conditions
will be read from
*in-file*in snapshot format [default: none - but is required]. **out=***out-file*- If given, results are written to
*out-file*in snapshot format [default: empty, no snapshot output file produced]. **potname=***name*- name of
file of
*potential(5NEMO)*descriptor [default:**harmonic**]. **potpars=***par-list*- List of parameters to the potential descriptor. The first parameter must
be the pattern speed in the x-y plane, although rotating frames of reference
are not yet supported. The remaining parameters are used by the
*_inipotential()*routine in the potential descriptor. [default: none - let them be defined by routine itself]. **potfile=***file-name*- name of an optional datafile to the potential descriptor. This might be an N-body snapshot or list of spline fit coefficients etc. [default: none].
**save=***state-file*- If given, the system
state will be saved in
*state-file*after each timestep. Useful for some recovery after system crashes. **freq=***integ-freq*- Inverse time-step, to be used with the
time integrator. [Default:
**64.0**(64 steps per unit time)]. **mode=***mode-number*- Integration modes, must be a number:

0=Euler (kick and drift at same time) 1=PC (Predictor Corrector) algorithm, 2=modified PC algorithm, 3=Runge-Kutta (JEB modified) 4=4th order Runge-Kutta. 5=leapfrog 6=modified Euler (kick done after drift) -1=special epicyclic motion, full motion) -2=special epicyclic motion, shows only epi motion

[default:**3].** **tstop=***stop-time*- Time to stop integration in N-body model units.
Default is
**2.0**. **freqout=***out-freq*- Frequency of major N-body data outputs. Default
is
**4.0**(4 frames per unit time). **minor_freqout=***out-freq*- Frequency of minor
diagnostic outputs. Default is
**32.0**(32 diagnostic measurements per unit time). **options=***option-string*- Miscellaneous control options, specified as
a comma-separated list of keywords. Currently recognized keywords are:
**reset_time**: when reading initial data, set*tnow*to zero;**new_tout**: when restarting, set new output times;**mass**,**phi**,**acc**: output mass, potential, acceleration data with major data outputs. **eta=***loss-fraction*- Fraction of random energy
that is dissipated [Default:
**0.0**]. **cell=***box-size*- Cell size in which dissipation
is performed after every timestep. Dissipation is current performed on a
cartesian grid, in which cells are square (2D)
or a cube (3D)
. [Default:
**0.1**]. **rmax=***max_box-size*- Maximum size of the "box" (actually cube) within which
dissipation is performed. If a negative number is given, the box is allow
to grow as large as is needed, though this may consume a lot of memory.
Default:
**-1**, i.e. box can grow indefinite. **fheat=***fheat*- The ratio of diffusion
angle to rms velocity dispersion in a cell. If
*fheat>0*, each time dissipation is applied, the rms velocity dispersion in a cell is computed, and a diffusion angle computed. The velocity vector of each particle is then rotated with a gaussian distributed value with rms*fheat*velsig*. This in effect gives a dissipation dependant heating source. See also**sigma=**, which is position independant. [Default:**0**]. **sigma=***angle*- Diffusion angle, gaussian distributed
with this sigma, about which each velocity vector is rotated after each
timestep. [Default:
**0**]. **seed=***random-seed*- Random number seed, only used when
diffusion (
**sigma=**) is used. 0 must be used to get the random seed from the time of the day. [Default:**0**]. **headline=**- Identifying text for this run. Default:
not used.

Since **cell** is a fixed number throughout the execution, is doesn’t
deal well with systems who’s lenght-scale changes, in particular, expanding
systems will allocate more and more space to hold the dissipation grid.

For each cell the relative position and velocity
for each particle within that cell is computed:

R = r - <r > i i i V = v - <v > i i iafter which the dimensionless viscosity parameter ’alpha’ controls the new velocity for each particle after a timestep:

< R x V > j j u = <v> - alpha R x ------------ + (1-alpha) V i i i i < R . R > j j

https://ui.adsabs.harvard.edu/abs/2004MNRAS.352.1199M (mode=-1)

5-feb-89V1.0 writtenPJT 3-apr-90V2.0 new potential(5), renamed keywordsPJT 9-jun-92V3.0 rotating pattern are now allowed - made rk4PJT 17-jun-92V3.1 fix energy conservation in ’dissipation’PJT 19-jun-92V3.2 added diffusionPJT 6-oct-92V4.0 added rmax=, made options= same as hackcode1PJT 24-oct-92V4.1 added fheat=PJT 5-mar-03V5.0 added mode=-1 to "integrate" orbits numerically on epicyclic orbitsPJT 6-jul-03(V5.1) compute guiding centerPJT/RPO 12-aug-09V5.1 added leapfrog and modified eulerPJT 2-jul-21V5.2 fheat added, but not implementedPJT