Hydrogen orbitals   

Matlab code 

% Movie of all hydrogen orbitals (n,l,m), up to nmax, determined by the user. 

% Adapted from code by Peter van Alem. See 

% https://ch.mathworks.com/matlabcentral/fileexchange/64274-hydrogen-orbitals 

 

close all

 

vid = VideoWriter('hydrogenOrbitals.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 30;

open(vid);

 

% Generate a sinusoidal ramp up and down of transparency when transitioning

% between orbitals 

x1 = (0:pi/15:14*pi/15);

fade1 = (-cos(x1)+1)/2; % Ramp up opacity

fade2 = ones(1,1); % One fully opaque image 

fade3 = (cos(x1+pi/15)+1)/2; % Ramp down opacity 

fade = [fade1 fade2 fade3]; % Concatenate 

 

% quantum numbers

% n

% 0 <= l < n

% -l <= m <= l

nmax = 4; % Maximum value of principle quantum number. Note that the number 

% of orbitals for a given n = n^2, and that the total number of orbitals 

% between n = 1 and nmax increases as nmax(nmax+1)(2nmax+1)/6, so for 

% nmax = 4, there are 30 orbitals, for nmax = 5, there are 55, etc, etc 

 

nOpaque = 20; % Number of fully opaque frames for a given orbital 

 

for n= 1:nmax

    for l = 0:n-1

        for m = -l:l

            for alphaVal = fade

                close all

                % plotting parameters

                probabilitydensity = 1.25e-5;

                a = 1;  % Bohr radius

                % angular part (Condon-Shortley)

                SphericalYlm = @(l,m,theta,phi)(-1)^m*sqrt((2*l+1)/(4*pi)* ...

                    factorial(l-abs(m))/factorial(l+abs(m)))* ...

                    AssociatedLegendre(l,m,cos(theta)).*exp(1i*m*phi);

                % real basis

                if (m < 0)

                    Y = @(l,m,theta,phi)sqrt(2)*(-1)^m* ...

                        imag(SphericalYlm(l,abs(m),theta,phi));

                elseif (m == 0)

                    Y = @(l,m,theta,phi)SphericalYlm(l,m,theta,phi);

                else

                    Y = @(l,m,theta,phi)sqrt(2)*(-1)^m*real(SphericalYlm(l,m,theta,phi));

                end

                % radial part

                R = @(n,l,r)sqrt((2/(a*n))^3*factorial(n-l-1)/(2*n*factorial(n+l))).* ...

                    exp(-r/(a*n)).*(2*r/(a*n)).^l*1/factorial(n-l-1+2*l+1) .* ...

                    AssociatedLaguerre(n-l-1,2*l+1,2*r/(a * n));

                % wave function

                psi = @(n,l,m,r,theta,phi)R(n,l,r).*Y(l,m,theta,phi);

                % setting the grid

                border = 32; 

                % "accuracy" slows program approximately with the cube of its value! 

                % Bump up to 500 only for n = 2, l = 0, due to artefact of

                % coloring of inner sphere. Just being super OCD here 

                if (n == 2) && (l ==0) 

                    accuracy = 500; % Set to 200 if you're not bothered by artefact 

                else

                    accuracy = 200; 

                end

                raster = linspace(-border,border,accuracy);

                [x, y, z] = ndgrid(raster,raster,raster);

                % conversion Cartesian to spherical coordinates

                r = sqrt(x.^2+y.^2+z.^2);

                theta = acos(z./r);

                phi = atan2(y,x);

                % plot orbital,  - and + wave function phase

                colors = sign(psi(n,l,m,r,theta,phi));

                isosurface(psi(n,l,m,r,theta,phi).^2,probabilitydensity,colors);

                alpha(0.8*alphaVal);

                colormap([0 0 1; 1 0.5 0])

                material shiny

                

                title(['$ n = ', num2str(n), ', l = ', num2str(l), ...

                    ', m = ', num2str(m) '$'], 'interpreter', 'latex', ...

                    'FontSize', 20)

                set(gcf,'color','w');

                set(gca,'CameraViewAngle',45,'Projection','perspective');

                camzoom(4)

                axis equal

                axis vis3d

                axis off % Comment out if you want axes shown. I find them unnecessary 

 

                % I find labels ugly so commented them out. You can

                % uncomment as you see fit 

                % xticklabels('');

                % yticklabels('');

                % zticklabels('');

                % xlabel('$x$', 'interpreter', 'latex', 'FontSize', 20)

                % ylabel('$y$', 'interpreter', 'latex', 'FontSize', 20)

                % zlabel('$z$', 'interpreter', 'latex', 'FontSize', 20)

                % Store the frame

                if (alphaVal == 1) % When full opaque, record nOpaque frames  

                    for ii = 1:nOpaque

                        frame = getframe(gcf);

                        writeVideo(vid,frame);

                    end

                else

                    frame = getframe(gcf);

                    writeVideo(vid,frame);

                end                

            end

        end

    end

end

% Output the movie as an mpg file

close(vid);

 

% Functions 

% _________

function Anm = AssociatedLaguerre(n,m,x)

Anm = 0;

for i = 0 : n

    Anm = Anm+factorial(m+n)*nchoosek(m+n,n-i)/factorial(i)*(-x).^i;

end 

end

 

function Alm = AssociatedLegendre(l,m,x)

Alm = 0;

for r = 0:floor(1/2 * l-1/2 * abs(m))

    Alm = Alm + (-1)^r * nchoosek(l-2*r,abs(m))*nchoosek(l,r) * ...

        nchoosek(2*l-2*r,l)*x.^(l-2*r-abs(m));

end

Alm = (1-x.^2).^(abs(m)/2).*(factorial(abs(m))/2^l*Alm);

end