.. _keplerian_system:

*******************************
Representing a Keplerian system
*******************************

To a good approximation the Earth--Sun system is Keplerian. The Earth moves on an elliptical orbit with semimajor axis :math:`a = 1~\text{AU}` and ellipticity :math:`e = 0.0167` in the gravitational central potential generated by the Sun, which has mass :math:`m = 1~\text{M}_{\odot}`. We may represent this system in Dyad using the class :class:`dyad.Orbit`.

.. doctest:: python

   >>> import dyad
   >>> m, a, e = 1., 1., 0.0167
   >>> orbit = dyad.Orbit(m, a, e)

The state (i.e. position and velocity) of the body is specified by its true anomaly, :math:`\theta`, which we can compute using the method :meth:`dyad.Orbit.state`. Suppose that :math:`\theta = 1`.

.. doctest:: python

   >>> theta = 1.
   >>> orbit.state(theta)
   array([ 0.53532139,  0.83371367,  0.        , -0.01447709,  0.00958295,
           0.        ])

This is given in the format :math:`(x, y, z, v_{x}, v_{y}, v_{z})`
where :math:`x`, :math:`y`, and :math:`z` are in units of
:math:`\mathrm{AU}` and :math:`v_{x}`, :math:`v_{y}`, and
:math:`v_{z}` are in units of :math:`\mathrm{AU}~\mathrm{d}^{-1}`. We
can also compute the orbital radius and speed of the body, which are
also in units of :math:`\mathrm{AU}` and
:math:`\mathrm{AU}~\mathrm{d}^{-1}`.

.. doctest:: python

   >>> orbit.radius(theta)
   0.9907812427855317
   >>> orbit.speed(theta)
   0.017361418667972937

An instance of the class :class:`dyad.Orbit` also has a number of
attributes that represent constants of the orbit. These include the arguments with which :class:`dyad.Orbit` was called.

.. doctest:: python

   >>> orbit.central_mass
   1.0
   >>> orbit.semimajor_axis
   1.0
   >>> orbit.eccentricity
   0.0167

They also include more interesting constants, such as the period (in :math:`\mathrm{d}`), periapsis (in :math:`\mathrm{AU}`) and apoapsis (in :math:`\mathrm{AU}`). A full list is available in the API documentation for :class:`dyad.Orbit`.

.. doctest:: python

   >>> orbit.period
   365.25689844815537
   >>> orbit.periapsis
   0.9833
   >>> orbit.apoapsis
   1.0167

If the perifocal frame differs from the observer's frame we can specify its orientation using the Euler angles :math:`\Omega`, :math:`i`, and :math:`\omega` (longitude of the ascending node, inclination, and argument of pericentre). Suppose that :math:`\Omega = 1`, :math:`i = 1`, and :math:`\omega = 1`.

.. doctest:: python

   >>> Omega, i, omega = 1., 1., 1.
   >>> orbit = dyad.Orbit(m, a, e, Omega, i, omega)

All methods and attributes are available as before.

We may not always wish to work with true anomaly but instead prefer mean anomaly or eccentric anomaly.
Dyad allows us to convert between these quantities.

.. doctest:: python

   >>> mu = dyad.mean_anomaly_from_true_anomaly(theta, e)
   >>> dyad.true_anomaly_from_mean_anomaly(mu, e)
   1.0
   >>> eta = dyad.eccentric_anomaly_from_true_anomaly(theta, e)
   >>> dyad.true_anomaly_from_eccentric_anomaly(eta, e)
   1.0
   >>> eta = dyad.eccentric_anomaly_from_mean_anomaly(mu, e)
   >>> dyad.mean_anomaly_from_eccentric_anomaly(eta, e)
   0.9720848452026836

To plot the evolution in time of any quantity first sample the mean anomaly
uniformly on the interval :math:`[0, 2\pi)` and convert it to true
anomaly. Consider, for example, the evolution of the orbital radius.

.. doctest:: python

   >>> import numpy as np
   >>> mu = np.linspace(0., 2.*np.pi)
   >>> theta = dyad.true_anomaly_from_mean_anomaly(mu, e)
   >>> r = orbit.radius(theta)

Plot this.

.. doctest:: python

   >>> import matplotlib.pyplot as plt
   >>> fig, ax = plt.subplots()
   >>> ax.plot(mu/(2.*np.pi), r)
   [<matplotlib.lines.Line2D object at 0x...>]
   >>> ax.set_xlabel(r"$t/\mathrm{yr}$")
   Text(0.5, 0, '$t/\\mathrm{yr}$')
   >>> ax.set_ylabel(r"$r/\mathrm{au}$")
   Text(0, 0.5, '$r/\\mathrm{au}$')
   >>> plt.show()

.. _radius:
.. figure:: ../figures/evolution_of_radius.jpg
   :figwidth: 75%
   :align: center

   The evolution of the Earth's radius over the course of a year.

.. _binary_system:

Delaunay elements
-----------------

Given the central mass we may convert any set of orbital elements, :math:`(a, e, \Omega, i, \omega, \theta)` to the equivalent set of Delaunay elements, :math:`(J_{1}, J_{2}, J_{3}, \Theta_{1}, \Theta_{2}, \Theta_{3})` and vice versa.
The functions are :func:`dyad.delaunay_elements_from_orbital_elements` and :func:`dyad.orbital_elements_from_delaunay_elements`.
See the API documentation for the definitions used by Dyad.

Delaunay elements are not fully defined for :math:`e = 0` or :math:`i = 0` (as we have here) so we must instead use modified Delaunay elements, :math:`(J_{\varpi}, J_{\Omega}, J_{\lambda}, \Theta_{\varpi}, \Theta_{\Omega}, \Theta_{\lambda})`.
The functions are :func:`dyad.modified_delaunay_elements_from_orbital_elements` and :func:`dyad.orbital_elements_from_modified_delaunay_elements`.
Again, see the API documentation for the definitions used by Dyad.

.. doctest:: python

   >>> m, a, e, Omega, i, omega, theta = 1., 1., 0.0167, 0., 0., 0., 1.
   >>> dyad.modified_delaunay_elements_from_orbital_elements(a, e, Omega, i, omega, theta, m)
   (2.398913957225284e-06, 0.0, 0.017202098944262438, -0.0, -0.0, 0.9720848452026836)
   >>> dyad.orbital_elements_from_modified_delaunay_elements(*_, m)
   (1.0, 0.01670000000000427, 0.0, 0.0, 0.0, 1.000000000000007)
   
In all threee cases (orbital elements, Delauany elements, and modified Delaunay elements) the first five elements define the orbit whilst the sixth defines the state. A list of the first five Delauanay elements is therefore available as a method of the :class:`dyad.Orbit` class.

.. doctest:: python

   >>> orbit.modified_delaunay_elements()
   array([ 2.39891396e-06,  7.90666244e-03,  1.72020989e-02, -2.00000000e+00,
	  -1.00000000e+00])

We can use the keyword ``theta`` to specify the true anomaly. In this case the function then returns :math:`\Theta_{\lambda}` also.

.. doctest:: python

   >>> orbit.modified_delaunay_elements(theta=1.)
   array([ 2.39891396e-06,  7.90666244e-03,  1.72020989e-02, -2.00000000e+00,
	  -1.00000000e+00,  2.97208485e+00])
   
****************************
Representing a binary system
****************************

To a good approximation the Alpha Centauri A--B system is an isolated binary. The two component stars have masses :math:`m_{A} = 1.0790~\text{M}_{\odot}` and :math:`m_{B} = 0.9092~\text{M}_{\odot}`. The relative orbit has a semimajor axis :math:`a = 23.182~\text{AU}` and an eccentricity :math:`e = 0.5195`. We may represent this system using the class :class:`dyad.TwoBody`.

.. doctest:: python

   >>> import dyad
   >>> m_A, m_B, a, e = 1.079, 0.9092, 23.182, 0.5195
   >>> binary = dyad.TwoBody(m_A, m_B, a, e)

The properties of the relative orbit are available as an instance of the class :class:`dyad.Orbit` with central mass :math:`m = m_{A} + m_{B}`. The class :class:`dyad.TwoBody` has all the attributes of the class :class:`dyad.Orbit` (see :ref:`keplerian_system`) plus three additional attributes: ``reduced_mass``, ``primary``, and ``secondary``. The attribute ``reduced_mass`` holds the reduced mass of the secondary star, :math:`\mu = m_{1}m_{2}/(m_{1} + m_{2})` while the attributes ``primary`` and ``secondary`` are themselves instances of the class :class:`Orbit` and represent the properties of the primary and secondary stars' orbits. The :class:`Orbit` instance ``primary`` has the additional attribute ``mass``, which holds the primary star mass. Similarly, the :class:`Orbit` instance ``secondary`` has the additional attribute ``mass``, which holds the secondary star mass.

.. doctest:: python

   >>> binary.reduced_mass
   0.49342460517050596
   >>> binary.primary
   <dyad._core.Orbit object at 0x...>
   >>> binary.secondary
   <dyad._core.Orbit object at 0x...>

For example, we may find the state of the secondary star in coordinates relative to the primary as follows.

.. doctest:: python

   >>> binary.state(1.)
   array([ 7.14066735e+00,  1.11209305e+01,  0.00000000e+00, -4.96110162e-03,
	   6.24832826e-03,  0.00000000e+00])

The properties of the primary and secondary stars are represented in the in the observer's frame, with the origin at the primary focus. They have masses as follows.

.. doctest:: python

   >>> binary.primary.mass
   1.079
   >>> binary.secondary.mass
   0.9092

We may wish to know the orbital state of the two stars when :math:`\theta = 1`.

.. doctest:: python

   >>> binary.primary.state(1.)
   array([ 3.26541331e+00,  5.08557992e+00,  0.00000000e+00, -2.26870214e-03,
           2.85734838e-03,  0.00000000e+00])
   >>> binary.secondary.state(1.)
   array([-3.87525403e+00, -6.03535056e+00, -0.00000000e+00,  2.69239948e-03,
	  -3.39097988e-03, -0.00000000e+00])

We might also like to know the system's other attributes. For example,
the periods and periapses of the two stars.

.. doctest:: python

   >>> binary.primary.period
   28913.10574161869
   >>> binary.secondary.period
   28913.105741618692
   >>> binary.primary.periapsis
   5.093820666532542
   >>> binary.secondary.periapsis
   6.045130333467458

Of course we might equivalently compute these last two values as follows.

.. doctest:: python

   >>> binary.primary.radius(0.)
   5.093820666532543
   >>> binary.secondary.radius(0.)
   6.045130333467459

Note that all properties of the primary and secondary orbits are
computed in the observer's frame.