where is the error?

Hi, running the following script:

% Create a data set:
x = rand(100,1)*4 - 2;
y = rand(100,1)*4 - 2;
S = x.*exp(-x.^2-y.^2) * 1000;
% Construct the interpolant:
F = TriScatteredInterp(x,y,S);
% 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);
% Define a set of ZI locations
ZI = -500:100:500;
% Find the node above 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
KImat = zeros(size(ZImat));
for i = 1:length(YI)
 for j = 1:length(XI)
  firstIndex = find(BW(i,j,:),1);
  if ~isempty(firstIndex)
  KImat(i,j) = firstIndex;
  end
 end
end
% Create a 3d grid where to compute a value at each node above the surface S
[X,Y,Z] = meshgrid(XI,YI,ZI);
for i = 1:length(YI)
for j = 1:length(XI)
 if KImat(i,j) == 1
  for i = 1:length(Y)
   for j = 1:length(X)
    for k = 1:length(Z)
     T_geotherm(i,j,k) = 18 + 0.003 .* (Z(i,j,k));
    end
   end
  end
 end
end
end

dislay this message: 'Attempted to access Z(1,1,12); index out of bounds because size(Z)=[41,17,11].'

where is the error?

Thanks, Gianluca

[gianluca]

Sorry, no message is displayed. But the result in not a 3d matrix where the values are computed at each node above the surface. I obtain a 1x1x11 matrix.

[Wayne King]

If I run this I do not get any outputs that are 1x1x11. Which variable do you get that is 1x1x11?

[gianluca]

Hi, I cleared the workspace and the message and the variable T_geotherm do not appear. OK, They were an old message and an old variable. Do you have any suggestion how to compute a value at each node of the 3d matrix [X,Y,Z] above the surface S? I hope to be right till the definition of the variable KImat. The last part of the script is wrong

[Wayne King]

The variable, T_geotherm, is not being generated because your if statement if KImat(i,j) == 1 is never true

Answer 1

[Jan]

Matlab code gets much cleaner and easier to debug, when vectorization is applied:

Ugly and due to a missing pre-allocation slow in addition:

for i = 1:length(Y)
  for j = 1:length(X)
    for k = 1:length(Z)
      T_geotherm(i,j,k) = 18 + 0.003 .* (Z(i,j,k));
    end
  end
end

Beautiful and fast:

T_geotherm = 18 + 0.003 .* Z;

Btw, using nested loops with the same counter is not really clean, but not an error:

for i = 1:length(YI)
  for j = 1:length(XI)
    if KImat(i,j) == 1
      for i = 1:length(Y)
        for j = 1:length(X)

It is strongly recommended to avoid such confusing constructions.

[gianluca]

OK, vectorization is better! But I need to define for which nodes and till which nodes along k index compute the variable T_geotherm. The logical matrix BW say me with a zero or an one if the nodes are above or below yhe surface S, KImat say me the first index where it realized. Then, how can I do this? foe example (but it do not work):

for i = 1:length(YI)
 for j = 1:length(XI)
  if KImat(i,j) >= 1
    T_geotherm = 18 + 0.003 .* Z;
  end
 end
end

[Jan]

index = (KImat >= 1);
T_geotherm(index) = 18 + 0.003 .* Z(index);

Answer 2

[Wayne King]

I can run this without error. I think you have some variable in your workspace that is conflicting with this and causing the error.

Can I suggest you first clear the workspace and then try to run this script?