4

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
  • 1
    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
12

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.

This adapts the python code written by Alex J. DeCaria, taken from the ipython notebook linked on this page.

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[1]
    lon_grid = meshed_grid[0]

    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}

enter image description here

  • 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
2

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

enter image description here

enter image description here enter image description here

Please check the Python code here

1

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;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

shading interp

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

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

degree = 1;
order = 0;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 1;
order = 1;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 2;
order = 0;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 2;
order = 1;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 2;
order = 2;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 3;
order = 0;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 3;
order = 1;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 3;
order = 2;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 3;
order = 3;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 4;
order = 0;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 4;
order = 1;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 4;
order = 2;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 4;
order = 3;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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

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

degree = 4;
order = 4;
amplitude = 0.5;
radius = 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)));
rho = radius + amplitude*yy/order;

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

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

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)
addpath(genpath([pwd '/../matlab2tikz']))

fig = figure(1);

matlab2tikz('SH1.tikz','height', '\figureheight', 'width', '\figurewidth','parseStrings',false, 'Floatformat', '%.4f');

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.