' h2alib.bas - jackord@kw.igs.net - revised 10 Sep 02 - Liberty Basic v3.01
' trapezoid and "vbar" integration algorithms for mass-spring system

' Initialize Window
  nomainwin
  button#1, "Plot", [plot], UL, 10, 10, 50, 20
  WindowWidth=510     ' pixel scale 0-500
  WindowHeight=389    ' pixel scale 0-360
  UpperLeftX=50: UpperLeftY=100
  open "Trapevoid and Vbar Algorithms" for graphics_nsb as #1
  #1 "trapclose [quit]"
  gosub [pin]

[waitHere]
  wait

[plot]
  gosub [pin]
  k=1: m=25: A=5: xa=0-4: xb=4: dx=(xb-xa)/200: x=xa
  z$="#.######"
  #1 "color blue"
  #1 "place "; 288+int(40*x); " "; 312-int(150/(k*(A*A-x*x)/m)^.5)
  for i=1 to 200
    x=xa+i*dx
    #1 "goto "; 288+int(40*x); " "; 312-int(150/(k*(A*A-x*x)/m)^.5)
  next i
  x=xa: dx=(xb-xa)/4
  #1 "color red"
  #1 "place "; 288+int(40*x); " 312"
  #1 "goto "; 288+int(40*x); " "; 312-int(150/(k*(A*A-x*x)/m)^.5)
  for i=1 to 6
    x=x+dx: x2=288+int(40*x): y2=312-int(150/(k*(A*A-x*x)/m)^.5)
    #1 "goto "; x2; " "; y2
    #1 "line "; x2; " 312 "; x2; " "; y2
    if i=2 then dx=dx/2
  next i
  #1 "color black"
  n=20: gosub [vbar]: s1=s: n=40: gosub [vbar]: s2=s: s3=s2+(s2-s1)/3
  #1 "place 100 32"
  #1 "\Vbar:   T(20) = "; using(z$, s1); "   T(40) = "; using(z$, s2); "   T = "; using(z$, s3)
  n=20: gosub [trap]: s1=s: n=40: gosub [trap]: s2=s: s3=s2+(s2-s1)/3
  #1 "place 100 48"
  #1 "\Trap:   T(20) = "; using(z$, s1); "   T(40) = "; using(z$, s2); "   T = "; using(z$, s3)
  s1=2*(m/k)^.5*atn(2/1.5)
  #1 "place 210 80"
  #1 "\Analytic:   T = "; using(z$, s1)
goto [waitHere]

[pin]
  #1 "cls": #1 "down": #1 "color black"
  #1 "place 88 12": #1 "box 488 312"
  #1 "place 85 332": #1 "\-5"
  #1 "place 485 332": #1 "\5"
  #1 "place 268 344": #1 "\x  (m)"
  #1 "place 72 318": #1 "\0"
  #1 "place 72 18": #1 "\2"
  #1 "place 16 93": #1 "\1/V(x)"
  #1 "place 24 243": #1 "\(s/m)"
return

[trap]
  s=0: dx=(xb-xa)/n
  for i=0 to n
    if i=0 or i=n then c=.5 else c=1
    x=xa+i*dx: s=s+c/(k*(A*A-x*x)/m)^.5
  next i
  s=s*dx
return

[vbar]
  s=0: x=xa: dx=(xb-xa)/n
  v1=(k*(A*A-x*x)/m)^.5
  for i=1 to n
    x=xa+i*dx: v2=(k*(A*A-x*x)/m)^.5
    s=s+2/(v1+v2): v1=v2
  next i
  s=s*dx
return

[quit]
  close #1
end