Error using Spherical harmonics: Array dimensions must match for binary array op.

I am trying to use the spherical harmonics to represnet a perturbation of a spherical object in an acoustic and fluid flow in 3D.

For example, the current code I have represnts a 3D sphere with a 2D perturbation so the perturbation is only in x and y:

Nx = 128;
Ny = 128; 
Nz = 128;

Lx =128; 
Ly = 128; 
Lz = 128; 

x  = (0:Nx-1)/Nx*2*pi;
y  = (0:Ny-1)/Ny*2*pi;
z   = (0:Nz-1)/Nz*2*pi;

dx = Lx/Nx;
dy = Ly/Ny;
dz = Lz/Nz;
%make mesh
[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) ;
%plot
figure
isosurface(X,Y,Z,Test);

Now, I would like to include the spherical harmonics to introduce a 3D perturbation to my sphere, so I did the following (from MATLAB):

%test case of m =2 and l =3
dx = pi/60;
col = 0:dx:pi;
az = 0:dx:2*pi;
[phi,theta] = meshgrid(az,col);
l = 3;
Plm = legendre(l,cos(theta));
m = 2;
if l ~= 0
    Plm = reshape(Plm(m+1,:,:),size(phi));
end
aa = (2*l+1)*factorial(l-m);
bb = 4*pi*factorial(l+m);
CC = sqrt(aa/bb);
Ylm = CC .*Plm .*exp(1i*m*phi);
size(d)

ans = 1×2

1 1

size(a)

ans = 1×2

1 1

size(p0)

ans = 1×3

128 128 128

size(Ylm)

ans = 1×2

61 121

size(p)

ans = 1×2

1 1

Test1 =d + a * exp((-1. * p0 .* Ylm).^p); %ERROR

Arrays have incompatible sizes for this operation.

%plot
figure
isosurface(X,Y,Z,Test1);

I get the error:

Array dimensions must match for binary array op.

Error in PlasmaCloudPert3D (line 100)
Test =d + a * exp((-1. * p0 .* Ylm).^p); 

The thing I am trying to get a plot that looks like:

Harmonium - An OpenGL spherical harmonic high resolution modeller (neolithicsphere.com)

Maybe my functional form is incorrect??

[Torsten]

See above. The sizes of p0 and Ylm are not compatible for that p0.*Ylm would be defined.

[Walter Roberson]

You cannot have p0.*Ylm when p0 is 1 x 3 and Ylm is 1 x 2. You also cannot have p0*Ylm when p0 is 1 x 3 and Ylm is 1 x 2.

If you were to transpose p0 you could get 3 x 2, and it is valid to take a 3 x 2 .^ with a 1 x 2, getting a 3 x 2 result. But a is 1 x 2 and you cannot * between a 1 x 2 and a 3 x 2 but you could use .* to get a 3 x 2 result. And you would be able to add the 1 x 2 d value to that, for an overall result of 3 x 2.

But... a 3 x 2 volume is not compatible with a 128 x 128 X Y Z

Or is Test1 intended to get down to a scalar ?

[Jamie Al]

I see the issue, mm well Test1 is suposed to have the same dimenstions as Test, but I guess I don't know how to represent Ylm to have dim of X*Y*Z.

[Walter Roberson]

Opps, I think I was looking a the sizes of sizes of the arrays, rather than at the sizes of the arrays !

[William Rose]

@Jamie Al, You wrote " I don't know how to represent Ylm to have dim of X*Y*Z". Why do you want Ylm to have dimensions of X*Y*Z? (I assume you mean 128^3.) Ylm() naturally has two dimensions: theta (polar angle, or colatitude, 0 to pi) and phi (azimuthal angle, 0 to 2pi).

See here for an example of how to compute spherical harmonics on the 2D grid (theta, phi), and plot the results as a nice surface in 3D. By the way, you will want to compute the surface values over the full range of angle [0,pi] and [0,2*pi], so that your surface does not have a hole at the south pole or a gap along the prime meridian. Your code has [0,pi) and [0,2pi).

If you are plotting the radius of a deformed sphere or droplet, then such 2D array is appropriate. You will not need to call isosurface(). The picture of a bumpy droplet which you shared suggests that you will use the spherical harmonics as a relatively small modulation to the droplet radius. Ylm() will go to zero for certain angles, ybut you do not want your radius to go to zero. So you do something like this:

radius=1+0.1*abs(Ylm);  
%array r has same size as Ylm; the radius is modulated with 10% of Ylm.

I assume that Ylm is a 2D array computed as in the example code pointed to above. Array radius will have same size as Ylm. The radius is modulated with 10% of Ylm.

Then you comvert the 2D grid locations of the surface to Cartesian coordinates with sph2cart(). Then you plot them with surf().

When using sph2cart(), note that the elevation angle is defined differently in sph2cart() than the polar angle theta in Matlab's legendre(). That is why, in the example cited above, they use pi/2-theta:

[Xm,Ym,Zm] = sph2cart(phi, pi/2-theta, abs(real(Ylm)));

If you follow the suggestion I gave above for the sphere radius, you would use

[Xm,Ym,Zm] = sph2cart(phi, pi/2-theta, abs(radius));

Here is an interesting posting about legendre() issues.

Answer 1

[William Rose]

@Jamie Al,

[edit: correct spelling errors]

The error is because Ylm has size 61x121, and p0 has size 128x128x128. Therefore you cannot multiply them together with ".*".

By the way, the top chunk of code generates an empty 3D plot. It is empty because you have omitted the final argument to isosurface(), which is "isovalue", i.e. the value in array Test whose surface you want to plot.

The values in Test are all 1's: The min and the max of Test is 1. Since Test=1 everywhere, isosurface() cannot make a surface from it, no matter what isovalue you specify.

The reason Test is all 1's is that is computed with

Test =d + a * exp(A1);

where A1 is a 3D array with 128^3 elements:

A1=(-1. * p0 .* (1 - b .* cos(c * theta))).^p ;

It turns out A1 has values that range from -21000 to -11000. Therefore exp(A1)=0 at all points, therefore Test=d=1 everywhere.

Since Test=1 everywhere, isosurface() cannot make a surface from it.

[Jamie Al]

Thanks, I think I might have forgot to add important lines in my code. Here is the old code and what the plot should look like (this generates a correct isosurface).

correction above code:

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;

New attempt using spherical harmonics:

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));

This returns error:

Error using griddedInterpolant
Grid arrays must have NDGRID structure.

[Walter Roberson]

@William Rose

isosurface(X, Y, Z, V) is defined (at least these days) as calculating a histogram of values to figure out an iso level.

[William Rose]

THank you @Walter Roberson.