% 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,



modeltrain_numbers = [1:3];
threshold_factors = [3.0:0.5:4];

for j=1:length(modeltrain_numbers)
    
modeltrain_number = modeltrain_numbers(j);
threshold_factor = threshold_factors(j); 

modeltraindir = cd;
inputdir = fullfile(fileparts(which(mfilename)),'..','InputData','Grids');
outputdir = fullfile(fileparts(which(mfilename)),'..','InputData',['Modeltrain_number_',num2str(modeltrain_number,'%02i')],['Storm_',num2str(1,'%02i')]);
paramdir = fullfile(fileparts(which(mfilename)),'..','XBeach');
xboutputdir = fullfile(fileparts(which(mfilename)),'..','Output');
recoverydir  = fullfile(fileparts(which(mfilename)),'..','Barrier_Island_Recovery');

XBRUN=1;

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

sealevelchangerate=0.64;
alpha=1;
Beta=0.05;
T=100;
add_sediment=00;

x= load(fullfile(inputdir,'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
InitialProfile=load(fullfile(inputdir,'z.dep'));
mkdir(fullfile(outputdir))
dlmwrite(fullfile(outputdir,'InitialProfile.txt'),InitialProfile);
initialprofile=InitialProfile(2,:);
D_low=max(initialprofile);

addpath(recoverydir)
[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

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

% derive how many hurricanes should be modelled and the jonswap spectrum
% are being maked
[storms recovery]=selectStorms(modeltrain_number,threshold_factor);
recoveryTimes=[(recovery.time(:)'/30),0 ];% months of recovery, no recovery after last hurricane

% Create the first param-file
outputdir_var = fullfile(fileparts(which(mfilename)),'..','InputData',['Modeltrain_number_',num2str(modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]);
paramdir_var = fullfile(fileparts(which(mfilename)),'..','InputData',['Modeltrain_number_',num2str(modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]);
boundaryConditions = fullfile(['jonswaplist_storm_',num2str(i,'%02i'),'.txt']);
surgeConditions=fullfile(['Scaled_Surge_Storm_',num2str(i,'%02i'),'.dat']);

[params_run]=CreateParam(1,crsection,ngridp,boundaryConditions,surgeConditions,paramdir_var);


% The modeltrain starts here
i=8;
while i <= size(storms,2) & D_low>Tide
    param=1;
    copyfile(fullfile(paramdir_var,['params_run_',num2str(i,'%02i'),'.txt']),fullfile(paramdir,'params.txt'));
    
    %hier wordt voor iedere run de bijbehorende param file gekopieerd naar
    %de input param file
    Max_WVHT_storm = max(storms{1,i}.Katrinas.waveheight);
    Hswash=(4*(D_low-0.45))/((4*0.83*(0.0035/sqrt(0.016))+0.8)+0.5);
    Hinundation=(0.56*(D_low-0.45)+3.72)/(0.56*(-0.02*(0.0035/sqrt(0.016))+0.14)+1);
    
    if Max_WVHT_storm<=Hswash
        Hurricane.regime='Swash';
    else if Max_WVHT_storm<=Hinundation
             Hurricane.regime='Overwash';
        else
            Hurricane.regime='Inundation';
        end
    end
    
    if XBRUN==1
        
        inputspectrums  = (fullfile(fileparts(which(mfilename)),'..','InputData',['Modeltrain_number_',num2str(modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]));
        copyfile(inputspectrums,fullfile(inputspectrums,'..','..','..','XBeach'))
        
        cd(fullfile(inputspectrums,'..','..','..','XBeach'))
        
        ! xbeach.exe
        XB=xb_getdimensions;
        xbeachrundir = fullfile(fileparts(which(mfilename)),'..','Runs');
        xbeachrunstormdir = fullfile(xbeachrundir,['Modeltrain_number_',num2str(modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]);
        mkdir(xbeachrunstormdir);
    end
    
    
    
    profile.new.what='zb';%['zb_run_',num2str(i,'%02i')];
    profile.new.readvar.all=xb_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
    temp = load(crsection);
    A(i,:) = temp(2,:); %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
    cd (recoverydir);
    
    
    %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);
    
    %Create new equilibriumprofile
    [IntersectA IntersectB intersections]=CalcintersectionsArea(A(i,:),B(i,:),MSLevel,wavebase,xas);

    [equilibriumprofile stormprofile]=CreateEquilibriumProfile4(i,OverwashV,initialprofile,A,B,wavebase,MSLevel,xas);
    allequilibriumprofiles(i,:)=equilibriumprofile;
    
    % recovery
    recoveryTime=ceil(recoveryTimes(i));
    
    if recoveryTime<=0
        recoveryTime=1;
    end
    
    [Continuity recoveryprofilet recoveryprofileall widthIsland VolIsland HeightIsland distance seatransl MaxHeight]=recoverprof1910(alpha,recoveryTime,stormprofile,equilibriumprofile,MSLevel,wavebase,Tide,theta,T,Beta,xas,sealevelchange,add_sediment);
      
    %write the new cross-section for the next run
    crsection=['z_run_',num2str(i,'%02i'),'.dep'];
    D_low=max(recoveryprofilet);
    paramdir_var = fullfile(fileparts(which(mfilename)),'..','InputData',['Modeltrain_number_',num2str(modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]);

    
    
    
    %% evaluation
    cd(fullfile(inputspectrums,'..','..','..','XBeach'))
    
    zbrun=['zb_run_',num2str(i,'%02i')];
    zb=xb_readvar([zbrun,'.dat'],XB);
    zsrun=['zs_run_',num2str(i,'%02i')];
    zs=xb_readvar([zsrun,'.dat'],XB);
    Hrun=['H_run_',num2str(i,'%02i')];
    H=xb_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.islandheight=HeightIsland;
    Hurricane.recovery.islandwidth=widthIsland;
    Hurricane.recovery.islandwvolume=VolIsland;
    Hurricane.recovery.islandroll_over=distance-XDistance;
    Hurricane.recovery.islandheight_max=MaxHeight;
    Hurricane.equilibriumprofile=equilibriumprofile;
    Analyses{1,i}=Hurricane;
    
    try
        delete('Storm*.inp');
        movefile('*.dat',xbeachrunstormdir);
        movefile('*.bcf',xbeachrunstormdir);
        movefile('*Storm*.inp',xbeachrunstormdir);
        movefile('*.txt',xbeachrunstormdir);
        movefile('*.dep',xbeachrunstormdir);
    catch
        delete('Storm*.inp');
        movefile('*.txt',xbeachrunstormdir);
        movefile('*.dep',xbeachrunstormdir);
    end
    
    cd(modeltraindir);
    
    mkdir(fullfile(xboutputdir,['Modeltrain_number_',num2str(modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]));
    cd(fullfile(xboutputdir,['Modeltrain_number_',num2str(modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]))
    RunProperties=Analyses{1,i};
    save('RunProperties','RunProperties');
    cd(modeltraindir);
    i=i+1;
    %create the new param-file
    while i <= size(storms,2) 
    paramdir_var = fullfile(fileparts(which(mfilename)),'..','InputData',['Modeltrain_number_',num2str(modeltrain_number,'%02i')],['Storm_',num2str(i,'%02i')]);
    boundaryConditions = fullfile(['jonswaplist_storm_',num2str(i,'%02i'),'.txt']);
    surgeConditions=fullfile(['Scaled_Surge_Storm_',num2str(i,'%02i'),'.dat']);
    [params_run]=CreateParam(i,crsection,ngridp,boundaryConditions,surgeConditions,paramdir_var);
    dlmwrite(fullfile(paramdir_var,crsection),[recoveryprofilet;recoveryprofilet;recoveryprofilet],'delimiter',' ','precision','%.10f');
    end
    cd(modeltraindir);
    
end
xboutputdir = fullfile(fileparts(which(mfilename)),'..','Output');
cd(fullfile(xboutputdir,['Modeltrain_number_',num2str(modeltrain_number,'%02i')]));
save('Analyses','Analyses')
cd(modeltraindir);
end
