我正在尝试使用球谐函数来表示三维声学和流体流动中球形物体的扰动。到目前为止，我已经能够在3D球体上使用2D微扰，然而，我想扩展一下。
当前代码表示具有2D扰动的3D球体，因此该扰动仅在x和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;

它返回等值面：

但是，我希望实现类似于此图的结果(也是z中的扰动)：

这是我尝试用球谐函数来再现上面的画面：

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

这会给我带来以下错误：

Error using griddedInterpolant
Grid arrays must have NDGRID structure.

问题出在我建立球谐函数的方式上吗？还是Test1的函数形式？？谢谢