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


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

% 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

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

varnames = {'H','zb','zs'};


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

[IntersectA IntersectB intersections]=CalcintersectionsArea(initialprofile,initialprofile,Tide,wavebase,xas);
XDistance=intersections.A.MSL(1,1);
% 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);

alpha_array=1;
Beta_array=[0,0.001, 0.005,0.01, 0.05, 0.1,0.5,1];
T_array=[0.01, 1, 6, 100];
sealevelchangerate_array=[0.29, 0.64, 0.99];
add_sediment_array=[-10,-6,-3,-0.3,0,0.3, 3, 6, 10];

CollectedData=zeros(nr_hurricane,length(alpha_array),length(Beta_array),length(T_array),length(sealevelchangerate_array),length(add_sediment_array),4,100);


[katrina steepness Hm0]=getKatrina();

% The modeltrain starts here
i=1;
while i <= nr_hurricane
    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);
    

    
    
    datestr(datevec(now),0)
    for mm =1:length(alpha_array)
        for nn =1:length(Beta_array);
            for oo =1:length(T_array);
                for pp =1:length(sealevelchangerate_array);
                    for qq =1:length(add_sediment_array);
                        try
                            recoveryTime=hurricanetrain.recoverytime(i,1);
                            alpha=alpha_array(1,mm);
                            Beta=Beta_array(1,nn);
                            T=T_array(1,oo);
                            sealevelchangerate=sealevelchangerate_array(1,pp);
                            sealevelchange=sealevelchangerate/1200; %per month
                            add_sediment=add_sediment_array(1,qq);
%                             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);
                            allequilibriumprofiles(i,:)=equilibriumprofile;
                            % recovery
                            
                            [Continuity recoveryprofilet recoveryprofileall widthIsland VolIsland HeightIsland distance seatransl MaxHeight]=recoverprof1910(alpha,recoveryTime,stormprofile,equilibriumprofile,MSLevel,wavebase,Tide,theta,T,Beta,xas,sealevelchange,add_sediment);
                            
                            
                            CollectedData(i,mm,nn,oo,pp,qq,1,:)=HeightIsland;
                            CollectedData(i,mm,nn,oo,pp,qq,2,:)=widthIsland;
                            CollectedData(i,mm,nn,oo,pp,qq,3,:)=VolIsland;
                            CollectedData(i,mm,nn,oo,pp,qq,4,:)=distance-XDistance;
                            CollectedData(i,mm,nn,oo,pp,qq,5,:)=MaxHeight;
                        catch
                            CollectedData(i,mm,nn,oo,pp,qq,1,:)=NaN;
                            CollectedData(i,mm,nn,oo,pp,qq,2,:)=NaN;
                            CollectedData(i,mm,nn,oo,pp,qq,3,:)=NaN;
                            CollectedData(i,mm,nn,oo,pp,qq,4,:)=NaN;
                            CollectedData(i,mm,nn,oo,pp,qq,5,:)=NaN;
                        end
                    end
                    
                end
            end
        end
    end
    datestr(datevec(now),0)
    % %     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;
    %
    
    %bestanden ordenen
    if XBRUN==1
    %     bestanden ordenen
    try
        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
    catch
        Moving_failed=1;
    end
    end
    
    
    i=i+1;
end

save('CollectedData','CollectedData')

