Step1
interp1が適切かと思われます。
図を参考にすると、x = 0:10:80; z = 0:10:150;とサイズが違うからです。
interp1でxを補間し、meshgridで n×mこの行列を作るのがよいと思われます。
--> meshgridについて
Step2
z軸中心に回転するのであればcylinder関数が使えます。
--> cylinder関数について
Step3
パラメータがないので、磁場強度を適当な関数で行っています。
まずは円筒を作ります。
clear,clc,close all;
% 最大半径40mm,中心(40,40),高さ150mmの円筒
xlen = 80;
z = 0:5:150;
% 半径1~40の円をつくる
r = 1:40;
for i = 1:length(r)
[xx,yy,zz] = cylinder(r(i));
X(i,:) = xx(1,:);
Y(i,:) = yy(1,:);
end
clear i r xx yy zz
% 中心(0,0)から(40,40)にシフト
X = X + 40;
Y = Y + 40;
次に磁場強度を作ります。
ポイントはmeshgridでf(x,z)の関数にすることです。
cnt = 1;
while cnt <= length(X(:,1))
x = X(cnt,:);
% point f(x,z)の関数にする
[xt,zt] = meshgrid(x,z);
rxt = 110.*(1-sin(2*pi*xt./xlen/2)); % xについてはここをinterp1
rzt = 250.*(1-sin(2*pi*zt./max(z)/2)); % zについてはここをinterp1
mesh=rzt.*rxt;
str(:,:,cnt) = mesh';
cnt = cnt+1;
end
最後に次元を合わせます。
ここでは、(θ,r) -> (x,y)にしているので、21×半径の数×高さの次元に直してます。
コードとしては正直汚いです。
Z = repmat(z,length(X(1,:)),1)';
Z = repmat(Z,1,1,length(X(:,1)));
X = X'; Y = Y';
X = permute(repmat(X,1,1,length(z)),[1 2 3]);
Y = permute(repmat(Y,1,1,length(z)),[1 2 3]);
Z = permute(Z,[2 3 1]);
str = permute(str,[1 3 2]);
clear x z xt zt rxt rzt mesh cnt
今回はsurfで円盤を積み重ねます。
x軸上＆z軸上の磁場強度を見る場合はplot関数を使うと可能です。
subplotを使うと並記できるのでお勧めです。
--> subplotについて
cnt = 1;
while cnt <= length(X(1,1,:))
surf(X(:,:,cnt),Y(:,:,cnt),Z(:,:,cnt),str(:,:,cnt)), hold on
xlabel('X[mm]'); ylabel('Y[mm]'); zlabel('Z[mm]');
cnt = cnt +1;
end
ylim([40 80]);
