M = 600000
T = 0.6
z = 0.10
p0 = 400000
t0, t1, t2, t3 = 0.0, 1.0, 3.0, 6.0

"""
Explanation: Duhamel Integral
Problem Data
End of explanation
"""


wn = 2*np.pi/T
wd = wn*np.sqrt(1-z**2)

"""
Explanation: Natural Frequency, Damped Frequency
End of explanation
"""


dt = 0.05
edt = np.exp(-z*wn*dt)
fac = dt/(2*M*wd)

"""
Explanation: Computation
Preliminaries
We chose a time step and we compute a number of constants of the integration procedure that depend on the time step
End of explanation
"""


t = dt*np.arange(1+int(t3/dt))

"""
Explanation: We initialize a time variable
End of explanation
"""


p = np.where(t<=t1, p0*(t-t0)/(t1-t0), np.where(t<t2, p0*(1-(t-t1)/(t2-t1)), 0))

s = np.sin(wd*t)
c = np.cos(wd*t)
sp = s*p
cp = c*p

plt.plot(t, p/1000)
plt.xlabel('Time/s')
plt.ylabel('Force/kN')
plt.xlim((t0,t3))
plt.grid();

"""
Explanation: We compute the load, the sines and the cosines of $\omega_D t$ and their products
End of explanation
"""


A, B, x = 0, 0, [0]

for i, _ in enumerate(t[1:], 1):
    A = A*edt+fac*(cp[i-1]*edt+cp[i])
    B = B*edt+fac*(sp[i-1]*edt+sp[i])
    x.append(A*s[i]-B*c[i])

"""
Explanation: The main (and only) loop in our code, we initialize A, B and a container for saving the deflections x,
then we compute the next values of A and B, the next value of x is eventually appended to the container.
End of explanation
"""


x = np.array(x)
k = M*wn**2
Dst = p/k

plt.plot(t, x*1000)
plt.plot(t, Dst*1000)
plt.xlabel('Time/s')
plt.ylabel('Deflection/mm')
plt.xlim((t0,t3))
plt.grid()
plt.show();

"""
Explanation: It is necessary to plot the response.
End of explanation
"""


