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

Question

Jamie Al 2022-11-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1844913-error-using-spherical-harmonics-array-dimensions-must-match-for-binary-array-op

评论： William Rose 2022-11-7

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??

5 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2022-11-7

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

William Rose 2022-11-7

在 MATLAB Online 中打开

@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 2022-11-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1844913-error-using-spherical-harmonics-array-dimensions-must-match-for-binary-array-op#answer_1092648

编辑：William Rose 2022-11-7

@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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Jamie Al 2022-11-7

在 MATLAB Online 中打开

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 2022-11-7

@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 2022-11-7

THank you @Walter Roberson.

请先登录，再进行评论。

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

回答（1 个）

3 个评论
显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

5 个评论 显示 3更早的评论隐藏 3更早的评论

回答（1 个）

3 个评论 显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论