ASTR615 Fall 2015 Problem Set #6

Due Mon Dec 7th, 2015

Generalize the integration package you wrote for PS5 so that it can solve the $N$-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: $m_i x_i y_i z_i \dot{x}_i \dot{y}_i \dot{z}_i$). 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. (Worth 85% of the grade) Check your code by solving the two-body problem for equal masses ($m_1 = m_2 = 1$) without softening: plot phase diagrams of $r$ vs. $v_r$, where $r$ is the separation and $v_r = (\mathbf{v} \mbox{\boldmath$\cdot$}\mathbf{r}) / r$ is the radial component of the relative velocity, for 100 loops around orbits of eccentricity $e = 0.5$ and $e = 0.9$, using both leapfrog and Runge-Kutta. Start at unit apoapse: $r = r_a = 1$ (you will need to derive the starting $v$ for yourself). By symmetry your initial conditions are then $\mathbf{r}_1 = (-1/2,0,0)$, $\mathbf{r}_2 = (1/2,0,0)$, $\mathbf{v}_1 = (0,-v/2,0)$, $\mathbf{v}_2 = (0,v/2,0)$, where we have chosen to start on the $x$-axis with the orbital angular momentum aligned with the positive $z$-axis. Use a timestep of 0.05 for $e = 0.5$ and 0.003 for $e = 0.9$ (plot approximately 1000 points for each case). Also plot the energy $E = \frac{1}{2} v_1^2 + \frac{1}{2} v_2^2 - 1/r$ vs. time and comment.

   You may find the following relations helpful for this problem:
   

   $$\frac{1}{a} = \frac{2}{r} - \frac{v^2}{GM}; \hspace{0.5in} r_a \equiv (1 + e)a .$$
   
   
   Here $a$ is the semi-major axis and $M = m_1 + m_2 = 2$ is the total mass. Assume $G \equiv 1$. You will also need Newton's version of Kepler's Third Law to get the orbital period.

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

2. (Worth 15% of the grade) 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 $N$ as you expect? Summarize the evolution of your system in a plot of some kind to hand in.

   Make a movie of the results.