3. Solving Ordinary Differential Equations

Consider bodies interacting through Newtonian gravity.  That is, the gravitational potential between two bodies of mass m_A and m_B that are separated by a distance vector r_AB= r_A - r_B is given by

$$ V(r_{AB}) = - \frac{Gm_Am_B}{|\vec r_{AB}|} $$

and the force on body A would be given by

$$ \vec F_A(r_{AB}) = - \frac{Gm_am_b}{|\vec r_{AB}|^3}\vec r_{AB} $$

The force on body B would be equal and opposite.

Here, G can be taken to be approximately 6.7 x 10^-11 m^3/kg/s^2.

There are three parts to this assignment:

1. Write a module that contains functions to compute the interaction of two bodies (potential and force), and to compute the interactions of three or more bodies.
2. Use this module together with the ODE solvers from the GNU Scientific Library to write a code that solves the equations of motion: 

$$\frac{d\vec p_A}{dt} = \vec F_A $$
$$   \frac{d\vec p_B}{dt} = \vec F_B $$

for the following parameters:
$$  m_A = 2 \times 10^{30} kg; m_B = 6 \times 10^{24} kg $$
and the following initial conditions:
$$\vec r_A(0)= (0,0,0),\vec r_B(0) = (1.48 \times 10^{11} m,0,0) $$
$$\vec p_A(0)= (0,0,0),\vec p_B(0) = (0, m_B \times 3\times 10^{4} m/s,0) $$

The code should write out the result either in ASCII or in NetCDF format (your choice).  Read the GSL docs and pick an integrator that does not require a "Jacobian". From the results, estimate the period P of the orbit this forms.

3. The orbit found in part 2 should be circular, and body A should be effectively stationary.  Thus, In a co-rotating coordinate system with A as its center and rotating with a period P, both bodies would seem to be fixed. An additional test particle C, i.e., a particle of very low mass, e.g. m_C=6000 kg, would, however, feel an additional centrifugal force 

   $$m_C\omega^2\vec r_{AC}$$, 

   where 

   $$\omega=2\pi/P$$

   We can say that the particle $C$ feels an effective potential

   $$V_{eff}(r_C) = V(r_{AB}) + V(r_{BC}) + V(r_{AC}) - \frac{m_C}{2}\omega^2|\vec r_{AC}|^2$$

   where you may assume A and B are not moving in the co-rotating frame. Write a program to find, using the root finding machinery of the GSL and the module from part 1, the values of r where the derivative of V_eff is zero (these are called the Lagrange points). You need to only consider situations where A, B and C are aligned in the co-rotating frame.

Use git and other best practices throughout this assignment. 

Submit all the code for these three parts, as well as the Makefile and the git log by Friday February 18, at midnight.