REM qm3bbc.bas - jackord@kw.igs.net - 3 Apr 07 - BBC Basic v5.50b
      REM quantum wave (Pulse B) incident on a variable thickness region ("film") where
      REM   the kinetic energy increases by a factor 1.44 followed by a region where
      REM   it increases again by a factor of 1.44
      REM the motion is generated either directly from the Schrodinger Equation
      REM   or from an expansion in terms if the stationary states of the system
      
      *FLOAT64
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize Window
      b1% = FN_button("P Dens", 4, 5, 60, 20, FN_setproc(PROCp1), 0)
      b2% = FN_button("R Comp", 68, 5, 60, 20, FN_setproc(PROCp2), 0)
      b3% = FN_button("Motion", 132, 5, 60, 20, FN_setproc(PROCp3), 0)
      c1% = FN_combobox("", 196, 5, 60, 120, 0, 3)
      c2% = FN_combobox("", 260, 5, 100, 120, 0, 3)
      WindowWidth=480
      WindowHeight=220
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "Quantum Pulse Incident on Film"
      SYS "SendMessage", c1%, &143, 0, "Weq"
      SYS "SendMessage", c1%, &143, 1, "Ser"
      SYS "SendMessage", c1%, &14E, 0, 0
      SYS "SendMessage", c2%, &143, 0, "Dfilm=0"
      SYS "SendMessage", c2%, &143, 1, "Dfilm=5"
      SYS "SendMessage", c2%, &143, 2, "Dfilm=10"
      SYS "SendMessage", c2%, &14E, 0, 0
      DIM yr(480): DIM yrr(480): DIM dyr(480): DIM vv(480)
      DIM yi(480): DIM yii(480): DIM dyi(480)
      DIM y(60, 480)
      DIM k(64): DIM b(60): DIM d(60): DIM w(60)
      n=480: nf=320: nk=4: cc=20: no%=104: nw%=96: kk=0: plt=0
      dt=1.25*(n/2*n/2)/(4*PI*cc*nf*nk*nk): p2=PI*PI/(n/2*n/2)
      kk=0: plt=0
      OFF: VDU 5
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      QUIT
      
      DEF PROCp1: kk=1: jt=0                         : REM P Dens
      
      DEF PROCp2: kk=2: jt=0                         : REM R Comp
      
      DEF PROCp3: IF kk>0 THEN plt=1 ELSE ENDPROC    : REM Motion
      
      IF jt=0 THEN
        SYS "SendMessage", c1%, &147, 0, 0 TO ch1%
        SYS "SendMessage", c2%, &147, 0, 0 TO ch2%
        FOR i=0 TO n                                 : REM Initialize
          yr(i)=0: yi(i)=0: yrr(i)=0: yii(i)=0: dyr(i)=0: dyi(i)=0
        NEXT i
        FOR i=0 TO nw%                               : REM Initial wave pulse
          phi=i*2*PI/nw%: yr(no%+i)=2/SQR(6*nw%)*(1-COS(phi))
          yi(no%+i)=-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%                         : REM 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))        : REM Initialize dy's
          dyi(i)=-dt*(yr(i-1)+yr(i+1)-2*yr(i))
        NEXT i
        e1=e1/p2
        dl=5*ch2%                                    : REM Film thickness
        FOR i=0 TO n                                 : REM Potential energy
          vv(i)=0
          IF i>=n/2 THEN vv(i)=-.44*e1
          IF i>=n/2+dl THEN vv(i)=-1.0736*e1
        NEXT i
        IF ch1%=1 THEN
          CLS: GCOL 0: MOVE 40, 340: PRINT "Busy"
          PROCfindk                                  : REM Find k(i)
          PROCfey                                    : REM Find y(i, j)
          FOR i=1 TO 60                              : REM 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              : REM Frequencies (including constants)
          NEXT i
        ENDIF
      ENDIF
      IF plt=1 THEN * REFRESH OFF
      tt=TIME
      WHILE jt=0 OR plt=1
        CLS: GCOL 0: IF kk=1 THEN LINE 0, 0, 960, 0
        MOVE 740, 420: PRINT "Frame "; jt; "/320"
        GCOL 12: MOVE 0, 200
        FOR i=1 TO n
          IF i=n/2 THEN GCOL 10
          IF i=n/2+dl THEN GCOL 9
          IF kk=1 THEN
            yy%=2*INT(6800*(yr(i)*yr(i)+yi(i)*yi(i)))
            IF yy%>0 THEN LINE 2*i, yy%, 2*i, 2
          ELSE
            DRAW 2*i, 200+INT(960*yr(i)+1)
          ENDIF
        NEXT i
        IF plt=1 THEN * REFRESH
        jt=jt+1
        IF ch1%=1 THEN
          FOR i=1 TO n                               : REM FFSS
            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
          FOR jj=1 TO cc                             : REM Accuracy loop
            FOR i=0 TO n                             : REM Look ahead dy/2
              yrr(i)=yr(i)+dyr(i)/2: yii(i)=yi(i)+dyi(i)/2
            NEXT i
            FOR i=1 TO n-1                           : REM 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)=-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
        ENDIF
        IF jt=nf+1 THEN
          plt=0: kk=0: s0=0
          FOR i=0 TO n/2-1
            s0=s0+yr(i)*yr(i)+yi(i)*yi(i)
          NEXT i
          GCOL 0: MOVE 8, 360
          PRINT "Reflection Probability = "; INT(10000*s0)/100; " %"
          MOVE 8, 320
          PRINT "Time "; (TIME-tt)/100; " s"
          * REFRESH ON
        ENDIF
      ENDWHILE
      ENDPROC
      
      DEF PROCf
      q=PI/240                                       : REM The upper right element
      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
      ENDPROC
      
      DEF PROCfindk
      cnt=0: e=.00001                                : REM Find the k-values
      k1=4: PROCf: a=fb                              : REM   for which the
      FOR i=41 TO 410                                : REM   matrix element = 0
        k1=.1*i: PROCf                               : REM Scan for zero crossings
        IF a*fb<0 THEN cnt=cnt+1: k(cnt)=k1          : REM   in the range of interest
        a=fb
      NEXT i
      FOR i=1 TO cnt                                 : REM Use a binary search
        a=k(i)-.1: b=k(i): k1=a: PROCf: fa=fb        : REM   to refine the values
        WHILE (b-a)>e
          k1=(a+b)/2: PROCf: IF fb*fa>0 THEN a=k1 ELSE b=k1
        ENDWHILE
        k(i)=(a+b)/2
      NEXT i
      k(0)=cnt
      ENDPROC
      
      DEF PROCfey
      c=PI/240: c=c*c: dx=.05: s=0                   : REM Find 60 statinary states
      FOR i=1 TO 60
        e=k(i)*k(i): a=-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=-(.44*e1+e)*c
          IF j=(4800+20*dl) THEN a=-(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    : REM Normalize them
      NEXT i
      ENDPROC