function result=getDuinAfslag(d,D,ws,WL_t,Hsig_t,Tp_t,n_d,calcMeth,outputDir,BOfactor,norm,Scen,flume)
%GETDUINAFSLAG   1-D dune erosion calculation on a transect
%
%   This is the McTools dune erosion routine. Several methods are available to the users. 
%
%   syntax:
%   result=GETDUINAFSLAG(d,D,ws,WL_t,Hsig_t,Tp_t,n_d,calcMeth,outputDir,BOfactor,norm,Scen,flume)
%
%   input: 
%       d           = input structure with McTools transect data
%       D           = grainsize D50 [m]   
%       Rekenpeil   = stormsurge level [m]
%       H0s         = significant wave height [m]
%       Tp          = peakperiod of the wave spectrum [s] 
%       calcMeth    = specifiec dune erosion calculation mathod (zie setCalcMeth)
%       BOfactor    = exaggeration factor (e.g.: 1.5 yields 50% more erosion)
%       norm        = indicates which norm frequency is used in the calculations
%       Scen        = indicates which climate scenario is used in the calculations '(min)', '(mid)' or '(max)'
%
%   example:
%       D           = 2.2900e-004;
%       ws          = -1;
%       WL_t        = 5.5;
%       Hsig_t      = 9.57;
%       Tp_t        = 16.0900;
%       n_d         = 1;
%       calcMeth    = 'VTV 2006 - nu';
%       outputDir   = 'D:\';
%       BOfactor    = NaN;
%       norm        = '1/100.000';
%       Scen        = '(-)';
%       flume       = 0;
%       d           = readTransectData('Jarkus Data','Noord-Holland','04000','2005');
%       result      = getDuinAfslag(d,D,ws,WL_t,Hsig_t,Tp_t,n_d,calcMeth,outputDir,BOfactor,norm,Scen,flume);
%       UCIT_plotDuinAfslag(result)
%
%   See also getDuinAfslagLZK, getDuinAfslagBO, setCalcMeth
% 

% --------------------------------------------------------------------------
% Copyright (c) WL|Delft Hydraulics 2004-2007 FOR INTERNAL USE ONLY 
% Version:      Version 1.2, June 2005 (Version 1.0, February 2004)
% By:           <M. van Koningsveld (email: mark.vankoningsveld@wldelft.nl)>                                                            
% --------------------------------------------------------------------------

% if ~isfield(d, 'jaar')
%     d.year=str2num(d.year);
% end
% if  ~isfield(d, 'kustvak')
%     d.areacode=str2num(d.areacode);
% end
% if  ~isfield(d, 'raai')
%     d.transectID=str2num(d.transectID);
% end

dold=d;  %dold is the old unmodified JARKUS data ... convenient to show the differences in the plots

guiH2=findobj('tag','UCIT_batchDuinafslagInput');
if ~isempty(guiH2)&& get(findobj(guiH2,'tag','pushbutton5'),'value')~=1;
    % If the batch window is still open prevent normal use of getduinafslag
  
    errordlg('You are still in batchmode!!','UCIT: Error message');
    return    
end 


% *** Versterkingsinvoer (landwaarts)
Versterking = str2double(getINIValue('mctoolbox.ini','Versterking'));  
if isempty(Versterking)
    setINIValue('mctoolbox.ini','Versterking',0);
end 

xVersterking = str2double(getINIValue('mctoolbox.ini','xVersterking'));  
if isempty(xVersterking)
    setINIValue('mctoolbox.ini','xVersterking',-71);
end 


% *** Consolidatie invoer (zeewaarts)
Consolideren = str2double(getINIValue('mctoolbox.ini','Consolideren'));  
if isempty(Consolideren)
    setINIValue('mctoolbox.ini','Consolideren',0);
end 

xConsolidatie = str2double(getINIValue('mctoolbox.ini','xConsolidatie'));  
if isempty(xConsolidatie)
    setINIValue('mctoolbox.ini','xConsolidatie',30);
end 

xMarge = str2double(getINIValue('mctoolbox.ini','xMarge'));  
if isempty(xMarge)
    setINIValue('mctoolbox.ini','xMarge',300);
end 


% *** Suppletie invoer (zeewaarts)
Suppleren = str2double(getINIValue('mctoolbox.ini','Suppleren'));  
if isempty(Suppleren)
    setINIValue('mctoolbox.ini','Consolideren',0);
end 

Begin = str2double(getINIValue('mctoolbox.ini','Begin'));  
if isempty(Begin)
    setINIValue('mctoolbox.ini','Begin',200);
end 

Lengte = str2double(getINIValue('mctoolbox.ini','Lengte'));  
if isempty(Lengte)
    setINIValue('mctoolbox.ini','Lengte',200);
end 

Hoogte = str2double(getINIValue('mctoolbox.ini','Hoogte'));  
if isempty(Hoogte)
    setINIValue('mctoolbox.ini','Hoogte',1);
end 

% Some important parameters are read from the ini file
GridSize   = str2double(getINIValue('mctoolbox.ini','GridSize'));                    
zDuinvoet  = str2double(getINIValue('mctoolbox.ini','zDuinvoet'));                  
Precisie   = str2double(getINIValue('mctoolbox.ini','IterPrecisie'));  % Nauwkeurigheid van de iteraties
optOptions = optimset('TolFun',Precisie);                           % Necessary input for the search routine                               

% Check if we have to run in batch mode
guiH=findobj('tag','UCIT_mainWin');
if strcmp(get(findobj(guiH,'tag','UCIT_batchCommand'),'checked'),'on')
    UCIT_batchGetDuinAfslag;
    result=[];
    return
end

if ~isempty(d)    % Profile data struct d must NOT be empty
%     try % In case of flume this field does not exist.
%         MLW=d.MLW;
%     end
    
    % Boundary conditions for the dune erosion calculation are determined
    % here if they are not already in the argument list of this function.
    % The batch mode for instance will not use this if-end section as it
    % also specifies the boundary conditions to be used.
    if nargin==1
        setINIValue('mctoolbox.ini','FrameTxt6','');
        
        % A dialog will prompt the user to select a method.
        warning off
        try
            [calcMeth,norm,Scen]=UCIT_DuinafslagSelector(str2num(d.areacode));
        catch
            errordlg('No input selected!!','UCIT: Error message');
            return
        end
        warning on
        
        % When the method is known the boundary conditions are extracted
        % from files that are indicated in the ini-file.
        RVW=getRandvoorwaarden(d.areacode,d.transectID,calcMeth);
        
        D=RVW.D;
        Rekenpeil=RVW.Rekenpeil;
        H0s=RVW.H0s;
        Tp=RVW.Tp;
        BOfactor=RVW.BOfactor;
        outputDir='d:\UCIT_Projects\';
        
        % If a different exceedance frequency of the design waterlevel is to be used. The boundary conditions are transformed accordingly here.
        if ~(strcmp(norm,'1/100.000')||strcmp(norm,'Det100')||strcmp(norm,'Det1000')||strcmp(norm,'Det10000'))&&(str2double(d.areacode)==7||str2double(d.areacode)==8||str2double(d.areacode)==9)
            [Rekenpeil,H0s] = getNewNorm(d,Rekenpeil,H0s,norm);
            if ~isnan(BOfactor)
                BOfactor=(100+round(-(BOfactor*100-100)*((log10(UCIT_convertNorm(norm))/2 + 1))))/100;
            end
        end
        
        % Next these boundary conditions are subjected to the selected climate scenarios
        [zDuinvoet,Rekenpeil,H0s,d] = UCIT_getScenario(calcMeth,Scen,d,zDuinvoet,Rekenpeil,H0s);
        
        % The final input is shown to the user one more time. Changes can
        % be made if necessary.
        warning off
        if ~strcmp(calcMeth,'Durosta - nu')&&~strcmp(calcMeth,'Durosta - 50 jaar')&&~strcmp(calcMeth,'Durosta - 100 jaar')&&~strcmp(calcMeth,'Durosta - 200 jaar')
            duinInput=inputdlg({'Rekenpeil [m]','Korreldiameter [m]','Significante golfhoogte [H0s in m]','Piek periode [Tp, in sec.]','Factor voor extra afslag','Directory voor resultaat files'},...
                ['Input: ' calcMeth ' ' Scen ' - ' norm],1,{num2str(Rekenpeil),num2str(D),num2str(H0s),num2str(Tp),num2str(BOfactor),'D:\UCIT_Projects\'});
            if isempty(duinInput)
                errordlg('No input selected!!','UCIT: Error message');
                return
            end
            
            warning on
            
            % The values that are in the input window will be used in the
            % further steps.
            Rekenpeil=str2double(duinInput{1});                                       
            D=str2doublenum(duinInput{2});                                               
            H0s=str2double(duinInput{3});                                             
            Tp=str2double(duinInput{4});                                              
            BOfactor=str2double(duinInput{5});
            outputDir=duinInput{6};
            ws = -1;              % negative value means it will be determined automatically later
%             n_d = 1;
%             flume = 0;

        else
            defaultStormperiod = 5;
            WL_t = [defaultStormperiod Rekenpeil];
            ws = -1;              % negative value means it will be determined automatically later
            Hsig_t = [defaultStormperiod H0s];                                             
            Tp_t = [defaultStormperiod Tp];     
            n_d = 1;
            flume = 0;
        end
    else
        % WL_t, Hsig_t, Tp_t are time series (used for DUROSTA and possible future rules
        % Here the maximum value of the time series is assumed to be valid for rules requiring constant input values
        if size(Hsig_t,2)>1
            H0s                 = max(Hsig_t(:,2));
            Rekenpeil           = max(WL_t(:,2));
            Tp                  = max(Tp_t(:,2));
        else
            H0s                 = Hsig_t;
            Rekenpeil           = WL_t;
            Tp                  = Tp_t;
            defaultStormperiod  = 5;
            WL_t                = [defaultStormperiod WL_t];
            ws                  = -1;              % negative value means it will be determined automatically later
            Hsig_t              = [defaultStormperiod Hsig_t];                                             
            Tp_t                = [defaultStormperiod Tp_t];     
 %           n_d = 1;
 %           flume = 0;
        end  
        
    end
    
    % When running in batchmode we still need to subject the selected boundary
    % conditions to the exeendance frequency and climate scenario functions.
    % We chose to put only the basic input in the batch file. This way we
    % can use the same batch input and play around with different
    % algorithms for the exceedance frequency and the climate scenarios.
    if ~isempty(findobj('tag','UCIT_batchDuinafslagInput')) 
        % if another exceedance frequency is used get the new conditions
        if ~(strcmp(norm,'1/100.000')||strcmp(norm,'Det100')||strcmp(norm,'Det1000')||strcmp(norm,'Det10000'))&&ismember(str2double(d.areacode),[7 8 9])
            if ~isnumeric(calcMeth)
                Temp_cM=setCalcMeth(calcMeth,1);
                if Temp_cM==36||Temp_cM==37||Temp_cM==38||Temp_cM==39
                    [Rekenpeil,H0s,Tp] = getNewNormVTV2006(d,Temp_cM+1,norm);
                else
                    [Rekenpeil,H0s]    = getNewNorm(d,Rekenpeil,H0s,norm);
                end
            end
        end
        
        % Subject these conditions to the climate scenario if necessary
        [zDuinvoet,Rekenpeil,H0s,d] = UCIT_getScenario(calcMeth,Scen,d,zDuinvoet,Rekenpeil,H0s);
    end
    
    % *** Start the calculation ***
    % In the next routines we shall use the basic modified input 'dmod' stored here, in stead of d, 
    % as d is polluted later due to all kinds of interpolation efforts. It
    % is important to have this statement at this location, as the implementation 
    % of the climate scenarios alters the profile.
    dmod=d; 
    xc=d.xi;
    zc=d.ze;
    
    startTime=now;  % The start time is recorded to assess the duration of the calculation.
    switch calcMeth
        case {'Leidraad Duinafslag - nu', 'Leidraad Duinafslag - 50 jaar', 'Leidraad Duinafslag - 100 jaar', 'Leidraad Duinafslag - 200 jaar'}
            % This is the dune erosion rule from the Leidraad Duinafslag (TAW 1984)
            
            % If seaward consolidation is turned on and the field xConsolidatie 
            % in UCIT_Options is set to -999, the automatic routine will start working. 
            % Otherwise the progam will use the user input value.
            Dold=D; % For this method D will always be equal to Dold. Dold has to be defined for output reasons!
            
            if (Consolideren==1 && xConsolidatie==-999)
                % The next two lines automatically determine the amount of sand that is needed to shift the afslaglijn seaward to a certain location 
                optOptions2=optimset('TolFun',1); % Settings for automatic dune strengthening routine
                if exist([outputDir 'ConsolidatieTargetsLZK.txt'])
                    ConsolidatieTargetsLZK=load([outputDir 'ConsolidatieTargetsLZK.txt']);
                    idGoal=ConsolidatieTargetsLZK(:,1)==dmod.transectID;
                    Goal=ceil(ConsolidatieTargetsLZK(idGoal,2)); % we maken er een mooi rond getal van
                    
                    id = max(dmod.ze>zDuinvoet);                                               % meest zeewaarts gelegen z>+3m
                    pn = polyfit(dmod.xi(id:id+1),dmod.ze(id:id+1),1);
                    xDuinvoet = (zDuinvoet-pn(2))/pn(1);   
                    xMarge=abs(Goal-xDuinvoet);
                else
                    Goal=-50;
                end
                %                 xConsolidatie=Goal+200;                     % We start with a large value and let the routine work its way down.
                xConsolidatie=200;                     % We start with a large value and let the routine work its way down.
                % getDuinafslagLZK_cons is identical to getDuinafslagLZK
                % ... only the position of the in and output argument has
                % been changed to enable automatic optimisation.
                xConsolidatie = ...
                    fminsearch('getDuinafslagLZK_cons',xConsolidatie,optOptions2,Scen,calcMeth,dmod,xc,zc,Rekenpeil,D,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions,'abs(BProfiel(1,idB2(1))-(varargin{1}))',Goal);
            end
            
            [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                getDuinafslagLZK(Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
            
            if (criterium>5 && Versterking ==0 && Consolideren==0 && Suppleren==0)
                disp(['No result could be obtained ... the automatic landward strengthening routine is attempted ...'])
                Versterking = 1;
                [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                    getDuinafslagLZK(Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
                if (isempty(BatchReport)&&criterium<=5)
                    BatchReport='No result could be obtained at first. A result was forced by automatically strengthening the dune at the landward side!';   
                end
            elseif (Consolideren==1 && ~isempty(Xc_cons) && ~isempty(Vol_cons))
                disp('The dune was strengthened on the seaward side!')  
                BatchReport='The dune was strengthened on the seaward side!';  
            end  
            
        case {'Beheerdersoordeel - nu', 'Beheerdersoordeel - 50 jaar', 'Beheerdersoordeel - 100 jaar', 'Beheerdersoordeel - 200 jaar'}
            % This is the dune erosion rule from the Crash Actions (also used in the Beheerdersoordeel Veiligheid 2003) 
            
            % If seaward consolidation is turned on and the field xConsolidatie 
            % in UCIT_Options is set to -999, the automatic routine will start working. 
            % Otherwise the progam will use the user input value.
            
            % Start with a brief input check
            if isnan(BOfactor)
                errordlg('Specify factor for extra erosion!','No value specified');
                return
            end
            
            Dold=D; % all routines will be started with Dold ... to avoid potential confusion
            try
                setINIValue('mctoolbox.ini','Dtry',num2str(D));
            end
            
            
            if (Consolideren==1 && xConsolidatie==-999)
                % The next two lines automatically determine the amount of sand that is needed to shift the afslaglijn seaward to a certain location 
                optOptions2=optimset('TolX',5,'TolFun',1); % Settings for automatic dune strengthening routine
                if exist([outputDir 'ConsolidatieTargetsBO.txt'])
                    ConsolidatieTargetsBO=load([outputDir 'ConsolidatieTargetsBO.txt']);
                    idGoal=find(ConsolidatieTargetsBO(:,1)==dmod.transectID);
                    Goal=ceil(ConsolidatieTargetsBO(idGoal,2)); % we maken er een mooi rond getal van
                    D=ConsolidatieTargetsBO(idGoal,3);          % Here we enter the preset grainsize
                    try
                        setINIValue('mctoolbox.ini','Dtry',num2str(D));
                    end
                    
                    id = find(dmod.ze>zDuinvoet, 1, 'last' );                                               % meest zeewaarts gelegen z>+3m
                    pn = polyfit(dmod.xi(id:id+1),dmod.ze(id:id+1),1);
                    xDuinvoet = (zDuinvoet-pn(2))/pn(1);   
                    xMarge=abs(Goal-xDuinvoet);
                else
                    Goal=-50;
                end
                % hier ergens maken dat ook de korreldiameter wordt
                % meegegeven
                
                xConsolidatie=200;                     % We start with a large value and let the routine work its way down.
                
                % getDuinafslagBO_cons is identical to getDuinafslagBO ... only the position of the in and output argument has been changed to enable automatic optimisation.
                %[xConsolidatie,dum,exitflag,fin] = ...
                %    fminsearch('getDuinafslagBO_cons',xConsolidatie,optOptions2,Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions,'abs(BProfiel(1,idB2(1))-(varargin{1}))',Goal);

                [xConsolidatie] = getConsolidationDistance(Goal,Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
                
            end
            
            if exist([outputDir 'Grainsize.txt'])
                Grainsize=load([outputDir 'Grainsize.txt']);
                idGrainsize=Grainsize(:,1)==dmod.transectID;
                D=Grainsize(idGrainsize,2);          % Here we enter the preset grainsize
                try
                    setINIValue('mctoolbox.ini','Dtry',num2str(D));
                end
            end

            [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,D,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                getDuinafslagBO(Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
            
            if (criterium>5 && Versterking ==0 && Consolideren==0 && Suppleren==0)
                disp('No result could be obtained ... the automatic landward strengthening routine is attempted ...')
                Versterking = 1;
                [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,D,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                    getDuinafslagBO(Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
                if (isempty(BatchReport)&&criterium<=5)
                    BatchReport='No result could be obtained at first. A result was forced by automatically strengthening the dune at the landward side!';   
                end
            end
            
        case {'VTV 2006 - nu', 'VTV 2006 - 50 jaar', 'VTV 2006 - 100 jaar', 'VTV 2006 - 200 jaar'}
            % This is the dune erosion rule developed for the VTV 2006 (includes peak period in the parabola formulation) 
            
            Dold=D; % For this method D will always be equal to Dold. Dold has to be defined for output reasons!
            %             dold=d;
            zDuinvoet=Rekenpeil;
            [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                getDuinafslag_VTV2006(Scen,calcMeth,dmod,xc,zc,Rekenpeil,D,ws,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
            
        case {'Koster1d - nu', 'Koster1d - 50 jaar', 'Koster1d - 100 jaar', 'Koster1d - 200 jaar'}
            zDuinvoet=Rekenpeil;
            Dold=D; % all routines will be started with Dold ... to avoid potential confusion
            %             dold=d;
            try
                setINIValue('mctoolbox.ini','Dtry',num2str(D));
            end
            
            s_op = (2*pi*H0s)/(9.81*Tp^2);
            if (s_op < 0.02)
                BOfactor = 1 + 0.40;
            elseif (s_op > 0.05)
                BOfactor = 1 - 0.36;
            else
                BOfactor = 1 + 266.6667*(s_op.^2)-44.13333*s_op+1.1807;
            end                
            disp(['Additional Erosion Koster-1D: ' num2str(BOfactor, '%8.2f')]);
            [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,D,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                UCIT2006_getDuinafslagKOS1d(Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,ws,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
            
            if (criterium>5 && Versterking ==0 && Consolideren==0 && Suppleren==0)
                disp('No result could be obtained ... the automatic landward strengthening routine is NOT attempted ...')
                %                 Versterking = 1;
                %                 [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,D,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                %                     UCIT2006_getDuinafslagBO(Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
                if (isempty(BatchReport)&&criterium<=5)
                    BatchReport='No result could be obtained at first. A result was forced by automatically strengthening the dune at the landward side!';   
                end
            end
            
            % Koster alternative 3
        case {'Koster3 - nu', 'Koster3 - 50 jaar', 'Koster3 - 100 jaar', 'Koster3 - 200 jaar'}
            Dold=D; % For this method D will always be equal to Dold. Dold has to be defined for output reasons!
            %             dold=d;
            zDuinvoet=Rekenpeil;
            if flume == 1 % Model conditions, Tp will be much smaller than 12 s 
%                 s_op = (2*pi*H0s)/(9.81*Tp^2); % Calculate wave steepness
                s_op_ref = 0.034;
                %if s_op < s_op_ref % Reference wave steepness
                Tp_ref = sqrt((2*pi*H0s)/(9.81*s_op_ref));
                w_coeff = Tp/Tp_ref;
                ws = ws/w_coeff;
                %end
            elseif (flume == 0 || flume == 2) % Model translated to prototype conditions
                %if Tp>12
                Tp_ref = 12;
                w_coeff = Tp/Tp_ref;
                ws = ws/w_coeff;
                %end
            end
            
            [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                UCIT2006_getDuinafslag(Scen,calcMeth,dmod,xc,zc,Rekenpeil,D,ws,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
            
        case {'DuinafslagRegel2006 - nu', 'DuinafslagRegel2006 - 50 jaar', 'DuinafslagRegel2006 - 100 jaar', 'DuinafslagRegel2006 - 200 jaar'}
            Dold=D; % For this method D will always be equal to Dold. Dold has to be defined for output reasons!
            %             dold=d;
            zDuinvoet=Rekenpeil;
            [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                UCIT2006_getDuinafslag(Scen,calcMeth,dmod,xc,zc,Rekenpeil,D,ws,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
            
        case {'WLConcept1A - nu', 'WLConcept1A - 50 jaar', 'WLConcept1A - 100 jaar', 'WLConcept1A - 200 jaar'}
            Dold=D; % For this method D will always be equal to Dold. Dold has to be defined for output reasons!
            %             dold=d;
            zDuinvoet=Rekenpeil;
            [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                UCIT2006_getDuinafslagWL1A(Scen,calcMeth,dmod,xc,zc,Rekenpeil,D,ws,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
            
        case {'WLConcept1B - nu', 'WLConcept1B - 50 jaar', 'WLConcept1B - 100 jaar', 'WLConcept1B - 200 jaar'}
            Dold=D; % For this method D will always be equal to Dold. Dold has to be defined for output reasons!
            %             dold=d;
            zDuinvoet=Rekenpeil;
            [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                UCIT2006_getDuinafslagWL1B(Scen,calcMeth,dmod,xc,zc,Rekenpeil,D,ws,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
            
        case {'WLConcept2A - nu', 'WLConcept2A - 50 jaar', 'WLConcept2A - 100 jaar', 'WLConcept2A - 200 jaar'}
            zDuinvoet=Rekenpeil;
            Dold=D; % all routines will be started with Dold ... to avoid potential confusion
            %             dold=d;
            try
                setINIValue('mctoolbox.ini','Dtry',num2str(D));
            end
            
            s_op = (2*pi*H0s)/(9.81*Tp^2);
            if (s_op < 0.015)
                BOfactor = 1 + 0.31;
            elseif (s_op > 0.04)
                BOfactor = 1 - 0.1;
            else
                BOfactor = 1 - 16.501*(s_op-.034);
            end                
            disp(['Additional Erosion WL-Concept 2A: ' num2str(BOfactor, '%8.2f')]);
            
            [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,D,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                UCIT2006_getDuinafslagKOS1d(Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,ws,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
            
            if (criterium>5 && Versterking ==0 && Consolideren==0 && Suppleren==0)
                disp('No result could be obtained ... the automatic landward strengthening routine is NOT attempted ...')
%                 Versterking = 1;
                %                 [d, Rekenpeil,Tp, H0s,xrp,id_AfslagZone,idB2,idA2, idG2,AfslagProfiel,AProfiel,BProfiel,GProfiel,GrensProfiel, xDuinvoet ,XVerschuiving,XVerschuivingP, XVerschuivingG,yTempJarkus,erz,erl,depl,depz,AVolume, BVolume, TVolume, GVolume, Er, Dep,id_yTempVerschil_neg,calcMeth,criterium,Scen,zDuinvoet,D,BatchReport,Xc_cons,Vol_cons,Xc_verst,Vol_verst] = ...
                %                     UCIT2006_getDuinafslagBO(Scen,calcMeth,dmod,xc,zc,Rekenpeil,Dold,H0s,Tp,GridSize,Versterking,xVersterking,Consolideren,xConsolidatie,xMarge,Suppleren,Begin,Lengte,Hoogte,zDuinvoet,BOfactor,optOptions);
                if (isempty(BatchReport)&&criterium<=5)
                    BatchReport='No result could be obtained at first. A result was forced by automatically strengthening the dune at the landward side!';   
                end
            end
            
        case {'Durosta - nu', 'Durosta - 50 jaar', 'Durosta - 100 jaar', 'Durosta - 200 jaar'}
            Dold=D; % For this method D will always be equal to Dold. Dold has to be defined for output reasons!
            BatchReport=[];
            Xc_cons=0;
            Vol_cons=0;
            Xc_verst=0;
            Vol_verst=0;
            zDuinvoet=Rekenpeil;
            x = (min(xc):GridSize:max(xc))';                                                % Interpolate to finer grid to obtain more accurate estimate
            warning off
            z = interp1(xc,zc,x);
            warning on
            id = find(z>zDuinvoet, 1, 'last' );                                                    % meest zeewaarts gelegen z>+3m
            pn = polyfit(x(id:id+1),z(id:id+1),1);
            xDuinvoet = (zDuinvoet-pn(2))/pn(1);                                            % Via lokale lineaire fit
            
            %**** Vul hier in wat je wilt doen voor zeewaartse consolidatie
            if Versterking == 1 %(landwaarts)
                HoogteGrensprofiel=Rekenpeil+0.12*Tp*sqrt(H0s);
                xVersterking=-500;                                                          % We take a large number here so that we always get results. Then chop off what is not needed at the end of this file.
                JarkusNew=getVersterkingsVolume(xVersterking, d, HoogteGrensprofiel+0.1);   % verschuiving profiel, d (jarkus), xduinvoet, zoekmarge 
                d.xi = JarkusNew.xi;
                d.ze = JarkusNew.ze;
                Xc_verst=xVersterking;
                Vol_verst=JarkusNew.VersVol;
            end
            
            if Consolideren == 1 %(zeewaarts)
                JarkusNew=getConsolidatieVolume(xConsolidatie,d,xDuinvoet,xMarge);     % verschuiving profiel, d (jarkus), xduinvoet, zoekmarge
                d.xi = JarkusNew.xi;
                d.ze = JarkusNew.ze;
                Xc_cons=xConsolidatie;
                Vol_cons=JarkusNew.ConsVol;
            end
            
            if Suppleren == 1
                JarkusNew=getConsolidatieVolume2(d,Begin,Lengte,Hoogte);     % verschuiving profiel, d (jarkus), xduinvoet, zoekmarge
                d.xi = JarkusNew.xi;
                d.ze = JarkusNew.ze;
            end
            
            [DurostaBaseDir,DurostaProjectName]=UCIT_createDurostaInputFile(d,WL_t,D,ws,Hsig_t,Tp_t,outputDir,n_d,flume);
            
            % now call durosta from its executable location and pass the input file as an argument.
            if flume >= 1
                c1=['!' which('durosta.exe') ' ' DurostaBaseDir DurostaProjectName '.dat'];
                c2=[DurostaBaseDir DurostaProjectName];
            else
                c1=['!' which('durosta.exe') ' ' DurostaBaseDir filesep DurostaProjectName filesep DurostaProjectName '.dat'];
                c2=[DurostaBaseDir filesep DurostaProjectName filesep DurostaProjectName];
            end
            evalin('base',c1);
            
            % handling of the resultfile
            fid=tekal('open',[c2  '.tek']);
            DurostaResults=tekal('read',fid,1);
            DurostaResults(:,1)=-DurostaResults(:,1);
            
            % plot preparation
            hf=figure(2005);clf;
            set(hf,'tag','plotDurosta');
            
            RaaiInformatie=['Raai: ' num2str(d.transectID) ' Jaar: ' num2str(d.year) ' Methode: ' num2str(calcMeth) ' ' num2str(Scen) ' Norm: ' norm ' per jaar'];
            set(hf,'Name', RaaiInformatie,'NumberTitle','Off','Units','normalized','Position', [0 0.5 0.5 0.41]);
            
            hold off;
            ph(1)=plot(dold.xi,dold.ze,':b');
            legendtext{1}='Measured JARKUS profile';
            hold on;
            ph(2)=plot(d.xi,d.ze,'b');
            legendtext{2}='Modified JARKUS profile';
            
            ph(3)=plot(DurostaResults(:,1),DurostaResults(:,3),'-or');
            set(ph(3),'markersize',2,'markerfacecolor','r');
            legendtext{3}='Unibest-DE post-storm estimation';
            
            axis([min(DurostaResults(:,1)) max(DurostaResults(:,1)) -25 25]); 
            set(gca,'XDir','reverse');
            
            lh=legend(ph,legendtext,2);
            set(lh,'fontsize',7);
            xlabel('Kustdwarse afstand [m tov RSP]');
            ylabel('Hoogteligging [m tov NAP]');
            title('Duinafslagberekening Unibest-DE');
            
            savePlot=0;
            if savePlot==1
                UCIT_preparePlot('plotDurosta')
                print(gcf,'-dpng','-r250',[c2 '.png']);     
            end
            
            VisualDune = str2double(getINIValue('mctoolbox.ini','VisualDune'));  
            if isempty(VisualDune)
                setINIValue('mctoolbox.ini','VisualDune','1');
            end 
            if VisualDune==1
                c1=['!' which('ani-de.exe') ' ' c2 '.daf'];
                evalin('base',c1);
            end
            
            criterium=10;
            
    end
    disp(' ')
    disp(['Derived in ' datestr(now-startTime,13)]);
    disp(' ')
    
    % ----------------------------------------------------------------------------------------------------------------------------------------
    % *** Store results in result structure
    if ~strcmp(calcMeth,'Durosta - nu')&&~strcmp(calcMeth,'Durosta - 50 jaar')&&~strcmp(calcMeth,'Durosta - 100 jaar')&&~strcmp(calcMeth,'Durosta - 200 jaar')
        if criterium>5
            [TVolume, xrp, AfslagProfiel, id_AfslagZone, idA2, AProfiel, BProfiel, XVerschuiving, XVerschuivingP, yTempJarkus, erz, erl, AVolume, Er, Dep, id_yTempVerschil_neg, depl, depz, Xhor1, Xhor, Yhor1, Yhor]=deal(NaN);
            disp(' ');
            disp('*** No result could be obtained for this scenario ***');
            disp(' ');
            disp(' ');
            
            result=[];    
        else
            % Store all relevant variables in the result structure
            disp(' ');
            disp('*** Test OK! ***');
            disp(' ');
            disp(' ');
            if isempty(BatchReport)
                BatchReport='Test OK!';
            end
            result.dold = dold;
            result.d = d;
            result.Rekenpeil = Rekenpeil;
            result.Tp = Tp;
            result.H0s = H0s;
            result.norm = norm;
            result.D = D;
            
            result.xrp = xrp;
            
            result.id_AfslagZone = id_AfslagZone;
            result.idB2 = idB2;
            result.idA2 = idA2;
            result.idG2 = idG2;
            
            result.AfslagProfiel=AfslagProfiel;
            result.AProfiel = AProfiel;
            result.BProfiel = BProfiel;
            result.GProfiel = GProfiel;
            result.GrensProfiel = GrensProfiel;
            
            
            result.xDuinvoet = xDuinvoet; 
            result.zDuinvoet = zDuinvoet;
            result.XVerschuiving = XVerschuiving;
            result.XVerschuivingP = XVerschuivingP;
            result.XVerschuivingG = XVerschuivingG;
            
            result.yTempJarkus = yTempJarkus;
            
            result.erz = erz;
            result.erl = erl;
            result.depl = depl;
            result.depz = depz;
            
            result.AVolume = AVolume;
            result.BVolume = BVolume;
            result.TVolume = TVolume;
            result.GVolume = GVolume;
            result.Er = Er;
            result.Dep = Dep;
            
            result.id_yTempVerschil_neg = id_yTempVerschil_neg;
            
            result.calcMeth = calcMeth;
            result.Scen = Scen;
            
            %UCIT_StoreResult(result);
        end
        
        % ----------------------------------------------------------------------------------------------------------------------------------------
        % *** Wegschrijven van enkele resulaten naar file
        
        % Voorbereiden enkele outputs
        xP=XVerschuiving+XVerschuivingP;
        try
            [xyRDP] = xyRSP2RD(d.xRD,d.yRD,d.GRAD,360,xP);
        catch
            xyRDP=[nan nan];
        end
        if isnan(xP)
            xR=nan;
            xyRDR=[nan nan];
        else
            try
                xR=BProfiel(1,idB2(1));
                [xyRDR] = xyRSP2RD(d.xRD,d.yRD,d.GRAD,360,xR);
            catch
                xR=nan;
                xyRDR=[nan nan];
            end
        end
        
        cM = setCalcMeth(calcMeth,1);
        
        % norm met een kleine letter is een andere dan Norm
        % norm gaat in de batch file
        norm=UCIT_convertNorm(norm);  % from  1/100.000 to 0.0001
        if ischar(Scen)
            Scen=UCIT_convertScen(Scen);
        end
        
        % ***  Wegschrijven van de belangrijkste resultaten
        % Jaar  Kustvak  Raai  Methode  Rekenpeil  H0s  Tp  Xp  XpParijsCoordinaatX  XpParijsCoordinaatY Xr  XrParijsCoordinaatX  XrParijsCoordinaatY AVolume  TVolume  E  norm  Scen 
        if ~exist([outputDir filesep 'batchresults.txt'])
            if ~exist(outputDir,'dir')
                mkdir(outputDir)
            end
            fid=fopen([outputDir filesep 'batchresults.txt'],'w');
            fprintf(fid,'%s\n','% Jaar  Kustvak  Raai  Methode  Rekenpeil  H0s  Tp  Xp  XpParijsCoordinaatX  XpParijsCoordinaatY Xr  XrParijsCoordinaatX  XrParijsCoordinaatY AVolume  TVolume  E  norm  Scen ');
            fclose(fid);
        end
        
        fid=fopen([outputDir filesep 'batchresults.txt'],'a');
        fprintf(fid,'%i   %i   %i   %i   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %8.6f    %i\n',str2double(d.year),str2double(d.areacode),str2double(d.transectID),cM,Rekenpeil,H0s,Tp,xP,xyRDP(1,1),xyRDP(1,2),xR,xyRDR(1,1),xyRDR(1,2),AVolume,TVolume,Er,norm,Scen);
        fclose(fid);
        
        % *** Wegschrijven redenering per resultaat
        if ~exist([outputDir filesep 'batchreports.txt'])
            fid=fopen([outputDir filesep 'batchreports.txt'],'w');
            fprintf(fid,'%s\n','% Jaar  Kustvak  Raai  Methode  Rekenpeil  H0s  Tp       norm        BOfactor      Dold         Dnew     Scen    Xp          Xr       Xc_cons   Vol_cons   Xc_verst  Vol_verst  BatchReport');
            fclose(fid);
        end
        
        % Jaar  Kustvak  Raai  Methode  Rekenpeil  H0s  Tp   norm  Dold  Dnew  Scen Xp Xr Xc_cons  Vol_cons  Xc_verst  Vol_verst  BatchReport
        fid=fopen([outputDir filesep 'batchreports.txt'],'a');
        fprintf(fid,'%i   %i   %i   %i   %8.2f   %8.2f   %8.2f   %8.6f   %8.8f   %8.8f   %8.8f   %i   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %8.2f   %s\n',str2double(d.year),str2double(d.areacode),str2double(d.transectID),cM,Rekenpeil,H0s,Tp,norm,BOfactor,Dold,D,Scen,xP,xR,Xc_cons,Vol_cons,Xc_verst,Vol_verst,BatchReport);
        fclose(fid);
    else
        result=[];    
        
    end
else
    result=[];    
end 

