% @(#) 6.728 Demo code -- Matlab Code % @(#) This code shows an animated gaussian wave packet in a SHO % @(#) potential. The width of the initial gaussian wave packet % @(#) can be changed to see "breathing" effects. Also, % @(#) and the variance of x are plotted vs. time to see % @(#) the classical result of sinusiodal motion. % 6.728 -- Simple Hamronic Oscillator % Mark Schweizer (mschweiz@mit.edu) clear all; % Clear workspace clc; % Clear the command window close all; % Get rid of old figures wavename=['Gaussian Wave']; % Used for titles in plots MOVIE = 1; % Make this zero to skip the % movie generation % Defines some standard constants used in quantum mechanics and % statistical mechanics. h = 6.62607E-34; % Planck constant [J*s] hbar = h/(2*pi); % Rationalized Planck constant [J*s] Me = 9.1094E-31; % Electron rest mass [kg] Mp = 1.672623E-27; % Proton rest mass [kg] Mn = 1.674929E-27; % Neutron rest mass [kg] c = 2.99792458E8; % Speed of light in vacuum [m/s] % 1 nm is the the first eigenfunction for this SHO. L = input('What is the gaussian width paramter, L, in nanometers? ')*1E-9; % Create our x vector dx=L/8; x = [-12*1E-9:dx:12*1E-9]'; % Define our initial position x0 = L*input('What is the initial displacement in units of L? '); % Create the column vector psi. psi = (L^2*2*pi)^(-1/4) * exp(-((x-x0) / (2*L)).^2); % Here we set up the q vector. "dq" is determined by the periodicity % of the DFT, as is the Max/Min q available. w0=.5*hbar/Me/(1e-9)^2 ; % Picked so ground state gaussian has L=1E-9; g = sqrt(Me*w0/hbar)*x; % Intermediate value. Morrison calls this Zeta. fprintf('Computing hermite polynomials ...\n'); h = herm(25); % Computes 1st 25 hermite polynomials % Compute our basis matrix fprintf('Computing basis matrix ...\n'); phi = zeros(length(x),size(h,1)); for j =[1:size(h,1)] A = (Me*w0/(hbar*pi))^(1/4)*1/sqrt( prod([1:(j-1)])*2^(j-1)); phi(:,j) = A*polyval(h(j,:),g).*(exp(-(g.^2)/2)); end; fprintf('Computing the time dependent wave function...\n'); a = dx * phi' * psi; % Create our energy and time matrices tau = 2*pi/w0; % Chosen to show one full periods of oscillations t = tau*[0:.02:2]'; % 10 steps per period E = w0*hbar*([0:(size(phi,2)-1)]+.5)'; eta = exp(-i* E/hbar*t.'); % Create our time dependent expansion coefficients a_t = (a*ones(size(t'))).*eta; % Find our time dependent psi psi_t = phi *a_t ; psi_t_pdf = conj(psi_t).*psi_t; fprintf('Computing Statistical Quantities...\n'); % Compute the statisitical quantitites of interest. mx = real(diag(psi_t'*diag(x)*psi_t))*dx; mx2=real(diag(psi_t'*diag(x.^2)*psi_t))*dx; msdx=sqrt(mx2-mx.^2); if MOVIE == 1 fprintf('Press any key to generate an animation.\n'); pause; fprintf('Generating the movie sequence...\n'); % Now we'll generate an animated movie ymax =1.1*max(max(psi_t_pdf)); pota =ymax / 12E-9^2; pot =pota*x.^2; frames=(length(t)-1)/4; figure; M = moviein(frames); for f=1:frames plot(x/1E-9,psi_t_pdf(:,f)); hold on; plot(x/1E-9,pot,'g');hold off; axis([-12 12 0 ymax]); M(:,f) = getframe; end delete(gcf); fprintf('Press any key to view the animation.\n'); pause; % Create a title axis to play the movie in figure; axis; mh=gca; % Save the axis handle for later plot(x/1E-9,psi_t_pdf(:,1)); % Setup for proper axis labels title([wavename,' Function in Time'],'FontSize',18); axis([-12 12 0 ymax]); xlabel('X [nm]','FontSize',18); % 00.10.09 ekd % sylabel('\18 {\symbol y^*y} [m^{-1}]'); ylabel('\Psi^* \Psi [m^{-1}]'); movie(mh,M,-4,6); % Play the movie 3 times end fprintf('Press any key to generate a x-t color plot.\n'); pause; % Now, we'll generate a two dimensional plot of x vs t, but the % color will indicate psi^2. The brighter the color, the higher the % probability! figure surf(x/1E-9,t/tau,psi_t_pdf'); view(2); brighten(.7); shading flat; title([wavename,' Function in Time'],'FontSize',18); xlabel('X [nm]','FontSize',18); % 00.10.09 ekd % sylabel('\18 t [{\symbol t}]'); ylabel('t [\tau]'); fprintf('Press any key to plot the statistical qunatitites.\n'); pause; figure plot(t/tau,mx); title(' vs. t','FontSize',18); ylabel('x [nm]','FontSize',18); % 00.10.09 ekd % sxlabel('\18 t [{\symbol t}]'); xlabel('t [\tau]'); %label changed by T. Orlando to read x^2, 29-Sept-00 figure plot(t/tau,msdx.^2); %stitle('\18 ({\symbol D} x)^2 vs. t'); title('(\Delta x)^2 vs. t'); ylabel('x^2 [nm^2]','FontSize',18); % 00.10.09 ekd % sxlabel('\18 t [{\symbol t}]'); xlabel('t [\tau]');