Spring through damping medium

from __future__ import division
from visual.graph import *

scene.y = 400 # moves render box 400 down.


# Constants
pi = 3.14159
L0 = .4 # equilibrium length
g = 0 # no grav
k = 17.6 # spring const
b = 1.8 # dat const

scene.center = vector(0,-.1,0)

# Objects
ceiling = box(pos=vector(0,0,0), size=(.3,0.005,0.005))
block=box(pos=vector(0,-0.1,0), size=(.02,0.02,0.02), color=color.yellow)
spring = cylinder(pos=ceiling.pos, axis=block.pos, radius=.005)

# Initial values
block.m = 6.8
block.v = vector(0,0,0) #block V_i = 0
block.p = block.m * block.v
block.pos = vector(0,-L0-0.2,0) # initial position 0.05m from equilibrium

deltat = .0001
t = 0
W = 0
displacement = 0
Kgraph = gcurve(color=color.cyan)
Ugraph = gcurve(color=color.yellow)
KplusUgraph = gcurve(color=color.red)
Wgraph = gcurve(color=color.green)

# Loop for repetitive calculations
scene.autoscale=0
while t < 10:
    Fnet = -(((block.pos-vector(0,-L0,0)) * k)) - ((block.p / block.m) * b) #force spring on block
    displacement = (mag(block.p) / block.m) * deltat
    block.p = block.p + Fnet * deltat # updates the momentum
    block.pos = block.pos + block.p / block.m * deltat # updates the position
    spring.axis = block.pos #updates the spring axis so it stays on the block
    t = t + deltat
    pmag = mag(block.p)

    K = (pmag**2)*.5/block.m #kinetic energy ((1/2)*|p|^2) / m = K
    U = ((mag(block.pos)-L0)**2)*.5*k #potential energy (|r| - L0)^2 * (1/2) * k = U
    W = W - displacement*(mag(block.p)/block.m)*b #work dr*(|p|/m)*b (ini)
    Kgraph.plot(pos=(t,K))
    Ugraph.plot(pos=(t,U))
    KplusUgraph.plot(pos=(t,K+U))
    Wgraph.plot(pos=(t,W))
    E = K + U
    #print E,'--------',U