' fosalib.bas - jackord@kw.igs.net - revised 3 Apr 05 - Liberty Basic v4.01
' identical masses hanging from springs one below the other
' they are released from rest with the upper spring unstretched and
'   Case 1: with the lower mass supported by the lower spring
'   Case 2: with the lower mass raised enough to generate the
'             low-frequency normal mode
'   Case 3: with the lower mass pulled down enough to generate the
'             high-frequency normal mode
' each spring (length 3*d) stretches d under the weight of one mass

' Initialize Window
  nomainwin
  button#1, "Case1", [c1], UL, 90, 5, 55, 20
  button#1, "Case2", [c2], UL, 155, 5, 55, 20
  button#1, "Case3", [c3], UL, 220, 5, 55, 20
  button#1, "Motion", [mm], UL, 285, 5, 55, 20
  WindowWidth=435     ' pixel scale 0-425
  WindowHeight=349    ' pixel scale 0-320
  UpperLeftX=100: UpperLeftY=100
  open "Two Oscillators" for graphics_nsb as #1
  #1 "trapclose [quit]"
  dim xs(121)
  pi=4*atn(1): d=25: dt=.1: plt=0
  for i=0 to 120                                       ' Define spring
    xs(i)=10+int(8*sin(i*pi/10))
  next i

[waitHere]
  wait

[c1]
  kk=1: jt=0
goto [plot]

[c2]
  kk=2: jt=0
goto [plot]

[c3]
  kk=3: jt=0
goto [plot]

[mm]
  if kk>0 then plt=1: goto [plot]
goto [waitHere]

[td]                                                   ' Find period difference
  t=0: ya=3*d: yb=yt: ta=0: tb=0
  aa=d-2*ya+yb: ab=4*d-yb+ya
  va=aa*dt/2: vb=ab*dt/2: vb1=vb                       ' Feynman half-step
  while ta=0 or tb=0
    t=t+dt: ya=ya+va*dt: yb=yb+vb*dt                   ' Step ya, yb
    aa=d-2*ya+yb: ab=4*d-yb+ya                         ' Newton's law
    va=va+aa*dt: vb=vb+ab*dt                           ' Step va, vb
    if ta=0 and va<0 then ta=t+dt/2-va/aa
    if tb=0 and vb*vb1<0 then tb=t+dt/2-vb/ab
  wend
  x=abs(ta-tb)
return

[plot]
  if jt=0 then
    #1 "backcolor white ; cls ; down ; rule xor ; color black"
    if kk=1 then yb=7*d
    if kk>1 then                                       ' Vary yb to match periods
      dd=d/4: yt=(4*kk-2)*d: gosub [td]
      while abs(dd)>.02
        yt=yt+dd: xo=x: gosub [td]
        if x>xo then dd=0-dd/2
      wend
      yb=yt+2*dd
      #1 "place 40 40": #1 "\Amplitude Ratio = "; using("##.###",(9-yb/d)/2)
    end if
    #1 "line 20 125 425 125 ; line 20 225 425 225"
    ya=3*d: aa=d-2*ya+yb: ab=4*d-yb+ya: nm=0
    va=aa*dt/2: vb=ab*dt/2                             ' Feynman half-step
  end if
  #1 "backcolor red"
  tt=time$("ms")
  while jt=0 or plt=1
    ya2=int(ya+.5): yb2=int(yb+.5)
    #1 "color red"
    if jt>0 then #1 "place 10 "; ya1: #1 "circlefilled 7 ; circle 7"
    #1 "place 10 "; ya2: #1 "circlefilled 7 ; circle 7"
    if jt>0 then #1 "place 10 "; yb1: #1 "circlefilled 7 ; circle 7"
    #1 "place 10 "; yb2: #1 "circlefilled 7 ; circle 7"
    #1 "color black"
    for j=1 to 120                                     ' Draw spring
      if jt>0 then #1 "line "; xs(j-1); " "; int(ya1*(j-1)/120); " "; xs(j); " "; int(ya1*j/120)
      #1 "line "; xs(j-1); " "; int(ya2*(j-1)/120); " "; xs(j); " "; int(ya2*j/120)
      if jt>0 then #1 "line "; xs(j-1); " "; int(ya1+(yb1-ya1)*(j-1)/120); " "; xs(j); " "; int(ya1+(yb1-ya1)*j/120)
      #1 "line "; xs(j-1); " "; int(ya2+(yb2-ya2)*(j-1)/120); " "; xs(j); " "; int(ya2+(yb2-ya2)*j/120)
    next j
    if jt>0 then
      #1 "color blue"
      #1 "line "; 19+jt; " "; ya1; " "; 20+jt; " "; ya2
      #1 "line "; 19+jt; " "; yb1; " "; 20+jt; " "; yb2
    end if
    ya1=ya2: yb1=yb2
    ya=ya+va*dt: yb=yb+vb*dt                           ' Step ya, yb
    aa=d-2*ya+yb: ab=4*d-yb+ya                         ' Newton's law
    va=va+aa*dt: vb=vb+ab*dt                           ' Step va, vb
    jt=jt+1
    if jt=406 then
      kk=0: plt=0: tt=time$("ms")-tt
      #1 "backcolor white ; color black ; place 240 40"
      #1 "\Run Time(ms)="; tt
    end if
    scan: #1 "discard"
  wend
goto [waitHere]

[quit]
  close #1
end