import pylab

delta_t = 1e-3 # 1 ms

t = 0.0
b = 13.0
M = 2.0
f = 30 # newton
g = 9.81

w = 0.0
theta = 0.0

vettore_theta = [ ]
vettore_w = [ ]
vettore_tempi = [ ]

while t < 40:

    theta_temp = theta + delta_t * w
    w_temp = -g * delta_t * theta + (1 - b*delta_t/M) * w \
		+ delta_t / M * f

    theta = theta_temp;
    w = w_temp;
    t = t + delta_t

    vettore_w.append(w)
    vettore_theta.append(theta)
    vettore_tempi.append(t)


pylab.figure(1)
pylab.plot(vettore_tempi, vettore_w, 'r-+', label='vel, w(t)')
pylab.xlabel('time')
pylab.legend()

pylab.figure(2)
pylab.plot(vettore_tempi, vettore_theta, 'b-+', 
		label='position, theta(t)')
pylab.xlabel('time')
pylab.legend()

pylab.show()