%% Comments
% Changed defenition of last_point_below_zs0 for the cases in which dune
% height drops below zs0 after the crest

%% Function

function [s s0] = compute_bedslope_berm_based_on_accretion_zone(x,z0,zend,zs0)
    
    plot_result = 0;    % Only use this [1] while finetuning this function (as it will generate a figure and pause at the end of this function)
      
    % Step 1: determine x_coordinate SSL
    if zend(end)<zend(1)
        s = 0;
        return
        %error('Assuming the dunes to be at the right side of the domain');
    end
%     last_point_below_zs0 = max(find(zte<zs0));
    last_point_below_zs0 = find(zend>zs0,1, 'first')-1; 
    x_SSL = interp1(zend(last_point_below_zs0:(last_point_below_zs0+1)),x(last_point_below_zs0:(last_point_below_zs0+1)),zs0);
    
    % Step 2: determine which cells to consider
    % Top point: SSL point
    z_to_consider = zend;
    z_to_consider(min(find(zend>zs0)):end) = NaN;    
    % Bottom point, end of the accretion zone
    % Here defined as the lowest bed level where still 50% of the maximum accretion occured (excluding first point)
    max_accretion = max(zend-z0);
    z_to_consider(1:(min(find(zend-z0>0.5*max_accretion)))) = NaN;
    
    % Step 3: determine the average slope of the remaining cells   
    % First interpolate to a uniform (tight) grid to ensure fair averaging
    new_x = x(min(find(~isnan(z_to_consider)))):min(diff(x))/10:x(max(find(~isnan(z_to_consider))));
    z_new = interp1(x,z_to_consider,new_x);
    slopes = (z_new(2:end)-z_new(1:end-1))./(new_x(2:end)-new_x(1:end-1));
    s = mean(slopes);
    
    % ADDITIONAL: bed slope on initial profile from the bottom of the
    % accretion zone to the original SSL intersect --> s0
    z_to_consider_for_s0 = z0;
    z_to_consider_for_s0(1:min(find(~isnan(z_to_consider)))-1) = NaN; % Remove cells left of deposition zone
    z_to_consider_for_s0(min(find(z0>zs0)):end) = NaN;
    new_x_for_s0 = x(min(find(~isnan(z_to_consider_for_s0)))):min(diff(x))/10:x(max(find(~isnan(z_to_consider_for_s0))));
    z_new_for_s0 = interp1(x,z_to_consider_for_s0,new_x_for_s0);
    slopes_for_s0 = (z_new_for_s0(2:end)-z_new_for_s0(1:end-1))./(new_x_for_s0(2:end)-new_x_for_s0(1:end-1));
    s0 = mean(slopes_for_s0);

    if plot_result == 1 
        fig = figure;
        plot(x,z0,'k.-');
        hold on
        plot(x,zend,'.-');
        hold on
        plot(x,z_to_consider,'.-r')
        plot(x,x.*0+zs0)
        plot(x_SSL,zs0,'ok')
        hold on
        plot([x_SSL x(max(find(~isnan(z_to_consider))))]  ,[zs0 z_to_consider(max(find(~isnan(z_to_consider))))],'k')
        title('Slope at t')
        
        
        fig = figure;
        plot(x,z0,'k.-');
        hold on
        plot(x,zend,'.-');
        hold on
        plot(x,z_to_consider_for_s0,'.-r')
        plot(x,x.*0+zs0)
        plot(x_SSL,zs0,'ok')
        hold on
        plot([x_SSL x(max(find(~isnan(z_to_consider_for_s0))))]  ,[zs0 z_to_consider_for_s0(max(find(~isnan(z_to_consider_for_s0))))],'k')
        title('Slope at t0')
        
        pause
    end
end