shell oscillator

#
# Set program constants
#
BobMass = .02
SpringConstant = .7
DampingBeta = 0
InitialStretch = .2
EquilLength = .5
ExpectedPeriod = 2*pi*sqrt(BobMass/SpringConstant)
#
scene2 = canvas(title='Damped Oscillator',caption='Animated Display',
     center=vector(0.50,0,0), background=color.white)
#
bob    = sphere(pos=vector(0,0,0),radius=0.05,color=color.red)
wall   = box(pos=vector(0,0,0),size=vector(0.05,.1,.1),color=color.blue)
spring = helix(pos=wall.pos,axis=bob.pos-wall.pos,radius=0.01,thickness=.004,coils=10,color=color.green)
#
# equilibrium position of the end of the spring.
#
length  = vector(EquilLength,0,0)
stretch = vector(InitialStretch,0,0)
#
# Set the initial position of the bob
#
bob.pos = wall.pos + length + stretch
bob.mom  = vector(0,0,0)
#
# Input Parameters needed in the program. Be sure to
# choose sensible values.
#
bob.mass = BobMass
spring.ks  = SpringConstant
beta = DampingBeta
#
print('Damped Mass on a Spring')
#
# Time step and total elapsed time
#
dt = 0.005
t  = 0.0
#
# Used to look for zero crossings to measure the period.
#
told = 0.0
xold = bob.pos.x
#
# Setup a graph window to plot things in
#
s='<b>Mass and Spring: Graph</b>'
#
# Move the mouse over the graph to explore its interactivity.
# Drag a rectangle in the graph to zoom. Examine the icons at the upper right.
# Click the "Reset axes" icon to restore. Drag along the bottom or left to pan.
#

graph(title=s, xtitle='Time',ytitle='Energy', xmas=10, ymax=.05, ymin=-.05, x=0, y = 500, width=500, height=300)
#

drawKE = gcurve(color=color.cyan,label='Kinetic Energy')
drawPE = gcurve(color=color.blue,label='Potential Energy')
drawTE = gcurve(color=color.magenta,label='Total Energy')
#
while(t<10):
    rate(100)
    t += dt
#
#    spring.stretch = block.pos-wall.pos
#
    spring.force = -SpringConstant*(bob.pos-length)
    bob.mom = bob.mom + spring.force*dt
    bob.pos = bob.pos + bob.mom*(dt/bob.mass)
    spring.axis = bob.pos-wall.pos

    KE = .5*((bob.mom.x)**2/bob.mass)
    PE = .5*SpringConstant*(bob.pos.x-length.x)**2
    TE = KE + PE

    drawKE.plot(pos=(t,KE))
    drawPE.plot(pos=(t,PE))
    drawTE.plot(pos=(t,TE))
#
# Check for a zero crossing
#
    xnew = bob.pos.x - wall.pos.x - length.x
    if xnew*xold <= 0:
        period = 2*(t - told)
        if told != 0:
            scene2.caption=('Expected period is',ExpectedPeriod,' actual period is',period,'.')
        told = t
    xold = xnew
#
# Plot the x-coordinate of the block as a function of time.
#
#
print("Expected Period (sec)",ExpectedPeriod)
print("Actual Period   (sec)",period)
print('All Done')