# 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. – Charles Staats Sep 22 '15 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? – Sanjay Sekaran Sep 22 '15 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. – Charles Staats Sep 22 '15 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. – Sanjay Sekaran Sep 29 '15 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? – Stefano Jun 6 '18 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 – Mark Wibrow Jun 29 '18 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] – PerroNoob Jul 7 at 20:45

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)