ASTR415/688C Spring 2005 Problem Set #6

Due May 12, 2005

Generalize the integration package you wrote for PS5 so that it can solve the -body problem in 3-D for particles of arbitrary mass using the PP method with fixed timesteps. The code should read the masses and initial positions and velocities of all the particles from a text file (use the following input format, one particle per line: ). Use #defines or command-line arguments to specify the integrator (either leapfrog or Runge-Kutta), the softening parameter, the timestep, the number of steps, the output frequency, and the name of the initial conditions file. The output format can be the same as the input format, with one file per output (use sprintf() to imbed the current step number to create unique output filenames)--this has the advantage that you can restart from any output step. However, for problem 1 you may want to have the option of appending outputs to a single file, with one line per output step and each line listing data for all particles in order.

1. Check your code by solving the two-body problem for equal masses () without softening: plot phase diagrams of vs. , where is the separation and is the radial component of the relative velocity, for 100 loops around orbits of eccentricity and , using both leapfrog and Runge-Kutta. Start at unit apoapse: (you will need to derive the starting for yourself). By symmetry your initial conditions are then , , , , where we have chosen to start on the -axis with the orbital angular momentum aligned with the positive -axis. Use a timestep of 0.05 for and 0.003 for (plot approximately 1000 points for each case). Also plot the energy vs. time and comment.

You may find the following relations helpful for this problem:

Here is the semi-major axis and is the total mass. Assume . You will also need Newton's version of Kepler's Third Law to get the orbital period.

BONUS: Compare your results with the analytical solution. You will need to solve Kepler's Equation.

2. Generate random initial conditions of your choosing for between 100 and 1000 particles, e.g. particles in a spherical region, particles in a disk, whatever you want. Use a range of masses. Integrate this system long enough for something interesting'' to happen. What should the timestep be? Do you need softening? Is energy conserved? Does the code scale with as you expect? What is the megaflop rate? Summarize the evolution of your system in a plot of some kind to hand in.

BONUS: Make a movie of the results.