How can I calculate the right obit period of planets with python?

Question

Here is my python code which can plot the 2D orbit of Earth, Mars and Mercury with Sun. However, I cannot calculate the proper orbit of Earth and others .

attempt: I use the $r_\max$, $r_\min$ of orbit, mass of Earth and $G$ constant as input. Using the formula of $T=(a^3/mu)^{0.5}$ to calculate the orbit period.

However,the calculated orbit is around $91667$ which is not around the theoretical value like $365$ days.

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, Circle
import numpy as np

#Set axes aspect to equal as orbits are almost circular; hence square is needed
ax = plt.figure(0).add_subplot(111, aspect='equal')

#Setting the title, axis labels, axis values and introducing a grid underlay
#Variable used so title can indicate user inputed date
plt.title('Inner Planetary Orbits at[user input date]')
plt.ylabel('x10^6 km')
plt.xlabel('x10^6 km')
ax.set_xlim(-300, 300)
ax.set_ylim(-300, 300)
plt.grid()

#Creating the point to represent the sun at the origin (not to scale), 
ax.scatter(0,0,s=200,color='y')
plt.annotate('Sun', xy=(5,-30))
plt.annotate('Mercury', xy=(30,-75))
plt.annotate('Earth', xy=(30,-180))
plt.annotate('Mars', xy=(30,-250))

#Implementing ellipse equations to generate the values needed to plot an ellipse
#Using only the planet's min (m) and max (M) distances from the sun
#Equations return '2a' (the ellipses width) and '2b' (the ellipses height)
def OrbitLength(M, m):
    a=(M+m)/2
    c=a-m
    e=c/a
    b=a*(1-e**2)**0.5
    
    return 2*a, 2*b

#This function uses the returned 2a and 2b for the ellipse function's variables
#Also generating the orbit offset (putting the sun at a focal point) using M and m
def PlanetOrbit(Name, M, m):
    w, h = OrbitLength(M, m)
    Xoffset= ((M+m)/2)-m
    Name = Ellipse(xy=((Xoffset),0), width=w, height=h, angle=0, linewidth=1, fill=False)
    ax.add_artist(Name)

    
    
from math import *

EPSILON = 1e-12
def solve_bisection(fn, xmin,xmax,epsilon = EPSILON):
  while True:
      xmid = (xmin + xmax) * 0.5
      if (xmax-xmin < epsilon):
        return xmid
      fn_mid = fn(xmid)
      fn_min = fn(xmin)
      if fn_min*fn_mid < 0:
          xmax = xmid
      else:
          xmin = xmid
    
'''
Found something similar at this gamedev question:
https://gamedev.stackexchange.com/questions/11116/kepler-orbit-get-position-on-the-orbit-over-time?newreg=e895c2a71651407d8e18915c38024d50

Equations taken from:
https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion#Position_as_a_function_of_time
'''
def SolveOrbit(rmax, rmin, t):
  # calculation precision
  epsilon = EPSILON
  # mass of the sun [kg]
  Msun = 1.9891e30
  # Newton's gravitational constant [N*m**2/kg**2]
  G = 6.6740831e-11
  # standard gravitational parameter
  mu = G*Msun
  # eccentricity
  eps = (rmax - rmin) / (rmax + rmin)
  # semi-latus rectum
  p = rmin * (1 + eps)
  # semi/half major axis
  a = p / (1 - eps**2)
  
  # period
  P = np.sqrt(a**3 / mu)
  print(p)
  # mean anomaly
  M = (t / P) % (2*pi)
  # eccentric anomaly
  def fn_E(E):
    return M - (E-eps*np.sin(E))
  E = solve_bisection(fn_E, 0, 2*pi)
  # true anomaly
  
  theta = 2*atan(sqrt((((1+eps)*tan(E/2)**2)/(1-eps))))
  # if we are at the second half of the orbit
  if (E > pi):
    theta = 2*pi - theta
  # heliocentric distance
  r = a * (1 - eps * cos(E))
  return theta, r

def DrawPlanet(name, rmax, rmin, t):
  SCALE = 1e9
  theta, r = SolveOrbit(rmax * SCALE, rmin * SCALE, t)
  x = -r * cos(theta) / SCALE
  y = r * sin(theta) / SCALE
  planet = Circle((x, y), 8)
  ax.add_artist(planet)


#These are the arguments taken from hyperphysics.phy-astr.gsu.edu/hbase/solar/soldata2.html
#They are the planet names, max and min distances, and their longitudinal angle
#Also included is Halley's Comet, used to show different scale  and eccentricity
PlanetOrbit('Mercury', 69.8, 46.0)
PlanetOrbit('Earth', 152.1, 147.1)
PlanetOrbit('Mars', 249.1, 206.7)

for i in range(0, 52):
  DrawPlanet('Earth', 152.1, 147.1, i/52 * 365.25 *60*60*24)
  
plt.show()

rmax=152.1*10**6
rmin=147.1*10**6
epsilon = EPSILON
# mass of the sun [kg]
Mearth = 5.97e24
# Newton's gravitational constant [N*m**2/kg**2]
G = 6.6740831e-11
# standard gravitational parameter
mu = G*Mearth
# eccentricity
eps = (rmax - rmin) / (rmax + rmin)
# semi-latus rectum
z= rmin * (1 + eps)
# semi/half major axis
a = z/ (1 - eps**2)

# period
print('period is',np.sqrt(a**3/mu))

Can anyone help me to calculate correctly?Thank you!

Answer 1

$H=(1/2)m\dot{r}^2+(1/2)mr^2\dot{\theta}^2+GM/r$

$L=mr^2\dot{\theta}$ by conservation of angular momentum.

$m\ddot{r}=mr\dot{\theta}^2-GM/r^2$

$m\ddot{r}/\dot{\theta}=L/r-GMm/L$

$-\frac{d^2}{d\theta^2}(\frac{L}{r})-\frac{L}{r}=-GMm/L$

$\frac{d^2}{d\theta^2}(\frac{1}{r})+\frac{1}{r}=GMm/L^2$

Solving the differential equation:

$\frac{1}{r}=c_1\cos\theta+c_2\sin \theta+GMm/L^2$, the equation of an ellipse with a focus at the origin.

WLOG $c_2=0$. Setting $\theta$ to 0 and $\pi$ to find perihelion and aphelion.

Aphelion and Perihelion happen on a straight line so one can solve for $c_1$ by letting $\theta$ be $0$ and $\pi$ and setting the equation equal to closest and farthest distance from the sun.

Given that $L=mr^2\frac{d\theta}{dt}\implies \tau=\int_0^{2\pi} \frac{mr^2d\theta}{L}=\int_0^{2\pi} \frac{m}{L}\frac{d\theta}{[(c_1\cos \theta)+GMm/L^2]^2}$ gives you the period.

Of course this breaks up a multibody system inaccurately into a series of two body problems which will necessarily remove some import long term dynamics. To account for that you'll need some non-analytical methods. I seem to recall a leapfrogging runge-kutta approach gives decent results.