Exercise 3
==========

Forward integration
-------------------

In ``class01``, we wrote a julia program ``orbit.jl`` that simulates a 
*two-dimensional* gravitating system. It describes a test mass 
:math:`\mathsf{m}` in orbit around a heavy mass :math:`\mathsf{M}`,
following the usual :math:`\mathsf{1/r^2}` law and also subject
to a velocity-dependent drag force:

.. math::

   \mathsf{m\frac{d^2\vec{r}}{dt^2}} 
   = \mathsf{m\ddot{\vec{r}}}
   = \mathsf{-\frac{GMm}{r^2}\hat{r} - C_0 m\vec{v}}
   = \mathsf{-\frac{GMm}{r^3}\vec{r} - C_0 m\vec{v}.}
   
Here, :math:`\mathsf{\vec{r}(t)} = \mathsf{x\hat{x} + y\hat{y}}` 
is the position vector, and :math:`\mathsf{\vec{v}(t)} = 
\mathsf{\dot{\vec{r}}} = \mathsf{\dot{x}\hat{x} + \dot{y}\hat{y}}
= \mathsf{v_x\hat{x} + v_y\hat{y}}` is the vector velocity. In units
where ``GM = 1`` and ``C_0 = 0.01``, the initial conditions are
:math:`\mathsf{\vec{r}(0)} = \mathsf{\hat{x}}` and 
:math:`\mathsf{\vec{v}(0)} = \mathsf{\hat{y}}` 
(corresponding to the planet starting at the "eastern" edge of the orbit 
and moving "north"). These are chosen to produce a near-circular orbit 
that slowly spirals in toward the central mass :math:`\mathsf{M}`.

In the code, we don't use vectors. Rather, we define individual 
floating-point variables defined for the :math:`\mathsf{x}` and
:math:`\mathsf{y}` components of the position (``x`` and ``y``), velocity (``vx`` and ``\mathsf{vy}``), and acceleration
(``ax`` and ``ay``). Breaking the 
vector differential equation into components gives

.. math::

    \begin{aligned}
    \mathsf{\dot{x}} &= \mathsf{v_x},\\
    \mathsf{\dot{y}} &= \mathsf{v_y},\\
    \mathsf{\ddot{x}} &= \mathsf{a_x} 
                       = \mathsf{-GM \frac{x}{(x^2+y^2)^{3/2}} - C_0 v_x}, \\
    \mathsf{\ddot{y}} &= \mathsf{a_y} 
                       = \mathsf{-GM \frac{y}{(x^2+y^2)^{3/2}} - C_0 v_y.}
    \end{aligned}

Integrating over small finite steps of duration :math:`\mathsf{\Delta t}`
produces

.. math::

    \begin{aligned} 
    \mathsf{\int_{t}^{t+\Delta t}\!dt'\,\dot{\vec{r}}(t')} 
    &= \mathsf{\vec{r}(t+\Delta t) - \vec{r}(t)},\\
    \mathsf{\int_{t}^{t+\Delta t}\!dt'\,\ddot{\vec{r}}(t')} 
    &= \mathsf{\dot{\vec{r}}(t+\Delta t) - \dot{\vec{r}}(t)} \\
    &= \mathsf{-\int_{t}^{t+\Delta t}\!dt'\,\biggl(\frac{GM\vec{r}(t')}{r(t')^3} + C_0 \vec{v}(t') \biggr)}\\
    &= \mathsf{-\biggl(\frac{GM\vec{r}(t^*)}{r(t^*)^3} + C_0 \vec{v}(t^*) \biggr)\Delta t},
    \end{aligned}

where the final equality is guaranteed to hold (by a theorem of 
calculus, very easy to justify geometrically) for some instant in 
time satisfying :math:`\mathsf{t < t^* < t + \Delta t}`.

Verlet algorithm
----------------

A quick and dirty (but surprisingly effective) approximation is to put 
:math:`\mathsf{t^* = t + \Delta t/2}` at the interval midpoint and to 
update the position in two half-steps; this leads to the famous 
`Verlet algorithm <https://en.wikipedia.org/wiki/Verlet_integration>`_,
viz.,

.. math::

    \begin{aligned} 
    \mathsf{x(t^*)} &:= \mathsf{x(t) + v_x(t)(\Delta t/2)}\\
    \mathsf{y(t^*)} &:= \mathsf{y(t) + v_y(t)(\Delta t/2)}\\
    \mathsf{v_x(t+\Delta t)} &:= \mathsf{v_x(t) - 
       \biggl(\frac{GMx(t^*)}{[x(t^*)^2+y(t^*)^2]^3} 
       + C_0 v_x(t) \biggr)\Delta t},\\
    \mathsf{v_y(t+\Delta t)} &= \mathsf{v_y(t) -
       \biggl(\frac{GMy(t^*)}{[x(t^*)^2+y(t^*)^2]^3} 
       + C_0 v_y(t) \biggr)\Delta t},\\
    \mathsf{x(t+\Delta t)} &= \mathsf{x(t^*) + v_x(t+\Delta) (\Delta t/2)}\\
    \mathsf{y(t+\Delta t)} &= \mathsf{y(t^*) + v_y(t+\Delta) (\Delta t/2)}\\
    \end{aligned}

Some things to try
------------------

+  Try generalizing the program from two to three dimensions, and
   give the test mass a small out-of-plane initial velocity, 
   :math:`\mathsf{z(0) = 0.1}`. Let's call that code ``orbit3D.jl``.

+  Try converting your julia program to python. Confirm that
   ``orbit3D.jl`` and the new ``orbit3D.py`` both produce the same
   trajectory.
   
   .. code-block:: console
   
       $ julia orbit3D.jl > traj1.dat
       $ python3 orbit3D.py > traj2.dat
       $ gnuplot
       > plot for [n=1:2] "traj".n.".dat" using 1:2 with lines

+  Notice that the stopping condition (``r < R0``) terminates
   the program when the test mass hits the surface of the "planet"
   (we model the large mass :math:`\mathsf{M}` as a sphere or radius 
   :math:`\mathsf{R_0}`). What if that stopping condition is removed? 
   In physical terms, what if :math:`\mathsf{M}` is a point mass with no 
   spatial extent? See what the outcome is and try to relate it to your
   understanding of the floating-point representation.