Variable layer thickness using loops

1 visualizzazione (ultimi 30 giorni)
Kathleen il 15 Set 2023
Modificato: Voss il 16 Set 2023
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);
0 CommentiMostra -1 commenti meno recentiNascondi -1 commenti meno recenti

Accedi per commentare.

Risposta accettata

Voss il 15 Set 2023
% 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')
8 CommentiMostra 7 commenti meno recentiNascondi 7 commenti meno recenti
Voss il 16 Set 2023
Modificato: Voss il 16 Set 2023
You're welcome! Any questions, let me know.

Accedi per commentare.

Più risposte (1)

Walter Roberson il 15 Set 2023
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
0 CommentiMostra -1 commenti meno recentiNascondi -1 commenti meno recenti

Accedi per commentare.

Categorie

Scopri di più su Orange in Help Center e File Exchange

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by