Error plotting spherical harmonics in MATLAB: Grid arrays must have NDGRID structure

Question

I am trying to use the spherical harmonics to represent a perturbation of a spherical object in an acoustic and fluid flow in 3D. So far, I have been able to use a 2D perturbation on a 3D spherical object, however, I would like to extend that.

The current code represents a 3D sphere with a 2D perturbation, so the perturbation is only in x and y:

clearvars; clc; close all;
Nx = 128;
Ny = 128;
Nz = 128; 
Lx =128; 
Ly = 128; 
Lz = 128; 
xi  = (0:Nx-1)/Nx*2*pi;
xi_x =  2*pi/Lx;
x =  xi/xi_x;
yi  = (0:Ny-1)/Ny*2*pi;
yi_y =  2*pi/Ly;
y = yi/yi_y;
zi  = (0:Nz-1)/Nz*2*pi;
zi_z =  2*pi/Lz;
z = zi/zi_z;
[X,Y,Z] = meshgrid(x,y,z);
A = 2*pi / Lx;
B = 2*pi / Ly;
C = 2*pi / Lz;
x0 = 64; 
y0 = 64;
z0 = 64; 
rx0 = 20; 
ry0 =  20; 
rz0 =  20;
p = 3;
b = 0.1; % pert amplitude 
c = 12; 
d = 1;
a = 4;
theta = atan2(Y -y0,  X-x0) - (pi/c);
p0 = ((X-x0) .* (X-x0)) /(rx0 * rx0) + ((Y-y0) .* (Y-y0))/(ry0 * ry0) + ((Z-z0) .* (Z-z0))/(rz0 * rz0);
Test =d + a * exp((-1. * p0 .* (1 - b .* cos(c * theta))).^p) ;
figure
isosurface(X,Y,Z,Test);
shading flat;
grid on;

Which returns the isosurface:

However, I would like to achieve something similar to this plot (perturbation in z as well):

This is my attempt using spherical harmonics to reproduce the above picture:

clearvars; clc; close all;
%in spherical coord
%calculate r
Nx = 128;  
Ny = 128;  
Nz = 128;
Lx =128; 
Ly = 128; 
Lz = 128; 
xi  = (0:Nx-1)/Nx*2*pi;
xi_x =  2*pi/Lx;
x =  xi/xi_x;
yi  = (0:Ny-1)/Ny*2*pi;
yi_y =  2*pi/Ly;
y = yi/yi_y;
zi  = (0:Nz-1)/Nz*2*pi;
zi_z =  2*pi/Lz;
z = zi/zi_z;
r = sqrt(x.^2 + y.^2 + z.^2); 
% Create the grid
delta = pi/127; 
%Taking for instance l=1, m=-1 you can generate this harmonic on a (azimuth, elevation) grid like this:
azimuths = 0 : delta : pi; 
elevations = 0 : 2*delta : 2*pi; 
[R, A, E] = ndgrid(r, azimuths, elevations); %A is phi and E is theta
H = 0.25 * sqrt(3/(2*pi)) .* exp(-1j*A) .* sin(E) .* cos(E);
%transform the grid back to cartesian grid like this:
%can also add some radial distortion to make things look nicer:
%the radial part depends on your domain
X = r .* cos(A) .* sin(E); 
Y = r .* sin(A) .* sin(E); 
Z = r .* cos(E); 
%parameters
x0 = 64;  
y0 = 64; 
z0 = 64;  
rx0 = 20; 
ry0 =  20; 
rz0 =  20;
p = 3;
b = 0.1; % pert amplitude 
%c = 12; 
d = 1;
a = 4;
p0 = ((X-x0) .* (X-x0)) /(rx0 * rx0) + ((Y-y0) .* (Y-y0))/(ry0 * ry0) + ((Z-z0) .* (Z-z0))/(rz0 * rz0);
Test1 =d + a * exp((-1. * p0 .*H).^p) ;
figure
isosurface(X,Y,Z,real(Test1)); %ERROR

This gives me the following error:

Error using griddedInterpolant
Grid arrays must have NDGRID structure.

Is the issue in the way I am setting up the spherical harmonics? or the functional form of Test1?? Thanks

Answer 1

I think the problem is that you've created a NDGrid in the second code. In the first code that worked you created a meshgrid.

isosurface expects a grid in meshgrid format and transforms it later into an NDGrid format with permute just for the usage of the griddedInterpolant function. By the input of an NDGrid this operation will fail.

Why did you switch to creating a NDGrid? Did you get the same error when using meshgrid?

EDIT

OK, new theory: I think the problem is the grid itsself. Read the documentation on ndgrid for more information but for short: A NDGrid format is a completely rectangular grid where all nodes are exclusivly surrounded by 90° angles. By spanning up a grid with ndgrid(r, azimuths, elevations) or meshgrid(r, azimuths, elevations) you are getting this rectangular grid but of course this grid is meaningless because r, azimuths and elevations represent spherical coordinates. By later converting R, A and E into the carthesian coordinates X, Y and Z you get a propper spherical grid but the grid isn't a NDGrid structure anymore since it's not rectangular anymore.

So you have to find a workaround to calculate with your spherical grid.

• Maybe you can use another function to visualize your data that dosn't take any grid format as input
• Or maybe you can try working with a rectangular cartesian grid by setting all values outside the sphere to 0. (You have to refine your grid accordingly to achieve a good estimation of the sphere)