REM mw6bbc.bas - jackord@kw.igs.net - revised 31 Dec 07 - BBC BASIC 5.50b
      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 the plot shows ONE (real or imaginary) component of the wavefunction
      
      *FLOAT64
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize Window
      button1% = FN_button("Component", 10, 5, 80, 20, FN_setproc(PROCp1), 0)
      button2% = FN_button("Resultant", 110, 5, 80, 20, FN_setproc(PROCp2), 0)
      WindowWidth=480                                : REM Scale 0-960
      WindowHeight=240                               : REM Scale 0-480
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "Quantum wave striking barrier"
      OFF
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      
      DEF PROCp1 : kk=1
      DEF PROCp2 : kk=2
      
      ORIGIN 0, 240
      jd=120: n=4: b=2*PI/n
      c=(EXP(-b)+EXP(b))/2                           : REM Amplitudes
      s=(EXP(-b)-EXP(b))/2                           : REM   and
      f=s*(n-1/n)/2                                  : REM     phases
      t=1/(c*c+f*f)^.5: ph1=ATN(f/c)
      r1=(n+1/n)*s*t/2: ph2=ph1-PI/2
      a2=(1+n*n)^.5*t/2: r2=a2
      ph4=ph1+ATN(n)
      ph3=ph1-ATN(n)
      a1=120: r1=a1*r1: a2=a1*a2: r2=a1*r2: t=a1*t
      b=PI/120: c=PI/60
      tt=TIME: del=1.5
      * REFRESH OFF
      FOR i=0 TO 960
        CLS: GCOL 0
        LINE 0, 0, 960, 0 : LINE 480, -240, 480, 240
        LINE 720, -240, 720, 240
        IF kk=1 THEN
          GCOL 9: MOVE 0, INT(a1*SIN(i*b))           : REM Plot components
          FOR j=1 TO 240
            DRAW 2*j, INT(a1*SIN(i*b-j*c))
          NEXT j
          GCOL 6: MOVE 0, INT(r1*SIN(i*b+ph2))
          FOR j=1 TO 240
            DRAW 2*j, INT(r1*SIN(i*b+j*c+ph2))
          NEXT j
          GCOL 13: MOVE 480, INT(a2*SIN(i*b+ph3)*EXP(-jd*c/n))
          FOR j=1 TO jd
            DRAW 480+2*j, INT(a2*SIN(i*b+ph3)*EXP((j-jd)*c/n))
          NEXT j
          GCOL 10: MOVE 480, INT(r2*SIN(i*b+ph4)*EXP(jd*c/n))
          FOR j=1 TO jd
            DRAW 480+2*j, INT(r2*SIN(i*b+ph4)*EXP((jd-j)*c/n))
          NEXT j
        ELSE
          GCOL 12: MOVE 0, INT(a1*SIN(i*b)+r1*SIN(i*b+ph2))
          FOR j=1 TO 240                             : REM Plot resultant
            DRAW 2*j, INT(a1*SIN(i*b-j*c)+r1*SIN(i*b+j*c+ph2))
          NEXT j
          FOR j=1 TO jd
            DRAW 480+2*j, INT(a2*SIN(i*b+ph3)*EXP((j-jd)*c/n)+r2*SIN(i*b+ph4)*EXP((jd-j)*c/n))
          NEXT j
        ENDIF
        GCOL 12: MOVE 720, INT(t*SIN(i*b+ph1))       : REM Plot transmitted wave
        FOR j=jd+1 TO 240
          DRAW 480+2*j, INT(t*SIN(i*b-(j-jd)*c+ph1))
        NEXT j
        tt=tt+del: WHILE tt>TIME: ENDWHILE
        * REFRESH
      NEXT i
      * REFRESH ON
      ENDPROC