REM qm5bbc.bas - jackord@kw.igs.net - revised 5 Mar 08 - BBC Basic v5.50b
      REM quantum harmonic oscillator
      REM   Pulse A from the quantum string placed between the ground state
      REM     turning points of the quantum oscillator
      REM   direct solution of the time-dependent Schrodinger Equation or
      REM     expansion in terms of stationary states (a 54-term series)
      REM classical mass/spring system motion given by Feynman algorithm
      REM differs from Java version: does not use expanded range for the
      REM   direct solution and hence has boundary reflections
      REM the series solution runs faster than the direct solution
      
      *FLOAT64
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize Window
      b1% = FN_button("Schr Eqn", 2, 2, 65, 20, FN_setproc(PROCp1), 0)
      b2% = FN_button("Series", 2, 27, 65, 20, FN_setproc(PROCp2), 0)
      WindowWidth=384
      WindowHeight=280
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "Quantum Oscillator"
      VDU 5
      n=24: n2=n*n: np=384: nc=np/2: nf=320: ns=54: sf=4000: sy=500: yo=80: xc=384
      c=160: dt=n2*PI/c/nf: du=PI/c/nf/n2: w=2*PI/nf: plt=0
      DIM y(54, 384): DIM xt%(120): DIM yt%(120)
      DIM yr(384): DIM yrr(384): DIM dyr(384)
      DIM yi(384): DIM yii(384): DIM dyi(384)
      DIM ar(54): DIM ai(54)
      ss=1: dx=0.1
      FOR k=1 TO ns                                  : REM Initialize y(k, i)...
        e=2*k-1: z=1+ss: y(k, nc)=z: v=1-ss
        s=z*z/2: a=-e*z/n2: v=v+a*dx/2
        x=0: jj=0: ii=0: tst=0: tp=n*SQR(e)
        WHILE tst=0
          jj=jj+1: x=x+dx: z=z+v*dx: a=(x*x/n2-e)*z/n2: v=v+a*dx
          IF x>tp AND v*z>0 THEN z=0: tst=1
          s=s+z*z
          IF jj=10 THEN
            jj=0: ii=ii+1: IF ii<=nc THEN y(k, nc+ii)=z
          ENDIF
        ENDWHILE
        s=SQR(2*s*dx)
        FOR i=0 TO nc                                : REM ...and normalize them
          y(k, nc+i)=y(k, nc+i)/s: y(k, nc-i)=ss*y(k, nc+i)
        NEXT i
        ss=0-ss
      NEXT k
      FOR i=0 TO 60                                  : REM Potential function
        yt%(60+i)=yo+8*i: yt%(60-i)=yt%(60+i)
        xt%(60+i)=2*nc+INT(2*n*SQR(4*i/15)): xt%(60-i)=2*np-xt%(60+i)
      NEXT i
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      QUIT
      
      DEF PROCp1: kk=1                               : REM Schr Eqn
      
      DEF PROCp2: kk=2                               : REM Series
      
      FOR jt=0 TO nf
        IF jt=0 THEN
          FOR i=0 TO n
            yr(i)=0: yrr(i)=0: dyr(i)=0
            yi(i)=0: yii(i)=0: dyi(i)=0
          NEXT i
          FOR i=nc-n TO nc+n                         : REM Init yr(i), yi(i)
            phi=(i-(nc-n))*PI/n: yr(i)=(1-COS(phi))/SQR(3*n)
            yi(i)=-yr(i)*SIN(phi): yr(i)=yr(i)*COS(phi)
          NEXT i
          FOR i=nc-n TO nc+n                         : REM Init dyr(i), dyi(i)
            dyr(i)=dt*(yi(i-1)+yi(i+1)-2*yi(i))-du*(i-nc)*(i-nc)*yi(i)
            dyi(i)=-dt*(yr(i-1)+yr(i+1)-2*yr(i))+du*(i-nc)*(i-nc)*yr(i)
          NEXT i
          enf=0                                      : REM Init enf
          FOR i=nc-n TO nc+n
            s1=yr(i)*(-n2*(yr(i-1)+yr(i+1)-2*yr(i))+(i-nc)*(i-nc)/n2*yr(i))
            s2=yi(i)*(-n2*(yi(i-1)+yi(i+1)-2*yi(i))+(i-nc)*(i-nc)/n2*yi(i))
            enf=enf+s1+s2
          NEXT i
          ens=0
          FOR k=1 TO ns-1 STEP 2                     : REM Init ar(k)
            ar(k)=0
            FOR i=nc-n TO nc+n
              ar(k)=ar(k)+yr(i)*y(k, i)
            NEXT i
            ens=ens+ar(k)*ar(k)*(2*k-1)
          NEXT k
          FOR k=2 TO ns STEP 2                       : REM Init ai(k)
            ai(k)=0
            FOR i=nc-n TO nc+n
              ai(k)=ai(k)+yi(i)*y(k, i)
            NEXT i
            ens=ens+ai(k)*ai(k)*(2*k-1)              : REM Init ens
          NEXT k
          IF kk=2 THEN energy=ens ELSE energy=enf
          xtp=2*n*SQR(energy): xs=xc: vx=xtp
          yen=yo+INT(30*energy)
          FOR jj=1 TO 2
            CLS
            IF kk=2 THEN
              FOR i=0 TO np                          : REM Series
                y(0, i)=0
                FOR k=1 TO ns
                  IF jj=1 THEN s1=ar(k) ELSE s1=ai(k)
                  y(0, i)=y(0, i)+s1*y(k, i)
                NEXT k
              NEXT i
              GCOL 10
              FOR k=1 TO ns                          : REM Components
                MOVE 0, 260
                FOR i=1 TO np
                  IF jj=1 THEN s1=ar(k)*y(k, i) ELSE s1=ai(k)*y(k, i)
                  DRAW 2*i, 260+2*INT(sy*s1+.5)
                NEXT i
              NEXT k
            ELSE
              FOR i=0 TO np
                IF jj=1 THEN s1=yr(i) ELSE s1=yi(i)
                y(0, i)=s1
              NEXT i
            ENDIF
            GCOL 12: MOVE 0, 260
            FOR i=1 TO np
              DRAW 2*i, 260+2*INT(sy*y(0, i))
            NEXT i
            MOVE 280-20*jj, 460
            IF jj=1 THEN PRINT "REAL COMPONENT" ELSE PRINT "IMAGINARY COMPONENT"
            WAIT 200
          NEXT jj
        ENDIF
        CLS: GCOL 0: LINE 0, 0, 768, 0: GCOL 12
        FOR i=1 TO np-1
          IF kk=2 THEN
            re=0: im=0                               : REM Series probability density
            FOR k=1 TO ns
              re=re+y(k, i)*(ar(k)*COS((k-1)*jt*w)-ai(k)*SIN((k-1)*jt*w))
              im=im+y(k, i)*(ar(k)*SIN((k-1)*jt*w)+ai(k)*COS((k-1)*jt*w))
            NEXT k
            yy%=yo+2*INT(sf*(re*re+im*im))
          ELSE
            yy%=yo+2*INT(sf*(yr(i)*yr(i)+yi(i)*yi(i)))  : REM W eqn probability density
          ENDIF
          IF yy%>yo THEN LINE 2*i, yy%, 2*i, yo+2
        NEXT i
        x2=INT(xs)-20: RECTANGLEFILL x2, 20, 40,40   : REM Plot mass...
        GCOL 0: MOVE 0, 40
        FOR i=2 TO x2 STEP 2                         : REM ...and spring
          y2=40+INT(20*SIN(16*PI*i/x2)): DRAW i, y2
        NEXT i
        GCOL 9: MOVE xt%(0), yt%(0)
        FOR i=1 TO 120                               : REM Plot potential and...
          DRAW xt%(i), yt%(i)
        NEXT i                                       : REM ...gs energy and turning points
        LINE 336, yo+30, 432, yo+30: LINE 432, yo+30, 432, yo: LINE 336, yo+30, 336, yo
        GCOL 10: LINE np-xtp, yen, np+xtp, yen
        LINE np+xtp, yen, np+xtp, yo
        LINE np-xtp, yen, np-xtp, yo
        GCOL 0: MOVE 10, yo+365: PRINT "Frm "; jt; "/320"
        MOVE 570, yen+10: PRINT "E/gs="; INT(1000*energy)/1000
        LINE 0, 80, 768, 80
        IF jt=0 THEN
          WAIT 200: tt=TIME: * REFRESH OFF
        ELSE
          * REFRESH
        ENDIF
        xs=xs+vx*w: vx=vx+(xc-xs)*w                  : REM Update mass on spring
        IF kk=1 THEN
          FOR jj=0 TO c
            FOR i=1 TO np-1                          : REM Project ahead
              yrr(i)=yr(i)+dyr(i)/2: yii(i)=yi(i)+dyi(i)/2
            NEXT i
            FOR i=1 TO np-1                          : REM Schrodinger equation
              dyr(i)=dt*(yii(i-1)+yii(i+1)-2*yii(i))-du*(i-nc)*(i-nc)*yii(i)
              yr(i)=yr(i)+dyr(i)
              dyi(i)=-dt*(yrr(i-1)+yrr(i+1)-2*yrr(i))+du*(i-nc)*(i-nc)*yrr(i)
              yi(i)=yi(i)+dyi(i)
            NEXT i
          NEXT jj
        ENDIF
      NEXT jt
      GCOL 0: MOVE 10, yo+330: PRINT "Time "; INT((TIME-tt)/100); " s"
      * REFRESH ON
      ENDPROC