def deriv(X, t):
    x, v = X.reshape(2, -1)
    acc = -GM * x * ((x**2).sum())**-1.5
    return np.hstack((v, acc))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

# https://physics.stackexchange.com/questions/348854/parker-solar-probe-passing-extremely-close-to-the-sun-what-relativistic-effects

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

c    = 2.9979E+08
a    = 57.8E+06 * 1000.
GM   = 1.327E+20
peri = 6.6E+06 * 1000.
apo  = 2*a - peri
vapo, vperi = [np.sqrt(GM*(2./r - 1./a)) for r in (apo, peri)]

T    = twopi * np.sqrt(a**3/GM)

X0   = np.hstack([apo, 0, 0, vapo])

time = np.linspace(0, T, 1000001)
days = time/(24.*3600.)

answer, info = ODEint(deriv, X0, time, full_output=True, rtol=1E-10)

x, v = answer.T.reshape(2, 2, -1)
r    = np.sqrt((x**2).sum(axis=0))

dfof = -(GM/c**2) * (2./r - 0.5/a)
dt   = time[1] - time[0]

delta_t = dt * dfof
DT      = delta_t.cumsum()

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    plt.plot(x[0]/1000., x[1]/1000.)
    plt.plot([0], [0], 'ok')
    plt.xlim(-0.1E+08, 1.1E+08)
    plt.ylim(-0.6E+08, 0.6E+08)
    plt.title('heliocentric orbit (km)', fontsize=14)
    plt.subplot(4, 1, 3)
    plt.plot(days, dfof)
    plt.title('delta f/f vs time (days)', fontsize=14)
    plt.subplot(4, 1, 4)
    plt.plot(days, DT)
    plt.title('cumulative time (sec) vs time (days)', fontsize=14)
    plt.show()