' nm7lib.bas - jackord@kw.igs.net - Revised 17 Feb 08 - Liberty Basic v4.02
' wave pulse incident on a variable thickness segment where the mass density
'   increases by a factor of 1.44 followed by a region where the mass density
'   increases again by a factor of 1.44
' the motion is generated directly from Newton's Law
'   or from a superposition of (non-orthogonal) normal mode displacements

  nomainwin
  button#1, "Pulse B", [p1], UL, 110, 5, 60, 20
  button#1, "Motion", [p2], UL, 180, 5, 60, 20
  dim a$(3)
  a$(0)="D(film)=0": a$(1)="D(film)=5": a$(2)="D(film)=10"
  combobox#1.c1, a$(, [waitHere], 10, 5, 90, 80
  checkbox#1.c2, "Series", [waitHere], [waitHere], 250, 5, 70, 20
  WindowWidth=498     ' Pixel Scale 0-480
  WindowHeight=258    ' Pixel Scale 0-220
  UpperLeftX=10: UpperLeftY=10
  open "Pulse Striking Interface" for graphics_nsb as #1
  #1 "trapclose [quit]"
  #1.c1 "selectindex 1"
  #1 "down ; rule xor"
  n=480: nm=21: ns=80: nf=320: nk=4: kk=0: dt=0.9: pi=4*atn(1): no=48: nw=96

[waitHere]
  wait

[p1]
  kk=1: jt=0: plt=0: goto [motion]
goto [waitHere]

[p2]
  if kk>0 then plt=1: goto [motion]
goto [waitHere]

[f]
  q=pi/n                                                    ' The upper right...
  km=1.2*k1: k2=1.2*km
  c1=cos(q*k1*n/2): s1=sin(q*k1*n/2)
  c=cos(q*km*dl): s=sin(km*q*dl)
  c2=cos(q*k2*(n/2-dl)): s2=sin(q*k2*(n/2-dl))
  fb=km*km*s*s1*s2-k1*c1*km*c*s2-k1*c1*k2*c2*s-km*c*k2*c2*s1
return                                                      ' ...matrix element

[findk]
  cnt=0: e=.00001                                           ' Find the k-values for
  k1=.5: gosub [f]: a=fb: k1=.6                              '   which the matrix element = 0
  while cnt<ns                                              '
    gosub [f]                                               ' Scan for zero crossings
    if a*fb<0 then cnt=cnt+1: k(cnt)=k1                     '   in the range of interest
    a=fb: k1=k1+.05
  wend
  for i=nm to ns                                            ' Use a binary search
    a=k(i)-.1: b=k(i): k1=a: gosub [f]: fa=fb               '   to refine the values
    while (b-a)>e
      k1=(a+b)/2: gosub [f]: if fb*fa>0 then a=k1 else b=k1
    wend
    k(i)=(a+b)/2
  next i
return

[fey]
  c=pi/n: c=c*c: dx=.05                                     ' Find modes nm to ns
  for i=nm to ns
    e=k(i)*k(i): a=0-e*c: s=0
    z=0: v=1: ii=0: jj=0
    for j=1 to 9600
      z=z+v*dx
      if j=4800 then a=1.44*a
      if j=(4800+20*dl) then a=1.44*a
      v=v+a*z*dx: s=s+z*z: jj=jj+1
      if jj=20 then jj=0: ii=ii+1: q(i, ii)=z
    next j
    s=sqr(s*dx)
    for j=1 to n: q(i, j)=q(i, j)/s: next j                 ' Normalize them
    scan
  next i
return

[motion]
  if jt=0 then
    redim y(n+1): redim v(n+1): redim yy(n+1): redim z(n+1)
    redim k(ns+1): redim w(ns+1): redim b(ns+1): redim d(ns+1)
    redim bb(ns+1): redim dd(ns+1): redim o(ns+1, ns+1): redim q(ns+1, n+1)
    e1=0: e2=0
    for i=no to no+nw                                       ' Init Pulse B
      q=(i-no)*2*pi/nw
      y(i)=40*(1-cos(q))
      v(i)=0-40*2*pi/nw*sin(q)
      v(i)=v(i)*cos(nk*q)+y(i)*nk*2*pi/nw*sin(nk*q)
      y(i)=y(i)*cos(nk*q)
      e1=e1+v(i)*v(i)
    next i
    #1 "cls"
    #1.c1 "selectionindex?": input#1.c1, ind
    dl=5*(ind-1)
    #1.c2 "value?": input#1.c2, z$
    if z$="set" then
      #1 "place 370 20 ; color black"
      #1 "\busy"
      gosub [findk]
      gosub [fey]
      for i=nm to ns
        w(i)=k(i)*pi/n                                      ' No dispersion
        for j=1 to n-1                                      ' Expansion coefficients
          bb(i)=bb(i)+y(j)*q(i, j)                          '   (first approximation)
          dd(i)=dd(i)+v(j)/w(i)*q(i,j)
        next j
      next i
      for i=nm to ns-1                                      ' Mode overlap integrals
        for j=i+1 to ns
          s1=0
          for jj=1 to n-1
            s1=s1+q(i, jj)*q(j, jj)
          next jj
          o(i, j)=s1: o(j, i)=s1
        next j
      next i
      for i=nm to ns                                        ' Expansion coefficients
        s1=0: s2=0                                          '   (second approximation)
        for j=nm to ns
          s1=s1+bb(j)*o(i, j)
          s2=s2+dd(j)*o(i, j)
        next j
        b(i)=bb(i)-s1
        d(i)=dd(i)-s2
      next i
      s1=0
      for i=nm to ns
        s1=s1+b(i)*q(i, no+nw/2)
      next i
      for i=nm to ns                                        ' Expansion coefficients
        b(i)=b(i)*80/s1: d(i)=d(i)*80/s1                    '   (final correction)
      next i
      for j=1 to n-1                                        ' Series expansion
        y(j)=0
        for i=nm to ns
          y(j)=y(j)+b(i)*q(i, j)
        next i
      next j
    else
      for i=no to no+nw                                     ' Feynman half-step
        v(i)=v(i)+(y(i-1)+y(i+1)-2*y(i))*dt/2
      next i
    end if
    yy(0)=110: z(0)=110: z(n)=110
  end if
  while jt=0 or plt=1                                       ' Main loop
    #1 "place 370 20 ; color black"
    #1 "\Frame "; jt; "/320"
    #1 "color blue"
    for i=1 to n                                            ' Plot string
      yy(i)=110-int(y(i)+.5)
      if i=n/2 then #1 "color green"
      if i=n/2+dl then #1 "color red"
      if jt>0 then #1 "line "; i-1; " "; z(i-1); " "; i; " "; z(i)
      #1 "line "; i-1; " "; yy(i-1); " "; i; " "; yy(i)
      z(i-1)=yy(i-1)
    next i
    jt=jt+1
    if z$="set" then
      for i=1 to n-1
        y(i)=0
        for j=nm to ns
          y(i)=y(i)+q(j, i)*(b(j)*cos(w(j)*jt)+d(j)*sin(w(j)*jt))
        next j
      next i
    else
      for i=1 to n-1                                        ' Feynman algorithm
        y(i)=y(i)+v(i)*dt
      next i
      m=1
      for i=1 to n-1
        if i=n/2 then m=1.44*m
        if i=n/2+dl then m=1.44*m
        v(i)=v(i)+(y(i-1)+y(i+1)-2*y(i))/m*dt
      next i
    end if
    scan
    if jt=nf+1 then
      plt=0: kk=0
      for i=1 to n/2-1
        if z$="set" then
          v(i)=0
          for j=nm to ns
            v(i)=v(i)+w(j)*q(j, i)*(d(j)*cos(w(j)*nf)-b(j)*sin(w(j)*nf))
          next j
        end if
        e2=e2+v(i)*v(i)
      next i
      #1 "place 10 45 ; color black"
      #1 "\Reflects "; using("#.#",100*e2/e1); " % of incident energy"
    end if
    #1 "discard"
  wend
goto [waitHere]

[quit]
  close #1
end