' mw7lib.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 the modulus of the wavefunction and the probability density

' Initialize Window
  nomainwin
  button#1, "Plot", [p1], UL, 20, 10, 50, 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]"
  #1 "cls ; down ; color black"
  #1 "line 240 0 240 240 ; line 360 0 360 240"

[waitHere]
  wait

[p1]
  goto [plt]

[plt]
  dim yy(121)
  pi=4*atn(1): dx=pi/20: 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=100: r1=a1*r1: a2=a1*a2: r2=a1*r2: t=a1*t
  for j=0 to 80
    yy(j)=a1*a1+r1*r1+2*a1*r1*cos(2*j*dx+phase2)
  next j
  for j=81 to 120
    a=a2*exp((j-120)*dx/n): r=r2*exp((120-j)*dx/n)
    yy(j)=a*a+r*r+2*a2*r2*cos(phase4-phase3)
  next j
  #1 "color blue ; place 0 "; 240-int(yy(0)^.5)
  for j=1 to 120
    #1 "goto "; 3*j; " "; 240-int(yy(j)^.5)
  next j
  #1 "goto 480 "; 240-int(t)
  #1 "place 380 60": #1 "\Probability"
  #1 "place 400 75": #1 "\Amplitude"
  #1 "color red ; place 0 "; 240-int(yy(0)/250)
  for j=1 to 120
    #1 "goto "; 3*j; " "; 240-int(yy(j)/250)
  next j
  #1 "goto 480 "; 240-int(t*t/250)
  #1 "place 380 120": #1 "\Probability"
  #1 "place 400 135": #1 "\Density"
goto [waitHere]

[quit]
  close #1
end