.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/classical/newtonian/plot_03_kepler_perturbations.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_classical_newtonian_plot_03_kepler_perturbations.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_classical_newtonian_plot_03_kepler_perturbations.py:


KeplerSystem's other perturbation: the power-law term
==========================================================

:doc:`plot_02_kepler_precession` shows apsidal precession from the
effective post-Newtonian :math:`c_\mathrm{pn}/r^3` correction.
:class:`~physicskit.classical.systems.newtonian.KeplerSystem` has a
second, independent perturbation family added to its radial potential:

.. math::

    V(r) = -\frac{k}{r} + \frac{c_\varepsilon}{r^{1+\varepsilon}} ,

an adjustable power-law addition to the pure :math:`1/r` potential. At
:math:`\varepsilon=0` this is degenerate -- it is just another
:math:`1/r` term, rescaling the effective coupling constant, so by
Bertrand's theorem the orbit stays a perfectly closed, non-precessing
ellipse. For any :math:`\varepsilon \neq 0`, Bertrand's theorem no
longer applies and the orbit precesses, faster for larger
:math:`\varepsilon` -- shown here for three exponents at the same
coefficient, all still exactly energy-conserving.

.. GENERATED FROM PYTHON SOURCE LINES 25-48

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from physicskit.classical.systems.newtonian import KeplerSystem

    fig, axes = plt.subplots(1, 3, figsize=(13, 4.3))
    for ax, eps in zip(axes, (0.5, 1.0, 2.0)):
        system = KeplerSystem.from_orbital_elements(a=1.0, e=0.4, eps=eps, c_eps=0.01)
        lrl0 = system.lrl_vector()
        result = system.integrate((0, 100), dt=1e-3, method="yoshida4")

        lrl_vecs = np.array([system.lrl_vector(q, p) for q, p in zip(result.q, result.p)])
        angle = np.arctan2(lrl_vecs[:, 1], lrl_vecs[:, 0]) - np.arctan2(lrl0[1], lrl0[0])
        drift = np.max(np.abs(result.energy - result.energy[0])) / abs(result.energy[0])
        print(f"eps={eps}: precession over t=100 = {np.unwrap(angle)[-1]:+.3f} rad, energy drift = {drift:.2e}")

        ax.plot(result.q[:, 0], result.q[:, 1], color="steelblue", lw=0.6)
        ax.plot(0, 0, "o", color="gold")
        ax.set_title(f"eps = {eps}")
        ax.set_aspect("equal")
    fig.suptitle(r"$V(r) = -k/r + c_\varepsilon / r^{1+\varepsilon}$: precession grows with $\varepsilon$", fontsize=13)
    fig.tight_layout(rect=[0, 0, 1, 0.92])




.. image-sg:: /api/gallery/classical/newtonian/images/sphx_glr_plot_03_kepler_perturbations_001.png
   :alt: $V(r) = -k/r + c_\varepsilon / r^{1+\varepsilon}$: precession grows with $\varepsilon$, eps = 0.5, eps = 1.0, eps = 2.0
   :srcset: /api/gallery/classical/newtonian/images/sphx_glr_plot_03_kepler_perturbations_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    eps=0.5: precession over t=100 = -0.362 rad, energy drift = 5.71e-12
    eps=1.0: precession over t=100 = -1.057 rad, energy drift = 5.71e-12
    eps=2.0: precession over t=100 = -3.237 rad, energy drift = 5.05e-12




.. GENERATED FROM PYTHON SOURCE LINES 49-54

See :doc:`plot_04_conservation_diagnostics` for a related point: even
with this perturbation switched on, the orbit's angular momentum -- a
*different* conserved quantity from energy -- stays exactly constant,
since ``V(r)`` depends only on ``r`` and the force therefore still
exerts no torque about the origin.

.. GENERATED FROM PYTHON SOURCE LINES 54-56

.. code-block:: Python


    plt.show()








.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 2.554 seconds)


.. _sphx_glr_download_api_gallery_classical_newtonian_plot_03_kepler_perturbations.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_03_kepler_perturbations.ipynb <plot_03_kepler_perturbations.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_03_kepler_perturbations.py <plot_03_kepler_perturbations.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_03_kepler_perturbations.zip <plot_03_kepler_perturbations.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_