' nm4lib.bas - jackord@kw.igs.net - revised 17 Feb 08 - Liberty Basic v4.02
' n-segment string clamped at both ends for n=8, 48, 96, 192, and 384
' time dependence from FFSS (only the 64 lowest modes) or directly from Newton's Law
' 'pluck' and 'pulse' initial displacements
' y(n/4) plotted versus time

' Initialize Window
  nomainwin
  button#1, "Pluck", [p1], UL, 80, 5, 60, 20
  button#1, "Pulse", [p2], UL, 150, 5, 60, 20
  button#1, "Motion", [p3], UL, 220, 5, 60, 20
  dim a$(6)
  a$(0)="n=8": a$(1)="n=48": a$(2)="n=96": a$(3)="n=192": a$(4)="n=384"
  combobox#1.c1, a$(, [waitHere], 290, 5, 70, 120
  checkbox#1.c2, "FFSS", [waitHere], [waitHere], 370, 5, 80, 20
  WindowWidth=543     ' Pixel Scale 0-525
  WindowHeight=298    ' Pixel Scale 0-260
  UpperLeftX=10: UpperLeftY=10
  open "FFSS for n-Segment String" for graphics_nsb as #1
  #1 "trapclose [quit] ; down"
  #1.c1 "selectindex 1"
  kk=0: plt=0: nf=768
  dim y(385): dim v(385): dim z(385): dim b(65)
  dim w(65): dim xx(385): dim yy(385)
  plt=0

[waitHere]
  wait

[p1]
  kk=1: jt=0: goto [motion]
[p2]
  kk=2: jt=0: goto [motion]
[p3]
  if kk>0 then plt=1: goto [motion]
goto [waitHere]

[motion]
  if jt=0 then
    #1.c1 "selectionindex?": input#1.c1, ind
    if ind=1 then n=8 else n=48*int(2^(ind-2))
    if n<64 then nm=n else nm=64
    redim y(n+1): redim v(n+1): redim b(nm+1): redim w(nm+1)
    redim z(n+1): redim xx(n+1): redim yy(n+1)
    #1 "cls ; rule over ; down ; backcolor blue"
    #1 "color black ; line 395 144 523 144"               ' y(n/4) plot axis
    #1 "line 4 124 4 164 ; line 389 124 389 164"
    for i=0 to n: xx(i)=4+384/n*i: next i
    dt=n/384: pi=4*atn(1)
    if kk=1 then                                          ' Init pluck
      for i=0 to n/2
        y(i)=224*i/n: y(n-i)=y(i)
      next i
    end if
    if kk=2 then                                          ' Init pulse
      for i=3*n/8 to n/2
        y(i)=896*(i-3*n/8)/n: y(n-i)=y(i)
      next i
    end if
    #1.c2 "value?": input#1.c2, z$                        ' FFSS
    if z$="set" then
      for i=1 to nm-1 step 2
        w(i)=2*sin(pi*i/2/n)
        for j=1 to n-1
          b(i)=b(i)+2*y(j)*sin(pi*i*j/n)/n
        next j
      next i
    else
      dt=dt/4
      for i=1 to n-1                                      ' Feynman half-step
        v(i)=v(i)+(y(i-1)+y(i+1)-2*y(i))*dt/2
      next i
    end if
    #1 "rule xor ; color blue"
    yy(0)=144: z(0)=144: z(n)=144
    if n<96 then #1 "place 4 144 ; circlefilled 5 ; circle 5 ; place 388 144 ; circlefilled 5 ; circle 5"
  end if
  tt=time$("ms")
  while jt=0 or plt=1                                     ' Main loop
    for i=1 to n                                          ' Plot y(i)
      yy(i)=144-int(y(i)+.5)
      if jt>0 then
        #1 "line "; xx(i-1); " "; z(i-1); " "; xx(i); " "; z(i)
        if n<96 then
          if i=n/4 then #1 "color red ; backcolor red"
          if i<n then #1 "circlefilled 5 ; circle 5"
          if i=n/4 then #1 "color blue ; backcolor blue"
        end if
      end if
      #1 "line "; xx(i-1); " "; yy(i-1); " "; xx(i); " "; yy(i)
      if n<96 then
        if i=n/4 then #1 "color red ; backcolor red"
        if i<n then #1 "circlefilled 5 ; circle 5"
        if i=n/4 then #1 "color blue ; backcolor blue"
      end if
      z(i-1)=yy(i-1)
    next i
    x2=395+int(jt/6): y2=yy(n/4)
    if jt>0 then                                          ' Plot y(n/4) vs t
      #1 "rule over ; color red ; line "; x1; " "; y1; " "; x2; " "; y2
      #1 "rule xor ; color blue ; backcolor blue ; discard"
    end if
    x1=x2: y1=y2
    jt=jt+1
    if z$="set" then                                      ' FFSS
      for i=1 to n/2
        y(i)=0
        for k=1 to nm-1 step 2
          y(i)=y(i)+b(k)*sin(pi*i*k/n)*cos(w(k)*jt*dt)
        next k
        y(n-i)=y(i)
      next i
    else                                                  ' Feynman algorithm
      for j=1 to 4
        for i=1 to n-1                                    ' Step y(i)
          y(i)=y(i)+v(i)*dt
        next i
        for i=1 to n-1                                    ' Step v(i) (F=ma)
          v(i)=v(i)+(y(i-1)+y(i+1)-2*y(i))*dt
        next i
      next j
    end if
    scan
    if jt=769 then
      plt=0: kk=0: tt=time$("ms")-tt
      #1 "color black ; backcolor white ; place 5 255"
      #1 "\Run time(ms)="; tt
    end if
    #1 "discard"
  wend
goto [waitHere]

[quit]
  close #1
end