Variable layer thickness using loops
1 次查看(过去 30 天)
显示 更早的评论
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);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1484107/image.png)
0 个评论
采纳的回答
Voss
2023-9-15
% 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
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
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')
8 个评论
更多回答(1 个)
Walter Roberson
2023-9-15
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Orange 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!