from visual import *

#================================
length = 16.0 # meters
springk = 5.0 # Newtons/meter
smass = 0.001 # kg/m
momentum_amplitude = 0.05
#================================

dt = 0.005 # sec time step
threshold = 0.1
g = []
for i in range(1):
  g.append(display(y = 80 + 200*i, width=600, height=170))

g[0].title="Wave1.py"

bands = [ curve( x = arange(-length/2.0,length/2.0,length/100.0), display=g[0], color=color.red, k = springk, mass = smass),
      ]
ampli = 0.0
for band in bands:
    band.radius = 0.05
    band.momentum = zeros((100,3),Float)
    band.momentum[:25,1] = sin(band.x[:25]*pi / 60.0)*momentum_amplitude   # half-wave pulse

t = 0.0
flag = 0
while 1:
    rate(100)
    for band in bands:
        # Keep endpoints fixed:
        band.momentum[0] = band.momentum[-1] = vector(0,0,0)

        # Integrate velocity:
        band.pos = band.pos + ((band.momentum/band.mass)*dt)

        if (max(band.y) > ampli):
          ampli = max(band.y)

        if ((flag == 0) and (band.pos[50,1]>threshold)):
          flag=1
        if (band.pos[50,1]<=threshold):
          flag=0

        # force[n] is the force on point n from point n+1 (to the right):
        force = band.k * (band.pos[1:] - band.pos[:-1])

        # all points but the last experience forces to the right:
        band.momentum[:-1] = band.momentum[:-1] + force * dt

        # all points but the first experience forces to the left:
        band.momentum[1:] = band.momentum[1:] - force * dt

        if (flag == 1):
          flag = 2
          threshold = 0.3*ampli
          print "Amplitude= ",ampli
          print "Time  (when amplitude goes over  ",threshold," m) = ",t," sec"
        t = t + dt