% Fig_1_2and3.m % Essential Semiconductor Laser Device Physics % this program calculates and plots the spherical harmonics using % Y(l,m)=((2l+1)(l-m)!/((l+m)!4pi))^0.5 P(l,m)exp(i*m*phi) % where exp(i*m*phi) = cos(m*phi) +i*sin(m*phi) clear all; clf; FS = 18; %label fontsize 18 FSN = 16; %number fontsize 16 LW = 2; %linewidth % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times New Roman'); set(0,'DefaultAxesFontSize', FSN); % Change default text fonts. set(0,'DefaultTextFontname', 'Times New Roman'); set(0,'DefaultTextFontSize', FSN); l=input('Input quantum number l (l>=0) = '); m=input('Input quantum number m (l>=|m|) = '); if (abs(m) > l) err='Wrong input l or m value, m can NOT be larger than l!' return; end rslv = 401; % the angle resolution for 3D drawing, it is 2*pi/(rslv-1) U = linspace(0, 2*pi, rslv); % Phi in the spherical coordinates V = linspace(0, 2*pi, rslv); % Theta in the spherical coordinates [u,v] = meshgrid(U,V); %===================================================== factor1=factorial(l-m); factor2=factorial(l+m); a=((2*l+1)*factor1/(4*pi*factor2)); pp=legendre(l,cos(v)); %the legendre polynomial if (l == 0) gs = a.*pp; else plm=pp(abs(m)+1,:); tt=1; for k1=1:rslv for k2=1:rslv gs(k2,k1) = a*(plm(tt)^2); tt=tt+1; end; end end gs=gs.^0.5; %|Y(l,m)| gsr=gs.*cos(m*u); %real part gsi=gs.*sin(m*u); %imaginary part % %====================================================== % |Y(l,m)| is modulus of spherical harmonic X=cos(u).*sin(v).*gs; Y=sin(u).*sin(v).*gs; Z=cos(v).*gs; xmax = 0; ymax = 0; zmax = 0; tmp=max(X,[],1); xmax=max(tmp); tmp=max(Y,[],1); ymax=max(tmp); tmp=max(Z,[],1); zmax=max(tmp); amax = max([xmax ymax zmax]); amaxn= amax*(-1.0); %******************** plot figures ************************ figure(1); surf(X,Y,Z,'FaceColor','m','EdgeColor','none') camlight left; lighting phong; axis equal; ttl=['Spherical harmonic |Y(',num2str(l),',',num2str(m),')|']; title(ttl); xlabel('X'); ylabel('Y'); zlabel('Z'); axis([amaxn amax amaxn amax amaxn amax]); %===================================================== % Re(Y(l,m)) is real part of spherical harmonic Xr=cos(u).*sin(v).*gsr; Yr=sin(u).*sin(v).*gsr; Zr=cos(v).*gsr; xrmax = 0; yrmax = 0; zrmax = 0; tmp=max(Xr,[],1); xrmax=max(tmp); tmp=max(Yr,[],1); yrmax=max(tmp); tmp=max(Zr,[],1); zrmax=max(tmp); armax = max([xrmax yrmax zrmax]); armaxn= armax*(-1.0); figure(2); surf(Xr,Yr,Zr,'FaceColor','m','EdgeColor','none') camlight left; lighting phong; axis equal; ttl=[' Real part of spherical harmonic Y(',num2str(l),',',num2str(m),')']; title(ttl); xlabel('X'); ylabel('Y'); zlabel('Z'); axis([armaxn armax armaxn armax armaxn armax]); if (m == 0) return; end %===================================================== % Im(Y(l,m)) is imaginary part of spherical harmonic Xi=cos(u).*sin(v).*gsi; Yi=sin(u).*sin(v).*gsi; Zi=cos(v).*gsi; ximax = 0; yimax = 0; zimax = 0; tmp=max(Xi,[],1); ximax=max(tmp); tmp=max(Yi,[],1); yimax=max(tmp); tmp=max(Zi,[],1); zimax=max(tmp); aimax = max([ximax yimax zimax]); aimaxn= aimax*(-1.0); figure(3); surf(Xi,Yi,Zi,'FaceColor','m','EdgeColor','none') camlight left; lighting phong; axis equal; ttl=['Imaginary part of spherical harmonic Y(',num2str(l),',',num2str(m),')']; title(ttl); xlabel('X'); ylabel('Y'); zlabel('Z'); axis([aimaxn aimax aimaxn aimax aimaxn aimax]);