clear all, clc
movie=0;

i=1;
j=1;

% OverwashV=355;
OverwashV=9.8895e-005
% initialprofile=load('InitialProfile.txt');
initialprofile=initialprofile(2,:);
% 
% load('Analyse')
% 
% 
% A=Analises{1,1}.pre_storm;%load('z_run_01.dep');
% % A=A(2,:);
% B=Analises{1,1}.post_storm;%load('z_run_02.dep');
% B=B(2,:);
wavebase=-7;
MSLevel=0;
xas=load('x.grd');
xas=xas(2,:);

% This function creates a new equilibriumprofile

% The initial profile (profile at t=0) is placed in such a way that the
% x position of the maximum elevation of the initial profile coinsist with
% the x position of the maximum elevation of the post-storm profile.

[row1 col1]=max(B(1,:));
maxBx=xas(col1);
[row2 col2]=max(initialprofile);
maxIx=xas(col2);
xtranslated=maxBx-maxIx;

[IntersectA IntersectB intersections]=CalcintersectionsArea(initialprofile,B(1,:),MSLevel,wavebase,xas);
[r1 c1]=find(xas>=intersections.A.wavebase(1,1),1);
x1=xas(1:c1-1);
x2=xas(c1:end)+xtranslated;
newxas=[x1,x2];
equilibriumprofile(j,:)=interp1(newxas,initialprofile,xas);
[IntersectA IntersectB intersections]=CalcintersectionsArea(equilibriumprofile(j,:),B(1,:),MSLevel,wavebase,xas);
[hhr hcc]=find(intersections.A.B(1,:)<=xas(col1));
[hr2 hc2]=max(intersections.A.B(2,1:hcc(end)));

% if the Volume is <0 the initial profile is moved to much to the right
% if the Volume is >0 the initial profile needs to be moved more to the
% right.
[Volume]=calcVolume(IntersectA, IntersectB,xas(c1-1),intersections.A.B(1,hc2));

j=1;
while Volume.total<0 & abs(Volume.total)>OverwashV(1)
    % translated to the left
    xtrans=j*-0.01+0.01;
    x1=xas(1:c1-1);
    x2=xas(c1:end)+xtrans;
    newxas=[x1,x2];
    equilibriumprofile(j+1,:)=interp1(newxas,equilibriumprofile(j,:),xas);
    zend=equilibriumprofile(j,((abs(equilibriumprofile(j+1,:))>=0)==1));
    equilibriumprofile(j+1,((abs(equilibriumprofile(j+1,:))>=0)==0))=zend(end);
    [IntersectA IntersectB intersections]=CalcintersectionsArea(equilibriumprofile(j+1,:),B(1,:),MSLevel,wavebase,xas);
    [hr1 hc1]=max(equilibriumprofile(j+1,:));
    [hhr hcc]=find(intersections.A.B(1,:)<=xas(hc1));
    [hr2 hc2]=max(intersections.A.B(2,1:hcc(end)));
    [Volume]=calcVolume(IntersectA, IntersectB,xas(c1-1),intersections.A.B(1,hc2));
    j=j+1;
end

while Volume.total>0 & abs(Volume.total)<OverwashV(1)
    % translate to the right
    xtrans=j*0.01-0.01;
    x1=xas(1:c1-1);
    x2=xas(c1:end)+xtrans;
    newxas=[x1,x2];
    equilibriumprofile(j+1,:)=interp1(newxas,equilibriumprofile(j,:),xas);
    zend=equilibriumprofile(j,((abs(equilibriumprofile(j+1,:))>=0)==1));
    equilibriumprofile(j+1,((abs(equilibriumprofile(j+1,:))>=0)==0))=zend(end);
    [IntersectA IntersectB intersections]=CalcintersectionsArea(equilibriumprofile(j+1,:),B(1,:),MSLevel,wavebase,xas);
    [hr1 hc1]=max(equilibriumprofile(j+1,:));
    [hhr hcc]=find(intersections.A.B(1,:)<=xas(hc1));
    [hr2 hc2]=max(intersections.A.B(2,1:hcc(end)));
    [Volume]=calcVolume(IntersectA, IntersectB,xas(c1-1),intersections.A.B(1,hc2));
    j=j+1;
end

XM=xas(hc1); %max;
[hhhr hccc]=find(xas>=intersections.A.B(1,hc2));
[hhhhr hcccc]=find(xas<=intersections.B.MSL(1,end));
xas_island=[xas(hccc(1):hc1),xas(hcccc(end))];
yas_island=[equilibriumprofile(j,(hccc(1):hc1)),intersections.B.MSL(2,end)];
island.eq.y=max(interp1(xas_island,yas_island,xas(hccc(1):hcccc(end))),B(1,(hccc(1):hcccc(end))));    

equilibriumprofile(j,(hccc(1):hcccc(end)))=island.eq.y;
equilibriumprofile(j,hcccc(end)+1:end)=B(1,hcccc(end)+1:end);
stormprofile=B(1,:);
if movie==0
figure()
end
if movie==1
    fig=figure;
    MAP=colormap;
    aviobj=avifile('EquilibriumProfile.avi','compression','none','fps',5,'COLORMAP',MAP,'quality',100);
end
for j=1:size(equilibriumprofile,1)
    plot(xas,A,'--r','LineWidth',2),hold on
    harea=area(xas,B);
    set(harea(1),'FaceColor',[215/255 240/255 255/255]);
    fill([xas,xas(end)+0.1],[B,B(1,1)],[255/255 255/255 210/255]);
    plot(xas,equilibriumprofile(j,:),'k','LineWidth',2)
    plot(xas,B,'--k','LineWidth',1)
    plot(xas,A,'--r','LineWidth',2)
    plot(xas,(ones(1,length(xas))*MSLevel),'b','LineWidth',2)
    axis([0 7499 -20 5])
    line(7499,[-20,5])
    
%      mkdir('afbeeldingen')
    cd([cd,'\afbeeldingen'])
    teller=([num2str(j,'%02i'),'.png']);
    saveas(gcf,teller)
    cd ..
        
    pause(0.01),hold off,
    if movie==1
        F=getframe(fig);
        aviobj=addframe(aviobj,F);
    end
end
if movie==1;
    close(fig);
    aviobj=close(aviobj);
end
