Compute values in the nodes of a 3d matrix

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

Answer 1

[Andrei Bobrov]

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

[gianluca]

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

Answer 2

[Walter Roberson]

Change

[XImat,YImat] = meshgrid(XI,YI);

to be

[XImat,YImat] = ndgrid(XI,YI);

[gianluca]

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.