' nm2lib - jackord@kw.igs.net - revised 3 Feb 08 - Liberty Basic v4.02
' motion of 8-segment string clamped at both ends
' time dependence from Newton's law
' "Pluck" and "Pulse" and Normal Mode initial displacements
' y(n/4) plotted versus time for Pluck and Pulse
' dispersion relation plotted for Normal Modes

' Initialize Window
  nomainwin
  button#1, "Pluck", [p1], UL, 120, 5, 60, 20
  button#1, "Pulse", [p2], UL, 190, 5, 60, 20
  button#1, "Motion", [p3], UL, 260, 5, 60, 20
  button#1, "Nmodes", [p4], UL, 330, 5, 60, 20
  WindowWidth=543     ' Pixel Scale 0-525
  WindowHeight=298    ' Pixel Scale 0-260
  UpperLeftX=10: UpperLeftY=10
  open "Newton's Law for 8-Segment String" for graphics_nsb as #1
  #1 "trapclose [quit]"
   kk=0: plt=0: n=8: nf=768

[waitHere]
  wait

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

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

[p3]
  if kk=1 or kk=2 then plt=1: goto [plot]
goto [waitHere]

[p4]
  redim y(9): redim v(9): redim xx(9): redim yy(9)
  #1 "cls ; rule over ; down ; color black ; backcolor white"
  #1 "line 393 144 521 144 ; line 393 16 393 144 ; line 393 16 401 16"
  #1 "place 405 21": #1 "\0.4 Hz"
  #1 "line 3 124 3 164 ; line 388 124 388 164"
  #1 "color blue ; backcolor blue"
  #1 "place 4 144 ; circlefilled 5 ; place 388 144 ; circlefilled 5"
  n=8: pi=4*atn(1): yy(0)=144: z(0)=144: z(n)=144
  for i=0 to 8: xx(i)=4+48*i: next i
  for mode=1 to n-1                                         ' Mode loop
    #1 "rule xor"                                           ' xor animation
    t=0: dt=0.1/mode/2: tst=1: plt=1: tf=1
    for i=0 to n
      y(i)=108*sin(pi*mode*i/n)
    next i
    #1 "color red"                                          ' Displacement curve
    y1=144
    for i=1 to 384
      y2=144-int(108*sin(pi*mode*i/384)+.5)
      #1 "line "; i+3; " "; y1; " "; i+4; " "; y2
      y1=y2
    next i
    #1 "color blue ; backcolor blue"
    for i=1 to n-1                                          ' Feynman half-step
      v(i)=(y(i-1)+y(i+1)-2*y(i))*dt/2
    next i
    while plt=1                                             ' Main loop
      for i=1 to n                                          ' Plot y(i)
        yy(i)=144-int(y(i)+.5)
        if t>0 then
          #1 "line "; xx(i-1); " "; z(i-1); " "; xx(i); " "; z(i)
          if i<n then #1 "circlefilled 6 ; circle 6"
        end if
        #1 "line "; xx(i-1); " "; yy(i-1); " "; xx(i); " "; yy(i)
        if i<n then #1 "circlefilled 6 ; circle 6"
        z(i-1)=yy(i-1)
      next i
      if t=0 then
        timer 1000, [pause1]: wait
        [pause1]
        timer 0
      end if
      t=t+dt                                                ' Step t
      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
      if tst=1 and y(1)<0 then
        tf=4*(t-y(1)/v(1)): tst=0
      end if
      scan
      if tst=0 and t>tf+dt/2 then plt=0
    wend
    #1 "place "; 393+16*mode; " "; 144-int(128/.4/tf+.5)    ' Plot frequency
    #1 "circlefilled 5 ; circle 5"
    #1 "backcolor white ; color black ; place "; 393+16*mode-4; " 160"
    #1 "\"; mode: #1 "discard ; color blue ; backcolor blue"
    timer 2000, [pause2]: wait
    [pause2]
    timer 0
    if mode<n-1 then                                        ' Clear display
      #1 "color red"
      y1=144
      for i=1 to 384
        y2=144-int(108*sin(pi*mode*i/384)+.5)
        #1 "line "; i+3; " "; y1; " "; i+4; " "; y2
        y1=y2
      next i
      #1 "color blue ; backcolor blue"
      for i=1 to n
        #1 "line "; xx(i-1); " "; z(i-1); " "; xx(i); " "; z(i)
        if i<n then #1 "circlefilled 6 ; circle 6"
      next i
    end if
  next mode
  #1 "color red ; place 393 144"                            ' Dispersion relation
  for i=1 to 128
    #1 "goto "; 393+i; " "; 144-int(128/.4/pi*sin(i*pi/256)+.5)
  next i
goto [waitHere]

[plot]
  if jt=0 then
    redim y(9): redim v(9): redim z(9): redim xx(9): redim yy(9)
    #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)=5+384/n*i: next i
    dt=n/768
    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
    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
    #1 "rule xor ; color blue"                              ' xor animation
    yy(0)=144: z(0)=144: z(n)=144
    if n<96 then
      #1 "place 4 144 ; circlefilled 6 ; circle 6"
      #1 "place 388 144 ; circlefilled 6 ; circle 6"
    end if
  end if
  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 6 ; circle 6"
          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 6 ; circle 6"
        if i=n/4 then #1 "color blue ; backcolor blue"
      end if
      z(i-1)=yy(i-1)
    next i
    x2=395+int(jt/12): 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"
    end if
    x1=x2: y1=y2
    jt=jt+1                                                 ' Step t
    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
    if jt=1537 then plt=0: kk=0
    scan: #1 "discard"
  wend
goto [waitHere]

[quit]
  close #1
end