REM qm1bbc.bas - jackord@kw.igs.net - 30 Mar 07 - BBC Basic v5.50b
      REM quantum rectilinear pluck and pulse wavefunctions
      REM   and a standing (Pulse S) or travelling (Pulse A) localized sinusoidal pulse
      REM   in a square well with time dependence calculated either directly
      REM   from the Schrodinger Equation or from FFSS expansion
      REM the FFSS expansion runs MUCH faster than the direct solution
      
      *FLOAT64
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize Window
      b1% = FN_button("Pluck", 4, 5, 60, 20, FN_setproc(PROCp1), 0)
      b2% = FN_button("Pulse", 68, 5, 60, 20, FN_setproc(PROCp2), 0)
      b3% = FN_button("PulseS", 132, 5, 60, 20, FN_setproc(PROCp3), 0)
      b4% = FN_button("PulseA", 196, 5, 60, 20, FN_setproc(PROCp4), 0)
      b5% = FN_button("Motion", 260, 5, 60, 20, FN_setproc(PROCp5), 0)
      c1% = FN_combobox("", 324, 5, 58, 120, 0, 3)
      WindowWidth=384
      WindowHeight=260
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "The Quantum String"
      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=192: nm=64: nf=320: c=600: dt=n*n/4/PI/c/nf: plt=0
      OFF: VDU 5
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      QUIT
      
      DEF PROCp1: kk=1: jt=0                         : REM Pluck
      
      DEF PROCp2: kk=2: jt=0                         : REM Pulse
      
      DEF PROCp3: kk=3: jt=0                         : REM PulseA
      
      DEF PROCp4: kk=4: jt=0                         : REM PulseB
      
      DEF PROCp5: 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
        IF kk=1 THEN
          FOR i=0 TO n/2
            yr(i)=30*i/n: yr(n-i)=yr(i)
          NEXT i
          energy=12/PI/PI
        ENDIF
        IF kk=2 THEN
          FOR i=3*n/8 TO n/2
            yr(i)=4*30/n*(i-3*n/8): yr(n-i)=yr(i)
          NEXT i
          energy=16*12/PI/PI
        ENDIF
        IF kk=3 OR kk=4 THEN
          FOR i=3*n/8 TO 5*n/8
            phi=(i-3*n/8)*PI*8/n: yr(i)=7.5*(1-COS(phi))
            IF kk=4 THEN
              yi(i)=-yr(i)*SIN(phi): yr(i)=yr(i)*COS(phi)
            ENDIF
          NEXT i
          energy=21.333: IF kk=4 THEN energy=64+energy
        ENDIF
        IF ch1%=0 THEN
          s=0: se=0
          FOR i=1 TO nm-1
            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
            s=s+b(i)*b(i)+d(i)*d(i): se=se+(b(i)*b(i)+d(i)*d(i))*i*i
          NEXT i
          energy=se/s
        ELSE
          FOR i=1 TO n-1
            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: LINE 0, 0, 768, 0
        MOVE 8, 450: PRINT "Energy = "; INT(1000*energy)/1000
        MOVE 540, 450: PRINT "Frame "; jt; "/320"
        GCOL 12
        FOR i=1 TO n-1
          yy%=2*INT(yr(i)*yr(i)+yi(i)*yi(i))
          IF yy%>0 THEN LINE 4*i, yy%, 4*i, 2
        NEXT i
        IF plt=1 THEN * REFRESH
        jt=jt+1
        IF kk<4 THEN
          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
            NEXT i
          ELSE
            FOR jj=1 TO c                            : REM Accuracy loop
              FOR i=1 TO n/2+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/2                         : REM Schrodinger Equation
                dyr(i)=dt*(yii(i-1)+yii(i+1)-2*yii(i)): dyr(n-i)=dyr(i)
                yr(i)=yr(i)+dyr(i): yr(n-i)=yr(i)
                dyi(i)=-dt*(yrr(i-1)+yrr(i+1)-2*yrr(i)): dyi(n-i)=dyi(i)
                yi(i)=yi(i)+dyi(i): yi(n-i)=yi(i)
              NEXT i
            NEXT jj
          ENDIF
        ELSE
          IF ch1%=0 THEN
            FOR i=1 TO n-1
              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 STEP 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 Project 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
        ENDIF
        IF jt=nf+1 THEN
          plt=0: kk=0
          GCOL 0: MOVE 540, 410: PRINT "Time "; (TIME-tt)/100; " s"
          * REFRESH ON
        ENDIF
      ENDWHILE
      ENDPROC