REM qm2bbc.bas - jackord@kw.igs.net - 1 Apr 07 - BBC Basic v5.50b
      REM quantum wave (Pulse B) in a square well with time dependence calculated
      REM   either directly from the Schrodinger Equation or from FFSS expansion
      
      *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, 58, 120, 0, 3)
      WindowWidth=320
      WindowHeight=220
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "Quantum Pulse"
      SYS "SendMessage", c1%, &143, 0, "FFSS"
      SYS "SendMessage", c1%, &143, 1, "Weq"
      SYS "SendMessage", c1%, &14E, 0, 0
      DIM yr(192): DIM yrr(192): DIM dyr(192)
      DIM yi(192): DIM yii(192): DIM dyi(192)
      DIM b(64): DIM d(64): DIM w(64)
      n=160: nm=64: nf=320: nk=4: c=8: dt=n*n/4/PI/c/nf/nk/nk
      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%
        FOR i=0 TO n: yr(i)=0: yi(i)=0: NEXT i       : REM Initialize
        FOR i=3*n/8 TO 5*n/8
          phi=(i-3*n/8)*PI/(n/8): yr(i)=(1-COS(phi))/2
          yi(i)=-yr(i)*SIN(phi*nk): yr(i)=yr(i)*COS(phi*nk)
        NEXT i
        IF ch1%=0 THEN
          FOR i=1 TO nm-1                            : REM FFSS
            b(i)=0: d(i)=0
            FOR j=1 TO n-1
              sn=SIN(PI*i*j/n)
              b(i)=b(i)+yr(j)*2/n*sn: d(i)=d(i)+yi(j)*2/n*sn
            NEXT j
            w(i)=i*i*2*PI/8/nf/16
          NEXT i
        ELSE
          FOR i=1 TO n-1                             : REM Look ahead
            dyr(i)=dt*(yi(i-1)+yi(i+1)-2*yi(i))
            dyi(i)=-dt*(yr(i-1)+yr(i+1)-2*yr(i))
          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, 640, 0
        MOVE 8, 380: PRINT "Frame "; jt; "/320"
        GCOL 12: MOVE 0, 200
        FOR i=1 TO n-1
          IF kk=1 THEN
            yy%=2*INT(160*(yr(i)*yr(i)+yi(i)*yi(i)))
            IF yy%>0 THEN LINE 4*i, yy%, 4*i, 2
          ELSE
            DRAW 4*i, 200+INT(160*yr(i)+1)
          ENDIF
        NEXT i
        IF plt=1 THEN * REFRESH
        jt=jt+1
        IF ch1%=0 THEN
          FOR i=1 TO n-1                           : REM FFSS
            yr(i)=0: yi(i)=0
            FOR k=1 TO nm-1 STEP 2
              sn=SIN(PI*i*k/n)
              yr(i)=yr(i)+sn*b(k)*COS(w(k)*jt)
              yi(i)=yi(i)+sn*b(k)*SIN(w(k)*jt)
            NEXT k
            FOR k=2 TO nm-2
              sn=SIN(PI*i*k/n)
              yr(i)=yr(i)-sn*d(k)*SIN(w(k)*jt)
              yi(i)=yi(i)+sn*d(k)*COS(w(k)*jt)
            NEXT k
          NEXT i
        ELSE
          FOR jj=1 TO c                            : REM Accuracy loop
            FOR i=1 TO n-1                         : 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))
              yr(i)=yr(i)+dyr(i)
              dyi(i)=-dt*(yrr(i-1)+yrr(i+1)-2*yrr(i))
              yi(i)=yi(i)+dyi(i)
            NEXT i
          NEXT jj
        ENDIF
        IF jt=nf+1 THEN
          plt=0: kk=0
          GCOL 0: MOVE 8, 340: PRINT "Time "; (TIME-tt)/100; " s"
          * REFRESH ON
        ENDIF
      ENDWHILE
      ENDPROC