Compute values in the nodes of a 3d matrix

1 view (last 30 days)
Hi,
Let S2 a surface below the surface S1, I would to compute values in the nodes below S2, starting from the values computed at the boundary between S1 and S2, using coefficients C2.
% Create the surfaces:
x = rand(100,1)*4 - 2;
y = rand(100,1)*4 - 2;
S1 = x.*exp(-x.^2-y.^2) * 1000;
S2 = S1 - 100;
% Create positive coefficients
C1 = 0.02 + (0.08 - 0.02).*rand(100,1);
C2 = 0.03 + (0.05 - 0.03).*rand(100,1);
% Construct the interpolant:
F1 = TriScatteredInterp(x,y,S1);
F2 = TriScatteredInterp(x,y,S2);
G1 = TriScatteredInterp(x,y,C1);
G2 = TriScatteredInterp(x,y,C2);
% Evaluate the interpolant at the locations [XI,YI].
XI = -2:0.25:2;
YI = -2:0.1:2;
[XImat,YImat] = ndgrid(XI,YI);
ZI_1 = F1(XImat,YImat);
ZI_2 = F2(XImat,YImat);
C_1 = G1(XImat,YImat);
C_2 = G2(XImat,YImat);
% Define a set of ZI locations
ZI = -500:100:500;
% Find the nodes below the surface S1 and S2
BW1 = bsxfun(@le, ZI_1, reshape(ZI,1,1,[]));
BW2 = bsxfun(@le, ZI_2, reshape(ZI,1,1,[]));
% Create a 3d grid where to compute the value in each node
[Z,Z,Z] = ndgrid(XI,YI,ZI);
% Compute value below S1 with coefficients C1
T = 18 + bsxfun(@times, C_1, Z .* BW1);
Whith the help of a file-exchange find_ndim, I've found the index where the surface S2 starts.
index = find_ndim(BW2,3,'first')
T = ???? + bsxfun(@time, C2, Z .* BW2);
Thanks for any suggestion, Gianluca
  5 Comments
gianluca
gianluca on 19 Sep 2012
hi, the values of T (dependent variable) change as function of Z (independent variable) from a initial value (T = 18) at the surface S1 following a linear trend defined by the coefficient C1 (gradient). When we pass below the surface S2 the linear trend changes from C1 to C2.
T = 18 + bsxfun(@times, C_1, Z .* BW1);
computes in all the nodes of the 3d matrix the variable T (also below S2). I would in a second step compute T only in the nodes below S2 rewriting the old values. In this second step the variable T starts from the values computed previously at surface S2.
index = find_ndim(BW2,3,'first')
Finds the first indexes of S2 in the third dimension, that is the indexes from which the surface S2 begins.
T2 = T(index)
don't works, an error message occurs. I expect that T2 is double and the final T double.
thanks
gianluca
gianluca on 19 Sep 2012
I expect that T2 is 17x41 double and the final T 17x41x11 double.

Sign in to comment.

Accepted Answer

Andrei Bobrov
Andrei Bobrov on 19 Sep 2012
d = cumsum(BW2(:,:,end:-1:1),3);
index = bsxfun(@eq,sum(BW2,3) - 1,d(:,:,end:-1:1) - cumsum(BW2,3));
T2 = T.*index + bsxfun(@times, C_2, Z .* BW2);
  1 Comment
gianluca
gianluca on 19 Sep 2012
Hi, you given me an useful tool with which try again. I would store the results in the same variable T.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!