' qm3lib.bas - jackord@kw.igs.net - revised 18 Feb 08 - Liberty Basic v4.02
' quantum wave pulse incident on a variable thickness region ("film") where
'   the kinetic energy increases by a factor 1.44 followed by a region where
'   it increases again by a factor of 1.44
' the motion is generated either directly from the Schrodinger Equation
'   or from an expansion in terms if the stationary states of the system


' Initialize Window
  nomainwin
  button#1, "P Dens", [p1], UL, 20, 5, 60, 20
  button#1, "R Comp", [p2], UL, 90, 5, 60, 20
  button#1, "Motion", [p3], UL, 160, 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], 230, 5, 105, 90
  checkbox#1.c2, "Series", [waitHere], [waitHere], 345, 5, 70, 20
  WindowWidth=498     ' Pixel Scale 0-480
  WindowHeight=258    ' Pixel Scale 0-220
  UpperLeftX=10: UpperLeftY=10
  open "Quantum Pulse Striking Inferface" for graphics_nsb as #1
  #1.c1 "selectindex 1"
  #1 "trapclose [quit]"
  #1 "rule xor ; down"
  n=480: nf=320: nk=4: cc=20: jt=0: pi=4*atn(1): no=104: nw=96
  dt=1.25*(n/2*n/2)/(4*pi*cc*nf*nk*nk): p2=pi*pi/(n/2*n/2)

[waitHere]
  wait

[p1]
  kk=1: jt=0: plt=0: goto [plot]                            ' P Dens

[p2]
  kk=2: jt=0: plt=0: goto [plot]                            ' R Comp

[p3]
  if kk>0 then                                              ' Motion
    plt=1: goto [plot]
  else
    goto [waitHere]
  end if

[f]
  q=pi/240                                                  ' The upper right...
  k=sqr(k1*k1+.44*e1): k2=sqr(k1*k1+1.076*e1)
  c1=cos(q*k1*240): s1=sin(q*k1*240)
  c=cos(q*k*dl): s=sin(k*q*dl)
  c2=cos(q*k2*(240-dl)): s2=sin(q*k2*(240-dl))
  fb=k*k*s*s1*s2-k1*c1*k*c*s2-k1*c1*k2*c2*s-k*c*k2*c2*s1    ' ...matrix element
return


[findk]
  cnt=0: e=.00001                                           ' Find the k-values
  k1=4: gosub [f]: a=fb                                     '   for which the
  for i=41 to 410                                           '   matrix element = 0
    k1=.1*i: gosub [f]                                      ' Scan for zero crossings
    if a*fb<0 then cnt=cnt+1: k(cnt)=k1                     '   in the range of interest
    a=fb
  next i
  for i=1 to cnt                                            ' 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
  k(0)=cnt
return

[fey]
  c=pi/240: c=c*c: dx=.05                                   ' Find 60 statinary states
  for i=1 to 60
    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=0-(.44*e1+e)*c
      if j=(4800+20*dl) then a=0-(1.0736*e1+e)*c
      v=v+a*z*dx: s=s+z*z: jj=jj+1
      if jj=20 then jj=0: ii=ii+1: y(i, ii)=z
    next j
    s=sqr(s*dx)
    for j=1 to n: y(i, j)=y(i, j)/s: next j                 ' Normalize them
  next i
return

[plot]
  if jt=0 then
    redim yr(n+1): redim yrr(n+1): redim dyr(n+1)
    redim yi(n+1): redim yii(n+1): redim dyi(n+1)
    redim z(n+1): redim yy(n+1): redim vv(n+1)
    redim k(64): redim b(61): redim d(61): redim w(61)
    redim y(61, n+1)
    #1 "cls"
    for i=0 to nw                                           ' Initial wave pulse
      phi=i*2*pi/nw: yr(no+i)=2/sqr(6*nw)*(1-cos(phi))
      yi(no+i)=0-yr(no+i)*sin(phi*nk): yr(no+i)=yr(no+i)*cos(phi*nk)
    next i
    e1=0
    for i=no to no+nw                                       ' Calculate e1
      e1=e1-yr(i)*(yr(i-1)+yr(i+1)-2*yr(i))-yi(i)*(yi(i-1)+yi(i+1)-2*yi(i))
      dyr(i)=dt*(yi(i-1)+yi(i+1)-2*yi(i))                   ' Initialize dy's
      dyi(i)=0-dt*(yr(i-1)+yr(i+1)-2*yr(i))
    next i
    e1=e1/p2
    for i=0 to n: yy(i)=220/kk: z(i)=yy(i): next i
    #1.c1 "selectionindex?": input#1.c1, ind
    dl=5*(ind-1)                                            ' Film thickness
    for i=1 to n-1                                          ' Potential energy
      vv(i)=0
      if i>=n/2 then vv(i)=0-.44*e1
      if i>=n/2+dl then vv(i)=0-1.0736*e1
    next i
    #1.c2 "value?": input#1.c2, z$
    if z$="set" then                                        ' Series
      #1 "place 20 45 ; color black"
      #1 "\Busy"
      gosub [findk]                                         ' Find k(i)
      gosub [fey]                                           ' Find y(i, j)
      for i=1 to 60                                         ' Expansion coefficients
        for j=1 to n-1
          b(i)=b(i)+yr(j)*y(i, j): d(i)=d(i)+yi(j)*y(i, j)
        next j
        w(i)=1.2*pi/64*k(i)*k(i)/nf                         ' Frequencies (including constants)
      next i
    end if
  end if
  tt=time$("ms")
  while jt=0 or plt=1                                       ' Frame loop
    #1 "place 20 45 ; color black"
    #1 "\Frame "; jt; "/320"
    #1 "color blue"
    for i=1 to n
      if i=n/2 then #1 "color green"
      if i=n/2+dl then #1 "color red"
      if kk=1 then                                          ' Prob density
        s=0: if abs(yr(i))>.001 then s=s+yr(i)*yr(i)        ' Underflow kluge
        if abs(yi(i))>.001 then s=s+yi(i)*yi(i)             '   ditto
        yy(i)=220-int(6800*s+.5)
        if yy(i)<z(i) then #1 "line "; i; " "; yy(i); " "; i; " "; z(i)
        if z(i)<yy(i) then #1 "line "; i; " "; z(i); " "; i; " "; yy(i)
      else                                                  ' Real Comp
        yy(i)=110-int(480*yr(i)+.5)
        if jt>0 then #1 "line "; i-1; " "; z(i-1); " "; i; " "; z(i)
        #1 "line "; i-1; " "; yy(i-1); " "; i; " "; yy(i)
      end if
      z(i-1)=yy(i-1)
    next i
    jt=jt+1
    if z$="set" then                                        ' Series
      for i=1 to n
        yr(i)=0: yi(i)=0
        for j=1 to 60
          yr(i)=yr(i)+y(j, i)*(b(j)*cos(w(j)*jt)-d(j)*sin(w(j)*jt))
          yi(i)=yi(i)+y(j, i)*(b(j)*sin(w(j)*jt)+d(j)*cos(w(j)*jt))
        next j
      next i
    else                                                    ' Direct solution
      if jt<nf/3 then nm=3*n/4 else nm=n                    ' Underflow kluge
      for jj=1 to cc                                        ' Accuracy loop
        for i=1 to nm-1                                     ' Project ahead dy/2
          yrr(i)=yr(i)+dyr(i)/2: yii(i)=yi(i)+dyi(i)/2
        next i
        for i=1 to nm-1                                     ' Schrodinger Equation
          dyr(i)=dt*(yii(i-1)+yii(i+1)-2*yii(i)-p2*vv(i)*yii(i))
          yr(i)=yr(i)+dyr(i)
          dyi(i)=0-dt*(yrr(i-1)+yrr(i+1)-2*yrr(i)-p2*vv(i)*yrr(i))
          yi(i)=yi(i)+dyi(i)
        next i
      next jj
    end if
    scan                                                    ' Exit if bored
    if jt=nf+1 then
      s0=0                                                  ' Calculate reflection
      for i=1 to n/2-1                                      '   probability
        s0=s0+yr(i)*yr(i)+yi(i)*yi(i)
      next i
      tt=time$("ms")-tt
      #1 "place 20 65 ; color black"
      #1 "\Reflection Probability = "; using("#.##",100*s0); " %"
      #1 "\Run Time = "; tt/1000; " s"
      plt=0: kk=0
    end if
    #1 "discard"
  wend
goto [waitHere]

[quit]
  close #1
end