データ点を補完して３次元マップを作成したい

Matlabを用いて下の図のような半円柱内の磁場分布を示す３次元マップを作成したいです。その過程の補完・マップ作成について教えていただきたいです。

高さをz方向として、原点を通るxz平面に関して格子点１cm間隔の磁場データを取得してあります。データはエクセルに記録し、位置x,y(=0),zとその地点の磁場強度Bを[x,0,z,B]として４× n の形で記録しています。（原点は円柱の中心にしています）

このデータを出発点に以下の操作をしたいです。

step1.xz平面の磁場データの補完（内挿）

　　　１cm間隔のデータを線形補完して滑らかな分布を示すデータにしたいです

　　　interp2を使用することを考えましたがどのように引数を設定すれば良いかわからず詰まっています。

step2.3次元への拡張

　　　原点からの距離と高さが等しい点は同じ磁場強度として扱うことで平面の磁場分布データを立体に拡張したいです

　　　(得られた平面のデータを、z軸（高さ方向の軸）を回転の軸としてまわすことで３次元に拡張するイメージです)

step3.マップ表示

　　　データから図のように色で磁場強度を示す半円柱の３次元マップを作成したいです

　　　またx軸上およびz軸上の磁場強度を表すグラフも一緒に作成したいです（まさに図のような形で）

どこか一部分だけでも構いませんのでご教授いただけると幸いです。

Answer 1

[Hernia Baby]

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]);

[kou sasai]

ご丁寧にありがとうございます。わかりやすく、初心者の私でも理解しやすかったです。while内について１つ質問です。f(x,z)の関数を作る=座標(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;

rxt、rztは仮で作成していただいた物ですがこれをどのように置き換えればいいのでしょうか？

x z B

0 0 2

1 0 3

2 0 6

3 0 9

4 0 13

という形で位置と磁場強度が保存されているのですが、置き換え方が分からないです。

[Hernia Baby]

xは等間隔と仮定した場合、reshape関数で磁場強度をmeshにしてinterp2をかけるとよさそうです。

reshapeまでの方法を書いておきます。

clc,clear;
x = repmat([1:20],1,5);
z = repelem([1:5],20);
M = randi(10,[1,100]);
table(x',z',M','VariableNames',{'x','z','B'})

ans = 100×3 table
    x     z    B
    __    _    _

     1    1    4
     2    1    4
     3    1    8
     4    1    5
     5    1    6
     6    1    2
     7    1    6
     8    1    9
     9    1    1
    10    1    5
    11    1    2
    12    1    4
    13    1    5
    14    1    9
    15    1    5
    16    1    9

meshの形にします

xu = unique(x);
zu = unique(z);

Mu = reshape(M,[],length(zu))

Mu = 20×5

4 3 9 6 8 4 10 2 5 2 8 5 6 9 7 5 7 4 9 9 6 4 5 2 5 2 1 4 4 1 6 4 6 1 1 9 6 5 2 4 1 1 5 5 4 5 2 6 8 5

surf(zu,xu,Mu,'FaceAlpha',0.5)

[kou sasai]

疑問全て解決しました。ご丁寧に対応してくださりありがとうございます。