音圧分布の計算とプロット

Question

現在私は点音源が作る音圧分布をプロットしようとしています。
そこで以下のようなプログラミングを書きました。

% 定数定義
f = 100;
w = 2*pi*f;
lambda = 340/f;
k = 2*pi/lambda;
p = 1.293;
c = 340;
A = 1/k^2*p*c;
r = 1;
theta = (2 * pi)/4;
% グリッド作成
[X,Y,Z] = meshgrid(-3:0.1:3);
soundposition_x = [];
soundposition_y = []; 
soundposition_z = [];
%主な原因と思われる場所
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%音源の位置計算
x1 = X - r*cos(theta);
%x1 = X; 本来r*cos(theta)はゼロになるため誤差を含めなけらばこれと同等のはず
y1 = Y - r*sin(theta);
z1 = Z;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%音圧計算
soundsorce1 = (1i*w*p*A/sqrt(x1.^2+y1.^2+z1.^2).*exp(-1i*k*sqrt(x1.^2+y1.^2+z1.^2)));
anser = abs(soundsorce1);
xslice = 0;
yslice = [];
zslice = 0;
% プロット部
slice(X,Y,Z,anser,xslice,yslice,zslice)
title('3次元空間の音圧分布')
xlabel('x[m]') 
ylabel('y[m]') 
zlabel('z[m]')
c = colorbar;
shading flat

なのですが、音源の位置を

x1 = X;
y1 = Y - r*sin(theta);
z1 = Z;
%プログラミング中ではr = 1,theta = 2*pi/4

とした時は音圧分布がうまくプロットできるのですが、

x1 = X - r*cos(theta);
y1 = Y - r*sin(theta);
z1 = Z;
%プログラミング中ではr = 1,theta = 2*pi/4

とした時に、うまくプロットできません。

本来r*cos(theta)=0なのでどちらのプログラミングも同じ結果が出るはずなのですが、matlab上では誤差の影響なのか6.1232e-17という計算結果になることはわかりました。
しかしこの誤差がどう悪さをしてプロットがうまくできないのかが解明できていません。
多少の誤差が生じたとしても、音源の位置が本来の場所からちょっとずれて音が広がっていく様子が二つ目の図と同じようにプロットされると思ったのですが、うまく点音源の音圧分布がプロットされず困っています。
もし原因がわかる方がいたらご教授よろしくお願いいたします。

※今回音圧の計算に用いたものは以下の式を変形したものです。
また今後比較的楽に音源の数を増やせるようにするためにcos,sinの計算を使用しています。

Answer 1

x1 = X;で計算すると、x1, y1, z1が同時に0になる場所があり、その場所のanserの値はInfになります。具体的にはanser(41,31,31)です。Infがあるときのsliceの動作の仕様は確認できませんでしたが、実際の動作を見る限り、うまく無視しているようです。

一方x1 = X - r*cos(theta);で計算すると、x1が0にならないため、anser(41,31,31)は、非常に大きいがInfではない値になっています。このため、その他の値は元とほぼ同じなのにも関わらず、相対的に非常に小さい値になってしまうため、プロットが変わっています。試しにanser(41,31,31)=Inf;とすると、同じプロットが得られます。

(41,31,31)は音源の場所だと思いますが、うまく無視する必要があると思います。

Comments

[sibazyun]

わかりやすく、細かいところまで説明ありがとうざいました！！infだった部分がただ非常に大きいだけの値に変わってしまっていたんですね、、。無視する処理を実装したいと思います。
ありがとうございました