#by Ruediger Kuersten, April 18, 2019 for teaching
# start this skript on Linux or MacOs by typing
# python blatt2_3.py
# for Windows users: google what to do
import scipy as sp
import scipy.integrate as integrate
import scipy.optimize as opt
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as matplotlib
import math

# Parameter
k=1.0
gamma=2.5
tau=1
lamda=1.5
m=1.0
F=1.0
omega=2.0*np.pi/tau

# Schrittweite
deltat=0.0001

# Anfangsbedingungen
x=[1.0]
v=[0.0]
t=[0.0]

# Langzeitverhalten als Fourierreihe
def x2(t):
    result=F/k*(1.0+1.0/tau/lamda*(np.exp(-lamda*tau)-1.0))
    modenumber=10
    for n in range(1, modenumber+1):
        omegan=n*omega
        an=2.0*(k*omegan+gamma*lamda*omegan-m*omegan**3)/((k-m*omegan**2)**2 + (gamma*omegan)**2)/( lamda**2 +omegan**2)*F/tau*(np.exp(-lamda*tau)-1.0)
        bn=2.0*(k*lamda-lamda*m*omegan**2-gamma*omegan**2)/((k-m*omegan**2)**2 + (gamma*omegan)**2)/( lamda**2 +omegan**2)*F/tau*(np.exp(-lamda*tau)-1.0)
        result=result+an*np.sin(omegan*t)
        result=result+bn*np.cos(omegan*t)
    return result


vecx=np.vectorize(x2)

# externe Kraft, sehr ineffizient implementiert
def f(t):
    tprime=t
    while (tprime>1.0):
        tprime=tprime-1.0
    return F*(1.0-np.exp(-lamda*tprime))

# Numerische Integration der Bewegungsgleichung (Eulerverfahren)
for i in range(1,500000):
    t.append(deltat*i)
    x.append(x[i-1] + deltat*v[i-1])
    v.append(v[i-1] + deltat*(  -gamma/m*v[i-1] -k/m*x[i-1] +1.0/m*f(t[i-1])  )   )


#Plot x(t)
plot1 = plt.subplot(111)
# x-Achse Darstellungsbereich
plot1.set_xlim(1, 50)
plt.plot( t, x, '-', color='r', mfc='none')
plt.plot( t, vecx(t), '-', color='b', mfc='none')
plt.ylabel(r'$x(t)$')
plt.xlabel(r'$t$')
#speichern als pdf?
#plt.savefig('xt.pdf')
plt.show()