function [equilibriumprofile stormprofile]=CreateEquilibriumProfile4(i,OverwashV,initialprofile,A,B,wavebase,MSLevel,xas)
% 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.

try

[row1 col1]=max(B(i,:));
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=interp1(newxas,initialprofile,xas);
[IntersectA IntersectB intersections]=CalcintersectionsArea(equilibriumprofile,B(i,:),MSLevel,wavebase,xas);
[hhr hcc]=find(intersections.A.B(1,:)<=xas(col1));%+100);



[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(i)
    % translated to the left
    xtrans=j*-0.1+0.1;
    x1=xas(1:c1-1);
    x2=xas(c1:end)+xtrans;
    newxas=[x1,x2];
    equilibriumprofile=interp1(newxas,equilibriumprofile,xas);
    zend=equilibriumprofile((abs(equilibriumprofile)>=0)==1);
    equilibriumprofile((abs(equilibriumprofile)>=0)==0)=zend(end);
    [IntersectA IntersectB intersections]=CalcintersectionsArea(equilibriumprofile,B(i,:),MSLevel,wavebase,xas);
    [hr1 hc1]=max(equilibriumprofile);
    [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(i)
    % translate to the right
    xtrans=j*0.1;
    x1=xas(1:c1-1);
    x2=xas(c1:end)+xtrans;
    newxas=[x1,x2];
    equilibriumprofile=interp1(newxas,equilibriumprofile,xas);
    zend=equilibriumprofile((abs(equilibriumprofile)>=0)==1);
    equilibriumprofile((abs(equilibriumprofile)>=0)==0)=zend(end);
    [IntersectA IntersectB intersections]=CalcintersectionsArea(equilibriumprofile,B(i,:),MSLevel,wavebase,xas);
    [hr1 hc1]=max(equilibriumprofile);
    [hhr hcc]=find(intersections.A.B(1,:)<=xas(hc1));%+100);
    [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

[IntersectA IntersectB intersections]=CalcintersectionsArea(equilibriumprofile,B(i,:),MSLevel,wavebase,xas);
[hr1 hc1]=max(equilibriumprofile);
[hhr hcc]=find(intersections.A.B(1,:)<=xas(hc1));
[hr2 hc2]=max(intersections.A.B(2,1:hcc(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(hccc(1):hc1),intersections.B.MSL(2,end)];
island.eq.y=max(interp1(xas_island,yas_island,xas(hccc(1):hcccc(end))),B(i,(hccc(1):hcccc(end))));    

equilibriumprofile(hccc(1):hcccc(end))=island.eq.y;

catch
    
% [IntersectA IntersectB intersections]=CalcintersectionsArea(initialprofile,B(i,:),MSLevel,wavebase,xas);
%     
% xtranslated=(intersections.B.wavebase(1,1)-intersections.A.wavebase(1,1));
% [r1 c1]=find(xas>=intersections.A.wavebase(1,1),1);
% x1=xas(1:c1-1);
% x2=xas(c1:end)+xtranslated;
% newxas=[x1,x2];

% no overwash therefore no recovery so 
equilibriumprofile=B(i,:);
end
stormprofile=B(i,:);