import numpy as np
import matplotlib.pyplot as plt

# konstanter
k = 1. # fjaerkonstant i kg/s**2
m = 1. # masse i kg

# egenfrekvens i Hz
omega_0 = np.sqrt( k/m )

# tidskritt i s
delta_t = 0.01

# antal tidsskritt
nn = 10000

# array for tid
t = np.empty( nn )
# array for utslaget
x = np.empty( nn )
# array for hastighet
v = np.empty( nn )

# initialbetingelser
t[0] = 0.
x[0] = 0. # initiale utslaget er 0
v[0] = 1. # initiale hastigheten er 1 m/s

for j in range( 1, nn ):

    # gaa et halvt tidsskritt
    vtemp = v[j-1] - omega_0**2*x[j-1]*delta_t/2.
    # posisjon
    xtemp = x[j-1] + v[j-1]*delta_t/2.

    # hastighet
    v[j] = v[j-1] - omega_0**2*xtemp*delta_t
    # posisjon
    x[j] = x[j-1] + vtemp*delta_t

    # gaa et tidskritt
    t[j] = t[j-1] + delta_t

# beregner potensielle energi
e_p = 0.5*k*x**2

# beregner kinetisk energi
e_k = 0.5*m*v**2

plt.figure()
plt.plot( t, x, 'r', t, v, 'g' )
plt.xlabel("tid")
plt.ylabel("utslag/hastighet")
plt.show()

plt.figure()
plt.plot( t, e_p, 'r', t, e_k, 'g', t, e_p+e_k, 'b' )
plt.xlabel("tid")
plt.ylabel("utslag/hastighet")
plt.show()