在matlab中绘制球谐函数

Posted 2023-02-21

【中文标题】在matlab中绘制球谐函数

【英文标题】：Plotting spherical harmonics in matlab

【发布时间】：2015-03-05 07:44:54

【问题描述】：

我想在笛卡尔坐标中绘制零阶球谐函数，但 Matlab 的输出不同，不是球面。

x11 = linspace(-1,1,100);
x22 = linspace(-1,1,100);
x33 = linspace(-1,1,100);

[x1 ,x2, x3]= meshgrid(x11,x22,x33);

G = (exp((-1)*(x1.^2 + x2.^2 + x3.^2)));
isosurface(G);

有人可以建议我错在哪里，如果可能的话，请告诉我如何绘制高阶球谐函数。谢谢你。

【问题讨论】：

我不太确定，但你不想用linspace 代替logspace吗？

我同意@knedlsepp，你应该使用linspace 来覆盖沿x、y、z 轴延伸[-1,1] 的立方体。此外，在isosurface 中，您必须指出要在其中绘制事物的isovalue，例如isosurface(x1, x2, x3, G, 0.5)（即绘制G == 0.5 的等值面）

对于更高的订单，我真心建议您购买蓝皮书 (Hansen) 作为必备参考。请参阅附录 A.1.2，了解直接在球坐标 （等值面） 中的公式和递推公式。不幸的是，我不允许在 matlab 中为您提供我的实现，无论如何周围有很多代码（例如在 fileexchange 上）。

会不会它看起来不像一个球体？有时绘图会扭曲外观，因为不同的暗淡在屏幕上的缩放比例不同。

@CitizenInsane 感谢您的回答。

【参考方案1】：

使用 MATLAB 的球谐函数

degree=6;
order=0;
grid=40;
radius=5;

% Create the grid
delta = pi/grid;
theta = 0 : delta : pi; % altitude
phi = 0 : 2*delta : 2*pi; % azimuth
[phi,theta] = meshgrid(phi,theta)

% Calculate the harmonic
Ymn = legendre(degree,cos(theta(:,1)))
Ymn = Ymn(order+1,:)';
yy = Ymn;
for kk = 2: size(theta,1);
    yy = [yy Ymn]
end;
yy = yy.*cos(order*phi);

order = max(max(abs(yy)));
rho = radius + 2*yy/order;

% Apply spherical coordinate equations
r = rho.*sin(theta);
x = r.*cos(phi);  % spherical coordinate equations
y = r.*sin(phi);
z = rho.*cos(theta);

% Plot the surface
clf
surf(x,y,z)
light
lighting phong
axis tight equal off
view(40,30)
camzoom(1.5)

这给出了：

【参考方案2】：

首先，我不建议您先评估立方体内所有可能坐标的球谐函数，然后使用isosurface 来绘制事物（我什至认为您将isosurface 误解为以某个常数切割数据半径（绝对不是这种情况，请参阅documentation）。

最好绘制球谐函数，是在球坐标(r, phi, theta) 中使用它们的公式。您可以找到一些模式here 的其中一些公式。公式仅针对角度部分提供，径向部分取决于您的领域。

以l=1、m=-1 为例，您可以像这样在（方位角、仰角）网格上生成这个谐波：

azimuths = linspace(0, 360, 361) * pi / 180;
elevations = linspace(0, 180, 181) * pi / 180;

[A, E] = ndgrid(azimuths, elevations);
H = 0.25 * sqrt(15/(2*pi)) .* exp(-1j*A) .* sin(E) .* cos(E);

然后您可以像这样将网格转换回笛卡尔网格：

X = cos(A) .* sin(E);
Y = sin(A) .* sin(E);
Z = cos(E);

您还可以添加一些径向扭曲以使事情看起来更好：

Data = abs(imag(H));
minData = min(Data(:));
maxData = max(Data(:));
Distord = (Data - minData)/(maxData-minData);

X = Distord .* cos(A) .* sin(E);
Y = Distord .* sin(A) .* sin(E);
Z = Distord .* cos(E);

surf(X, Y, Z, Data);
shading flat;

这给出了：

【参考方案3】：

将axis equal 添加到您的代码底部，这将导致您的坐标区长度相等（以像素为单位）。现在你的表面看起来像一个球体！

【讨论】：

我用过axis equal。但我看起来只是一些等值的球体。

@Wouster Kuijsters 谢谢。我已经将axis equal 用于某些等值，它看起来像球体，但不适用于其他值。可能是情节的形状也取决于我使用的等值。

@YugeshEgala 当然，所有低于exp(-1) 的等值都不会看起来像一个球体，而是它的一部分（因为它们超出了半径为 1.0 的球体，即立方体的一侧）......并且您不会在exp(-3) 下方绘制任何内容...目前还不清楚您的期望是什么？