import numpy as np
import matplotlib.pyplot as plt


x0 = 2
v0 = 0
dt = 0.5
tMax = 2000


dtl = 10.**(-np.arange(1,5))
feil = np.zeros(len(dtl))
feilE = np.zeros(len(dtl))
feilM = np.zeros(len(dtl))


k = 0

for dt in dtl:
    dt2 = dt*dt/2
    n = int(tMax/dt)
    t = np.zeros(n)
    a = np.zeros(n)
    v = np.zeros(n)
    x = np.zeros(n)
    e = np.zeros(n)
    ae = np.zeros(n)
    ve = np.zeros(n)
    xe = np.zeros(n)
    ee = np.zeros(n)  
    aem = np.zeros(n)
    vem = np.zeros(n)
    xem = np.zeros(n)
    eem = np.zeros(n)  
    x[0] = x0
    v[0] = v0
    xe[0] = x0
    ve[0] = v0
    xem[0] = x0
    vem[0] = v0
    t[0] = 0
    e[0] = 0.5*v[0]**2 + 1-np.cos(x[0])
    for i in range(n-1):
        t[i+1] = t[i] +dt
        a[i] = -x[i]
        v[i+1] = v[i] + a[i]*dt
        x[i+1] = x[i] + v[i+1]*dt
        ae[i] = -xe[i]
        ve[i+1] = ve[i] + ae[i]*dt
        xe[i+1] = xe[i] + ve[i]*dt
        vem[i+1] = vem[i] -xem[i]*dt - vem[i]*dt2
        xem[i+1] = xem[i] + vem[i]*dt - xem[i]*dt2
        e[i+1] = 0.5*v[i+1]**2 + 1-np.cos(x[i+1])
        ee[i+1] = 0.5*ve[i+1]**2 + 1-np.cos(xe[i+1])

    feil[k] = x[-1]-x0*np.cos(t[-1])
    feilE[k] = xe[-1]-x0*np.cos(t[-1])
    feilM[k] = xem[-1]-x0*np.cos(t[-1])
    print(dt,k,feil[k],feilE[k],feilM[k])
    k = k+1
     
plt.loglog(dtl,abs(feil),'o', label='Euler Cromer')
plt.loglog(dtl,abs(feilE),'o', label='Eulers metode')
plt.loglog(dtl,abs(feilM),'o', label='Euler midtpunkt')
plt.loglog(dtl,dtl*50)
plt.legend(loc='upper left')
plt.xlabel(r'$\Delta t$')
plt.ylabel(r'$F(\Delta t)$')
plt.ylim(10**-6,10**5)
plt.savefig('HOVarDt2000.pdf')
#fig, ax = plt.subplots(nrows=2)
#ax[0].plot(t,x0*np.cos(t))
#ax[0].plot(t,x)
#ax[0].plot(t,xe)
#ax[1].plot(t,e)
#ax[1].plot(t,ee)


plt.show()