% function [Analises]=modeltrain(nr_hurricane,MSLevel,wavebase,Tide,sealevelchangerate,alpha,Beta,theta,T,add_sediment)
% This function contains the modeltrain
% nr_hurricanes is the number of the hurricanes
% MSLevel=MSLevel
% wavebase = kind of a closure depth
% Tide = reference level for the width of the island
% sealevelchangerate = rate of change per century
% alpha = vegetation zone factor
% Beta = max. rate of growth due to for instance vegetation
% theta = thetarule 0 - 0.5 - 1
% T = recovery time scale
% add_sediment = add sediment to the foreshore >0

% clear all,clc,
tic
addpath D:\pruis\XBeach\MatlabScripts\benodigd
addpath D:\pruis\XBeach\MatlabScripts\benodigd\Visualisation
addpath ([cd, '\jonswap_files'])

prompt = {'alpha value ', 'Beta value ',...
'sealevelchangerate ','T ','added Sediment'};
dlg_title = 'Recovery parameters';
num_lines = 1;
def = {'1','0.05','0.64','6','0'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
alpha=str2double(answer{1});
Beta=str2double(answer{2});
sealevelchangerate=str2double(answer{3});
T=str2double(answer{4});
add_sediment=str2double(answer{5});

XBRUN=0;
nr_hurricane=1;
MSLevel=0;
wavebase=-7;
Tide=0.37;
theta=1;
D_low=1;

x= load('x.grd');
xas=x(1,:);
scale = xas;
ngridp=size(x,2);

varnames = {'H','zb','zs'};
sealevelchange=sealevelchangerate/1200; %per month

%initial profile is needed as an equilibrium shape
copyfile('z.dep','InitialProfile.txt');
InitialProfile=load('InitialProfile.txt')';
initialprofile=InitialProfile(:,2)';

% take into account =1 don't take into account=0

%Om XB te beginnen zijn er een aantal zaken nodig
% - params = initeel = params_run_01.txt                                    % dit is een afwijking van de normale procedure!!!
% - depfile = initeel = z_run_00.dep

copyfile('z.dep','z_run_00.dep');
crsection='z_run_00';
fid = fopen('z_run_00.dep');
z_initial=textscan(fid,'%f32',452,'delimiter','');
fclose(fid);



[katrina steepness Hm0]=getKatrina();

% The modeltrain starts here
i=1;
while i <= nr_hurricane & D_low>Tide
    param=1;
    
    %hier wordt voor iedere run de bijbehorende param file gekopieerd naar
    %de input param file
    copyfile(['params_run_',num2str(i,'%02i'),'.txt'],'params.txt');
    
    [Hswash Hinundation scaling hurricanetrain profile]=calcRegime(crsection,i,param,Hm0,ngridp);
    
    [regime]=scaleKatrina(katrina,i,scaling,steepness,Hm0);
    if XBRUN==1
        ! xbeach.exe
    end
    
    XB=getdimensions;
    
    profile.new.what='zb';%['zb_run_',num2str(i,'%02i')];
    profile.new.readvar.all=readvar([profile.new.what,'.dat'],XB);
    profile.new.readvar.end = squeeze(profile.new.readvar.all(:,:,end))';
    profile.end=profile.new.readvar.end(1,:);
    
    %put this in a .txt file for control
    
    %% Profiles
    
    %een profiel voor iedere run wegschrijven voor latere controle
    % A=pre-storm
    A(i,:)=profile.initial'; %crsection.dep
    
    % post-storm
    B(i,:)=profile.end;
    
    % MSL
    MSL=(ones(1,(size(profile.end,2)))*MSLevel)';
    
    
    
    %% dat bestanden voor de run wegschrijven
    
    for j = 1:length(varnames)
        varn{j} = [varnames{j},'_run_',num2str(i,'%02i')];
        copyfile([varnames{j},'.dat'],[varn{j},'.dat']);
    end
    
    
    %% recovery
    
    %calculate intersections
    [IntersectA IntersectB intersections]=CalcintersectionsArea(A(i,:),B(i,:),MSLevel,wavebase,xas);
    
    % calculate overwash volume
    % highest intersection point between pre- and post-storm profiles
    [x y]=max(intersections.A.B(2,:));
    x1=intersections.A.B(1,y);
    x2=xas(end);
    [Volume]=calcVolume(IntersectA, IntersectB,x1,x2);
    
    if Volume.accretion<0
        OverwashV(i)=0;
    end
    if Volume.accretion>=0
        OverwashV(i)=Volume.accretion;
    end
    
    
    x1=0;
    x2=xas(end);
    [Volume]=calcVolume(IntersectA, IntersectB,x1,x2);
    
    XBeach.erosion=Volume.erosion;
    XBeach.accretion=Volume.accretion;
    XBeach.continuity=Volume.total;
    
    Volumetranslation(i)=sum(OverwashV);
    
%     if add_sediment<0
%         OverwashV=OverwashV+abs(add_sediment);
%     end
    %Create new equilibriumprofile
    [IntersectA IntersectB intersections]=CalcintersectionsArea(A(i,:),B(i,:),MSLevel,wavebase,xas);
%     mIP(1)=0;
        [equilibriumprofile stormprofile]=CreateEquilibriumProfile4(i,OverwashV,initialprofile,A,B,wavebase,MSLevel,xas);

%     [equilibriumprofile stormprofile translated]=CreateEquilibriumProfile3(i,Volumetranslation,mIP,initialprofile,A,B,wavebase,MSLevel,xas);
    %     [equilibriumprofile stormprofile]=CreateEquilibriumProfile2(i,Volumetranslation,initialprofile,A(i,:),B(i,:),intersections,wavebase,MSLevel,xas);
    
    
    
    allequilibriumprofiles(i,:)=equilibriumprofile;
    % recovery
%     recoveryTime=hurricanetrain.recoverytime(i,1);
    recoveryTime=100;
    [Continuity recoveryprofilet recoveryprofileall widthIsland VolIsland HeightIsland distance seatransl]=recoverprof1910(alpha,recoveryTime,stormprofile,equilibriumprofile,MSLevel,wavebase,Tide,theta,T,Beta,xas,sealevelchange,add_sediment);
    
%     mIP(i+1)=translated+seatransl;
    
    %write the new cross-section for the next run
    crsection=['z_run_',num2str(i,'%02i')];
    D_low=max(recoveryprofilet);
    dlmwrite([crsection,'.dep'],[recoveryprofilet;recoveryprofilet;recoveryprofilet],'delimiter',' ','precision','%.10f');
    
    
    %create the new param-file
    [params_run]=CreateParam(i,crsection,ngridp);
    
    
    
    %% evaluation
    
    zbrun=['zb_run_',num2str(i,'%02i')];
    zb=readvar([zbrun,'.dat'],XB);
    zsrun=['zs_run_',num2str(i,'%02i')];
    zs=readvar([zsrun,'.dat'],XB);
    Hrun=['H_run_',num2str(i,'%02i')];
    H=readvar([Hrun,'.dat'],XB);
    
    
    
    Hurricane.pre_storm=A(i,:);
    Hurricane.post_storm=B(i,:);
    Hurricane.XBeach.Volumes=XBeach;
    Hurricane.XBeach.hydro.zb=zb;
    Hurricane.XBeach.hydro.zs=zs;
    Hurricane.XBeach.hydro.H=H;
    Hurricane.recovery.profiles=recoveryprofileall;
    Hurricane.recovery.continuity=Continuity.total;
    Hurricane.recovery.islandwidth=widthIsland;
    Hurricane.equilibriumprofile=equilibriumprofile;
    Analises{1,i}=Hurricane;
    
    if XBRUN==1
    %     bestanden ordenen
        mkdir(['Folder_',num2str(i,'%02i')])
        to_be_moved=regime.input.files.list;
        varn{j}
        for j=1:length(to_be_moved)
            movefile(to_be_moved{j},[['Folder_',num2str(i,'%02i'),'/'],to_be_moved{j}])
            movefile(['q_',to_be_moved{j}],[['Folder_',num2str(i,'%02i'),'/'],['q_',to_be_moved{j}]])
            movefile(['E_',to_be_moved{j}],[['Folder_',num2str(i,'%02i'),'/'],['E_',to_be_moved{j}]])
    
        end
    end
    
    
    i=i+1;
end

save('Analyse','Analises')
toc

movie=0;
XB=getdimensions;
% load('allequilibriumprofiles.txt')
% load('checkcontinuity')
load('Analyse')
nr_hurricanes=1;
j=1;
for j=1:nr_hurricanes
    zb_simpl=squeeze(Analises{1,j}.XBeach.hydro.zb(:,2,:))';
    zs_simpl=squeeze(Analises{1,j}.XBeach.hydro.zs(:,2,:))';
    H_simpl=squeeze(Analises{1,j}.XBeach.hydro.H(:,2,:))';
    
    recoveryphase=Analises{1,j}.recovery.profiles;
    allequilibriumprofiles(j,:)=Analises{1,j}.equilibriumprofile;
    
    
    blanco=zeros(size(Analises{1,j}.recovery.profiles,1),size(Analises{1,j}.recovery.profiles,2));
    
    profile=[zb_simpl; recoveryphase];
    wave=[H_simpl; blanco];
    surface=[zs_simpl; blanco];
    
    Matrix{1,j}=profile;
    WaveM{1,j}=wave;
    SurfaceM{1,j}=surface;
    
    j+1;
end


%j is nummer of hurricane
%r is row number
if movie==0
    figure()
end
if movie==1
    fig=figure;
    MAP=colormap;
    aviobj=avifile('9timesOverwash.avi','compression','Cinepak','fps',5,'COLORMAP',MAP,'quality',100);
end
for j=1:nr_hurricanes
    for r=87:size(Matrix{1,j},1)
        plot(xas,(SurfaceM{1,j}(r,:)+WaveM{1,j}(r,:)),'b','LineWidth',2);hold on
%         harea=area(xas,[(SurfaceM{1,j}(r,:)'+WaveM{1,j}(r,:)'),Matrix{1,j}(r,:)'])
        hulpMatrix=(0-Matrix{1,j}(r,:));
%         hulpMatrix(Matrix{1,j}(r,:)>=0)=0;
        nieuweMatrix=[Matrix{1,j}(r,:)',hulpMatrix'+(SurfaceM{1,j}(r,:)'+WaveM{1,j}(r,:)')];
        harea=area(xas,nieuweMatrix);
        set(harea(1),'FaceColor',[255/255 255/255 255/255]);
        set(harea(2),'FaceColor',[215/255 240/255 255/255]);
        fill([xas,xas(end)+0.1],[Matrix{1,j}(r,:),Matrix{1,j}(r,1)],[255/255 255/255 210/255]);
        plot(xas,Matrix{1,j}(r,:),'k','LineWidth',2); 
%         plot(xas,SurfaceM{1,j}(r,:),'r');
        if r<size(Matrix{1,j},1)
            plot(xas,Matrix{1,1}(1,:),'--r','LineWidth',2);
        else 
            sealevellowered=Matrix{1,1}(1,:)-sealevelchange*r
            plot(xas,sealevellowered,'--r','LineWidth',2);
        end
        if r>size(zb_simpl,1)
            if r<size(Matrix{1,j},1)
                plot(xas,allequilibriumprofiles(j,:),':k','LineWidth',2);
            else
                sealevelloweredEQ=Matrix{1,1}(1,:)-sealevelchange*r
                plot(xas,sealevelloweredEQ,':k','LineWidth',2);
            end
        end
%         if r<=size(H_simpl,1)
%             title(['hurricane ',num2str(j),'-',num2str(r)]);
%         end
%         if r>size(H_simpl,1)
%             title(['recovery ',num2str(j),'-',num2str(r-size(H_simpl,1))])
%         end
        if r<=size(H_simpl,1)
            text(2850,7,['hurricane ',num2str(j,'%02i'),'-',num2str(r)],'FontSize',14,'FontWeight','bold');
        end
        if r>size(H_simpl,1)
            text(2850,7,[' recovery ',num2str(j,'%02i'),'-',num2str(r-size(H_simpl,1))],'FontSize',14,'FontWeight','bold');
        end
        
%         text(3000,-7.2,['continuity XBeach phase      ',num2str(Analises{1,j}.XBeach.Volumes.continuity)],'FontSize',11);
%         text(3000,-8,['continuity recovery phase   ',num2str(Analises{1,j}.recovery.continuity)],'FontSize',11);
%         if r>size(H_simpl,1)
%         text(3000,-8.8,['width of the island recovery phase   ',num2str(Analises{1,j}.recovery.islandwidth(1,r-size(H_simpl,1)))],'FontSize',11);
%         end
        axis([0 7400 -20 10]);
        pause(0.1), hold off

        teller=([num2str(j,'%02i'),'-',num2str(r,'%02i'),'.png']);
%         saveas(fig,teller)
        saveas(gcf,teller)
        if movie==1
            F=getframe(fig);
            aviobj=addframe(aviobj,F);
            
        end
    end
end
if movie==1;
    close(fig);
    aviobj=close(aviobj);
end