Array indices must be positive integers or logical values.
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
This code plots an idealized cross section of a delta (geologic feature). It works when the shoreline (SH) transgression (lines 71-74) is commented out, however this doesn't allow the shoreline to transgress (which doesn't happen in real life). The lines that allow for shoreline transgression work only when line 55 is commented (Z).
Any thoughts on how I can remedy this error message? My original idea is to have an initial delta present at the start of the model run, but I am not sure how to write that into the code. Perhaps I need to adjust where the initial shoreline is to do that?
Thank you!
%% Input parameter values
Rab=0.3; %Sediment input
A=0.2; %sea level amplitude
B=1; %sea level frequency
%% setting up the domain
nx=1000; %length of delta
p = 2; %starting position
pos=p;
dx=.01;
%% Time vector
dt=.00005;
Tmax=10;
t=0:dt:Tmax;
nt=length(t);
%% Initialize time vectors
sc=zeros(1,nt); % cycle Shoreline position
rc=zeros(1,nt); % cycle Alluvial-basement position
sm=zeros(1,nt); % Fixed SL Shoreline position
rm=zeros(1,nt); % Fixed SL Alluvial-basement position
Zp=zeros(1,nt); % SL position
Vc=zeros(1,nt); % Volume
%% Initialize space vectors
x=zeros(1,nx);
%cycle values
Fc=zeros(1,nx); %sediment flux
Hc=zeros(1,nx); %Enthalpy
Hnewc=zeros(1,nx);
hc=zeros(1,nx); %height above current sea-level
% To define the dimension of F(flux)
for i=1:p
x(i)=(i-(p-0.5))*dx;
hc(i) = -x(i);
hm(i) = -x(i);
end
for i=p:nx
x(i)=(i-(p-0.5))*dx;
end
Z = 0;
Z0=1;
Fc(1) = Rab;
%Cycle time loop
for j=1:nt
%% Sea Level
Z = Z0 + A*sin(j*dt*B);
Zp(j) = Z;
sc(j)=(pos-p)*dx - Hc(pos)*dx/(-Z); %shoreline position
for i=p:nx-1
Fc(i) = min((hc(i) - hc(i+1))/dx, Hc(i)*dx/dt+Fc(i-1));
%% Exner equation
Hnewc(i) = Hc(i) + dt/dx*(Fc(i-1)-Fc(i));
%tracks SH regression
if Hnewc(pos)> Z
pos = pos + 1;
end
%tracks SH transregression
while Hc(pos-1) < Z
pos = pos-1;
end
Hc(i)=Hnewc(i);
hc(i) = max(Hc(i)-Z,0);
end
%% Counter
if mod(j,.5/dt) == 0
disp(['t=', num2str(j*dt)]);
end
%Checks if trajectories hit boundary of domain
if j>1000 %skips roughness at first few steps
if pos == nx-1 || rc(j)/dx >= p || abs(rc(j) - rc(j-1)) >= 100*dx
disp('end of domain reached')
break
end
end
end
plot(t,sc, 'b')
hold on
plot(t,Zp)
hold on
1 Kommentar
Walter Roberson
am 22 Feb. 2020
Line 66 you use Hc(pos-1) but pos can be 1. You do not know that it is certain that Hc(pos-1)<Z will become false so you might decrease pos even though you are at the boundary
Antworten (1)
Vimal Rathod
am 24 Feb. 2020
Hi,
Using the following lines decreases the pos value to less than 1, Thus the MATLAB throws an error. From your given code, all the Hc values have 0 value and Z value is 1 (Thus the Error). Commenting the following loop makes sure the error doesn't occur and if you comment the line 55 (Where Z is initialized) Z will always have a value of zero (Thus the loop won't run) which is same as commenting the loop.
%tracks SH transregression
while Hc(pos-1) < Z
pos = pos-1;
end
You could put a certain boundary check that to make sure the pos value doesn't fall behind the boundary.
You could also debug your code, check the values by using breakpoints and verifying if your values are appropriate.
Refer to the following link to know more about debugging using breakpoints.
0 Kommentare
Siehe auch
Kategorien
Mehr zu Geology finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!