REM qm5b.bas - jackord@kw.igs.net - revised 12 Mar 2014 - BBC Basic v5.94a REM Speed enhancement by richard@rtrussell.co.uk (the author of BBC Basic) 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 40-term series) REM with dispersion relation plotted REM Plot either probability density (Yr*Yr+Yi*Yi) or variables Yr and Yi REM Classical mass/spring system motion given by Feynman algorithm *FLOAT64 INSTALL @lib\$+"WINLIB5" REM Initialize Window b1% = FN_button("Motion", 2, 2, 70, 25, FN_setproc(PROCp1), 0) c1% = FN_combobox("", 2, 32, 70, 20, 0, 3) c2% = FN_combobox("", 2, 62, 70, 20, 0, 3) WindowWidth=384 WindowHeight=280 VDU 23,22,WindowWidth;WindowHeight;8,15,16,128 SYS "SetWindowText", @hwnd%, "Quantum Oscillator" SYS "SendMessage", c1%, &143, 0, "Sch Eq" SYS "SendMessage", c1%, &143, 1, "Series" SYS "SendMessage", c1%, &14E, 0, 0 SYS "SendMessage", c2%, &143, 0, "Pr Den" SYS "SendMessage", c2%, &143, 1, "Yr Yi" SYS "SendMessage", c2%, &14E, 0, 0 OFF: VDU 5 n=24: n2=n*n: ns=40 nq=576: nc=nq/2: np=384: nn=np/2: d=nc-nn nf=320: c=160: 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(ns, np), zr%(np), zi%(np), xt%(120), yt%(120) DIM yr(nq), yrr(nq), dyr(nq), dr(nq) DIM yi(nq), yii(nq), dyi(nq), di(nq) DIM ar(ns), ai(ns), au(ns), av(ns), iy(ns), re(ns), im(ns), ev(ns) FOR i=1 TO nq-1 dr(i)=-2*dt-du*(i-nc)*(i-nc) di(i)= 2*dt+du*(i-nc)*(i-nc) NEXT i PROCpotential PROCstates REPEAT WAIT 1 UNTIL FALSE QUIT DEF PROCp1 SYS "SendMessage", c1%, &147, 0, 0 TO ch1% SYS "SendMessage", c2%, &147, 0, 0 TO ch2% FOR jt=0 TO nf IF jt=0 THEN FOR i=0 TO nq 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: s3=0 REM Init ens, s3 FOR k=1 TO ns 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-d) NEXT i ens=ens+ar(k)*ar(k)*ev(k): s3=s3+ar(k)*ar(k) 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-d) NEXT i ens=ens+ai(k)*ai(k)*ev(k): s3=s3+ai(k)*ai(k) NEXT k IF ch1%=1 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 ch1%=1 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, 320 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, 320+2*INT(sy*s1+.5) NEXT i NEXT k ELSE FOR i=0 TO np IF jj=1 THEN s1=yr(i+d) ELSE s1=yi(i+d) y(0, i)=s1 NEXT i ENDIF GCOL 12: MOVE 0, 320 FOR i=1 TO np DRAW 2*i, 320+2*INT(sy*y(0, i)) NEXT i MOVE 280-20*jj, 500 IF jj=1 THEN PRINT "REAL COMPONENT" ELSE PRINT "IMAGINARY COMPONENT" WAIT 200 NEXT jj ENDIF FOR k=1 TO ns t=(k-1)*jt*w au(k)=ar(k)*COS(t)-ai(k)*SIN(t) av(k)=ar(k)*SIN(t)+ai(k)*COS(t) NEXT k CLS: GCOL 0: LINE 0, 0, 768, 0: GCOL 12 FOR i=0 TO np IF ch1%=1 THEN FOR K%=1 TO ns: iy(K%)=y(K%, i): NEXT REM Fast code for series solution re()=iy()*au(): re=SUM(re()) im()=iy()*av(): im=SUM(im()) yy%=yo+2*INT(sf*(re*re+im*im)) zr%(i)=320+2*INT(sy*re): zi%(i)=320+2*INT(sy*im) ELSE yy%=yo+2*INT(sf*(yr(i+d)*yr(i+d)+yi(i+d)*yi(i+d))) zr%(i)=320+2*INT(sy*yr(i+d)): zi%(i)=320+2*INT(sy*yi(i+d)) ENDIF IF (ch2%=0) AND (yy%>yo) THEN LINE 2*i, yy%, 2*i, yo+2 NEXT i IF ch2%=1 THEN MOVE 0, zr%(0) FOR i=1 TO np DRAW 2*i, zr%(i) NEXT i GCOL 9: MOVE 0, zi%(0) FOR i=1 TO np DRAW 2*i, zi%(i) NEXT i ENDIF GCOL 12 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 580, 540: 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 ch1%=0 THEN FOR jj=0 TO c yrr()=dyr()/2+yr(): yii()=dyi()/2+yi() FOR I%=1 TO nq-1 REM Schrodinger equation dyr(I%)=yii(I%-1)+yii(I%+1) dyi(I%)=yrr(I%-1)+yrr(I%+1) NEXT dyr()*=dt: dyi()*=-dt REM Fast code for Schrodinger Eqn dyr()+=dr()*yii(): dyi()+=di()*yrr() yr()+=dyr(): yi()+=dyi() NEXT jj ENDIF NEXT jt GCOL 0: MOVE 10, yo+330: PRINT "Time "; INT((TIME-tt)/100); " s" IF ch1%=1 THEN GCOL 15: RECTANGLEFILL 520, 10, 240 REM Plot Dispersion Curve GCOL 0: RECTANGLE 520, 10, 240 MOVE 550, 220: PRINT "S"; ns; "="; INT(s3*100000)/1000; "%" MOVE 528, 246: PRINT "Dispersion Rel" GCOL 9 FOR i=1 TO 30: y1=INT(500*(ar(i)*ar(i)+ai(i)*ai(i))) IF y1>0 THEN RECTANGLEFILL 520+8*i-2, 10, 6, 2*y1 NEXT i GCOL 12: MOVE 520, 10 FOR i=1 TO 30 DRAW 520+8*i, 10+INT(4*ev(i)) NEXT i ENDIF * REFRESH ON ENDPROC DEF PROCstates e=0: de=.00001: b=.5: dev=1.0: ss=1 REM Find ev(k) FOR k=1 TO ns cf=b+dev: e=b: PROCfey: yb=y: e=cf: PROCfey: yc=y WHILE yb*yc>0 b=cf: cf=cf+dev: e=b: PROCfey: yb=y: e=cf: PROCfey: yc=y ENDWHILE WHILE cf-b>de e=(b+cf)/2: PROCfey: IF y*yb>0 THEN b=e ELSE cf=e ENDWHILE ev(k)=e: b=e: dev=e-ev(k-1): ss=-ss NEXT k ss=1: dx=0.1 FOR k=1 TO ns REM Initialize y(k, i)... e=ev(k): z=1+ss: y(k, nn)=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<=nn THEN y(k, nn+ii)=z ENDIF ENDWHILE s=SQR(2*s*dx) FOR i=0 TO nn REM ...and normalize them y(k, nn+i)=y(k, nn+i)/s: y(k, nn-i)=ss*y(k, nn+i) NEXT i ss=-ss NEXT k ENDPROC DEF PROCfey REM Schr Eqn for test wave function i=0: x=0: dx=.002: y=1+ss: v=1-ss REM Inputs e, returns y a=e*e*(x*x-1)*y: v=v+a*dx/2 WHILE x<1 OR v*y<0 i=i+1: x=x+dx: y=y+v*dx a=e*e*(x*x-1)*y: v=v+a*dx ENDWHILE ENDPROC DEF PROCpotential FOR i=0 TO 60 REM Potential function yt%(60+i)=yo+8*i: yt%(60-i)=yt%(60+i) xt%(60+i)=2*nn+INT(2*n*SQR(4*i/15)): xt%(60-i)=2*np-xt%(60+i) NEXT i ENDPROC