' mw6lib.bas - jackord@kw.igs.net - revised 13 Mar 08 - Liberty Basic v4.02
' wave incident (k) on a barrier (ik/n) wave transmitted (k)
'   (n can have any value, but is set equal to 4 in the example)
' the barrier width equals the incident wavelength
' plot shows ONE (real or imaginary) component of the wavefunction

' Initialize Window
  nomainwin
  button#1, "Component", [p1], UL, 20, 10, 80, 20
  button#1, "Resultant", [p2], UL, 120, 10, 80, 20
  WindowWidth=498     ' Scale 0-480
  WindowHeight=278    ' Scale 0-240
  UpperLeftX=10: UpperLeftY=10
  open "Barrier Tunneling" for graphics_nsb as #1
  #1 "trapclose [quit]"

[waitHere]
  wait

[p1]
  kk=1: goto [motion]
[p2]
  kk=2: goto [motion]

[motion]
  #1 "cls ; down ; color black"
  #1 "line 0 120 480 120 ; line 240 0 240 240 ; line 360 0 360 240"
  #1 "rule XOR"
  dim yi(81): dim yr(81): dim ys(81)
  dim ya(41): dim yb(41): dim yc(41): dim yt(41)
  dim zi(81): dim zr(81): dim zs(81)
  dim za(41): dim zb(41): dim zc(41): dim zt(41)
  pi=4*atn(1): dx=pi/20: dt=dx/2: n=4: b=2*pi/n
  c=(exp(0-b)+exp(b))/2                                ' Amplitudes
  s=(exp(0-b)-exp(b))/2                                '   and
  f=s*(n-1/n)/2                                        '     phases
  t=1/(c*c+f*f)^.5: phase1=atn(f/c)
  r1=(n+1/n)*s*t/2: phase2=phase1-pi/2
  a2=(1+n*n)^.5*t/2: r2=a2
  phase4=phase1+atn(n)
  phase3=phase1-atn(n)
  a1=60: r1=a1*r1: a2=a1*a2: r2=a1*r2: t=a1*t
  for i=0 to 320
    tt=i*dt
    for j=0 to 80
      yi(j)=120-int(a1*sin(tt-j*dx))
      yr(j)=120-int(r1*sin(tt+j*dx+phase2))
      ys(j)=yi(j)+yr(j)-120
    next j
    for j=0 to 40
      ya(j)=120-int(a2*sin(tt+phase3)*exp((j-40)*dx/n))
      yb(j)=120-int(r2*sin(tt+phase4)*exp((40-j)*dx/n))
      yc(j)=ya(j)+yb(j)-120
    next j
    for j=0 to 40
      yt(j)=120-int(t*sin(tt-j*dx+phase1))
    next j
    if kk=1 then
      #1 "color darkpink"
      for j=1 to 80
        if i>0 then #1 "line "; 3*j-3; " "; zi(j-1); " "; 3*j; " "; zi(j)
        #1 "line "; 3*j-3; " "; yi(j-1); " "; 3*j; " "; yi(j)
      next j
      #1 "color red"
      for j=1 to 80
        if i>0 then #1 "line "; 3*j-3; " "; zr(j-1); " "; 3*j; " "; zr(j)
        #1 "line "; 3*j-3; " "; yr(j-1); " "; 3*j; " "; yr(j)
      next j
      #1 "color yellow"
      for j=1 to 40
        if i>0 then #1 "line "; 240+3*j-3; " "; za(j-1); " "; 240+3*j; " "; za(j)
        #1 "line "; 240+3*j-3; " "; ya(j-1); " "; 240+3*j; " "; ya(j)
      next j
      #1 "color green"
      for j=1 to 40
        if i>0 then #1 "line "; 240+3*j-3; " "; zb(j-1); " "; 240+3*j; " "; zb(j)
        #1 "line "; 240+3*j-3; " "; yb(j-1); " "; 240+3*j; " "; yb(j)
      next j
    else
      #1 "color blue"
      for j=1 to 80
        if i>0 then #1 "line "; 3*j-3; " "; zs(j-1); " "; 3*j; " "; zs(j)
        #1 "line "; 3*j-3; " "; ys(j-1); " "; 3*j; " "; ys(j)
      next j
      for j=1 to 40
        if i>0 then #1 "line "; 240+3*j-3; " "; zc(j-1); " "; 240+3*j; " "; zc(j)
        #1 "line "; 240+3*j-3; " "; yc(j-1); " "; 240+3*j; " "; yc(j)
      next j
    end if
    #1 "color blue"
    for j=1 to 40
      if i>0 then #1 "line "; 360+3*j-3; " "; zt(j-1); " "; 360+3*j; " "; zt(j)
      #1 "line "; 360+3*j-3; " "; yt(j-1); " "; 360+3*j; " "; yt(j)
    next j
    for j=0 to 80
      zi(j)=yi(j): zr(j)=yr(j): zs(j)=ys(j)
    next j
    for j=0 to 40
      za(j)=ya(j): zb(j)=yb(j): zc(j)=yc(j): zt(j)=yt(j)
    next j
  next i
goto [waitHere]

[quit]
  close #1
end