function [xb_coords_new,Vdune,Vcut] = make_artificial_dune(xb_coords,crest_pos,z_dune,width,z_excavation,duneslope_sea,duneslope_land,cut)

%function make_artificial_dune(xb_coords,crest_pos,z_dune,width,z_excavation,cut)
%
%Function that simulates the effects of scraping an artificial dune 

warning('off')
if nargin == 2
    z_dune = 3; %Default value of artificial dune height
    width = 2; %Default value of artificial dune crestwidth
    z_excavation = 0; %Default value of minimum excavation depth
    cut = 1;
    duneslope_sea = 0.25;
    duneslope_land = 0.75;
elseif nargin == 3
    width = 2; %Default value of artificial dune crestwidth
    z_excavation = 0; %Default value of minimum excavation depth
    cut = 1;
    duneslope_sea = 0.25;
    duneslope_land = 0.75;
elseif nargin == 4
    z_excavation = 0; %Default value of minimum excavation depth
    cut = 1;
    duneslope_sea = 0.25;
    duneslope_land = 0.75;
elseif nargin == 5
    cut = 1;
    duneslope_sea = 0.25;
    duneslope_land = 0.75;   
end


x = xb_coords(:,1); %XBeach alongshore coordinates
z = xb_coords(:,2);
Icrest = interp1(x,1:length(x),crest_pos,'nearest'); %Indice to round crest_pos to actual position on grid
crest_pos = x(Icrest); %Round round crest_pos to actual position on grid
z_new = xb_coords(:,2);  %z_new is what the new dune profile will be
res = x(Icrest)-x(Icrest-1); %Grid resolution around dune crest
width = ceil(width/res)*res; %Round up width to grid resolution
I = find(x>=crest_pos-width&x<=crest_pos); %Indices of dune crest
J = find(z(I)<z_dune);
z_new(I(J)) = z_dune; %Change dune crest to new height if it is lower than existing height

alpha_dune_sea(1) = duneslope_sea;
alpha_dune_sea(2) = z_new(min(I))-alpha_dune_sea(1)*x(min(I));
alpha_dune_land(1) = -duneslope_land; %Negative because it slopes in the other direction
alpha_dune_land(2) = z_new(max(I))-alpha_dune_land(1)*x(max(I));
z_dummy = polyval(alpha_dune_sea,x);
I1 = find(z_dummy-z<0&z<z_dune);
I1 = max(I1);
z_dummy = polyval(alpha_dune_land,x);
I2 = find(z_dummy-z>0);
I2 = max(I2);
if I2==max(I);
    I2 = I2+1;
end

%I1 = min(I)-ceil(z_dune*2/res); %Assume slope of dune on seaward side is approximately 1:2
%I2 = max(I)+ceil(0.5*z_dune/res); %Assume slope of dune on land side is steeper than 1:2
%I2 = max(I)+1;
alpha_dune_sea = polyfit([x(I1) x(min(I))],[z_new(I1) z_new(min(I))],1); %Polyfit for seaward dune slope
alpha_dune_land = polyfit([x(max(I)) x(I2)],[z_new(max(I)) z_new(I2)],1); %Polyfit for landward dune slope
z_new(I1:min(I)) = polyval(alpha_dune_sea,x(I1:min(I)));
z_new(max(I):I2) = polyval(alpha_dune_land,x(max(I):I2));
Inew = I1:I2; %Include indices at the base of the new dune
Vdune = trapz(x(Inew),z_new(Inew))-trapz(x(Inew),z(Inew)); %Calculate volume of new dune

if cut == 0
    Vcut = 0; %Do not excavate anything
else   
    Ifoot = min(Inew); %Indice for dune foot
    z_foot = z_new(Ifoot); %Elevation at dune foot
    x_excavation_dummy = find_xintercept(x,z_new-z_excavation);
    II = find(x_excavation_dummy<x(Ifoot));
    x_excavation_dummy = max(x_excavation_dummy(II));
    Iexcavation = min(find(x>x_excavation_dummy)); %Indice of where excavation starts
    Viniz = trapz(x(Iexcavation:Ifoot),z_new(Iexcavation:Ifoot));
    Isub = Iexcavation:Ifoot;
    smax = (z_foot-z_new(Iexcavation))./(x(Ifoot)-x(Iexcavation)); %Maximum slope to cut out
    s = 0:0.0001:smax;
    for i = 1:length(s)
        alpha(1) = s(i); %Zlinear = alpha(1)*x + alpha(2)
        alpha(2) = z_foot-s(i)*x(Ifoot); %alpha(2) = z_footIexcavation - alpha(1)*xpiede
        Zlinear = polyval(alpha,x(Isub)); %Cut line from dune foot with slope s
        J = find(z_new(Isub)>Zlinear); %Indices of all elevations greater than cut line
        z_new(Isub(J)) = Zlinear(J); %Change all existing elevations greater than cut line to the cut line
        Vcut = Viniz - trapz(x(Isub),z_new(Isub)); %Calculate volume of cut
        if Vcut < Vdune
            continue
        else
            break
        end
    end
    
    if Vcut<Vdune %If more volume still needs to be cut
        dmax = z_foot-z_new(Iexcavation);
        b = z_foot-alpha_dune_sea(1)*x(Ifoot);%y = mx+b, b = y-mx
        Zlinear2 = polyval(alpha_dune_sea,x(Isub)); %Doesn't change through the loop
        for d = 0.005:0.005:dmax
            xpivot = ((z_foot-d)-b)/alpha_dune_sea(1); %x = (y-b)/m
            J = find(z_new(Isub)>Zlinear2&x(Isub)>=xpivot);
            z_new(Isub(J)) = Zlinear2(J);
            alpha1 = polyfit([x(Iexcavation) xpivot],[z_new(Iexcavation) (z_foot-d)],1);
            Zlinear1 = polyval(alpha1,x(Isub));
            J = find(z_new(Isub)>Zlinear1&x(Isub)<=xpivot);
            z_new(Isub(J)) = Zlinear1(J);
            Vcut = Viniz - trapz(x(Isub),z_new(Isub));
            if Vcut < Vdune
                continue
            else
                break
            end
        end
    end
end

xb_coords_new = [x z_new];
