.. Created by Adam Cunnningham on Fri Jun 3 2016.

**The Magnetic Pendulum**
=========================

A pendulum is viewed from above, suspended from a pivot at the origin. Magnets at positions :math:`(x_i, y_i)` attract another magnet at the end of the pendulum. The position of the pendulum magnet is :math:`(x, y)`, and the state vector of the system is:

.. math::

    \vec{X} = \begin{bmatrix} x \\ y \\ \dot{x} \\ \dot{y} \end{bmatrix}

Three kinds of force act on the pendulum magnet: gravity, friction, and magnetic attraction.

Gravity
-------

Under the assumption that the horizontal displacement of the pendulum is small, the restoring force of gravity towards the origin is directly proportional to distance from the origin. The parameters governing gravity can be combined into a single coefficient *C*.

Friction
--------

For low velocities we make the simplifying assumption that frictional forces act opposite to and directly proportional to velocity. The parameters governing friction can be combined into a single coefficient *R*.

Magnetic Attraction
-------------------

The plane and pendulum magnets can be modelled as *magnetic dipoles* with magnetic moments :math:`\mathbf{m}` for the pendulum and :math:`\mathbf{m}_i` for plane magnet *i*. We make the simplifying assumption that the pendulum is sufficiently long that the height of the magnet over the plane is a constant *d* and that the dipoles are vertically aligned with each other.

There is a `standard formula <https://en.wikipedia.org/wiki/Force_between_magnets#Magnetic_dipole-dipole_interaction>`_ that we can use for the force between two magnetic dipoles separated by a distance vector :math:`\mathbf{r}`. To make use of this formula, we first express :math:`\mathbf{m}`, :math:`\mathbf{m}_i` and :math:`\mathbf{r}` in vector component form:

- :math:`\mathbf{m} = m\mathbf{k}`, where :math:`m = |\mathbf{m}|`
- :math:`\mathbf{m}_i = m_i\mathbf{k}`, where :math:`m_i = |\mathbf{m}_i|`
- :math:`\mathbf{r} = (x - x_i)\mathbf{i} + (y - y_i)\mathbf{j} + d\mathbf{k}`

Plugging these component vector representations into the standard formula yields the force *F* exerted by the *i* th plane magnet on the pendulum magnet:

.. math::
   :nowrap:

   \begin{eqnarray*}
    F(\mathbf{r}, \mathbf{m}_i, \mathbf{m}) & = & \frac{3\mu _0}{4\pi r^5}\left[(\mathbf{m}_i \cdot \mathbf{r} )\mathbf{m} + (\mathbf{m} \cdot \mathbf{r} )\mathbf{m}_i +(\mathbf{m}_i \cdot \mathbf{m}) \mathbf{r} - \frac{5(\mathbf{m}_i \cdot \mathbf{r} )(\mathbf{m} \cdot \mathbf{r})}{r^2} \mathbf{r}\right] \\
   & = & \frac{3\mu _0}{4\pi r^5}\left[m_i d m\mathbf{k} + m d m_i\mathbf{k} + m_i m ((x - x_i)\mathbf{i} + (y - y_i)\mathbf{j} + d\mathbf{k}) - 5 \frac{m_i m d^2}{r^2}((x - x_i)\mathbf{i} + (y - y_i)\mathbf{j} + h\mathbf{k})\right] \\
   & = & \frac{3m_i m \mu _0}{4\pi r^5}\left[(x - x_i)\left(1 - \frac{5 d^2}{r^2}\right)\mathbf{i} + (y - y_i)\left(1 - \frac{5 d^2}{r^2}\right)\mathbf{j} + \left(3d - \frac{5 d^3}{r^2}\right) \mathbf{k}\right]
   \end{eqnarray*}

The :math:`\mathbf{i}` and :math:`\mathbf{j}` components of the force *F* can then be used to directly calculate the :math:`x` and :math:`y` components of acceleration of the pendulum due to the plane magnets. For simplicity, the system is scaled so that the contribution to the acceleration of the constant term :math:`\frac{3m_i m \mu _0}{4\pi}` = 1.

Adding up the different forces, the differential equation governing the evolution of the system is:

.. math::

   \dfrac{d\vec{X}}{dt} = \begin{bmatrix} \dot{x} \\ \dot{y} \\ -R\dot{x} - Cx + \sum_i \frac{x_i - x}{r^5}\left(1 - \frac{5 d^2}{r^2}\right) \\ -R\dot{y} - Cy + \sum_i \frac{y_i - y}{r^5}\left(1 - \frac{5 d^2}{r^2}\right) \end{bmatrix}

Report Description
------------------

The report consists of the following exercises:

**Exercise 1**. Simulate the dynamics of this system using the improved Euler method, and plot the path of the pendulum from several different starting positions.

**Exercise 2**. Repeat step 1 using the same starting position with several different values for R, C, and d. Plot the path the pendulum takes on the same graph for comparison.

**Exercise 3**. Use the Matplotlib animation module to animate the path taken by the pendulum from one initial starting position.

**Exercise 4**. Assign a color to each magnet. Consider each pixel in the plane as a different starting point for the pendulum, and color the pixel according to the magnet the pendulum ends up over.

**Exercise 5**. Repeat step 4 using different values for R, and comment on the different kinds of image generated.

.. note:: The simulation is more accurate with a smaller step-size *h*, but takes longer to run. Explain your reasoning in choosing a particular value for *h*, and how you know that it is small enough.