Variable layer thickness using loops

Good afternoon,

I am trying to make the graphs as attached. However I keep running into error code:

Error using /

Matrix dimensions must agree.

Error in HW4_inverse3_KVW (line 46)

jj = ceil( sd / dz); %jj is layer that sd(i) resides within

dz is variable layer thickness dz = 1 2 3 (m).

How can I make a code for the variable layer thickness? I do not understand what is going on.

Thank you

% Make VSP G-mat for arbitrary sensor depths and non-uniform layer thicknesses

clear; close all; clc

set(0,'DefaultAxesLineWidth',2,'DefaultLineLineWidth',3)

set(0,'DefaultAxesXminorTick','on','DefaultAxesYminorTick','on')

set(0,'DefaultAxesFontSize',12,'DefaultAxesFontWeight','bold')

M =  8

N  = 3

dz = [ 1 : 1 : N ]                   %variable layer thickness (m)

zt  = cumsum( dz -dz(1) )            %layer tops (m)

zb = zt + dz                         %layer bottoms (m)

zmax = zb(end)                       %depth of model domain (m)

rng(6)

sd = zmax * sort( rand(1,M) )        %sensor depths (m)


figure(1)
plot(ones(1,M), sd,'bsq')
yline( [zt zb(end)], 'r-', 'linewidth',2 )
ylim([0 zmax]); set(gca,'ydir','reverse','Xtick', [ ] )
title('sensors and layers'); ylabel('depth (m)')
title('Sensors and Layers');
ylabel('Depth (m)')

G0 = zeros( M, N );

fprintf('======== single loop simplest method =========\n\n')

%% single loop simplest method (least lines)
for i = 1: M

    jj = ceil( sd(i) / dz);           %jj is layer that sd(i) resides within
    
    G0( i, jj ) = sd(i) - zt( jj );   %assign segment length in jj layer
    
    G0( i, 1: jj - 1) = dz;          %works for jj=1 because 1:0 is not executed
    
    fprintf('i = %d  jj = %d \n', i, jj )
end

fprintf('\n======== single loop logical-vec method =========\n\n')

figure(2)
subplot (2,2,1)
imagesc(G0); colorbar; colormap('cool(20)');  axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G0')

G1 = zeros(M,N);

Answer 1

[Voss]

% Make VSP G-mat for arbitrary sensor depths and non-uniform layer thicknesses

clear; close all; clc

set(0,'DefaultAxesLineWidth',2,'DefaultLineLineWidth',3)

set(0,'DefaultAxesXminorTick','on','DefaultAxesYminorTick','on')

set(0,'DefaultAxesFontSize',12,'DefaultAxesFontWeight','bold')

M =  8

M = 8

N  = 3

N = 3

dz = [ 1 : 1 : N ]                   %variable layer thickness (m)

dz = 1×3

1 2 3

zt  = cumsum( dz -dz(1) )            %layer tops (m)

zt = 1×3

0 1 3

zb = zt + dz                         %layer bottoms (m)

zb = 1×3

1 3 6

zmax = zb(end)                       %depth of model domain (m)

zmax = 6

rng(6)

sd = zmax * sort( rand(1,M) )        %sensor depths (m)

sd = 1×8

0.2502 0.6459 1.9919 2.5128 3.1789 3.5703 4.9274 5.3572

figure(1)
plot(ones(1,M), sd,'bsq')
yline( [zt zb(end)], 'r-', 'linewidth',2 )
ylim([0 zmax]); set(gca,'ydir','reverse','Xtick', [ ] )
title('sensors and layers'); ylabel('depth (m)')
title('Sensors and Layers');
ylabel('Depth (m)')

G0 = zeros( M, N );

fprintf('======== single loop simplest method =========\n\n')

======== single loop simplest method =========

%% single loop simplest method (least lines)
for i = 1: M

    % jj = ceil( sd(i) / dz )           %jj is layer that sd(i) resides within
    jj = find(zt <= sd(i), 1, 'last');
    
    G0( i, jj ) = sd(i) - zt( jj );   %assign segment length in jj layer
    
    % G0( i, 1: jj - 1) = dz;          %works for jj=1 because 1:0 is not executed
    G0( i, 1: jj - 1) = dz( 1: jj - 1);
    
    fprintf('i = %d  jj = %d \n', i, jj )
end

i = 1  jj = 1 
i = 2  jj = 1 
i = 3  jj = 2 
i = 4  jj = 2 
i = 5  jj = 3 
i = 6  jj = 3 
i = 7  jj = 3 
i = 8  jj = 3 

fprintf('\n======== single loop logical-vec method =========\n\n')

======== single loop logical-vec method =========

figure(2)
% subplot (2,2,1)
imagesc(G0); colorbar; colormap('cool(20)');  axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G0')

[Voss]

What's the error? And what line does it happen on?

You can see the code runs here, since it produces the second plot above.

[Walter Roberson]

Right at the very top of Voss's posting there are some tool icons to link the posting or flag it, and so on.

On the line immediately below that in Voss's posting, to the right, it says "Ran in R2023b"

On the right hand side of the line below the "Ran in" there is a gray box followed by the word "Copy". Click that word "Copy". Now go over to your MATLAB editor and use your system Paste facility (for example control-V for Windows) to get a copy of Voss's code in the editor on your system.

[Kathleen]

could you possible tell me what is wrong with the followinging lines? They all should be getting the same outcome

%% vectorized single loop logical-vec method
for i = 1 : M
    
    ip   =  sd(i) - zt >= 0     % logical: layers crossed by i'th ray == true (1)
    
    nl   = sum( ip );           % num layers crossed by ith ray
    
    blrl = sd(i) - zt(nl);      % bottom layer ray length
    
    G1( i ,ip )  = [  ones(1,nl-1)*dz   blrl ];
    
    fprintf('i= %3d nl= %3d sd = %.4f blrl= %.4f\n', i,nl,sd(i), blrl )
end


subplot(2,2,2)
imagesc(G1); colorbar; colormap('cool(20)');  axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G1')

fprintf('\n===========double loop method ======\n\n')


G2 = zeros(M,N);

%% double loop method
for i = 1 : M
    
    for j = 1 : N
        
        d = sd( i ) - zt (j );        %distance between i'th sensor and j'th layer top
        
        if  d <= dz & d >= 0          %sensor within j'th layer
            
            G2( i , j ) = d;
            
        elseif  sd( i ) > zt( j )     %sensor crosses j'th layer
            
            G2( i , j ) = dz;
            
        end                            %if
        
        fprintf('i=%3d j=%3d  sd=%.2f zt=%.2f d= %+.2f G(i,j)=%.2f\n', ...
            i,j,sd(i),zt(j), d, G2(i,j) )
        
    end  %inner loop
    fprintf('\n')
end      %outer loop

[Voss]

The fundamental problem, I think, is that these code snippets treat dz as a scalar, when in fact it is a vector.

[Kathleen]

yeah.... It is the code that our professor gave us to start with, so not sure why he wants us to do it that way. Nor my classmates nor I understand it and non of us know how to solve it :(

[Voss]

See below for changes required. The method for calculating G1 required indexing into dz rather than multiplying dz by a vector of ones (and assigning elements 1 through nl of row i on the left-hand side), and the method for calculating G2 just required replacing dz by dz(j).

% Make VSP G-mat for arbitrary sensor depths and non-uniform layer thicknesses

clear; close all; clc

set(0,'DefaultAxesLineWidth',2,'DefaultLineLineWidth',3)

set(0,'DefaultAxesXminorTick','on','DefaultAxesYminorTick','on')

set(0,'DefaultAxesFontSize',12,'DefaultAxesFontWeight','bold')

M =  8

M = 8

N  = 3

N = 3

dz = [ 1 : 1 : N ]                   %variable layer thickness (m)

dz = 1×3

1 2 3

zt  = cumsum( dz -dz(1) )            %layer tops (m)

zt = 1×3

0 1 3

zb = zt + dz                         %layer bottoms (m)

zb = 1×3

1 3 6

zmax = zb(end)                       %depth of model domain (m)

zmax = 6

rng(6)

sd = zmax * sort( rand(1,M) )        %sensor depths (m)

sd = 1×8

0.2502 0.6459 1.9919 2.5128 3.1789 3.5703 4.9274 5.3572

figure(1)
plot(ones(1,M), sd,'bsq')
yline( [zt zb(end)], 'r-', 'linewidth',2 )
ylim([0 zmax]); set(gca,'ydir','reverse','Xtick', [ ] )
title('sensors and layers'); ylabel('depth (m)')
title('Sensors and Layers');
ylabel('Depth (m)')

%% single loop simplest method (fewest lines)

fprintf('======== single loop simplest method =========\n\n')

======== single loop simplest method =========

G0 = zeros( M, N );

for i = 1: M

    % jj = ceil( sd(i) / dz )           %jj is layer that sd(i) resides within
    jj = find(zt <= sd(i), 1, 'last');
    
    G0( i, jj ) = sd(i) - zt( jj );   %assign segment length in jj layer
    
    % G0( i, 1: jj - 1) = dz;          %works for jj=1 because 1:0 is not executed
    G0( i, 1: jj - 1) = dz( 1: jj - 1);
    
    fprintf('i = %d  jj = %d \n', i, jj )
end

i = 1  jj = 1 
i = 2  jj = 1 
i = 3  jj = 2 
i = 4  jj = 2 
i = 5  jj = 3 
i = 6  jj = 3 
i = 7  jj = 3 
i = 8  jj = 3 

figure(2)
subplot (2,2,1)
imagesc(G0); colorbar; colormap('cool(20)');  axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G0')


%% vectorized single loop logical-vec method

fprintf('\n======== single loop logical-vec method =========\n\n')

======== single loop logical-vec method =========

G1 = zeros(M,N);

for i = 1 : M
    
    ip   =  sd(i) - zt >= 0;    % logical: layers crossed by i'th ray == true (1)
    
    nl   = sum( ip );           % num layers crossed by ith ray
    
    blrl = sd(i) - zt(nl);      % bottom layer ray length
    
    % G1( i ,ip )  = [  ones(1,nl-1)*dz   blrl ];
    G1( i , 1:nl )  = [ dz(1:nl-1)  blrl ];
    
    fprintf('i= %3d nl= %3d sd = %.4f blrl= %.4f\n', i,nl,sd(i), blrl )
end

i=   1 nl=   1 sd = 0.2502 blrl= 0.2502
i=   2 nl=   1 sd = 0.6459 blrl= 0.6459
i=   3 nl=   2 sd = 1.9919 blrl= 0.9919
i=   4 nl=   2 sd = 2.5128 blrl= 1.5128
i=   5 nl=   3 sd = 3.1789 blrl= 0.1789
i=   6 nl=   3 sd = 3.5703 blrl= 0.5703
i=   7 nl=   3 sd = 4.9274 blrl= 1.9274
i=   8 nl=   3 sd = 5.3572 blrl= 2.3572

subplot(2,2,2)
imagesc(G1); colorbar; colormap('cool(20)');  axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G1')


%% double loop method

fprintf('\n===========double loop method ======\n\n')

===========double loop method ======

G2 = zeros(M,N);

for i = 1 : M
    
    for j = 1 : N
        
        d = sd( i ) - zt (j );        %distance between i'th sensor and j'th layer top
        
        if  d <= dz(j) & d >= 0          %sensor within j'th layer
            
            G2( i , j ) = d;
            
        elseif  sd( i ) > zt( j )     %sensor crosses j'th layer
            
            G2( i , j ) = dz(j);
            
        end                            %if
        
        fprintf('i=%3d j=%3d  sd=%.2f zt=%.2f d= %+.2f G(i,j)=%.2f\n', ...
            i,j,sd(i),zt(j), d, G2(i,j) )
        
    end  %inner loop
    fprintf('\n')
end      %outer loop

i=  1 j=  1  sd=0.25 zt=0.00 d= +0.25 G(i,j)=0.25
i=  1 j=  2  sd=0.25 zt=1.00 d= -0.75 G(i,j)=0.00
i=  1 j=  3  sd=0.25 zt=3.00 d= -2.75 G(i,j)=0.00
i=  2 j=  1  sd=0.65 zt=0.00 d= +0.65 G(i,j)=0.65
i=  2 j=  2  sd=0.65 zt=1.00 d= -0.35 G(i,j)=0.00
i=  2 j=  3  sd=0.65 zt=3.00 d= -2.35 G(i,j)=0.00
i=  3 j=  1  sd=1.99 zt=0.00 d= +1.99 G(i,j)=1.00
i=  3 j=  2  sd=1.99 zt=1.00 d= +0.99 G(i,j)=0.99
i=  3 j=  3  sd=1.99 zt=3.00 d= -1.01 G(i,j)=0.00
i=  4 j=  1  sd=2.51 zt=0.00 d= +2.51 G(i,j)=1.00
i=  4 j=  2  sd=2.51 zt=1.00 d= +1.51 G(i,j)=1.51
i=  4 j=  3  sd=2.51 zt=3.00 d= -0.49 G(i,j)=0.00
i=  5 j=  1  sd=3.18 zt=0.00 d= +3.18 G(i,j)=1.00
i=  5 j=  2  sd=3.18 zt=1.00 d= +2.18 G(i,j)=2.00
i=  5 j=  3  sd=3.18 zt=3.00 d= +0.18 G(i,j)=0.18
i=  6 j=  1  sd=3.57 zt=0.00 d= +3.57 G(i,j)=1.00
i=  6 j=  2  sd=3.57 zt=1.00 d= +2.57 G(i,j)=2.00
i=  6 j=  3  sd=3.57 zt=3.00 d= +0.57 G(i,j)=0.57
i=  7 j=  1  sd=4.93 zt=0.00 d= +4.93 G(i,j)=1.00
i=  7 j=  2  sd=4.93 zt=1.00 d= +3.93 G(i,j)=2.00
i=  7 j=  3  sd=4.93 zt=3.00 d= +1.93 G(i,j)=1.93
i=  8 j=  1  sd=5.36 zt=0.00 d= +5.36 G(i,j)=1.00
i=  8 j=  2  sd=5.36 zt=1.00 d= +4.36 G(i,j)=2.00
i=  8 j=  3  sd=5.36 zt=3.00 d= +2.36 G(i,j)=2.36

subplot(2,2,3)
imagesc(G2); colorbar; colormap('cool(20)');  axis tight
set(gca,'Xtick',[1:1:N],'Ytick',[1:1:M])
title('G2')

[Kathleen]

Thank you so much!! I shed many tears over this assignment but I think I understand it now. Thank you so much!

[Voss]

You're welcome! Any questions, let me know.

Answer 2

[Walter Roberson]

sd(i) / dz

sd(i) is a scalar but dz is a row vector.

In MATLAB, the / operator is mrdivide, / matrix right divide. A/B is similar to A * pinv(B) where * is algebraic matrix multiplication. For the / operator to work properly, the number of columns in the numerator and denominator must be the same . You have a scalar numerator, so one column, and your denominiator is a row vector with N columns. Unless N is 1, then the / operator would fail.

If you wanted [sd(i) / dz(1), sd(i) / dz(2), sd(i) / dz(3)] and so on, a collection of individual divisions, then you need to use rdivide, ./

sd(i) ./ dz