PROGRAM TDSE ! motion of a wavepacket incident on a potential DIM Re(0:1100),Im(0:1100),Imold(0:1100) CALL parameters(x0,k0,width,V0,a,xmin,xmax,n,dx,dx2,dt) CALL initial_packet(Re(),Im(),x0,k0,width,xmin,n,dx,dt) CALL set_up_windows(xmin,xmax,V0,#1,#2,#3,#4) CALL draw_potential(V0,a,xmin,xmax,#4) DO CALL evolve(Re(),Im(),Imold(),t,V0,a,dx,dx2,dt,xmin,n) CALL plots(Re(),Im(),Imold(),t,xmin,n,dx,#1,#2,#3,#4) LOOP until key input END SUB parameters(x0,k0,width,V0,a,xmin,xmax,n,dx,dx2,dt) LET x0 = -15 ! initial position of packet LET width = 1 ! width of wave packet in x space LET k0 = 2 ! group velocity of packet LET xmax = 20 ! maximum value of x LET xmin = -xmax ! minimum value of x LET V0 = 2 LET a = 1 ! half-width of potential well LET dx = 0.4 ! grid size LET dx2 = dx*dx LET n = (xmax - xmin)/dx LET dt = 0.1 ! time interval END SUB SUB initial_packet(Re(),Im(),x0,k0,width,xmin,n,dx,dt) ! initial Gaussian wavepacket LET delta2 = width*width LET A = (2*pi*delta2)^(-0.25) LET b = 0.5*k0*dt FOR i = 1 to n LET x = xmin + (i-1)*dx LET arg = 0.25*(x - x0)^2/delta2 LET e = exp(-arg) LET Re(i) = A*cos(k0*(x - x0))*e ! value at t = 0 LET Im(i) = A*sin(k0*(x - x0 - 0.5*b))*e ! value at t = dt/2 NEXT i END SUB SUB set_up_windows(xmin,xmax,V0,#1,#2,#3,#4) OPEN #1: screen 0,1,0.75,1.0 ! probability SET WINDOW xmin,xmax,-0.1,0.5 PLOT LINES: xmin,0;xmax,0 OPEN #2: screen 0,1,0.5,0.75 ! real part SET WINDOW xmin,xmax,-1,1 PLOT LINES: xmin,0;xmax,0 OPEN #3: screen 0,1,0.25,0.5 ! imaginary part SET WINDOW xmin,xmax,-1,1 PLOT LINES: xmin,0;xmax,0 OPEN #4: screen 0,1,0,0.22 ! potential SET WINDOW xmin,xmax,-V0,V0 END SUB SUB draw_potential(V0,a,xmin,xmax,#4) DECLARE DEF V WINDOW #4 SET COLOR "red" PLOT LINES: xmin,0;a,0;a,0;a,V0;a,V0;xmax,V0 SET COLOR "black/white" PRINT "total probability = "; END SUB SUB evolve(Re(),Im(),Imold(),t,V0,a,dx,dx2,dt,xmin,n) DECLARE DEF V FOR i = 1 to n LET x = xmin + (i-1)*dx LET HIm = V(x,V0,a)*Im(i) -0.5*(Im(i+1) -2*Im(i) +Im(i-1))/dx2 ! real part defined at multiples of dt LET Re(i) = Re(i) + HIm*dt NEXT i FOR i = 1 to n LET x = xmin + (i-1)*dx ! dt/2 earlier than real part LET Imold(i) = Im(i) LET HRe = V(x,V0,a)*Re(i) -0.5*(Re(i+1) -2*Re(i) +Re(i-1))/dx2 ! dt/2 later than real part LET Im(i) = Im(i) - HRe*dt NEXT i LET t = t + dt ! time of real part END SUB SUB plots(Re(),Im(),Imold(),t,xmin,n,dx,#1,#2,#3,#4) WINDOW #2 CLEAR PLOT LINES: xmin,0;xmax,0 FOR i = 1 to n LET x = xmin + (i - 1)*dx PLOT x,Re(i); NEXT i WINDOW #3 CLEAR PLOT LINES: xmin,0;xmax,0 FOR i = 1 to n LET x = xmin + (i - 1)*dx PLOT x,Im(i); NEXT i WINDOW #1 CLEAR PLOT LINES: xmin,0;xmax,0 LET Psum = 0 FOR i = 1 to n LET x = xmin + (i - 1)*dx LET P = Re(i)*Re(i) + Im(i)*Imold(i) LET Psum = Psum + P*dx PLOT x,P; NEXT i WINDOW #4 SET CURSOR 1,20 PRINT " "; ! erase previous probability SET CURSOR 1,20 PRINT Psum ! total probability should be unity END SUB FUNCTION V(x,V0,a) ! step potential IF x > a then LET V = V0 ELSE LET V = 0 END IF END DEF