Orbit games#

We consider energy plots and orbital solutions in polar coordinates for the general potential energy

\(\begin{align} U(r) = k r^n \end{align}\)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# Change the common font size
font_size = 14
plt.rcParams.update({'font.size': font_size})
class Orbit:
    """
    Potentials and associated differential equations for central force motion
    with the potential U(r) = k r^n.
    """
    
    def __init__(self, ang_mom, n, k=1, mu=1):
        self.ang_mom = ang_mom
        self.n = n
        self.k = k
        self.mu = mu
    
    def U(self, r):
        """Potential energy of the form U = kr^n."""
        return self.k * r**self.n
    
    def Ucf(self, r):
        """Centrifugal potential energy"""
        return self.ang_mom**2 / (2. * self.mu * r**2)
    
    def Ueff(self, r):
        """Effective potential energy"""
        return self.U(r) + self.Ucf(r)
    
    def U_deriv(self, r):
        """dU/dr"""
        return self.n * self.k * r**(self.n - 1)
        
    def Ucf_deriv(self, r):
        """dU_cf/dr"""
        return -2. * self.ang_mom**2 / (2. * self.mu * r**3)
        
    def Ueff_deriv(self, r):
        """dU_eff/dr"""
        return self.U_deriv(r) + self.Ucf_deriv(r)
        
    def dy_dt(self, t, y):
        """
        This function returns the right-hand side of the diffeq: 
        [dr/dt d^2r/dt^2 dphi/dt]
        
        Parameters
        ----------
        t : float
            time 
        y : float
            3-component vector with y[0] = r(t), y[1] = dr/dt, y[2] = dphi/dt
            
        """
        return [ y[1], 
                -1./self.mu * self.Ueff_deriv(y[0]), 
                self.ang_mom / (self.mu * y[0]**2) ]
    
    
    def solve_ode(self, t_pts, r_0, r_dot_0, phi_0, 
                  abserr=1.0e-8, relerr=1.0e-8):
        """
        Solve the ODE given initial conditions.
        For now use odeint, but we have the option to switch.
        Specify smaller abserr and relerr to get more precision.
        """
        y = [r_0, r_dot_0, phi_0]  
        solution = solve_ivp(self.dy_dt, (t_pts[0], t_pts[-1]), 
                             y, t_eval=t_pts, 
                             atol=abserr, rtol=relerr)
        r, r_dot, phi = solution.y
        return r, r_dot, phi
    
    def energy(self, t_pts, r, r_dot):
        """Evaluate the energy as a function of time"""
        return (self.mu/2.) * r_dot**2 + self.Ueff(r)
def start_stop_indices(t_pts, plot_start, plot_stop):
    start_index = (np.fabs(t_pts-plot_start)).argmin()  # index in t_pts array 
    stop_index = (np.fabs(t_pts-plot_stop)).argmin()  # index in t_pts array 
    return start_index, stop_index

Gravity \(n = -1\)#

n = -1
k = -3. 
ang_mom = 1. 
o1 = Orbit(ang_mom, n=n, k=k, mu=1)

fig_2 = plt.figure(figsize=(7,5))
ax_2 = fig_2.add_subplot(1,1,1)

r_pts = np.linspace(0.001, 3., 200)
U_pts = o1.U(r_pts)
Ucf_pts = o1.Ucf(r_pts)
Ueff_pts = o1.Ueff(r_pts)

ax_2.plot(r_pts, U_pts, linestyle='dashed', color='blue', label='U(r)')
ax_2.plot(r_pts, Ucf_pts, linestyle='dotted', color='green', label='Ucf(r)')
ax_2.plot(r_pts, Ueff_pts, linestyle='solid', color='red', label='Ueff(r)')

ax_2.set_xlim(0., 3.)
ax_2.set_ylim(-5., 5.)
ax_2.set_xlabel('r')
ax_2.set_ylabel('U(r)')
ax_2.set_title(f'$n = {n},\ \ k = {k},\ \  l = {ang_mom}$')
ax_2.legend(loc='upper center')

ax_2.axhline(0., color='black', alpha=0.3)

fig_2.tight_layout()

fig_2.savefig('Gravitation_orbit_1.png')
../../_images/Orbit_games_6_0.png

Orbit (time dependence)#

We’ll directly solve the equations for r(t) and phi(t).

# Plotting time 
t_start = 0.
t_end = 10.
delta_t = 0.01

t_pts = np.arange(t_start, t_end+delta_t, delta_t)  

# Initial conditions
r_0 = 1.
r_dot_0 = 0.
phi_0 = 0.0
r_pts, r_dot_pts, phi_pts = o1.solve_ode(t_pts, r_0, r_dot_0, phi_0)

c = o1.ang_mom**2 / (np.abs(o1.k) * o1.mu)
epsilon = c / r_0 - 1.
energy_0 = o1.mu/2. * r_dot_0**2 + o1.Ueff(r_0)
print(f'energy = {energy_0:.2f}')
print(f'eccentricity = {epsilon:.2f}')
energy = -2.50
eccentricity = -0.67
fig_4 = plt.figure(figsize=(8,8))

overall_title = 'Gravitational orbit:  ' + \
                rf' $n = {o1.n},$' + \
                rf' $k = {o1.k:.1f},$' + \
                rf' $l = {o1.ang_mom:.1f},$' + \
                rf' $r_0 = {r_0:.1f},$' + \
                rf' $\dot r_0 = {r_dot_0:.2f},$' + \
                rf' $\phi_0 = {phi_0:.2f}$' + \
                '\n'     # \n means a new line (adds some space here)
fig_4.suptitle(overall_title, va='baseline')

ax_4a = fig_4.add_subplot(2,2,1)
ax_4a.plot(t_pts, r_pts, color='black')
ax_4a.set_xlabel(r'$t$')
ax_4a.set_ylabel(r'$r$')
ax_4a.set_title('Time dependence of radius')

ax_4b = fig_4.add_subplot(2,2,2)
ax_4b.plot(t_pts, phi_pts/(2.*np.pi), color='black')
ax_4b.plot(t_pts, phi_pts/(2.*np.pi)%1, color='red')
ax_4b.set_xlabel(r'$t$')
ax_4b.set_ylabel(r'$\phi/2\pi$')
ax_4b.set_title(r'Time dependence of $\phi$')

ax_4c = fig_4.add_subplot(2,2,3)
ax_4c.plot(r_pts*np.cos(phi_pts), r_pts*np.sin(phi_pts), color='black')
ax_4c.set_xlabel(r'$x$')
ax_4c.set_ylabel(r'$y$')
ax_4c.set_aspect(1)
ax_4c.set_title('Cartesian plot')

ax_4d = fig_4.add_subplot(2,2,4, polar=True)
ax_4d.plot(phi_pts, r_pts, color='black')
ax_4d.set_title('Polar plot', pad=20.)


fig_4.tight_layout()
../../_images/Orbit_games_10_0.png

Energy versus time#

E_tot_pts = o1.energy(t_pts, r_pts, r_dot_pts)
E_tot_0 = E_tot_pts[0]
E_tot_rel_pts = np.abs((E_tot_pts - E_tot_0)/E_tot_0)

print(f'    t        E_tot        rel. error')
for t, E_tot, E_tot_rel in zip(t_pts, E_tot_pts, E_tot_rel_pts):
    print(f'{t:8.5f}  {E_tot:8.5e}  {E_tot_rel:12.4e}')
    t        E_tot        rel. error
 0.00000  -2.50000e+00    0.0000e+00
 0.01000  -2.50000e+00    3.9790e-14
 0.02000  -2.50000e+00    5.4987e-11
 0.03000  -2.50000e+00    2.2135e-10
 0.04000  -2.50000e+00    2.8173e-10
 0.05000  -2.50000e+00    1.1803e-10
 0.06000  -2.50000e+00    1.1083e-10
 0.07000  -2.50000e+00    1.6178e-10
 0.08000  -2.50000e+00    5.9066e-12
 0.09000  -2.50000e+00    2.3637e-10
 0.10000  -2.50000e+00    2.2980e-10
 0.11000  -2.50000e+00    5.2911e-11
 0.12000  -2.50000e+00    3.2821e-10
 0.13000  -2.50000e+00    3.1704e-10
 0.14000  -2.50000e+00    5.3087e-11
 0.15000  -2.50000e+00    3.2847e-10
 0.16000  -2.50000e+00    1.5683e-10
 0.17000  -2.50000e+00    3.1309e-10
 0.18000  -2.50000e+00    6.0824e-10
 0.19000  -2.50000e+00    4.6531e-10
 0.20000  -2.50000e+00    2.7574e-10
 0.21000  -2.50000e+00    5.8478e-10
 0.22000  -2.50000e+00    6.9095e-11
 0.23000  -2.50000e+00    7.3493e-10
 0.24000  -2.50000e+00    1.0180e-09
 0.25000  -2.50000e+00    6.7889e-10
 0.26000  -2.50000e+00    7.4282e-10
 0.27000  -2.50000e+00    1.1964e-09
 0.28000  -2.50000e+00    1.0926e-10
 0.29000  -2.50000e+00    1.3503e-09
 0.30000  -2.50000e+00    1.7085e-09
 0.31000  -2.50000e+00    1.2202e-09
 0.32000  -2.50000e+00    1.1001e-09
 0.33000  -2.50000e+00    2.2014e-09
 0.34000  -2.50000e+00    7.6467e-10
 0.35000  -2.50000e+00    1.6446e-09
 0.36000  -2.50000e+00    2.6799e-09
 0.37000  -2.50000e+00    2.2847e-09
 0.38000  -2.50000e+00    1.1544e-09
 0.39000  -2.50000e+00    6.3315e-10
 0.40000  -2.50000e+00    1.7217e-09
 0.41000  -2.50000e+00    1.2599e-09
 0.42000  -2.50000e+00    9.7891e-10
 0.43000  -2.50000e+00    3.4940e-09
 0.44000  -2.50000e+00    6.4154e-09
 0.45000  -2.50000e+00    6.4049e-09
 0.46000  -2.50000e+00    5.9422e-11
 0.47000  -2.50000e+00    1.9396e-09
 0.48000  -2.50000e+00    3.9348e-09
 0.49000  -2.50000e+00    1.3715e-08
 0.50000  -2.50000e+00    7.6172e-09
 0.51000  -2.50000e+00    7.5057e-09
 0.52000  -2.50000e+00    1.5841e-09
 0.53000  -2.50000e+00    1.9661e-08
 0.54000  -2.50000e+00    1.0743e-08
 0.55000  -2.50000e+00    1.5008e-08
 0.56000  -2.50000e+00    4.3178e-09
 0.57000  -2.50000e+00    3.3344e-08
 0.58000  -2.50000e+00    1.6230e-08
 0.59000  -2.50000e+00    2.8139e-10
 0.60000  -2.50000e+00    4.6405e-08
 0.61000  -2.50000e+00    2.5538e-08
 0.62000  -2.50000e+00    3.5538e-09
 0.63000  -2.50000e+00    4.7292e-08
 0.64000  -2.50000e+00    4.2732e-08
 0.65000  -2.50000e+00    7.0327e-08
 0.66000  -2.50000e+00    5.2806e-08
 0.67000  -2.50000e+00    6.6673e-08
 0.68000  -2.50000e+00    2.7369e-08
 0.69000  -2.50000e+00    1.8511e-10
 0.70000  -2.50000e+00    6.3437e-08
 0.71000  -2.50000e+00    1.0228e-08
 0.72000  -2.50000e+00    5.7222e-08
 0.73000  -2.50000e+00    4.6488e-09
 0.74000  -2.50000e+00    1.0201e-07
 0.75000  -2.50000e+00    5.5432e-08
 0.76000  -2.50000e+00    4.1167e-08
 0.77000  -2.50000e+00    2.3906e-08
 0.78000  -2.50000e+00    2.6859e-08
 0.79000  -2.50000e+00    2.4018e-08
 0.80000  -2.50000e+00    1.6807e-08
 0.81000  -2.50000e+00    2.4489e-08
 0.82000  -2.50000e+00    3.1727e-08
 0.83000  -2.50000e+00    2.6533e-08
 0.84000  -2.50000e+00    2.5021e-08
 0.85000  -2.50000e+00    2.4827e-08
 0.86000  -2.50000e+00    2.4990e-08
 0.87000  -2.50000e+00    2.2530e-08
 0.88000  -2.50000e+00    1.8485e-08
 0.89000  -2.50000e+00    1.4529e-08
 0.90000  -2.50000e+00    1.4165e-08
 0.91000  -2.50000e+00    1.8196e-08
 0.92000  -2.50000e+00    9.4591e-09
 0.93000  -2.50000e+00    1.6221e-08
 0.94000  -2.50000e+00    2.8593e-08
 0.95000  -2.50000e+00    2.4612e-08
 0.96000  -2.50000e+00    8.1724e-09
 0.97000  -2.50000e+00    5.2331e-09
 0.98000  -2.50000e+00    1.5965e-07
 0.99000  -2.50000e+00    9.7842e-08
 1.00000  -2.50000e+00    1.5077e-07
 1.01000  -2.50000e+00    6.3371e-08
 1.02000  -2.50000e+00    9.2148e-08
 1.03000  -2.50000e+00    1.1760e-08
 1.04000  -2.50000e+00    3.8406e-08
 1.05000  -2.50000e+00    5.1112e-08
 1.06000  -2.50000e+00    7.8091e-09
 1.07000  -2.50000e+00    7.8441e-08
 1.08000  -2.50000e+00    1.0661e-08
 1.09000  -2.50000e+00    2.6395e-08
 1.10000  -2.50000e+00    6.9259e-08
 1.11000  -2.50000e+00    1.8172e-08
 1.12000  -2.50000e+00    2.3060e-08
 1.13000  -2.50000e+00    5.4972e-08
 1.14000  -2.50000e+00    4.3295e-08
 1.15000  -2.50000e+00    1.3891e-08
 1.16000  -2.50000e+00    2.9865e-08
 1.17000  -2.50000e+00    4.8707e-08
 1.18000  -2.50000e+00    3.9716e-08
 1.19000  -2.50000e+00    2.1374e-08
 1.20000  -2.50000e+00    2.8416e-08
 1.21000  -2.50000e+00    3.5244e-08
 1.22000  -2.50000e+00    3.8388e-08
 1.23000  -2.50000e+00    3.1432e-08
 1.24000  -2.50000e+00    2.6794e-08
 1.25000  -2.50000e+00    3.0628e-08
 1.26000  -2.50000e+00    3.1624e-08
 1.27000  -2.50000e+00    3.1109e-08
 1.28000  -2.50000e+00    3.0274e-08
 1.29000  -2.50000e+00    3.0224e-08
 1.30000  -2.50000e+00    3.1105e-08
 1.31000  -2.50000e+00    3.1598e-08
 1.32000  -2.50000e+00    2.9588e-08
 1.33000  -2.50000e+00    2.8757e-08
 1.34000  -2.50000e+00    2.9968e-08
 1.35000  -2.50000e+00    3.1428e-08
 1.36000  -2.50000e+00    3.1666e-08
 1.37000  -2.50000e+00    3.1202e-08
 1.38000  -2.50000e+00    2.9435e-08
 1.39000  -2.50000e+00    2.9077e-08
 1.40000  -2.50000e+00    3.0222e-08
 1.41000  -2.50000e+00    3.1305e-08
 1.42000  -2.50000e+00    3.1346e-08
 1.43000  -2.50000e+00    3.1013e-08
 1.44000  -2.50000e+00    3.0006e-08
 1.45000  -2.50000e+00    2.9737e-08
 1.46000  -2.50000e+00    3.0309e-08
 1.47000  -2.50000e+00    3.0871e-08
 1.48000  -2.50000e+00    3.0913e-08
 1.49000  -2.50000e+00    3.0759e-08
 1.50000  -2.50000e+00    3.0347e-08
 1.51000  -2.50000e+00    3.0184e-08
 1.52000  -2.50000e+00    3.0392e-08
 1.53000  -2.50000e+00    3.0570e-08
 1.54000  -2.50000e+00    3.0530e-08
 1.55000  -2.50000e+00    3.0443e-08
 1.56000  -2.50000e+00    3.0335e-08
 1.57000  -2.50000e+00    3.0314e-08
 1.58000  -2.50000e+00    3.0423e-08
 1.59000  -2.50000e+00    3.0425e-08
 1.60000  -2.50000e+00    3.0274e-08
 1.61000  -2.50000e+00    3.0160e-08
 1.62000  -2.50000e+00    3.0138e-08
 1.63000  -2.50000e+00    3.0212e-08
 1.64000  -2.50000e+00    3.0352e-08
 1.65000  -2.50000e+00    3.0330e-08
 1.66000  -2.50000e+00    3.0116e-08
 1.67000  -2.50000e+00    2.9924e-08
 1.68000  -2.50000e+00    2.9902e-08
 1.69000  -2.50000e+00    3.0020e-08
 1.70000  -2.50000e+00    3.0211e-08
 1.71000  -2.50000e+00    3.0217e-08
 1.72000  -2.50000e+00    2.9977e-08
 1.73000  -2.50000e+00    2.9718e-08
 1.74000  -2.50000e+00    2.9680e-08
 1.75000  -2.50000e+00    2.9859e-08
 1.76000  -2.50000e+00    3.0107e-08
 1.77000  -2.50000e+00    3.0113e-08
 1.78000  -2.50000e+00    2.9818e-08
 1.79000  -2.50000e+00    2.9498e-08
 1.80000  -2.50000e+00    2.9462e-08
 1.81000  -2.50000e+00    2.9805e-08
 1.82000  -2.50000e+00    3.0166e-08
 1.83000  -2.50000e+00    3.0116e-08
 1.84000  -2.50000e+00    2.9657e-08
 1.85000  -2.50000e+00    2.9217e-08
 1.86000  -2.50000e+00    2.9153e-08
 1.87000  -2.50000e+00    2.9797e-08
 1.88000  -2.50000e+00    3.0416e-08
 1.89000  -2.50000e+00    3.0252e-08
 1.90000  -2.50000e+00    2.9440e-08
 1.91000  -2.50000e+00    2.8764e-08
 1.92000  -2.50000e+00    2.8719e-08
 1.93000  -2.50000e+00    2.9682e-08
 1.94000  -2.50000e+00    3.0977e-08
 1.95000  -2.50000e+00    3.0806e-08
 1.96000  -2.50000e+00    2.9335e-08
 1.97000  -2.50000e+00    2.8027e-08
 1.98000  -2.50000e+00    2.8011e-08
 1.99000  -2.50000e+00    2.8931e-08
 2.00000  -2.50000e+00    3.1367e-08
 2.01000  -2.50000e+00    3.1871e-08
 2.02000  -2.50000e+00    2.9921e-08
 2.03000  -2.50000e+00    2.7504e-08
 2.04000  -2.50000e+00    2.6852e-08
 2.05000  -2.50000e+00    2.7363e-08
 2.06000  -2.50000e+00    2.9128e-08
 2.07000  -2.50000e+00    3.0888e-08
 2.08000  -2.50000e+00    3.1161e-08
 2.09000  -2.50000e+00    2.9821e-08
 2.10000  -2.50000e+00    2.7572e-08
 2.11000  -2.50000e+00    2.5998e-08
 2.12000  -2.50000e+00    1.9258e-08
 2.13000  -2.50000e+00    1.6918e-08
 2.14000  -2.50000e+00    2.8844e-08
 2.15000  -2.50000e+00    4.2448e-08
 2.16000  -2.50000e+00    3.9541e-08
 2.17000  -2.50000e+00    2.4564e-08
 2.18000  -2.50000e+00    6.7603e-09
 2.19000  -2.50000e+00    1.0100e-08
 2.20000  -2.50000e+00    4.1866e-08
 2.21000  -2.50000e+00    4.2528e-08
 2.22000  -2.50000e+00    1.9076e-08
 2.23000  -2.50000e+00    1.1391e-08
 2.24000  -2.50000e+00    3.7500e-08
 2.25000  -2.50000e+00    5.2163e-08
 2.26000  -2.50000e+00    8.3625e-09
 2.27000  -2.50000e+00    1.1108e-08
 2.28000  -2.50000e+00    6.8170e-08
 2.29000  -2.50000e+00    2.3860e-08
 2.30000  -2.50000e+00    2.9712e-08
 2.31000  -2.50000e+00    7.9435e-08
 2.32000  -2.50000e+00    4.4674e-09
 2.33000  -2.50000e+00    8.6412e-09
 2.34000  -2.50000e+00    5.1929e-08
 2.35000  -2.50000e+00    6.6709e-08
 2.36000  -2.50000e+00    1.0055e-07
 2.37000  -2.50000e+00    1.9473e-07
 2.38000  -2.50000e+00    2.1942e-07
 2.39000  -2.50000e+00    2.2803e-07
 2.40000  -2.50000e+00    2.8652e-07
 2.41000  -2.50000e+00    2.1659e-07
 2.42000  -2.50000e+00    1.8433e-07
 2.43000  -2.50000e+00    8.7791e-08
 2.44000  -2.50000e+00    5.3790e-08
 2.45000  -2.50000e+00    6.2019e-08
 2.46000  -2.50000e+00    5.0856e-08
 2.47000  -2.50000e+00    5.8734e-08
 2.48000  -2.50000e+00    5.7665e-08
 2.49000  -2.50000e+00    5.5052e-08
 2.50000  -2.50000e+00    6.4470e-08
 2.51000  -2.50000e+00    5.3802e-08
 2.52000  -2.50000e+00    5.4336e-08
 2.53000  -2.50000e+00    5.3307e-08
 2.54000  -2.50000e+00    5.2329e-08
 2.55000  -2.50000e+00    4.9606e-08
 2.56000  -2.50000e+00    5.0846e-08
 2.57000  -2.50000e+00    4.6301e-08
 2.58000  -2.50000e+00    5.1003e-08
 2.59000  -2.50000e+00    4.4954e-08
 2.60000  -2.50000e+00    4.8803e-08
 2.61000  -2.50000e+00    3.7095e-08
 2.62000  -2.50000e+00    3.2212e-08
 2.63000  -2.50000e+00    3.1163e-08
 2.64000  -2.50000e+00    5.2749e-09
 2.65000  -2.50000e+00    4.4275e-08
 2.66000  -2.50000e+00    1.7910e-07
 2.67000  -2.50000e+00    5.6989e-08
 2.68000  -2.50000e+00    2.1647e-07
 2.69000  -2.50000e+00    4.6915e-08
 2.70000  -2.50000e+00    1.5710e-07
 2.71000  -2.50000e+00    1.2379e-08
 2.72000  -2.50000e+00    1.0907e-07
 2.73000  -2.50000e+00    3.5037e-08
 2.74000  -2.50000e+00    5.4542e-08
 2.75000  -2.50000e+00    1.0177e-07
 2.76000  -2.50000e+00    2.3182e-08
 2.77000  -2.50000e+00    7.8895e-08
 2.78000  -2.50000e+00    8.1401e-08
 2.79000  -2.50000e+00    3.2958e-08
 2.80000  -2.50000e+00    6.8648e-08
 2.81000  -2.50000e+00    8.9304e-08
 2.82000  -2.50000e+00    4.4057e-08
 2.83000  -2.50000e+00    5.4867e-08
 2.84000  -2.50000e+00    7.8225e-08
 2.85000  -2.50000e+00    7.6573e-08
 2.86000  -2.50000e+00    4.9471e-08
 2.87000  -2.50000e+00    5.4450e-08
 2.88000  -2.50000e+00    6.5270e-08
 2.89000  -2.50000e+00    7.1986e-08
 2.90000  -2.50000e+00    6.1695e-08
 2.91000  -2.50000e+00    5.3705e-08
 2.92000  -2.50000e+00    5.9251e-08
 2.93000  -2.50000e+00    6.1535e-08
 2.94000  -2.50000e+00    6.1963e-08
 2.95000  -2.50000e+00    5.9979e-08
 2.96000  -2.50000e+00    5.8666e-08
 2.97000  -2.50000e+00    5.9844e-08
 2.98000  -2.50000e+00    6.0678e-08
 2.99000  -2.50000e+00    5.8797e-08
 3.00000  -2.50000e+00    5.8188e-08
 3.01000  -2.50000e+00    5.9335e-08
 3.02000  -2.50000e+00    6.0607e-08
 3.03000  -2.50000e+00    6.0797e-08
 3.04000  -2.50000e+00    6.0092e-08
 3.05000  -2.50000e+00    5.8217e-08
 3.06000  -2.50000e+00    5.8232e-08
 3.07000  -2.50000e+00    5.9659e-08
 3.08000  -2.50000e+00    6.0689e-08
 3.09000  -2.50000e+00    6.0535e-08
 3.10000  -2.50000e+00    5.9999e-08
 3.11000  -2.50000e+00    5.8824e-08
 3.12000  -2.50000e+00    5.8810e-08
 3.13000  -2.50000e+00    5.9632e-08
 3.14000  -2.50000e+00    6.0227e-08
 3.15000  -2.50000e+00    6.0155e-08
 3.16000  -2.50000e+00    5.9930e-08
 3.17000  -2.50000e+00    5.9353e-08
 3.18000  -2.50000e+00    5.9278e-08
 3.19000  -2.50000e+00    5.9602e-08
 3.20000  -2.50000e+00    5.9824e-08
 3.21000  -2.50000e+00    5.9771e-08
 3.22000  -2.50000e+00    5.9690e-08
 3.23000  -2.50000e+00    5.9493e-08
 3.24000  -2.50000e+00    5.9485e-08
 3.25000  -2.50000e+00    5.9607e-08
 3.26000  -2.50000e+00    5.9605e-08
 3.27000  -2.50000e+00    5.9470e-08
 3.28000  -2.50000e+00    5.9398e-08
 3.29000  -2.50000e+00    5.9363e-08
 3.30000  -2.50000e+00    5.9446e-08
 3.31000  -2.50000e+00    5.9557e-08
 3.32000  -2.50000e+00    5.9491e-08
 3.33000  -2.50000e+00    5.9275e-08
 3.34000  -2.50000e+00    5.9138e-08
 3.35000  -2.50000e+00    5.9138e-08
 3.36000  -2.50000e+00    5.9276e-08
 3.37000  -2.50000e+00    5.9436e-08
 3.38000  -2.50000e+00    5.9385e-08
 3.39000  -2.50000e+00    5.9127e-08
 3.40000  -2.50000e+00    5.8915e-08
 3.41000  -2.50000e+00    5.8914e-08
 3.42000  -2.50000e+00    5.9098e-08
 3.43000  -2.50000e+00    5.9298e-08
 3.44000  -2.50000e+00    5.9247e-08
 3.45000  -2.50000e+00    5.8950e-08
 3.46000  -2.50000e+00    5.8698e-08
 3.47000  -2.50000e+00    5.8732e-08
 3.48000  -2.50000e+00    5.9059e-08
 3.49000  -2.50000e+00    5.9332e-08
 3.50000  -2.50000e+00    5.9215e-08
 3.51000  -2.50000e+00    5.8777e-08
 3.52000  -2.50000e+00    5.8432e-08
 3.53000  -2.50000e+00    5.8489e-08
 3.54000  -2.50000e+00    5.9097e-08
 3.55000  -2.50000e+00    5.9512e-08
 3.56000  -2.50000e+00    5.9242e-08
 3.57000  -2.50000e+00    5.8522e-08
 3.58000  -2.50000e+00    5.8038e-08
 3.59000  -2.50000e+00    5.8071e-08
 3.60000  -2.50000e+00    5.9177e-08
 3.61000  -2.50000e+00    6.0024e-08
 3.62000  -2.50000e+00    5.9550e-08
 3.63000  -2.50000e+00    5.8224e-08
 3.64000  -2.50000e+00    5.7348e-08
 3.65000  -2.50000e+00    5.7481e-08
 3.66000  -2.50000e+00    5.8866e-08
 3.67000  -2.50000e+00    6.0762e-08
 3.68000  -2.50000e+00    6.0474e-08
 3.69000  -2.50000e+00    5.8331e-08
 3.70000  -2.50000e+00    5.6469e-08
 3.71000  -2.50000e+00    5.6480e-08
 3.72000  -2.50000e+00    5.7371e-08
 3.73000  -2.50000e+00    5.9931e-08
 3.74000  -2.50000e+00    6.0976e-08
 3.75000  -2.50000e+00    5.9631e-08
 3.76000  -2.50000e+00    5.7265e-08
 3.77000  -2.50000e+00    5.5761e-08
 3.78000  -2.50000e+00    5.5360e-08
 3.79000  -2.50000e+00    5.2059e-08
 3.80000  -2.50000e+00    5.3567e-08
 3.81000  -2.50000e+00    6.0945e-08
 3.82000  -2.50000e+00    6.6258e-08
 3.83000  -2.50000e+00    6.1504e-08
 3.84000  -2.50000e+00    5.3859e-08
 3.85000  -2.50000e+00    3.8428e-08
 3.86000  -2.50000e+00    4.3543e-08
 3.87000  -2.50000e+00    6.8924e-08
 3.88000  -2.50000e+00    7.1845e-08
 3.89000  -2.50000e+00    5.3010e-08
 3.90000  -2.50000e+00    2.1556e-08
 3.91000  -2.50000e+00    4.9470e-08
 3.92000  -2.50000e+00    8.7058e-08
 3.93000  -2.50000e+00    5.4437e-08
 3.94000  -2.50000e+00    8.8719e-09
 3.95000  -2.50000e+00    6.1144e-08
 3.96000  -2.50000e+00    8.6738e-08
 3.97000  -2.50000e+00    2.2088e-08
 3.98000  -2.50000e+00    3.5921e-08
 3.99000  -2.50000e+00    9.8443e-08
 4.00000  -2.50000e+00    1.6172e-09
 4.01000  -2.50000e+00    7.3203e-08
 4.02000  -2.50000e+00    6.1181e-08
 4.03000  -2.50000e+00    1.6670e-08
 4.04000  -2.50000e+00    1.1506e-07
 4.05000  -2.50000e+00    1.5661e-07
 4.06000  -2.50000e+00    2.3915e-07
 4.07000  -2.50000e+00    1.5732e-07
 4.08000  -2.50000e+00    2.7625e-07
 4.09000  -2.50000e+00    2.4448e-07
 4.10000  -2.50000e+00    4.0674e-07
 4.11000  -2.50000e+00    1.0520e-08
 4.12000  -2.50000e+00    4.2229e-08
 4.13000  -2.50000e+00    6.3992e-08
 4.14000  -2.50000e+00    7.5868e-08
 4.15000  -2.50000e+00    8.5191e-08
 4.16000  -2.50000e+00    1.0191e-07
 4.17000  -2.50000e+00    8.7854e-08
 4.18000  -2.50000e+00    9.0760e-08
 4.19000  -2.50000e+00    9.3205e-08
 4.20000  -2.50000e+00    8.9486e-08
 4.21000  -2.50000e+00    8.7530e-08
 4.22000  -2.50000e+00    8.7159e-08
 4.23000  -2.50000e+00    8.7231e-08
 4.24000  -2.50000e+00    8.6003e-08
 4.25000  -2.50000e+00    7.8471e-08
 4.26000  -2.50000e+00    8.2938e-08
 4.27000  -2.50000e+00    8.2454e-08
 4.28000  -2.50000e+00    7.7062e-08
 4.29000  -2.50000e+00    7.9631e-08
 4.30000  -2.50000e+00    9.3738e-08
 4.31000  -2.50000e+00    9.0892e-08
 4.32000  -2.50000e+00    1.0634e-07
 4.33000  -2.50000e+00    1.5764e-07
 4.34000  -2.50000e+00    7.7304e-08
 4.35000  -2.50000e+00    2.8403e-08
 4.36000  -2.50000e+00    2.0395e-07
 4.37000  -2.50000e+00    1.2710e-08
 4.38000  -2.50000e+00    1.5767e-07
 4.39000  -2.50000e+00    4.1623e-08
 4.40000  -2.50000e+00    1.0782e-07
 4.41000  -2.50000e+00    1.0091e-07
 4.42000  -2.50000e+00    8.0572e-08
 4.43000  -2.50000e+00    1.4757e-07
 4.44000  -2.50000e+00    5.4399e-08
 4.45000  -2.50000e+00    1.0339e-07
 4.46000  -2.50000e+00    1.2279e-07
 4.47000  -2.50000e+00    6.3985e-08
 4.48000  -2.50000e+00    1.0017e-07
 4.49000  -2.50000e+00    1.2513e-07
 4.50000  -2.50000e+00    7.6355e-08
 4.51000  -2.50000e+00    8.8859e-08
 4.52000  -2.50000e+00    1.1490e-07
 4.53000  -2.50000e+00    1.0932e-07
 4.54000  -2.50000e+00    8.0455e-08
 4.55000  -2.50000e+00    8.9902e-08
 4.56000  -2.50000e+00    1.0220e-07
 4.57000  -2.50000e+00    1.0692e-07
 4.58000  -2.50000e+00    9.2982e-08
 4.59000  -2.50000e+00    8.6409e-08
 4.60000  -2.50000e+00    9.3949e-08
 4.61000  -2.50000e+00    9.6351e-08
 4.62000  -2.50000e+00    9.6384e-08
 4.63000  -2.50000e+00    9.3235e-08
 4.64000  -2.50000e+00    9.1823e-08
 4.65000  -2.50000e+00    9.3836e-08
 4.66000  -2.50000e+00    9.3956e-08
 4.67000  -2.50000e+00    9.2318e-08
 4.68000  -2.50000e+00    9.2200e-08
 4.69000  -2.50000e+00    9.3378e-08
 4.70000  -2.50000e+00    9.4353e-08
 4.71000  -2.50000e+00    9.4412e-08
 4.72000  -2.50000e+00    9.3044e-08
 4.73000  -2.50000e+00    9.1569e-08
 4.74000  -2.50000e+00    9.2254e-08
 4.75000  -2.50000e+00    9.3808e-08
 4.76000  -2.50000e+00    9.4466e-08
 4.77000  -2.50000e+00    9.4117e-08
 4.78000  -2.50000e+00    9.3149e-08
 4.79000  -2.50000e+00    9.2212e-08
 4.80000  -2.50000e+00    9.2657e-08
 4.81000  -2.50000e+00    9.3609e-08
 4.82000  -2.50000e+00    9.4007e-08
 4.83000  -2.50000e+00    9.3781e-08
 4.84000  -2.50000e+00    9.3372e-08
 4.85000  -2.50000e+00    9.2847e-08
 4.86000  -2.50000e+00    9.2990e-08
 4.87000  -2.50000e+00    9.3392e-08
 4.88000  -2.50000e+00    9.3549e-08
 4.89000  -2.50000e+00    9.3442e-08
 4.90000  -2.50000e+00    9.3320e-08
 4.91000  -2.50000e+00    9.3112e-08
 4.92000  -2.50000e+00    9.3176e-08
 4.93000  -2.50000e+00    9.3297e-08
 4.94000  -2.50000e+00    9.3257e-08
 4.95000  -2.50000e+00    9.3131e-08
 4.96000  -2.50000e+00    9.3092e-08
 4.97000  -2.50000e+00    9.3059e-08
 4.98000  -2.50000e+00    9.3170e-08
 4.99000  -2.50000e+00    9.3233e-08
 5.00000  -2.50000e+00    9.3103e-08
 5.01000  -2.50000e+00    9.2898e-08
 5.02000  -2.50000e+00    9.2832e-08
 5.03000  -2.50000e+00    9.2862e-08
 5.04000  -2.50000e+00    9.3029e-08
 5.05000  -2.50000e+00    9.3129e-08
 5.06000  -2.50000e+00    9.2990e-08
 5.07000  -2.50000e+00    9.2721e-08
 5.08000  -2.50000e+00    9.2591e-08
 5.09000  -2.50000e+00    9.2647e-08
 5.10000  -2.50000e+00    9.2861e-08
 5.11000  -2.50000e+00    9.2995e-08
 5.12000  -2.50000e+00    9.2848e-08
 5.13000  -2.50000e+00    9.2534e-08
 5.14000  -2.50000e+00    9.2369e-08
 5.15000  -2.50000e+00    9.2496e-08
 5.16000  -2.50000e+00    9.2830e-08
 5.17000  -2.50000e+00    9.2992e-08
 5.18000  -2.50000e+00    9.2760e-08
 5.19000  -2.50000e+00    9.2327e-08
 5.20000  -2.50000e+00    9.2102e-08
 5.21000  -2.50000e+00    9.2343e-08
 5.22000  -2.50000e+00    9.2921e-08
 5.23000  -2.50000e+00    9.3120e-08
 5.24000  -2.50000e+00    9.2689e-08
 5.25000  -2.50000e+00    9.2024e-08
 5.26000  -2.50000e+00    9.1737e-08
 5.27000  -2.50000e+00    9.2060e-08
 5.28000  -2.50000e+00    9.3178e-08
 5.29000  -2.50000e+00    9.3566e-08
 5.30000  -2.50000e+00    9.2765e-08
 5.31000  -2.50000e+00    9.1562e-08
 5.32000  -2.50000e+00    9.1107e-08
 5.33000  -2.50000e+00    9.1393e-08
 5.34000  -2.50000e+00    9.3328e-08
 5.35000  -2.50000e+00    9.4466e-08
 5.36000  -2.50000e+00    9.3390e-08
 5.37000  -2.50000e+00    9.1218e-08
 5.38000  -2.50000e+00    9.0076e-08
 5.39000  -2.50000e+00    9.0499e-08
 5.40000  -2.50000e+00    9.2218e-08
 5.41000  -2.50000e+00    9.4541e-08
 5.42000  -2.50000e+00    9.4416e-08
 5.43000  -2.50000e+00    9.2165e-08
 5.44000  -2.50000e+00    8.9924e-08
 5.45000  -2.50000e+00    8.9358e-08
 5.46000  -2.50000e+00    8.8716e-08
 5.47000  -2.50000e+00    8.7797e-08
 5.48000  -2.50000e+00    9.0964e-08
 5.49000  -2.50000e+00    9.6036e-08
 5.50000  -2.50000e+00    9.7003e-08
 5.51000  -2.50000e+00    9.0949e-08
 5.52000  -2.50000e+00    8.4711e-08
 5.53000  -2.50000e+00    7.2830e-08
 5.54000  -2.50000e+00    8.4680e-08
 5.55000  -2.50000e+00    1.0515e-07
 5.56000  -2.50000e+00    1.0019e-07
 5.57000  -2.50000e+00    8.3412e-08
 5.58000  -2.50000e+00    5.6544e-08
 5.59000  -2.50000e+00    8.9380e-08
 5.60000  -2.50000e+00    1.1918e-07
 5.61000  -2.50000e+00    8.7445e-08
 5.62000  -2.50000e+00    4.4202e-08
 5.63000  -2.50000e+00    9.3330e-08
 5.64000  -2.50000e+00    1.2191e-07
 5.65000  -2.50000e+00    6.4392e-08
 5.66000  -2.50000e+00    5.8550e-08
 5.67000  -2.50000e+00    1.3726e-07
 5.68000  -2.50000e+00    5.8069e-08
 5.69000  -2.50000e+00    7.5382e-08
 5.70000  -2.50000e+00    1.2039e-07
 5.71000  -2.50000e+00    3.6755e-10
 5.72000  -2.50000e+00    1.6513e-07
 5.73000  -2.50000e+00    1.4659e-09
 5.74000  -2.50000e+00    1.7512e-07
 5.75000  -2.50000e+00    6.3895e-08
 5.76000  -2.50000e+00    2.4858e-07
 5.77000  -2.50000e+00    8.9148e-08
 5.78000  -2.50000e+00    3.4532e-07
 5.79000  -2.50000e+00    2.2874e-08
 5.80000  -2.50000e+00    1.1491e-07
 5.81000  -2.50000e+00    1.3583e-07
 5.82000  -2.50000e+00    1.1605e-07
 5.83000  -2.50000e+00    1.1138e-07
 5.84000  -2.50000e+00    1.1303e-07
 5.85000  -2.50000e+00    1.3341e-07
 5.86000  -2.50000e+00    1.2960e-07
 5.87000  -2.50000e+00    1.2762e-07
 5.88000  -2.50000e+00    1.1719e-07
 5.89000  -2.50000e+00    1.1728e-07
 5.90000  -2.50000e+00    1.1649e-07
 5.91000  -2.50000e+00    1.1549e-07
 5.92000  -2.50000e+00    1.1317e-07
 5.93000  -2.50000e+00    1.1349e-07
 5.94000  -2.50000e+00    1.1427e-07
 5.95000  -2.50000e+00    1.1540e-07
 5.96000  -2.50000e+00    1.1241e-07
 5.97000  -2.50000e+00    1.0673e-07
 5.98000  -2.50000e+00    1.0697e-07
 5.99000  -2.50000e+00    1.2400e-07
 6.00000  -2.50000e+00    1.4378e-07
 6.01000  -2.50000e+00    1.5294e-07
 6.02000  -2.50000e+00    1.0446e-07
 6.03000  -2.50000e+00    5.0984e-08
 6.04000  -2.50000e+00    2.8048e-07
 6.05000  -2.50000e+00    8.1322e-09
 6.06000  -2.50000e+00    2.4580e-07
 6.07000  -2.50000e+00    2.9638e-08
 6.08000  -2.50000e+00    1.8993e-07
 6.09000  -2.50000e+00    7.6227e-08
 6.10000  -2.50000e+00    1.3326e-07
 6.11000  -2.50000e+00    1.4204e-07
 6.12000  -2.50000e+00    9.9084e-08
 6.13000  -2.50000e+00    1.7043e-07
 6.14000  -2.50000e+00    1.0587e-07
 6.15000  -2.50000e+00    1.1656e-07
 6.16000  -2.50000e+00    1.6205e-07
 6.17000  -2.50000e+00    1.1655e-07
 6.18000  -2.50000e+00    1.0937e-07
 6.19000  -2.50000e+00    1.4201e-07
 6.20000  -2.50000e+00    1.4410e-07
 6.21000  -2.50000e+00    1.0769e-07
 6.22000  -2.50000e+00    1.1911e-07
 6.23000  -2.50000e+00    1.3636e-07
 6.24000  -2.50000e+00    1.3902e-07
 6.25000  -2.50000e+00    1.1834e-07
 6.26000  -2.50000e+00    1.1499e-07
 6.27000  -2.50000e+00    1.2379e-07
 6.28000  -2.50000e+00    1.2971e-07
 6.29000  -2.50000e+00    1.2720e-07
 6.30000  -2.50000e+00    1.2069e-07
 6.31000  -2.50000e+00    1.2030e-07
 6.32000  -2.50000e+00    1.2384e-07
 6.33000  -2.50000e+00    1.2325e-07
 6.34000  -2.50000e+00    1.2238e-07
 6.35000  -2.50000e+00    1.2225e-07
 6.36000  -2.50000e+00    1.2294e-07
 6.37000  -2.50000e+00    1.2372e-07
 6.38000  -2.50000e+00    1.2336e-07
 6.39000  -2.50000e+00    1.2127e-07
 6.40000  -2.50000e+00    1.2113e-07
 6.41000  -2.50000e+00    1.2266e-07
 6.42000  -2.50000e+00    1.2386e-07
 6.43000  -2.50000e+00    1.2376e-07
 6.44000  -2.50000e+00    1.2302e-07
 6.45000  -2.50000e+00    1.2150e-07
 6.46000  -2.50000e+00    1.2159e-07
 6.47000  -2.50000e+00    1.2276e-07
 6.48000  -2.50000e+00    1.2356e-07
 6.49000  -2.50000e+00    1.2342e-07
 6.50000  -2.50000e+00    1.2303e-07
 6.51000  -2.50000e+00    1.2216e-07
 6.52000  -2.50000e+00    1.2213e-07
 6.53000  -2.50000e+00    1.2269e-07
 6.54000  -2.50000e+00    1.2309e-07
 6.55000  -2.50000e+00    1.2303e-07
 6.56000  -2.50000e+00    1.2289e-07
 6.57000  -2.50000e+00    1.2252e-07
 6.58000  -2.50000e+00    1.2248e-07
 6.59000  -2.50000e+00    1.2268e-07
 6.60000  -2.50000e+00    1.2277e-07
 6.61000  -2.50000e+00    1.2268e-07
 6.62000  -2.50000e+00    1.2262e-07
 6.63000  -2.50000e+00    1.2252e-07
 6.64000  -2.50000e+00    1.2256e-07
 6.65000  -2.50000e+00    1.2266e-07
 6.66000  -2.50000e+00    1.2261e-07
 6.67000  -2.50000e+00    1.2243e-07
 6.68000  -2.50000e+00    1.2234e-07
 6.69000  -2.50000e+00    1.2233e-07
 6.70000  -2.50000e+00    1.2244e-07
 6.71000  -2.50000e+00    1.2257e-07
 6.72000  -2.50000e+00    1.2250e-07
 6.73000  -2.50000e+00    1.2226e-07
 6.74000  -2.50000e+00    1.2209e-07
 6.75000  -2.50000e+00    1.2210e-07
 6.76000  -2.50000e+00    1.2226e-07
 6.77000  -2.50000e+00    1.2244e-07
 6.78000  -2.50000e+00    1.2238e-07
 6.79000  -2.50000e+00    1.2211e-07
 6.80000  -2.50000e+00    1.2188e-07
 6.81000  -2.50000e+00    1.2189e-07
 6.82000  -2.50000e+00    1.2213e-07
 6.83000  -2.50000e+00    1.2236e-07
 6.84000  -2.50000e+00    1.2229e-07
 6.85000  -2.50000e+00    1.2194e-07
 6.86000  -2.50000e+00    1.2164e-07
 6.87000  -2.50000e+00    1.2169e-07
 6.88000  -2.50000e+00    1.2212e-07
 6.89000  -2.50000e+00    1.2245e-07
 6.90000  -2.50000e+00    1.2227e-07
 6.91000  -2.50000e+00    1.2173e-07
 6.92000  -2.50000e+00    1.2133e-07
 6.93000  -2.50000e+00    1.2137e-07
 6.94000  -2.50000e+00    1.2219e-07
 6.95000  -2.50000e+00    1.2276e-07
 6.96000  -2.50000e+00    1.2240e-07
 6.97000  -2.50000e+00    1.2144e-07
 6.98000  -2.50000e+00    1.2082e-07
 6.99000  -2.50000e+00    1.2086e-07
 7.00000  -2.50000e+00    1.2217e-07
 7.01000  -2.50000e+00    1.2344e-07
 7.02000  -2.50000e+00    1.2298e-07
 7.03000  -2.50000e+00    1.2125e-07
 7.04000  -2.50000e+00    1.1997e-07
 7.05000  -2.50000e+00    1.2014e-07
 7.06000  -2.50000e+00    1.2136e-07
 7.07000  -2.50000e+00    1.2384e-07
 7.08000  -2.50000e+00    1.2406e-07
 7.09000  -2.50000e+00    1.2187e-07
 7.10000  -2.50000e+00    1.1947e-07
 7.11000  -2.50000e+00    1.1896e-07
 7.12000  -2.50000e+00    1.1940e-07
 7.13000  -2.50000e+00    1.2064e-07
 7.14000  -2.50000e+00    1.2239e-07
 7.15000  -2.50000e+00    1.2351e-07
 7.16000  -2.50000e+00    1.2285e-07
 7.17000  -2.50000e+00    1.2009e-07
 7.18000  -2.50000e+00    1.1780e-07
 7.19000  -2.50000e+00    1.0612e-07
 7.20000  -2.50000e+00    1.0526e-07
 7.21000  -2.50000e+00    1.2398e-07
 7.22000  -2.50000e+00    1.4105e-07
 7.23000  -2.50000e+00    1.3263e-07
 7.24000  -2.50000e+00    1.1556e-07
 7.25000  -2.50000e+00    9.1176e-08
 7.26000  -2.50000e+00    1.0774e-07
 7.27000  -2.50000e+00    1.4297e-07
 7.28000  -2.50000e+00    1.2694e-07
 7.29000  -2.50000e+00    9.2943e-08
 7.30000  -2.50000e+00    8.8830e-08
 7.31000  -2.50000e+00    1.5125e-07
 7.32000  -2.50000e+00    1.2321e-07
 7.33000  -2.50000e+00    6.6419e-08
 7.34000  -2.50000e+00    1.2701e-07
 7.35000  -2.50000e+00    1.4498e-07
 7.36000  -2.50000e+00    5.9838e-08
 7.37000  -2.50000e+00    1.3789e-07
 7.38000  -2.50000e+00    1.2705e-07
 7.39000  -2.50000e+00    4.0647e-08
 7.40000  -2.50000e+00    1.8728e-07
 7.41000  -2.50000e+00    2.6776e-08
 7.42000  -2.50000e+00    2.0024e-07
 7.43000  -2.50000e+00    2.1272e-08
 7.44000  -2.50000e+00    2.4798e-07
 7.45000  -2.50000e+00    1.6370e-08
 7.46000  -2.50000e+00    2.8476e-07
 7.47000  -2.50000e+00    6.5484e-09
 7.48000  -2.50000e+00    2.7239e-07
 7.49000  -2.50000e+00    1.6836e-07
 7.50000  -2.50000e+00    1.4812e-07
 7.51000  -2.50000e+00    1.5368e-07
 7.52000  -2.50000e+00    1.4240e-07
 7.53000  -2.50000e+00    1.4178e-07
 7.54000  -2.50000e+00    1.6274e-07
 7.55000  -2.50000e+00    1.6037e-07
 7.56000  -2.50000e+00    1.4526e-07
 7.57000  -2.50000e+00    1.4874e-07
 7.58000  -2.50000e+00    1.4564e-07
 7.59000  -2.50000e+00    1.4500e-07
 7.60000  -2.50000e+00    1.4499e-07
 7.61000  -2.50000e+00    1.4545e-07
 7.62000  -2.50000e+00    1.3646e-07
 7.63000  -2.50000e+00    1.4560e-07
 7.64000  -2.50000e+00    1.4146e-07
 7.65000  -2.50000e+00    1.3857e-07
 7.66000  -2.50000e+00    1.4004e-07
 7.67000  -2.50000e+00    1.2486e-07
 7.68000  -2.50000e+00    1.2677e-07
 7.69000  -2.50000e+00    1.2313e-07
 7.70000  -2.50000e+00    9.1658e-08
 7.71000  -2.50000e+00    1.6326e-07
 7.72000  -2.50000e+00    2.3484e-07
 7.73000  -2.50000e+00    3.1348e-08
 7.74000  -2.50000e+00    2.7408e-07
 7.75000  -2.50000e+00    5.0500e-08
 7.76000  -2.50000e+00    2.1463e-07
 7.77000  -2.50000e+00    1.0502e-07
 7.78000  -2.50000e+00    1.6096e-07
 7.79000  -2.50000e+00    1.6837e-07
 7.80000  -2.50000e+00    1.3046e-07
 7.81000  -2.50000e+00    2.0332e-07
 7.82000  -2.50000e+00    1.2417e-07
 7.83000  -2.50000e+00    1.4711e-07
 7.84000  -2.50000e+00    1.9205e-07
 7.85000  -2.50000e+00    1.3156e-07
 7.86000  -2.50000e+00    1.4685e-07
 7.87000  -2.50000e+00    1.8211e-07
 7.88000  -2.50000e+00    1.5862e-07
 7.89000  -2.50000e+00    1.3283e-07
 7.90000  -2.50000e+00    1.5210e-07
 7.91000  -2.50000e+00    1.7387e-07
 7.92000  -2.50000e+00    1.5787e-07
 7.93000  -2.50000e+00    1.3936e-07
 7.94000  -2.50000e+00    1.5031e-07
 7.95000  -2.50000e+00    1.5824e-07
 7.96000  -2.50000e+00    1.6023e-07
 7.97000  -2.50000e+00    1.5100e-07
 7.98000  -2.50000e+00    1.4646e-07
 7.99000  -2.50000e+00    1.5159e-07
 8.00000  -2.50000e+00    1.5248e-07
 8.01000  -2.50000e+00    1.5207e-07
 8.02000  -2.50000e+00    1.5093e-07
 8.03000  -2.50000e+00    1.5068e-07
 8.04000  -2.50000e+00    1.5171e-07
 8.05000  -2.50000e+00    1.5212e-07
 8.06000  -2.50000e+00    1.5007e-07
 8.07000  -2.50000e+00    1.4948e-07
 8.08000  -2.50000e+00    1.5078e-07
 8.09000  -2.50000e+00    1.5214e-07
 8.10000  -2.50000e+00    1.5227e-07
 8.11000  -2.50000e+00    1.5165e-07
 8.12000  -2.50000e+00    1.4987e-07
 8.13000  -2.50000e+00    1.4975e-07
 8.14000  -2.50000e+00    1.5103e-07
 8.15000  -2.50000e+00    1.5204e-07
 8.16000  -2.50000e+00    1.5196e-07
 8.17000  -2.50000e+00    1.5155e-07
 8.18000  -2.50000e+00    1.5049e-07
 8.19000  -2.50000e+00    1.5036e-07
 8.20000  -2.50000e+00    1.5103e-07
 8.21000  -2.50000e+00    1.5158e-07
 8.22000  -2.50000e+00    1.5156e-07
 8.23000  -2.50000e+00    1.5139e-07
 8.24000  -2.50000e+00    1.5091e-07
 8.25000  -2.50000e+00    1.5080e-07
 8.26000  -2.50000e+00    1.5105e-07
 8.27000  -2.50000e+00    1.5123e-07
 8.28000  -2.50000e+00    1.5118e-07
 8.29000  -2.50000e+00    1.5110e-07
 8.30000  -2.50000e+00    1.5096e-07
 8.31000  -2.50000e+00    1.5096e-07
 8.32000  -2.50000e+00    1.5107e-07
 8.33000  -2.50000e+00    1.5106e-07
 8.34000  -2.50000e+00    1.5090e-07
 8.35000  -2.50000e+00    1.5081e-07
 8.36000  -2.50000e+00    1.5079e-07
 8.37000  -2.50000e+00    1.5087e-07
 8.38000  -2.50000e+00    1.5100e-07
 8.39000  -2.50000e+00    1.5096e-07
 8.40000  -2.50000e+00    1.5073e-07
 8.41000  -2.50000e+00    1.5056e-07
 8.42000  -2.50000e+00    1.5055e-07
 8.43000  -2.50000e+00    1.5069e-07
 8.44000  -2.50000e+00    1.5087e-07
 8.45000  -2.50000e+00    1.5084e-07
 8.46000  -2.50000e+00    1.5059e-07
 8.47000  -2.50000e+00    1.5035e-07
 8.48000  -2.50000e+00    1.5033e-07
 8.49000  -2.50000e+00    1.5052e-07
 8.50000  -2.50000e+00    1.5075e-07
 8.51000  -2.50000e+00    1.5073e-07
 8.52000  -2.50000e+00    1.5042e-07
 8.53000  -2.50000e+00    1.5013e-07
 8.54000  -2.50000e+00    1.5013e-07
 8.55000  -2.50000e+00    1.5048e-07
 8.56000  -2.50000e+00    1.5081e-07
 8.57000  -2.50000e+00    1.5072e-07
 8.58000  -2.50000e+00    1.5025e-07
 8.59000  -2.50000e+00    1.4985e-07
 8.60000  -2.50000e+00    1.4984e-07
 8.61000  -2.50000e+00    1.5050e-07
 8.62000  -2.50000e+00    1.5104e-07
 8.63000  -2.50000e+00    1.5081e-07
 8.64000  -2.50000e+00    1.5001e-07
 8.65000  -2.50000e+00    1.4941e-07
 8.66000  -2.50000e+00    1.4938e-07
 8.67000  -2.50000e+00    1.5047e-07
 8.68000  -2.50000e+00    1.5161e-07
 8.69000  -2.50000e+00    1.5129e-07
 8.70000  -2.50000e+00    1.4983e-07
 8.71000  -2.50000e+00    1.4866e-07
 8.72000  -2.50000e+00    1.4874e-07
 8.73000  -2.50000e+00    1.4984e-07
 8.74000  -2.50000e+00    1.5214e-07
 8.75000  -2.50000e+00    1.5235e-07
 8.76000  -2.50000e+00    1.5027e-07
 8.77000  -2.50000e+00    1.4799e-07
 8.78000  -2.50000e+00    1.4759e-07
 8.79000  -2.50000e+00    1.4819e-07
 8.80000  -2.50000e+00    1.5028e-07
 8.81000  -2.50000e+00    1.5189e-07
 8.82000  -2.50000e+00    1.5165e-07
 8.83000  -2.50000e+00    1.4989e-07
 8.84000  -2.50000e+00    1.4777e-07
 8.85000  -2.50000e+00    1.4671e-07
 8.86000  -2.50000e+00    1.4054e-07
 8.87000  -2.50000e+00    1.3985e-07
 8.88000  -2.50000e+00    1.5087e-07
 8.89000  -2.50000e+00    1.6185e-07
 8.90000  -2.50000e+00    1.5778e-07
 8.91000  -2.50000e+00    1.4506e-07
 8.92000  -2.50000e+00    1.2677e-07
 8.93000  -2.50000e+00    1.3102e-07
 8.94000  -2.50000e+00    1.6266e-07
 8.95000  -2.50000e+00    1.6495e-07
 8.96000  -2.50000e+00    1.4215e-07
 8.97000  -2.50000e+00    1.0815e-07
 8.98000  -2.50000e+00    1.5137e-07
 8.99000  -2.50000e+00    1.7817e-07
 9.00000  -2.50000e+00    1.3941e-07
 9.01000  -2.50000e+00    9.9571e-08
 9.02000  -2.50000e+00    1.8156e-07
 9.03000  -2.50000e+00    1.5103e-07
 9.04000  -2.50000e+00    8.2699e-08
 9.05000  -2.50000e+00    1.8719e-07
 9.06000  -2.50000e+00    1.4565e-07
 9.07000  -2.50000e+00    8.5522e-08
 9.08000  -2.50000e+00    2.0735e-07
 9.09000  -2.50000e+00    5.6883e-08
 9.10000  -2.50000e+00    2.2349e-07
 9.11000  -2.50000e+00    2.9129e-08
 9.12000  -2.50000e+00    2.3874e-07
 9.13000  -2.50000e+00    1.1519e-07
 9.14000  -2.50000e+00    2.0399e-07
 9.15000  -2.50000e+00    8.6206e-08
 9.16000  -2.50000e+00    3.3959e-07
 9.17000  -2.50000e+00    1.1163e-07
 9.18000  -2.50000e+00    1.6757e-07
 9.19000  -2.50000e+00    1.7018e-07
 9.20000  -2.50000e+00    1.7389e-07
 9.21000  -2.50000e+00    1.7951e-07
 9.22000  -2.50000e+00    1.9799e-07
 9.23000  -2.50000e+00    1.9276e-07
 9.24000  -2.50000e+00    1.9569e-07
 9.25000  -2.50000e+00    1.8118e-07
 9.26000  -2.50000e+00    1.8433e-07
 9.27000  -2.50000e+00    1.8299e-07
 9.28000  -2.50000e+00    1.8173e-07
 9.29000  -2.50000e+00    1.8068e-07
 9.30000  -2.50000e+00    1.8183e-07
 9.31000  -2.50000e+00    1.7588e-07
 9.32000  -2.50000e+00    1.8231e-07
 9.33000  -2.50000e+00    1.7908e-07
 9.34000  -2.50000e+00    1.7216e-07
 9.35000  -2.50000e+00    1.7486e-07
 9.36000  -2.50000e+00    1.9173e-07
 9.37000  -2.50000e+00    1.9592e-07
 9.38000  -2.50000e+00    2.2244e-07
 9.39000  -2.50000e+00    2.3829e-07
 9.40000  -2.50000e+00    1.3532e-07
 9.41000  -2.50000e+00    2.5405e-07
 9.42000  -2.50000e+00    1.7617e-07
 9.43000  -2.50000e+00    2.1315e-07
 9.44000  -2.50000e+00    1.7739e-07
 9.45000  -2.50000e+00    2.3203e-07
 9.46000  -2.50000e+00    1.6262e-07
 9.47000  -2.50000e+00    2.5312e-07
 9.48000  -2.50000e+00    1.3991e-07
 9.49000  -2.50000e+00    2.2077e-07
 9.50000  -2.50000e+00    1.9227e-07
 9.51000  -2.50000e+00    1.7129e-07
 9.52000  -2.50000e+00    2.2810e-07
 9.53000  -2.50000e+00    1.8622e-07
 9.54000  -2.50000e+00    1.7280e-07
 9.55000  -2.50000e+00    2.1159e-07
 9.56000  -2.50000e+00    2.0804e-07
 9.57000  -2.50000e+00    1.6912e-07
 9.58000  -2.50000e+00    1.8728e-07
 9.59000  -2.50000e+00    2.1053e-07
 9.60000  -2.50000e+00    2.0073e-07
 9.61000  -2.50000e+00    1.7702e-07
 9.62000  -2.50000e+00    1.8561e-07
 9.63000  -2.50000e+00    1.9497e-07
 9.64000  -2.50000e+00    1.9950e-07
 9.65000  -2.50000e+00    1.8972e-07
 9.66000  -2.50000e+00    1.8320e-07
 9.67000  -2.50000e+00    1.8846e-07
 9.68000  -2.50000e+00    1.9015e-07
 9.69000  -2.50000e+00    1.9016e-07
 9.70000  -2.50000e+00    1.8866e-07
 9.71000  -2.50000e+00    1.8789e-07
 9.72000  -2.50000e+00    1.8893e-07
 9.73000  -2.50000e+00    1.8966e-07
 9.74000  -2.50000e+00    1.8769e-07
 9.75000  -2.50000e+00    1.8698e-07
 9.76000  -2.50000e+00    1.8817e-07
 9.77000  -2.50000e+00    1.8954e-07
 9.78000  -2.50000e+00    1.8975e-07
 9.79000  -2.50000e+00    1.8917e-07
 9.80000  -2.50000e+00    1.8731e-07
 9.81000  -2.50000e+00    1.8714e-07
 9.82000  -2.50000e+00    1.8846e-07
 9.83000  -2.50000e+00    1.8954e-07
 9.84000  -2.50000e+00    1.8947e-07
 9.85000  -2.50000e+00    1.8904e-07
 9.86000  -2.50000e+00    1.8791e-07
 9.87000  -2.50000e+00    1.8776e-07
 9.88000  -2.50000e+00    1.8848e-07
 9.89000  -2.50000e+00    1.8908e-07
 9.90000  -2.50000e+00    1.8907e-07
 9.91000  -2.50000e+00    1.8888e-07
 9.92000  -2.50000e+00    1.8836e-07
 9.93000  -2.50000e+00    1.8823e-07
 9.94000  -2.50000e+00    1.8851e-07
 9.95000  -2.50000e+00    1.8872e-07
 9.96000  -2.50000e+00    1.8868e-07
 9.97000  -2.50000e+00    1.8859e-07
 9.98000  -2.50000e+00    1.8858e-07
 9.99000  -2.50000e+00    1.8859e-07
10.00000  -2.50000e+00    1.8859e-07
fig_5 = plt.figure(figsize=(6,6))

overall_title = 'Gravitational orbit:  ' + \
                rf' $n = {o1.n},$' + \
                rf' $k = {o1.k:.1f},$' + \
                rf' $l = {o1.ang_mom:.1f},$' + \
                rf' $r_0 = {r_0:.1f},$' + \
                rf' $\dot r_0 = {r_dot_0:.2f},$' + \
                rf' $\phi_0 = {phi_0:.2f}$' + \
                '\n'     # \n means a new line (adds some space here)
fig_5.suptitle(overall_title, va='baseline')

ax_5a = fig_5.add_subplot(1,1,1)
ax_5a.semilogy(t_pts, np.abs(E_tot_pts), color='black', label=r'$E(t)$')
ax_5a.semilogy(t_pts, E_tot_rel_pts, color='red', label=r'$\Delta E(t)$')
ax_5a.set_xlabel(r'$t$')
ax_5a.set_ylabel(r'Energy')
ax_5a.set_title('Change in energy with time')
ax_5a.legend(loc='lower right')

fig_5.tight_layout()
../../_images/Orbit_games_13_0.png

Verifying the virial theorem numerically#

First a quick calculation with the current orbit:

t_range = t_pts[-1] - t_pts[0]
U_pts = o1.U(r_pts)
T_pts = o1.energy(t_pts, r_pts, r_dot_pts) - U_pts
from scipy.integrate import simps
T_avg = simps(T_pts, t_pts, delta_t) / t_range
U_avg = simps(U_pts, t_pts, delta_t) / t_range
print(f'<T> = {T_avg:.6f}, <U> = {U_avg:.6f}, n<U>/2 = {o1.n*U_avg/2:.6f}')
<T> = 2.522990, <U> = -5.022990, n<U>/2 = 2.511495

Now a more complete study with different \(n\) values and varying integration times:

from scipy.integrate import simps

# For simplicity, set the angular momentum to 1 and k to +2 or -2 with n*k > 0 
n_values = (-1, 2, 3, 4, 7)   # different n values to test
t_end_values = (10., 100., 1000.)  # integration max values to test
k_magnitude = 2.

# Time start and spacing
t_start = 0.
delta_t = 0.001

# Initial conditions
r_0 = 1.
r_dot_0 = 0.
phi_0 = 0.0
ang_mom = 1.

print(f' n    t_end     <T>        <U>        n<U>/2     |<T>-n<U>/2|/<T>')
for n in n_values:
    k = np.sign(n) * k_magnitude   # ensure that n*k > 0
    o_now = Orbit(ang_mom, n=n, k=k, mu=1.)  # current orbit 
    for t_end in t_end_values:
        t_pts = np.arange(t_start, t_end+delta_t, delta_t)  
        r_pts, r_dot_pts, phi_pts = o_now.solve_ode(t_pts, r_0, 
                                                    r_dot_0, phi_0)

        t_range = t_pts[-1] - t_pts[0]
        # Generate the integrands to be averaged
        U_pts = o_now.U(r_pts)
        T_pts = o_now.energy(t_pts, r_pts, r_dot_pts) - U_pts
        # Calculate the time averages from 0 to t_end
        T_avg = simps(T_pts, t_pts, delta_t) / t_range
        U_avg = simps(U_pts, t_pts, delta_t) / t_range

        rel_error = np.fabs((T_avg - n*U_avg/2.) / T_avg)
        print(f'{n:2.0f}   {t_end:5.0f}   {T_avg:9.6f}  {U_avg:9.6f} ',
              f' {n*U_avg/2.:9.6f}       {rel_error:9.3e}')
 n    t_end     <T>        <U>        n<U>/2     |<T>-n<U>/2|/<T>
-1      10    1.468564  -2.968564    1.484282       1.070e-02
-1     100    1.494275  -2.994275    1.497137       1.915e-03
-1    1000    1.499960  -2.999956    1.499978       1.197e-05
 2      10    1.236029   1.263971    1.263971       2.261e-02
 2     100    1.251595   1.248404    1.248404       2.550e-03
 2    1000    1.250122   1.249865    1.249865       2.050e-04
 3      10    1.512491   0.987509    1.481264       2.065e-02
 3     100    1.498214   1.001784    1.502676       2.978e-03
 3    1000    1.500149   0.999830    1.499745       2.695e-04
 4      10    1.665613   0.834386    1.668773       1.897e-03
 4     100    1.668070   0.831928    1.663856       2.526e-03
 4    1000    1.666828   0.833147    1.666295       3.196e-04
 7      10    1.954318   0.545682    1.909887       2.273e-02
 7     100    1.943101   0.556897    1.949141       3.109e-03
 7    1000    1.944343   0.555637    1.944731       1.994e-04