# Drawing spherical harmonic density plots on the surface of a sphere in tikz/pgfplots

I would like to create spherical harmonic plots in tikz/pgfplots similar to the one below, except without the axes and with labels above every image: It would be good if I could specify different combinations of the degree l and the order m of the associated Legendre polynomial, as well as the viewing angle with the poles of the sphere demarcated.

One possibility is plotting the spherical harmonics in MATLAB, output them to matrices and use MATLAB2Tikz to transform them to tikz files.

• A simple way? Unlikely, unless by some chance gnuplot can compute spherical harmonics. Sep 22, 2015 at 2:59
• Perhaps simple is asking for a bit too much. Is there any way of doing this? Maybe even plotting them in MATLAB and doing a MATLAB2Tikz conversion? Sep 22, 2015 at 3:22
• If you can get Matlab to plot the spherical harmonics (as rho = function of (phi,theta)) in spherical coordinates and output the result as an array/table (preferably to a file), there's probably a way to get pgfplots to read them in as a color map for a sphere. That's likely your best bet if you want to be able to manipulate them using tikz. Sep 22, 2015 at 3:30

No idea if this could be done "natively" with asymptote, pstricks or TikZ, or even by calculating all the data in an external program and then plotting it with pgfplots. I went for doing everything in python and simply including the resulting image.

On the python side this requires numpy, scipy, matplotlib and mpl_tookits.basemap, and on the (pdf) latex side it must be compiled with -shell-escape. Most (but not all) options are parametrized with keys that can be called from latex.

\documentclass[tikz,border=5]{standalone}
\usepackage{filecontents}
\begin{filecontents*}{shpl.py}
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import scipy.special as sp
def plot(filename, m, n, longitude=0, latitude=0, inches=(1,1),
cmap='RdYlBu', points=100):

figure, ax = plt.subplots(1,1)
figure.set_size_inches(*inches)

lon = np.linspace(0, 2*np.pi, points)
lat = np.linspace(-np.pi / 2, np.pi / 2, points)
colat = lat + np.pi / 2
d = np.zeros((len(lon), len(colat)), dtype=np.complex64)

meshed_grid = np.meshgrid(lon, lat)
lat_grid = meshed_grid
lon_grid = meshed_grid

mp = Basemap(projection='ortho', lat_0=latitude, lon_0=longitude, ax=ax)
mp.drawmapboundary()
mp.drawmeridians(np.arange(0, 360, 30))
mp.drawparallels(np.arange(-90, 90, 30))

for j, yy in enumerate(colat):
for i, xx in enumerate(lon):
d[i,j] = sp.sph_harm(m, n, xx, yy)

drm = np.round(np.real(d) / np.max(np.real(d)), 2)
x, y = mp(np.degrees(lon_grid), np.degrees(lat_grid))
mp.pcolor(x, y, np.transpose(drm), cmap=cmap)

figure.savefig(filename, transparent=True)
\end{filecontents*}

\newif\ifshpoverwrite
\tikzset{%
spherical harmonics/.cd,
overwrite/.is if=shpoverwrite,
file/.store in=\shpfilename,
m/.store in=\shpm,
n/.store in=\shpn,
longitude/.store in=\shplongitude,
latitude/.store in=\shplatitude,
cmap/.store in=\shpcmap,
points/.store in=\shppoints,
inches/.store in=\shpinches,
longitude=0, latitude=0,
cmap=RdYlBu,  points=100, inches={(1,1)}
}
\def\sphericalharmonicplot#1{%
\tikzset{spherical harmonics/.cd,#1}%
\edef\pythoncommand{python -c "import shpl;
shpl.plot('\shpfilename', \shpm, \shpn,
latitude=\shplatitude, longitude=\shplongitude,
cmap='\shpcmap', points=\shppoints, inches=\shpinches)"}%
\ifshpoverwrite
\immediate\write18{\pythoncommand}%
\else
\IfFileExists{\shpfilename}{}{\immediate\write18{\pythoncommand}}%
\fi%
\includegraphics{\shpfilename}%
}
\begin{document}
\begin{tikzpicture}[x=1in,y=1in]
\foreach \m/\n [count=\i from 0] in {0/1, 0/2, 0/3, 1/1, 1/2, 1/3,
2/2, 2/3, 3/6, 4/5, 5/7, 6/10}
\node [label=270:{$m=\m,\,n=\n$}] at ({floor(\i/3)*1.5}, {-mod(\i,3)*1.5})
{\sphericalharmonicplot{file=sph\i.png, m=\m, n=\n,
longitude=-100, latitude=30}};
\end{tikzpicture}
\end{document} • Your solution looks great. However, I had managed to use MATLAB to come up with my own solution. Thank you for your efforts, Mark. Sep 29, 2015 at 2:33
• @Mark Wibrow - Your solution looks great and I would like to make my own ones with the python code you shared. However, I am having problems for importing the "basemap" from the "mpl_toolkits" module in Python 3. It seems that it no longer exist. How did you overcome that? Jun 6, 2018 at 16:06
• @Stefano Almost 3 years ago, I was using Python 2. Python 3 has some incompatibility issues with the basemap stuff. One possible solution is this stackoverflow answer Jun 29, 2018 at 11:51

I could not replicate @Mark solution in Python 3 due to the missing "basemap" in the current Matplotlib module.

# I finally found a neat solution which outputs 3 types of plots Please check the Python code here

• I checked your Python code and it seems the z axis is inverted. If π is in the -z direction (usual convention), then you have to define colat = (lat + np.pi / 2)[::-1] Jul 7, 2020 at 20:45

Here's a solution in pgfplots (a little weirdness in the coordinate names...):

\documentclass{standalone}
\usepackage{pgfplots}

\begin{document}
\begin{tikzpicture}
\begin{axis}
[
view={140}{30},
xmin=-1.1,xmax=1.2,ymin=-1.1,ymax=1.2,zmin=-1.25,zmax=1.3,
width=12cm,height=12cm,
axis equal,
axis lines=center,
ticks = none,
colormap/cool,
xlabel={$x$},
ylabel={$y$},
zlabel={$z$},
zlabel style={at={(zticklabel* cs:1)},anchor=south,},
ylabel style={at={(yticklabel* cs:1)},anchor=north west,},
xlabel style={at={(xticklabel* cs:1)},anchor=east,}
]
[
domain=0:180,   samples=37,   % polar angle
y domain=0:360, samples y=37, % azimuthal angle
surf,
z buffer=sort,
opacity=0.85,
%       note weirdness: acos(z/r) = polar angle, atan2(y,x) = azimuth in degrees
%        point meta={(z/sqrt(x*x+y*y+z*z))} % zonal harmonic
point meta={sin(acos(z/sqrt(x*x+y*y+z*z)))^2*cos(acos(z/sqrt(x*x+y*y+z*z)))*cos(2*atan2(y,x))*sqrt(105/32*pi)}% tesseral harmonic (2,3)
] (
{sin(x)*cos(y)}, % x coord
{sin(x)*sin(y)}, % y coord
{cos(x)}         % z coord
);
% hacking the opacity for axes
\draw[black, ] (0,1,0) -- (0,1.4,0);
\draw[black, ] (1,0,0) -- (1.4,0,0);
\draw[black, ] (0,0,1) -- (0,0,1.2);
\end{axis}
\end{tikzpicture}

\end{document}


Here's a tesseral harmonic (2,3) I had managed to use MATLAB's spherical harmonics function and use matlab2tikz to convert it into a tikz image. Solution adapted from http://au.mathworks.com/help/matlab/examples/animating-a-surface.html?refresh=true

MATLAB code:

theta = 0:pi/40:pi;                   % polar angle
phi = 0:pi/20:2*pi;                   % azimuth angle

[phi,theta] = meshgrid(phi,theta);    % define the grid

degree = 0;
order = 0;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,1)
surf(x,y,z, rho);
title('$\ell=0, m=0$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 1;
order = 0;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,6)
surf(x,y,z, rho);
title('$\ell=1, m=0$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 1;
order = 1;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,7)
surf(x,y,z, rho);
title('$\ell=1, m=\pm 1$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 2;
order = 0;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,11)
surf(x,y,z, rho);
title('$\ell=2, m=0$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 2;
order = 1;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,12)
surf(x,y,z, rho);
title('$\ell=2, m=\pm 1$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 2;
order = 2;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,13)
surf(x,y,z, rho);
title('$\ell=2, m=\pm 2$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 3;
order = 0;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,16)
surf(x,y,z, rho);
title('$\ell=3, m=0$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 3;
order = 1;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,17)
surf(x,y,z, rho);
title('$\ell=3, m=\pm 1$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 3;
order = 2;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,18)
surf(x,y,z, rho);
title('$\ell=3, m=\pm 2$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 3;
order = 3;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,19)
surf(x,y,z, rho);
title('$\ell=3, m=\pm 3$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 4;
order = 0;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,21)
surf(x,y,z, rho);
title('$\ell=4, m=0$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 4;
order = 1;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,22)
surf(x,y,z, rho);
title('$\ell=4, m=\pm 1$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 4;
order = 2;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,23)
surf(x,y,z, rho);
title('$\ell=4, m=\pm 2$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 4;
order = 3;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,24)
surf(x,y,z, rho);
title('$\ell=4, m=\pm 3$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degree = 4;
order = 4;
amplitude = 0.5;

Ymn = legendre(degree,cos(theta(:,1)));
Ymn = Ymn(order+1,:)';
yy = Ymn;

for kk = 2: size(theta,1)
yy = [yy Ymn];
end

yy = yy.*cos(order*phi);

order = max(max(abs(yy)));

r = radius.*sin(theta);    % convert to Cartesian coordinates
x = r.*cos(phi);
y = r.*sin(phi);

subplot(5,5,25)
surf(x,y,z, rho);
title('$\ell=4, m=\pm 4$')

axis equal off      % set axis equal and remove axis
view(0,30)         % set viewpoint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

map = makeColorMap([0.2 0.2 0.6],[1.0 0.99 0.72],[0.8 0.25 0.33],80);
colormap(map);
cd(Figures)