White paper for ring.f   --  Steven Janowiecki, 2-Dec-2007

ring -- makes a ring galaxy by colliding a satellite with a parent disk galaxy

General Description

This program takes an initial disk galaxy as an input, and collides a smaller satellite galaxy with it. The potentials for both galaxies are specified, as are their total masses and starting positions and velocities. The positions and velocities of the galaxy centers and each particle are outputted periodically, as the simulation evolves.


Inputs:

This code requires an input galaxy file, with 6 columns - position and velocity, in x,y,z components. Each line of the input file is considered an equal mass particle. (In subsequent output files, the first two lines are reserved for the parent and satellite galaxy centers, and the particles follow.) The standard input file is generated by the program "md", to make the disk. 

There is a collection of input parameters at the beginning of the code that need to be set before it is compiled and executed. After specifying the file where the initial galaxy particles are, the user needs to give the number of particles in that galaxy (npart) and details about the time steps (step size, total number of steps, how often to write output files along the way). The parameters for the galaxy potentials must also be set (r0, v0, cutr), and are explained in detail, below. The mass ratio is set with (satgal) and is M_sat/M_gal. Finally the initial positions and velocities for the main and satellite galaxies are set with the arrays wgal and wsat, both 6 elements long (x,y,z,v_x,v_y,v_z).


Potentials and Accelerations:

The main galaxy is assigned a potential that is a composite of a logarithmic and Keplerian potential, with a cutoff radius to transition between the two (set by cutr). Inside this radius, the potential is:

phi_inner(r) = v0^2 ln ( (r/r0)^2 )

And outside the cutoff it is:

phi_outer(r) = GM/r + C

C is a constant to make the potentials continuous at the cutoff radius (in our application, the value of the constant is irrelevant since the important quantity is the gradient of the potential, and the constant drops out).

The satellite galaxy is assigned a purely Keplerian potential, and treated as a point mass.

The component-by-component accelerations are calculated from these potentials as:

Logarithmic: a_i = -2 v0^2 i / r^2
Keplerian: a_i = -G M_tot i / r^3

For each time step, the amount a particle gets moved is determined using a 4th order Runge-Kutta integrator, from Numerical Recipes. The above equations for accelerations are used to calculate how far to move each particle at each step.

In the code, when the Runge-Kutta integrator is called, the command requires you to tell it which potentials to take into consideration. The variable 'pp' indicates which potential, as follows:
pp = 1 : use only the main galaxy's potential
pp = 2 : use only the satellite's potential
pp = 3 : use both the galaxy and satellite potentials

The Runge-Kutta integrator further requires the current position of the galaxy and satellite centers, as well as their masses, to properly compute the acceleration. Finally, it takes the value of v0 and cutr, in order to scale the potential and properly construct the composite potential for the main galaxy, respectively.

DERIVS is a small subroutine which the RK4 integrator relies on, and is where the aforementioned physics is applied (the accelerations in the potentials). The calculation of derivatives depends on where the particle is and which potentials you are using (see parameter pp).


Time Steps and Overall Calculations:

For each time step, the following is a flowchart of the progression of events:

1) A snapshot file of all the particles is written to disk. See "Output" below.

2) The main galaxy center (center of its potential) is stepped with respect to the satellite galaxy. That is, they are both treated as point masses and the main one is allowed to move closer to the satellite.

3) The satellite galaxy center is stepped with respect to the main galaxy. Again, they are both treated as point masses, and the satellite is allowed to move as a result of the force from the main galaxy.

4) Each particle is stepped for this time step. Each particle feels the potentials of both the main and satellite galaxy.


Output

At the beginning of each time step, if the modulo of the number of the current step and the variable "outrate" is 0, then a snapshot of every particle's position and velocity (in x,y,z) is written to a file. These files are automatically numbered, starting with the initial conditions at file0000, and continuing to file0001, or whatever the stepsize and outrate demand. These files will have one line per particle, and the first and second lines will be the main and satellite galaxies' positions and velocities, respectively.


See also: md - a program to initialise a rotating exponential disk of particles

Author:
Steven Janowiecki: spj4@cwru.edu