6

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.

Improve this question
3
  • A simple way? Unlikely, unless by some chance gnuplot can compute spherical harmonics.
    – Charles Staats
    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?
    – Sanjay Sekaran
    Sep 22, 2015 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, 2015 at 3:30

4 Answers 4

Reset to default
15

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

Improve this answer
3
  • 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, 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?
    – Stefano
    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
    – Mark Wibrow
    Jun 29, 2018 at 11:51
3

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

Improve this answer
1
  • 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, 2020 at 20:45
3

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

\documentclass{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}

\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,}
    ]
    \addplot3
      [
        domain=0:180,   samples=37,   % polar angle
        y domain=0:360, samples y=37, % azimuthal angle
        surf,
        z buffer=sort,
        shader=faceted interp, 
        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) enter image description here

Improve this answer
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');
Improve this answer

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .