McGill University
School of Computer Science

Computer Science COMP 199 (Winter term)
Excursions in Computing Science

%MATLABpak08cIII: circWave.m, spherHarm.m and supporting progs, atom.m and
% supporting progs.
%
%	spherHarm.m			spherHarmEval.m
%
%	atom.m				radialAtomEval.m
%	atom.m				spherHarmEval.m
%	atom.m				radialAtomNodes.m

% function circWave(m)		THM				090114,090123
% plot exp(i*m*phi) and extended to circle
function circWave(m)
phi = 0:pi/20/m:2*pi;
for k = m:m+4
 subplot(4,5,k-m+1)
 polar(phi,cos(k*phi),'r')
 title(['real: cos(' int2str(k) '\phi)'])
 subplot(4,5,k-m+6)
 polar(phi,1 + 0.5*cos(k*phi),'r')
 hold on
 polar(phi,cos(0*phi))
 title(['1 + 0.5cos(' int2str(k) '\phi)'])
 subplot(4,5,k-m+11)
 polar(phi,sin(k*phi),'r')
 title(['imag: sin(' int2str(k) '\phi)'])
 subplot(4,5,k-m+16)
 polar(phi,1 + 0.5*sin(k*phi),'r')
 hold on
 polar(phi,cos(0*phi))
 title(['1 + 0.5sin(' int2str(k) '\phi)'])
 hold off
end %for

% function spherHarm(lm)			THM			081224
% Six plots of spherical harmonics:
% Top: real, imaginary and squared components on phi-theta plane
% Bottom: real, imaginary and squared components colour coded on sphere
% Plot Y_l^mY*_l^m on sphere 0:2pi * 0:pi
% NB 090108 instead of using the ..Node.. programs to show nodes, it seems best
% to use this program, setting colormap(hsv(2)) or colormap(copper(2))from 
% colormap('default')
% Look for n-l-1 radial 0-nodes, l nodes on spherical harmonics (Herzberg p.41)
% Uses	SH = spherHarmEval(lm,theta,phi)
function spherHarm(lm)
[theta,phi] = meshgrid(0:pi/15:pi,0:2*pi/15:2*pi);
 SH = spherHarmEval(lm,theta,phi);
 subplot(2,3,1)
 surfc(theta,phi,real(SH))
%colorbar
%colormap(winter)
 xlabel('theta'), ylabel('phi'), zlabel(lm)
 title('real')
%hold on
%figure(2)
 subplot(2,3,2)		% 2nd figure replaces winter colormap by autumn
 surfc(theta,phi,imag(SH))
%colorbar
%colormap(autumn)
 xlabel('theta'), ylabel('phi'), zlabel(lm)
 title('imag')
 subplot(2,3,3)
 surfc(theta,phi,SH.*conj(SH))
 xlabel('theta'), ylabel('phi'), zlabel(lm)
 title('prob')
%figure(3)
[X Y Z] = sphere(20);
theta = acos(Z); phi = atan2(Y,X);
SH = spherHarmEval(lm,theta,phi);
subplot(2,3,4)
surf(X,Y,Z,real(SH))
%colorbar
xlabel('x'), ylabel('y'), zlabel('z')
title([lm ' real'])
subplot(2,3,5)
surf(X,Y,Z,imag(SH))
%colorbar
xlabel('x'), ylabel('y'), zlabel('z')
title([lm ' imag'])
subplot(2,3,6)
surf(X,Y,Z,SH.*conj(SH))
colorbar
xlabel('x'), ylabel('y'), zlabel('z')
title([lm ' prob'])

% function SH = spherHarmEval(lm,theta,phi)		THM		081225
% evaluate spherical harmonic selected by string lm (e.g., '0,0, '1,-1')
% at given angles theta, phi
% Used by spherHarm.m, scatterRadialTheta.m, atom.m
function SH = spherHarmEval(lm,theta,phi)
switch(lm)
 case {'0,0'}
  SH = (1/2)*sqrt(1/pi).*cos(0*theta).*cos(0*phi);
 case {'1,-1'}
  SH = (1/2)*sqrt(3/2/pi).*sin(theta).*exp(-i.*phi);
 case {'1,0'}
  SH = (1/2)*sqrt(3/pi).*cos(theta).*cos(0*phi);
 case {'1,1'}
  SH = -(1/2)*sqrt(3/2/pi).*sin(theta).*exp(i.*phi);
 case {'2,-2'}
  SH = -(1/4)*sqrt(15/2/pi).*sin(theta).^2.*exp(-2*i.*phi);
 case {'2,-1'}
  SH = (1/2)*sqrt(15/2/pi).*sin(theta).*cos(theta).*exp(-i.*phi);
 case {'2,0'}
  SH = (1/4)*sqrt(5/pi).*(3*cos(theta).^2 - 1).*cos(0*phi);
 case {'2,1'}
  SH = -(1/2)*sqrt(15/2/pi).*sin(theta).*cos(theta).*exp(i.*phi);
 case {'2,2'}
  SH = (1/4)*sqrt(15/2/pi).*sin(theta).^2.*exp(2*i.*phi);
 case {'3,-3'}
  SH = (1/8)*sqrt(35/pi)*sin(theta).^3.*exp(-3*i*phi);
 case {'3,-2'}
  SH = (1/4)*sqrt(105/2/pi)*sin(theta).^2.*cos(theta).*exp(-2*i*phi);
 case {'3,-1'}
  SH = (1/8)*sqrt(21/pi)*sin(theta).*(5*cos(theta).^2 - 1).*exp(-i*phi);
 case {'3,0'}
  SH = (1/4)*sqrt(7/pi).*(5*cos(theta).^3 - 3*cos(theta)).*cos(0*phi);
 case {'3,1'}
  SH = -(1/8)*sqrt(21/pi)*sin(theta).*(5*cos(theta).^2 - 1).*exp(i*phi);
 case {'3,2'}
  SH = (1/4)*sqrt(105/2/pi)*sin(theta).^2.*cos(theta).*exp(2*i*phi);
 case {'3,3'}
  SH = -(1/8)*sqrt(35/pi)*sin(theta).^3.*exp(3*i*phi);
 case {'4,-4'}
  SH = (3/16)*sqrt(35/2/pi)*sin(theta).^4.*exp(-4*i*phi);
 case {'4,-3'}
  SH = (3/8)*sqrt(35/pi)*sin(theta)^3.*cos(theta).*exp(-3*i*phi);
 case {'4,-2'}
  SH = (3/8)*sqrt(5/2/pi)*sin(theta).^2.*(7*cos(theta).^2 - 1).*exp(-2*i*phi);
 case {'4,-1'}
  SH = (3/8)*sqrt(5/pi)*sin(theta).*(7*cos(theta).^3 - 3*cos(theta)).*exp(-i*phi);
 case {'4,0'}
  SH = (3/16)*sqrt(1/pi)*(35*cos(theta).^4 - 30*cos(theta).^2 +3).*cos(0*phi);
 case {'4,1'}
  SH = -(3/8)*sqrt(5/pi)*sin(theta).*(7*cos(theta).^3 - 3*cos(theta)).*exp(i*phi);
 case {'4,2'}
  SH = (3/8)*sqrt(5/2/pi)*sin(theta).^2.*(7*cos(theta).^2 - 1).*exp(2*i*phi);
 case {'4,3'}
  SH = -(3/8)*sqrt(35/pi)*sin(theta)^3.*cos(theta).*exp(3*i*phi);
 case {'4,4'}
  SH = (3/16)*sqrt(35/2/pi)*sin(theta).^4.*exp(4*i*phi);
 case {'5,-5'}
  SH = (3/32)*sqrt(77/pi)*sin(theta).^5.*exp(-5*i*phi);
 case {'5,-4'}
  SH = (3/16)*sqrt(385/2/pi)*sin(theta).^4.*cos(theta).*exp(-4*i*phi);
 case {'5,-3'}
  SH = (1/32)*sqrt(385/pi)*sin(theta).^3.*(9*cos(theta).^2 - 1).*exp(-3*i*phi);
 case {'5,-2'}
  SH = (1/8)*sqrt(1155/2/pi)*sin(theta).^2.*(3*cos(theta).^3 - cos(theta)).*exp(-2*i*phi);
 case {'5,-1'}
  SH = (1/16)*sqrt(165/2/pi)*sin(theta).*(21*cos(theta).^4 - 14*cos(theta).^2 + 1).*exp(-i*phi);
 case {'5,0'}
  SH = (1/16)*sqrt(11/pi)*(63*cos(theta).^5 - 70*cos(theta).^3 + 15*cos(theta)).*cos(0*phi);
 case {'5,1'}
  SH = -(1/16)*sqrt(165/2/pi)*sin(theta).*(21*cos(theta).^4 - 14*cos(theta).^2 + 1).*exp(i*phi);
 case {'5,2'}
  SH = (1/8)*sqrt(1155/2/pi)*sin(theta).^2.*(3*cos(theta).^3 - cos(theta)).*exp(2*i*phi);
 case {'5,3'}
  SH = -(1/32)*sqrt(385/pi)*sin(theta).^3.*(9*cos(theta).^2 - 1).*exp(3*i*phi);
 case {'5,4'}
  SH = (3/16)*sqrt(385/2/pi)*sin(theta).^4.*cos(theta).*exp(4*i*phi);
 case {'5,5'}
  SH = -(3/32)*sqrt(77/pi)*sin(theta).^5.*exp(5*i*phi);
 otherwise disp(['spherHarmEval: unknown case ' lm])
end


% function atom(n,l,m)		THM			090120
% WAS function radialThetaNscatt(n,l,m)
% four plots involving radal and theta components for inverse-square
% centrally symmetric field (see Landau & Lifshitz pp.121--4, spherical harm)
% NB r=1 is r = hbar^2/mu/e^2 = 0.0529nm, Bohr radius, for 1 electron
% Look for n-l-1 radial 0-nodes, l nodes on spherical harmonics (Herzberg p.41)
% 1) R vs r; also R^2 and (rR)^2 for integrated effect
% 2) probability vs radius r and azimuth t (theta)
% 3) "cloud" from many atoms on r-t plane
% 4) "cloud" from many atoms in 3D (using uniform distribution for phi)
% Uses	R = radialAtomEval(n,l,r)
% 	SH = spherHarmEval(lm,theta,phi)
%	nodes = radialAtomNodes(n,l)
function atom(n,l,m)
r = 0.1:0.1*n:10*n;
R = radialAtomEval(n,l,r);
subplot(2,2,1)
if l==0 & n==1	% (pity we can't use plotyy() for 3 plots)
 plot(r,R,'b',r,R.^2,'g',r,(r.*R).^2,'r')
 legend('R','R^2','(rR)^2')
else	% enlarge R^2
 plot(r,R,'b',r,10*R.^2,'g',r,(r.*R).^2,'r')
 legend('R','10R^2','(rR)^2')
end
xlabel('r')
title(['Radial part of ' int2str(n) ',' int2str(l) ' wavefunction in 1/r potential'])
% mesh part
[r,t] = meshgrid( 0.1:0.5*n:10*n, 0:pi/20:pi );
RT2 = radialAtomEval(n,l,r).^2.*spherHarmEval([int2str(l) ',' int2str(m)],t,0.0).^2;
RT2max = max(max(RT2));
subplot(2,2,3)
meshc(r,t,RT2)
%axis([0 10*n 0 pi -RT2max/2 RT2max])
xlabel('r'), ylabel('\theta')
title(['Probability vs r, \theta for n,l,m =' int2str(n) ',' int2str(l) ',' int2str(m) ' in 1/r potential'])
% scattering part
% Following Ross p.263 to generate random variable X from uniform U:
%	generate uniform random U on [0,1];
%	if U < Sum pj = 0.01*j^2 j = 1:k then X = xk
sizeU = 1000;
Ur = rand(1,sizeU);	% 1-by-sizeU array
Ut = rand(1,sizeU);	% 1-by-sizeU array
r = 0.1:0.1*n:10*n;
%r = 0.1:2*n:10*n;
R1 = radialAtomEval(n,l,r).^2;
R = R1/sum(sum(R1));	% probabilities should sum to 1
theta = 0:pi/25:pi;
%theta = 0:pi/5:pi;
T1 = spherHarmEval([int2str(l) ',' int2str(m)],theta,0.0).^2;
T = T1/sum(sum(T1));	% probabilities should sum to 1
for k = 1:sizeU
  Rindx = cumSumPos(Ur(k),R);
  rad(k) = r(Rindx);
  Tindx = cumSumPos(Ut(k),T);
  thet(k) = theta(Tindx);
end %for n
S = 5; %size
switch n-l-1				% set colour according to shell
 case 0, C = 'b';			% all blue
 case 1
  nodes = radialAtomNodes(n,l);
  for k = 1:sizeU
    if rad(k) < nodes(1), C(k,:) = [0 0 1];	% 1st shell blue
    else C(k,:) = [0 1 0]; end % if		% 2nd shell green
  end %for k
 case 2
  nodes = radialAtomNodes(n,l) ;
  for k = 1:sizeU
    if rad(k) < nodes(1), C(k,:) = [0 0 1];	% 1st shell blue
    elseif rad(k) < nodes(2), C(k,:) = [0 1 0];	% 2nd shell green
    else C(k,:) = [1 0 0]; end % if		% 3rd shell red
  end %for k
end % switch n-l-1
subplot(2,2,2)
scatter(rad,thet,S,C)
xlabel('r'), ylabel('\theta')
title(['r-\theta plot of ' num2str(sizeU) ' atoms, n,l,m=' int2str(n) ',' int2str(l) ',' int2str(m)])
phi = rand(1,sizeU)*2*pi;	% 1-by-sizeU array
X = rad.*sin(thet).*cos(phi);
Y = rad.*sin(thet).*sin(phi);
Z = rad.*cos(thet);
subplot(2,2,4)
scatter3(X,Y,Z,S,C)
xlabel('x'), ylabel('y'), zlabel('z')
axis equal
title(['x-y-z plot of ' num2str(sizeU) ' atoms, n,l,m=' int2str(n) ',' int2str(l) ',' int2str(m)])

% function R = radialAtomAval(n,l)		THM			090106
% find radial dependence on wavenumber k and angular momentum l of motion in
% inverse-square centrally symmetric field (see Landau & Lifshitz pp.121--4)
% NB r=1 is r = hbar^2/mu/e^2 = 0.0529nm, Bohr radius, for 1 electron
% Look for n-l-1 radial 0-nodes, l nodes on spherical harmonics (Herzberg p.41)
% Used by scatterRadialTheta.m, atom.m
function R = radialAtomEval(n,l,r)
switch n
 case 1
 switch l
  case 0
  R = 2*exp(-r);
  otherwise disp(['radialAtomEval: n=1 cannot have l=' num2str(l)])
 end %switch l
 case 2
 switch l
  case 0
  R = (1/sqrt(2))*(1-r/2).*exp(-r/2);
  case 1
  R = (1/2/sqrt(6))*r.*exp(-r/2);
  otherwise disp(['radialAtomEval: n=2 cannot have l=' num2str(l)])
 end %switch l
 case 3
 switch l
  case 0
  R = (2/3/sqrt(3))*(1 - 2*r/3 + 2*r.^2/27).*exp(-r/3);
  case 1
  R = (8/27/sqrt(6))*r.*(1 - r/6).*exp(-r/3);
  case 2
  R = (4/81/sqrt(30))*r.^2.*exp(-r/3);
  otherwise disp(['radialAtomEval: n=3 cannot have l=' num2str(l)])
 end %switch l
 otherwise disp(['radialAtomEval: unknown case n=' num2str(n)])
end

% function nodes = radialAtomNodes(n,l)		THM			090120
% find zeros of radial dependence on wavenumber k and angular momentum l of
% motion in inverse-square centrally symmetric field (see Landau & Lifshitz
% pp.121--4). There are n-l-1 zeros, so needed only for l