% 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]);
