% function [sealevelprofile sealevelprofiletrain seatransl]=SeaLevelRise1510(seaprof_single,MSLevel,wavebase,sealevelchange,xas,recoveryTime)
% function [sealevelprofile sealevelprofiletrain seatransl stormprofile equilibriumprofile]=SeaLevelRise(seaprof_single,MSLevel,wavebase,sealevelchange,xas,recoveryTime,equilibriumprofile,stormprofile)
% This function takes into account the sealevel rise and creates a new
% profile in such a way that the amount of sediment 'eroded' on the
% foreshore is equal to the amount of accretion of the back side of the
% island. The input parameters needed for this function are listed below;
%   seaprof_single      = pre- sealevelrise profile
%   MSLevel             = Mean Sea Level
%   wavebase            = wavebase/closure depth
%   sealevelchange      = monthly sealevel change rate
%   xas                 = x-points (needed to calculate the area)
%   recoveryTime        = The recovery time after each hurricane
% recoveryTime=1;
% The initial pre-sealevel profile
sealevelchange=6/1200;
load('Analyse.mat')
seaprof_single=Analises{1,1}.recovery.profiles(end-1,:);
endvalue=seaprof_single;
MSL.total=MSLevel*ones(1,length(xas));
% The initial pre-sealevel profile is raised with the same amount as the
% sealevelchange = sealevelrate(century) / 1200 (months) *
% recoveryTime(months)
sealevelprof=seaprof_single+sealevelchange*recoveryTime;

% Find the wavebase intersection, below this intersection the profile
% doesn't change.
[IntersectA IntersectB intersections]=CalcintersectionsArea(seaprof_single,sealevelprof,MSLevel,wavebase,xas);
[rd cd]=find(xas>intersections.A.wavebase(1,1),1);
sealevelprof(1:(cd-1))=sealevelprof(1:(cd-1))-sealevelchange*recoveryTime;

x1=intersections.A.wavebase(1,1);
x2=xas(end);

SeaVol=-100;
j=1;
while SeaVol<0
    
    xtrans=j*1.1;
    xas1=xas(1:cd-1);
    xas2=xas(cd:end)+xtrans;
    xascombi=[xas1(1:cd-1),xas2];
    
    sealevelprof(j+1,:)=interp1(xascombi,sealevelprof(j,:),xascombi);
    sealevelprofx(j+1,:)=interp1(xascombi,sealevelprof(j,:),xas);
    [rowf columnf]=find(abs(sealevelprofx(j+1,:)-sealevelprofx(j+1,end))>0.01);
    sealevelprofx(j+1,columnf(end)+1:end)=endvalue(columnf(end)+1:end);
    
    [IntersectA IntersectB intersections]=CalcintersectionsArea(seaprof_single,sealevelprofx(j+1,:),MSLevel,wavebase,xas);
    
    x1=intersections.A.wavebase(1,1)-xtrans;
%     x2=intersections.B.MSL(1,2);
%     x2=intersections.A.MSL(1,2)+xtrans;
    x2=xas(end);
    [Volume]=calcVolume(IntersectA, IntersectB,x1,x2);
    
    SeaVol=Volume.total; %zou 0 moeten zijn.
    
    [rowsea columnsea]=find(xas>=x2,1);
    sealevelprofx(j+1,columnsea:end)=seaprof_single(1,columnsea:end);

    j=j+1;
end

sealevelprofile=sealevelprofx(j,:);
sealevelprofiletrain=sealevelprofx(j,:)-sealevelchange*recoveryTime;
sealevelprofiletrain(columnf(end)+1:end)=endvalue(columnf(end)+1:end);
seatransl=xtrans;

k=j;
figure()
for j=1:3:k
    plot(xas,sealevelprofx(j,:),'k','LineWidth',2),hold on
    text(500,3,' Visualisation of Sealevel Rise 6 m / Century ','FontSize',14,'FontWeight','bold');
    hulpMatrix=(0-(sealevelprofx(j,:)));
    nieuweMatrix=[(sealevelprofx(j,:))',hulpMatrix'];
    harea=area(xas,nieuweMatrix);
    set(harea(1),'FaceColor',[215/255 240/255 255/255]);
    set(harea(2),'FaceColor',[215/255 240/255 255/255]);
    fill([xas,xas(end)+0.1],[(sealevelprofx(j,:)),(sealevelprofx(j,1)-sealevelchange*recoveryTime)],[255/255 255/255 210/255],'LineWidth',2);
    plot(xas,seaprof_single,':r','LineWidth',2)
    plot(xas,MSL.total,'b','LineWidth',2)
    axis([0000 7000 -7.5 5]);
    teller=([num2str(j,'%02i'),'.png']);
    saveas(gcf,teller)
    pause(0.1),hold off;
end