% 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,
for principle=3:3
    [hurricanetrain]=getHurricanetrain();
    cd([cd,'\',hurricanetrain.intensity{principle,1}]);
    
    addpath ([cd, '\jonswap_files'])
    addpath D:\pruis\XBeach\MatlabScripts\benodigd
    addpath D:\pruis\XBeach\MatlabScripts\benodigd\Visualisation
    
    

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

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

    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);

        [Continuity recoveryprofilet recoveryprofileall widthIsland VolIsland HeightIsland distance seatransl MaxHeight]=recoverprof(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
        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('Analises')
    cd ..
end
