function [Continuity recoveryprofilet recoveryprofileall widthIsland VolIsland HeightIsland distance seatransl MaxHeight]=recoverprof1910(alpha,recoveryTime,stormprofile,equilibriumprofile,MSLevel,wavebase,Tide,theta,T,Beta,xas,sealevelchange,add_sediment)
inclEquilibriumIsland=0;
MSL.total=MSLevel*ones(1,length(xas));

%% Definition of the different sections

[IntersectA IntersectB intersections]=CalcintersectionsArea(equilibriumprofile,stormprofile,MSLevel,wavebase,xas);

[r1 c1]=find(xas>=intersections.B.wavebase(1,1),1);
[hr1 hc1]=max(equilibriumprofile);
[hhr hcc]=find(intersections.A.B(1,:)<=xas(hc1));
[hr2 hc2]=max(intersections.A.B(2,1:hcc(end)));%>=intersections.A.MSL(1,1),1);
[r2 c2]=find(xas>=intersections.A.B(1,hc2),1);
j=1;
while c2<=c1
    [r2 c2]=find(xas>=intersections.A.B(1,hc2-1+j),1);
    j=j+1;
end

% [r2 c2]=find(xas>=intersections.B.MSL(1,1),1);
[r3 c3]=find(xas>=intersections.B.MSL(1,2),1);
[xy xz]=find((intersections.A.B(1,:)>=intersections.A.MSL(1,2)),1);
x.A.B.bay=c3;


equilibrium.deep.y=stormprofile(1:(c1-1));
equilibrium.deep.x=xas(1:(c1-1));
equilibrium.gulf.y=equilibriumprofile(c1:(c2-1));

lastone=find(equilibrium.gulf.y>intersections.A.MSL(1,1),1)-1;
equilibrium.gulf.y_part1=equilibrium.gulf.y(1:lastone);
equilibrium.gulf.y_part2=equilibrium.gulf.y((lastone+1):end);

equilibrium.gulf.x=xas(c1:(c2-1));
equilibrium.island.y=equilibriumprofile(c2:(c3-1));
equilibrium.island.x=xas(c2:(c3-1));
equilibrium.bay.y=equilibriumprofile(c3:length(xas));
equilibrium.bay.x=xas(c3:length(xas));
equilibrium.sediment.x=xas(x.A.B.bay:length(xas));
MSL.sediment=MSL.total(x.A.B.bay:length(xas));

%% Equilibrium profile

dh_add_sediment=(add_sediment/12)/(xas((c2-1))-xas(c1)); %layer of sediment added per month

%post-storm profile
storm.gulf=stormprofile(c1:(c2-1));
storm.island=stormprofile(c2:(c3-1));
storm.bay=stormprofile(c3:length(xas));
storm.sediment=stormprofile(x.A.B.bay:length(xas));

recovery.gulf(1,:)=storm.gulf;
recovery.gulf_part1(1,:)=storm.gulf(1:lastone);
recovery.gulf_part2(1,:)=storm.gulf((lastone+1):end);
recovery.island(1,:)=storm.island;
recovery.bay(1,:)=storm.bay;
recovery.sediment(1,:)=storm.sediment;
recovery.gulfnegativesediment(1,:)=storm.gulf;

%recovery in time
for t=1:recoveryTime
    
    % the gulf-side of the island recovery depend on the distance
    % out of equilibrium. This is calculated by making use of the function
    % thetarule2.
    
    %recovery dh/dt
    if dh_add_sediment<0
%         if (c2-1)<=intersections.A.MSL(1,1)
        [h]=thetarule2(theta,equilibrium.gulf.y+dh_add_sediment*t,1,T,1,recovery.gulf(t,:));
        recovery.gulf(t+1,:)=h(2,:)+dh_add_sediment;
%         else
%         [h]=thetarule2(theta,equilibrium.gulf.y_part1+dh_add_sediment*t,1,T,1,recovery.gulf_part1(t,:));
%         recovery.gulf_part1(t+1,:)=h(2,:)+dh_add_sediment;
%         
%         [h]=thetarule2(theta,equilibrium.gulf.y_part2,1,T,1,recovery.gulf_part2(t,:));
%         recovery.gulf_part2(t+1,:)=h(2,:);
%         
%         recovery.gulf(t+1,:)=recovery.gulf_part1(t+1,:)+recovery.gulf_part2(t+1,:);
%         end
    else
        [h]=thetarule2(theta,equilibrium.gulf.y,1,T,1,recovery.gulf(t,:));
        recovery.gulf(t+1,:)=h(2,:)+dh_add_sediment;
    end
      
    % due to this recovery process of the foreshore (gulfside) sediment is
    % coming available for the 'island' part of the profile. The amount of
    % available sediment is calculated each timestep.
    
    [IntersectA IntersectB intersections]=CalcintersectionsArea(recovery.gulf(t,:),recovery.gulf(t+1,:)-dh_add_sediment,MSLevel,wavebase,equilibrium.gulf.x);
    [Volume]=calcVolume(IntersectA, IntersectB,equilibrium.gulf.x(1),equilibrium.gulf.x(end));
    
    Sediment.gulf(t)=Volume.total;%+(add_sediment/12);
    
    % The recovery of the island depends on the sediment gathered from the
    % foreshore (Sediment.gulf(t)), the distance out of equilibrium and a
    % factor Beta (i.e vegetation a.o.)
    
    %depending on the available sediment received from the foreshore
    Sediment.island(t)=alpha*Sediment.gulf(t);
    Sediment.available(t)=Sediment.gulf(t);
   
    % The amount of sediment is spread as a layer of sediment over the
    % whole width of the island
    deltahisland=Sediment.island(t)/((equilibrium.island.x(end)-equilibrium.island.x(1)));
    if deltahisland<0
        deltahisland=0;
    end
    % If all of this gathered sediment can settle the next profile will
    % develop
    recovery.island(t+1,:)=recovery.island(t,:)+deltahisland;
    
    % The amount of sediment is limited by the equilibriumprofile (the
    % island isn't growing higher than the equilibriumprofile)
    if inclEquilibriumIsland==1;
        [h]=thetarule2(theta,equilibrium.island.y,1,T,1,recovery.island(t,:));
        recoverydt.island(t+1,:)=h(2,:);
        recovery.island(t+1,:)=min(recovery.island(t+1,:),recoverydt.island(t+1,:));
    end
    recovery.island(t+1,:)=min(equilibrium.island.y,recovery.island(t+1,:));
    
    % The difference between the volume at t+1 and t is the amount of
    % sediment that settle at the 'island'. By looking at the difference
    % with the amount that was available the rest volume can be calculated.
    restvol=Sediment.island(t)-abs(abs(trapz(equilibrium.island.x,(recovery.island(t+1,:)-recovery.island(t,:)))));%+sum(Sediment.island));
%     if restvol<0
%         restvol=0;
%     end
    % restz is the sum of difference between the recovery of the island and
    % the equilibrium of the island. The reason for this is ....
    restz=sum(abs(recovery.island(t+1,:)-equilibrium.island.y));
    % while there is still some sediment available the profile recovers.
    while restvol>0.01 && restz>1
        deltahisland=deltahisland+restvol/((equilibrium.island.x(end)-equilibrium.island.x(1)));
        recovery.island(t+1,:)=recovery.island(t,:)+deltahisland;
        recovery.island(t+1,:)=min(equilibrium.island.y,recovery.island(t+1,:));
        restvol=Sediment.island(t)-abs(abs(trapz(equilibrium.island.x,(recovery.island(t+1,:)-recovery.island(t,:)))));%+Sediment.island(t));
        restz=sum(abs(recovery.island(t+1,:)-equilibrium.island.y));
    end
    
    % Beta is included to limit the speed of recovery, this can be due to a
    % lack of vegetation is the early months after a hurricane for
    % instance.
    recovery.vegetation(t+1,:)=min(recovery.island(t,:)+Beta,recovery.island(t+1,:));
    
    Sediment.left.Sum(1)=restvol;
    Sediment.left.Sum(t+1)=abs(trapz(equilibrium.island.x,(recovery.island(t+1,:)-recovery.vegetation(t+1,:)))+Sediment.left.Sum(t)+abs((1-alpha)*Sediment.available(t)))+restvol;
    recovery.island(t+1,:)=recovery.vegetation(t+1,:);
    
    % the recovery of the bay depends on the difference of the amount of
    % sediment coming from the foreshore and settle on the island. If the
    % amount of sediment that settle on the island is smaller than the
    % amount of sediment from the foreshore this difference will be
    % available for the bay to develop. This will be done in such a way
    % that the bay side is expanding towards the mainland at a hight equal
    % to the MSL
    
    % the difference in height between the bay stormprofile and the MSL.
    % (This hight will be filled up)
    difh=MSL.sediment-storm.sediment;
    
    % calculate the area at each x-point from the bay side.
    for j=1:length(difh)-1
        Opp(j)=difh(1,j)*(equilibrium.sediment.x(j+1)-equilibrium.sediment.x(j));
    end
    
    % add the areas of each seperata delta-x/delta-h part.
    for j=1:length(Opp)
        AreaBay(j)=sum(Opp(1:j));
    end
    % If the sediment that is left from the foreshore is smaller than the
    % firs delta-x/delta-h area. Nothing will happen at the bay-side at the
    % profile at the bay remains the storm profile.
    if Sediment.left.Sum(t)<=AreaBay(1)
        recovery.sediment(t,:)=storm.sediment(1,:);
    end
    % if the amount sediment left from the foreshore can fit in the area of
    % the bay the bay side will be filled up to the MSL level. The last
    % part of the sediment left, which is insufficient to fill until MSL is
    % filled up in such a way that the 'bucket' is empty.
    if Sediment.left.Sum(t)>AreaBay(1) & Sediment.left.Sum(t)<=AreaBay(end)
        [xxx yyy]=find((AreaBay<(Sediment.left.Sum(t))));
        recovery.sediment(t,(1:yyy(end)))=0;
        leftover=(Sediment.left.Sum(t)-AreaBay(yyy(end)));
        hleftover=leftover/((equilibrium.sediment.x(yyy(end)+1)-equilibrium.sediment.x(yyy(end))));
        recovery.sediment(t,(yyy(end)+1:yyy(end)+2))=storm.sediment(1,(yyy(end)+1:yyy(end)+2))+hleftover;
        recovery.sediment(t,yyy(end)+3:end)=storm.sediment(1,yyy(end)+3:end);
        %         recovery.sediment(t,:)=storm.sediment(1,:);
    end
    % The amount of sediment is higher than the amount that can settle in
    % the bay. Think about this!!
    if Sediment.left.Sum(t)>AreaBay(end)
        [xxx yyy]=find((AreaBay<(Sediment.left.Sum(t))));
        recovery.sediment(t,:)=0;
        leftover=(Sediment.left.Sum(t)-AreaBay(end));
        hleftover=leftover/((equilibrium.sediment.x(end))-equilibrium.sediment.x(1));
        recovery.sediment(t,:)=0+hleftover;
        %         recovery.sediment(t,yyy(end)+3:end)=storm.sediment(1,yyy(end)+3:end);
        %         recovery.sediment(t,:)=storm.sediment(1,:);
    end
    
    % At this section the different parts of the profile is composed to one
    % profile.
    recoveryprofile.x=[equilibrium.deep.x,equilibrium.gulf.x,equilibrium.island.x,equilibrium.sediment.x];
    recoveryprofile.y(t,:)=[equilibrium.deep.y,recovery.gulf(t,:),recovery.island(t,:),recovery.sediment(t,:)];
%     if add_sediment<0
%             recoveryprofile.y(t,:)=[equilibrium.deep.y,recovery.gulfnegativesediment(t,:),recovery.island(t,:),recovery.sediment(t,:)];
%     end
    recoveryprofile.t(t,:)=interp1(recoveryprofile.x,recoveryprofile.y(t,:),xas);
  
    Supratidel=Tide*ones(1,length(recoveryprofile.t(t,:)));%*recoveryTime;
    if max(recoveryprofile.t(t,:))>Supratidel
        [intersections.recoveryprofile.tide]=calcIntersections(Supratidel,recoveryprofile.t(t,:),xas);
        widthIsland(t)=intersections.recoveryprofile.tide(1,end)-intersections.recoveryprofile.tide(1,1);
        [IntersectA IntersectB intersections]=CalcintersectionsArea(recoveryprofile.t(t,:),Supratidel,MSLevel,wavebase,xas);
        [intersections.recoveryprofile.tide]=calcIntersections(Supratidel,recoveryprofile.t(t,:),xas);
        [IslandArea]=calcVolume(IntersectA, IntersectB,intersections.recoveryprofile.tide(1,1),intersections.recoveryprofile.tide(1,end));
        distance(t)=intersections.recoveryprofile.tide(1,1);
        VolIsland(t)=IslandArea.erosion;
        HeightIsland(t)=VolIsland(t)/widthIsland(t);
        MaxHeight(t)=max(recoveryprofile.t(t,:));
    end
    if max(recoveryprofile.t(t,:))<=Supratidel
        [intersections.recoveryprofile.tide]=calcIntersections(Supratidel,recoveryprofile.t(t,:),xas);
        widthIsland(t)=0;
        distance(t)=0;
        VolIsland(t)=0;
        HeightIsland(t)=0;
        MaxHeight(t)=max(recoveryprofile.t(t,:));
    end
    
    [IntersectA IntersectB intersections]=CalcintersectionsArea(equilibriumprofile,stormprofile,MSLevel,wavebase,xas);

    [r1 c1]=find(xas>=intersections.B.wavebase(1,1),1);
    [hr1 hc1]=max(equilibriumprofile);
    [hhr hcc]=find(intersections.A.B(1,:)<=xas(hc1));
    [hr2 hc2]=max(intersections.A.B(2,1:hcc(end)));%>=intersections.A.MSL(1,1),1);
    [r2 c2]=find(xas>=intersections.A.B(1,hc2),1);
    j=1;
    while c2<=c1
        [r2 c2]=find(xas>=intersections.A.B(1,hc2-1+j),1);
        j=j+1;
    end

    [r3 c3]=find(xas>=intersections.B.MSL(1,2),1);
    [xy xz]=find((intersections.A.B(1,:)>=intersections.A.MSL(1,2)),1);
    
    
    equilibrium.deep.y=stormprofile(1:(c1-1));
    equilibrium.deep.x=xas(1:(c1-1));
    equilibrium.gulf.y=equilibriumprofile(c1:(c2-1));
    equilibrium.gulf.x=xas(c1:(c2-1));
    equilibrium.island.y=equilibriumprofile(c2:(c3-1));
    equilibrium.island.x=xas(c2:(c3-1));
    equilibrium.bay.y=equilibriumprofile(c3:length(xas));
    equilibrium.bay.x=xas(c3:length(xas));
    equilibrium.sediment.x=xas(x.A.B.bay:length(xas));
    MSL.sediment=MSL.total(x.A.B.bay:length(xas));
    
    %% Equilibrium profile
    
    %post-storm profile
    storm.gulf=stormprofile(c1:(c2-1));
    storm.island=stormprofile(c2:(c3-1));
    storm.bay=stormprofile(c3:length(xas));
    storm.sediment=stormprofile(x.A.B.bay:length(xas));

    t=t+1;
end
t=t-1;

%% influence sealevelrise
if sealevelchange>0
    seaprof_single=recoveryprofile.t(t,:);
    [sealevelprofile sealevelprofiletrain seatransl]=SeaLevelRise(seaprof_single,MSLevel,wavebase,sealevelchange,xas,recoveryTime); %[Bsealevel]=SeaLevelRise(InitialProfile,MSLevel,wavebasedepth,sealevelchange)
    recoveryprofile.t(t,:)=sealevelprofiletrain;
end
if sealevelchange==0
    seatransl=0;
end

%check for continuity
[IntersectA IntersectB intersections]=CalcintersectionsArea(recoveryprofile.t(t,:),recoveryprofile.t(1,:),MSLevel,wavebase,xas);
[Continuity]=calcVolume(IntersectA, IntersectB,xas(1),xas(end));%3738.45,3874.34)

% The profile is lowered with the amount of the sealevelchange, this is
% done in order to keep the waterlevels in XBeach contstant.
if sealevelchange>0
    recoveryprofile.t(t,:)=sealevelprofiletrain;
end
recoveryprofilet=recoveryprofile.t(t,:);
recoveryprofileall=recoveryprofile.t;
if max(recoveryprofile.t(t,:))>Supratidel
    [intersections.recoveryprofile.tide]=calcIntersections(Supratidel,recoveryprofile.t(t,:),xas);
    widthIsland(t)=intersections.recoveryprofile.tide(1,end)-intersections.recoveryprofile.tide(1,1);
    [IntersectA IntersectB intersections]=CalcintersectionsArea(recoveryprofile.t(t,:),Supratidel,MSLevel,wavebase,xas);
    [intersections.recoveryprofile.tide]=calcIntersections(Supratidel,recoveryprofile.t(t,:),xas);
    [IslandArea]=calcVolume(IntersectA, IntersectB,intersections.recoveryprofile.tide(1,1),intersections.recoveryprofile.tide(1,end));
    distance(t)=intersections.recoveryprofile.tide(1,1);
    VolIsland(t)=IslandArea.erosion;
    HeightIsland(t)=VolIsland(t)/widthIsland(t);
    MaxHeight(t)=max(recoveryprofile.t(t,:));
end
if max(recoveryprofile.t(t,:))<=Supratidel
    [intersections.recoveryprofile.tide]=calcIntersections(Supratidel,recoveryprofile.t(t,:),xas);
    widthIsland(t)=0;
    distance(t)=0;
    VolIsland(t)=0;
    HeightIsland(t)=0;
    MaxHeight(t)=max(recoveryprofile.t(t,:));
end
