arnold orbit

arnold-orbit.py

from functools import partial

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp


def U(r):
    return 0.5 * r ** 3


def V(r):
    return U(r) + 0.5 / r ** 2


def dphi(r, E):
    return 1 / r ** 2 * (1 / np.sqrt(2 * (E - V(r))))


r_min = 0.50842208660526406332007805917088124001667877021803
r_max = 1.5286429152095940020185327716896748970412268021517
E = 2

fun = partial(dphi, E=E)
r_eval = np.arange(r_min, r_max, 0.001)
sol = solve_ivp(lambda r_, phi: fun(r_), [r_min, r_max], [0], t_eval=r_eval)
phi = sol.y.flatten()
r = sol.t

ax = plt.subplot(111, projection="polar")
for i in range(26):
    ax.plot(phi, r, color="k")
    r = r[::-1]
    phi = 2 * phi[-1] - phi[::-1]

circle_phi = np.arange(0, 2 * np.pi, 0.01)
ax.plot(circle_phi, r_min * np.ones(circle_phi.shape), color="k")
ax.plot(circle_phi, r_max * np.ones(circle_phi.shape), color="k")
ax.axis("off")
plt.show()