MATLABのTips 3次元空間に結果を表示する方法

2021.07.03

はじめに

MATLABって3次元空間に結果を表示する関数がないんですよね。。。。
“surf”や”plot3″、”scatter3″があるじゃんって思う方もいるかもしれませんが、これらの関数は2次元空間のz軸上に結果を表示しているか、3次元空間上に散布図や線を表示しているにすぎません。

「3次元空間に結果を表示する」というのは、例えば(x1,y1,z1)で値がp1、(x2,y2,z2)で値がp2……といった結果を図示することを意味します。

もっと具体的には、3次元モデルの有限要素解析の結果を表示するような関数がありません。

これって機械系の人には非常に不便なんですよね。。。。。
そんなことを感じまして、学生時代に3次元空間上に結果を表示するプログラムを作成しました。

ちなみに、学生時代も現在もMATLAB2007を使用しているので、全バージョンで動くと思います。
（MATLAB2006はWindows XPまでしか動かないので、もう使ってる人いないですよね？）

3次元空間での結果表示の考え方

（MATLABコードは一番下に載せているので、コードだけ知りたいという方は読み飛ばしてください。）

まず、考え方を説明します。
FEMソフトのように3次元空間に結果を表示させたい場合、(x1,y1,z1)に値p1を表示することになります。つまり、4変数を表示させるということになります。

そのため”surf”は使えません。”surf”はXY平面上に値をZ軸方向に表示しているので、3変数しか扱えません。

学生当時に候補に挙がった関数は”plot3″です。plot3に(x1,y1,z1)の3変数を3次元空間にプロットして、値を色で表現することで4変数を扱うことができます。

コードの検証

せっかくなので、過去の記事を参考に、3次元空間上に多数の点音源を設けたときの音圧分布を求めたいと思います。

MATLABで学ぶ音響工学　球面波2
2次元での球面波のプログラム方法を紹介しています。
mech-eng-uni.com
2021.03.15

音源を下記のように10個設けました。

% % % 音源位置　適当に10個ぐらい
input=[
    -10 0 0
    -10 -5 -5
    -10 -5 +5
    -10 +5 -5
    -10 +5 +5
    +10 0 0
    +10 -5 -5
    +10 -5 +5
    +10 +5 -5
    +10 +5 +5
    ];

X,Y,Zで-10～10の範囲で音圧分布を求めると下図のようになりました。

音圧分布がおかしいような気もしますが、3次元空間に結果を表示するという目標は達成できました。

MATLABのコード

学生時代に作成したコードなので、使い勝手のいいように各自で編集してください。

functionファイル

function plot3color(Z1,color_lenge,geometory,width,linestyle)
% geometory(x,y,z)上にZ1を表示する
% なお、色はZ1の最小値～最大値で表示する仕様になっている。変更したければ10行目を修正ください。
% 

colorrgb=color_lenge; %色を何分割で表現するか
J=jet(colorrgb);
RANK=zeros(1,length(Z1));
for nz=1:length(Z1)
Z=[Z1(nz) linspace(min(Z1),max(Z1),color_lenge-1)];
rank=tiedrank(Z);
RANK(nz)=rank(1);
end

for n=1:length(Z1)
    RGB=J( ceil(RANK(n)),: );
    plot3(geometory(n,1),geometory(n,2),geometory(n,3),linestyle,'Color',RGB,'linewidth',width)
    hold on
end  
 
end

 

実行ファイル

clear all;close all;clc

% % % 3次元空間
x=-10:1:10;
y=-10:1:10;
z=-10:0.25:10;

% % % 音源位置　適当に10個ぐらい
input=[
    -10 0 0
    -10 -5 -5
    -10 -5 +5
    -10 +5 -5
    -10 +5 +5
    +10 0 0
    +10 -5 -5
    +10 -5 +5
    +10 +5 -5
    +10 +5 +5
    ];

        
rho=1.293;
A=1;
c=340; %音速
freq=100; %周波数
w=2*pi*freq;
k=w/c;



P=zeros(length(x)*length(y)*length(z),1);
geo=zeros(length(x)*length(y)*length(z),3);
count=0;
for ii1=1:length(x)
    for ii2=1:length(y)
        for ii3=1:length(z)
            count=count+1;
            geo(count,:)=[x(ii1) y(ii2) z(ii3)];
            for ii4=1:size(input,1)
                r=sqrt( (x(ii1)-input(ii4,1))^2+(y(ii2)-input(ii4,2))^2+(z(ii3)-input(ii4,3)) );
                if r<0.1
                    r=0.1; %r=0の時にP(ii1,ii2)が無限大になるのを防ぐ目的
                end
                P(count,1)=P(count,1)+j*w*rho*A./(4*pi*r).*exp(j*(k*r));                
            end
        end
    end
end

figure
plot3color(10*log10(abs(P)),100,geo,1,'o')
axis equal
xlabel('X')
ylabel('Y')
zlabel('Z')
view([45 45])

 

 

コメント