Compute values in the nodes of a 3d matrix

3 views (last 30 days)
Hi,
I've to compute values in the nodes of a 3d matrix. The coefficients vary along i and j index (X and Y direction). They vary also along the k index (Z direction) as function of the node position (above or below a surface).
% Create the surface:
x = rand(100,1)*4 - 2;
y = rand(100,1)*4 - 2;
S = x.*exp(-x.^2-y.^2) * 1000;
% Create positive coefficients
coeff = 0.02 + (0.08 - 0.02).*rand(100,1);
% Construct the interpolant:
F = TriScatteredInterp(x,y,S);
G = TriScatteredInterp(x,y,coeff)
% Evaluate the interpolant at the locations [XI,YI].
XI = -2:0.25:2;
YI = -2:0.1:2;
[XImat,YImat] = meshgrid(XI,YI);
ZImat = F(XImat,YImat);
C = G(XImat,YImat)
% Define a set of ZI locations
ZI = -500:100:500;
% Find the nodes below the surface S
BW = false(length(YI),length(XI),length(ZI));
for i = 1:length(YI)
for j = 1:length(XI)
BW(i,j,:) = ZImat(i,j) > ZI;
end
end
% Create a 3d grid where to compute the value in each node below the surface S
[X,Y,Z] = meshgrid(XI,YI,ZI);
index = BW >= 1;
T_geotherm = zeros(size(X));
% If I use a single coefficient (e.g. 0.003), this statement works
T_geotherm(index) = 18 + 0.003 .* Z(index);
But I would use the coefficients defined in C, with a command of the type:
for i = 1:length(XI)
for j = 1:length(YI)
T_geotherm(index) = 18 + C(i,j) .* Z(index);
end
end
Error message: Attempted to access C(1,18); index out of bounds because size(C)=[41,17].
Thanks for any suggestion, Gianluca

Accepted Answer

Andrei Bobrov
Andrei Bobrov on 18 Sep 2012
Please, try this is code:
% Create the surface:
x = rand(100,1)*4 - 2;
y = rand(100,1)*4 - 2;
S = x.*exp(-x.^2-y.^2) * 1000;
% Create positive coefficients
coeff = 0.02 + (0.08 - 0.02).*rand(100,1);
% Construct the interpolant:
F = TriScatteredInterp(x,y,S);
G = TriScatteredInterp(x,y,coeff);
% Evaluate the interpolant at the locations [XI,YI].
XI = -2:0.25:2;
YI = -2:0.1:2;
[XImat,YImat] = ndgrid(XI,YI);
ZImat = F(XImat,YImat);
C = G(XImat,YImat);
% Define a set of ZI locations
ZI = -500:100:500;
BW = bsxfun(@gt,ZImat,reshape(ZI,1,1,[]));
[Z,Z,Z] = ndgrid(XI,YI,ZI);
T_geotherm = 18 + bsxfun(@times,C,Z.*BW);
  1 Comment
gianluca
gianluca on 18 Sep 2012
Hi Andrei, your code works, the values are computed in the T_geotherm matrix. I try to understand how the code works...
BW = bsxfun(@gt,ZImat,reshape(ZI,1,1,[]));
[Z,Z,Z] = ndgrid(XI,YI,ZI);
T_geotherm = 18 + bsxfun(@times,C,Z.*BW);
these three lines are unclear for me. If you had the time to give me some explanations, I would be very grateful.
Thanks very much

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 17 Sep 2012
Change
[XImat,YImat] = meshgrid(XI,YI);
to be
[XImat,YImat] = ndgrid(XI,YI);
  1 Comment
gianluca
gianluca on 18 Sep 2012
thanks, but it don't works. I've changed in
[XImat,YImat] = ndgrid(XI,YI);
With the command
for i = 1:length(XI)
for j = 1:length(YI)
T_geotherm(index) = 18 + C(i,j) .* Z(index);
end
end
I get the same results: in the T_geotherm matrix where I would compute the values there are NaN.

Sign in to comment.

Categories

Find more on Resizing and Reshaping Matrices in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!