REM qm4bbc.bas - jackord@kw.igs.net - 4 Apr 07 - BBC Basic v5.50b
      REM superposition of the two lowest energy states of a double well (an
      REM   infinite square well with a central barrier)
      
      *FLOAT64
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize Window
      b1% = FN_button("P Dens", 4, 5, 60, 20, FN_setproc(PROCp1), 0)
      b2% = FN_button("Yr  Yi", 68, 5, 60, 20, FN_setproc(PROCp2), 0)
      b3% = FN_button("Motion", 132, 5, 60, 20, FN_setproc(PROCp3), 0)
      WindowWidth=400
      WindowHeight=300
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "Double Well"
      DIM ya(400): DIM yb(400)
      n=400: c=PI/n: c=c*c: ea=3.85297: eb=4.07503
      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 Yr Yi
      
      DEF PROCp3: IF kk>0 THEN plt=1 ELSE ENDPROC    : REM Motion
      
      IF jt=0 THEN
        FOR i=0 TO n: ya(i)=0: yb(i)=0: NEXT i       : REM Initialize
        e=ea: v=0.9: sa=0
        FOR i=1 TO n-1
          ya(i)=ya(i-1)+v: sa=sa+ya(i)*ya(i)
          IF i>191 AND i<209 THEN a=(200-e)*c ELSE a=-e*c
          v=v+a*ya(i)
        NEXT i
        e=eb: v=0.9: sb=0
        FOR i=1 TO n-1
          yb(i)=yb(i-1)+v: sb=sb+yb(i)*yb(i)
          IF i>191 AND i<209 THEN a=(200-e)*c ELSE a=-e*c
          v=v+a*yb(i)
        NEXT i
        sb=SQR(sa/sb)
        FOR i=1 TO n-1: yb(i)=sb*yb(i): NEXT i
        IF kk=1 THEN nf=320 ELSE nf=960
        a=2*PI/(eb-ea)/nf
      ENDIF
      IF plt=1 THEN * REFRESH OFF
      ENDIF
      WHILE jt=0 OR plt=1
      CLS: GCOL 0
      MOVE 580, 580: PRINT "Frame "; jt; "/"; nf
      GCOL 10: RECTANGLE FILL 384, 0, 20, 600
      IF kk=1 THEN
        GCOL 0: LINE 0, 0, 800, 0: GCOL 12
        FOR i=1 TO n-1                               : REM Probability density
          yr=ya(i)*COS(jt*ea*a)+yb(i)*COS(jt*eb*a)
          yi=ya(i)*SIN(jt*ea*a)+yb(i)*SIN(jt*eb*a)
          yy%=2*INT(.019*(yr*yr+yi*yi))
          IF yy%>0 THEN LINE 2*i, yy%, 2*i, 2
        NEXT i
      ELSE
        GCOL 12: MOVE 0, 300
        FOR i=1 TO n                                 : REM Yr
          yr=ya(i)*COS(jt*ea*a)+yb(i)*COS(jt*eb*a)
          DRAW 2*i, 300+2*INT(yr)
        NEXT i
        GCOL 9: MOVE 0, 300
        FOR i=1 TO n                                 : REM Yi
          yi=ya(i)*SIN(jt*ea*a)+yb(i)*SIN(jt*eb*a)
          DRAW 2*i, 300+2*INT(yi)
        NEXT i
      ENDIF
      IF plt=1 THEN * REFRESH
      jt=jt+1
      IF jt=nf+1 THEN plt=0: kk=0: * REFRESH ON
      ENDWHILE
      ENDPROC