REM mw7bbc.bas - jackord@kw.igs.net - 23 Apr 07 - BBC BASIC (RTR)
      REM wave incident (k) on a barrier (ik/n) wave transmitted (k)
      REM   (n can have any value, but is set equal to 4 in the example)
      REM the barrier width equals the incident wavelength
      REM plot shows the modulus of the wavefunction and the probability density
      
      *FLOAT 64
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize Window
      button% = FN_button("Plot", 20, 10, 50, 20, FN_setproc(PROCp1), 0)
      WindowWidth=480     : REM Scale 0-480
      WindowHeight=240    : REM Scale 0-240
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "Barrier Tunneling"
      CLS : VDU 5
      GCOL 0
      LINE 480, 0, 480, 480 : LINE 720, 0, 720, 480
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      
      DEF PROCp1
      DIM yy(121)
      dx=PI/20: n=4: b=2*PI/n
      c=(EXP(0-b)+EXP(b))/2                         : REM Amplitudes
      s=(EXP(0-b)-EXP(b))/2                         : REM   and
      f=s*(n-1/n)/2                                 : REM     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
      GCOL 12 : MOVE 0, 2*INT(yy(0)^.5)
      FOR j=1 TO 120
        DRAW 6*j, 2*INT(yy(j)^.5)
      NEXT j
      DRAW 960, 2*INT(t)
      MOVE 760, 360 : PRINT "Probability";
      MOVE 800, 330 : PRINT "Amplitude";
      GCOL 9 : MOVE 0, 2*INT(yy(0)/250)
      FOR j=1 TO 120
        DRAW 6*j, 2*INT(yy(j)/250)
      NEXT j
      DRAW 960, 2*INT(t*t/250)
      MOVE 760, 240 : PRINT "Probability";
      MOVE 800, 210 : PRINT "Density";
      ENDPROC