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

1 view (last 30 days)
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:
Maybe my functional form is incorrect??
  5 Comments
Walter Roberson
Walter Roberson on 7 Nov 2022
Opps, I think I was looking a the sizes of sizes of the arrays, rather than at the sizes of the arrays !
William Rose
William Rose on 7 Nov 2022
@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));

Sign in to comment.

Answers (1)

William Rose
William Rose on 7 Nov 2022
Edited: William Rose on 7 Nov 2022
[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 Comments
Walter Roberson
Walter Roberson on 7 Nov 2022
isosurface(X, Y, Z, V) is defined (at least these days) as calculating a histogram of values to figure out an iso level.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!