Array indices must be positive integers or logical values.

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

[Walter Roberson]

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

Answer 1

[Vimal Rathod]

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.

https://www.mathworks.com/help/matlab/matlab_prog/debugging-process-and-features.html