function data = csb_analyse_utc(varargin)
%CSB_ANALYSE_UTC  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = csb_analyse_utc(varargin)
%
%   Input: For <keyword,value> pairs call csb_analyse_utc() without arguments.
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   csb_analyse_utc
%
%   See also

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2013 Deltares
%       Kees den Heijer
%
%       Kees.denHeijer@deltares.nl
%
%       P.O. Box 177
%       2600 MH  DELFT
%       The Netherlands
%       Rotterdamseweg 185
%       2629 HD  DELFT
%
%   This library is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this library.  If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://www.OpenEarth.eu">OpenEarthTools</a>.
% OpenEarthTools is an online collaboration to share and manage data and
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute
% your own tools.

%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
% Created: 03 Sep 2013
% Created with Matlab version: 7.14.0.739 (R2012a)

% $Id: $
% $Date: $
% $Author: $
% $Revision: $
% $HeadURL: $
% $Keywords: $

%%
OPT = struct(...
    'maindir', 'd:\runs2_Kees\',...
    'target', 87947544,... .daf file maller than target is crashed somewhere
    'iniProf', 'suppleties\iniProf.bot' ...
    );
% return defaults (aka introspection)
% if nargin==0;
%     varargout = {OPT};
%     return
% end
% overwrite defaults with user arguments
OPT = setproperty(OPT, varargin);
%% code
% find all result files
D = dir2(OPT.maindir,...
    'file_incl', 'utc001.daf$',...
    'no_dirs', true);
% calculate relative file size wrt a fully successful run
prog = [D.bytes] / OPT.target;

% read base profile (without nourishment)
iniproffname = fullfile(OPT.maindir, OPT.iniProf);
[x00, z00] = readBOT(iniproffname);
idx = true(size(D));

% calculate indicators (before nourishment)
x_df00 = getdunefoot(x00, z00);
x_mkl00 = getmkl(x00, z00);

data = struct(...
    'Time_horizon', {},...
    'Nourishment_Type', {},...
    'Nourishment_Volume', {},...
    'Dune_foot_change', {},...
    'MKL_change', {},...
    'irun', {} ...
    );

for i = 1:length(D)
    irun = regexp(D(i).pathname, 'run\d+', 'match');
    % read initial profile (including nourishment)
    botfile = fullfile(D(i).pathname, sprintf('%s.bot', irun{:}));
    [x0, z0] = readBOT(botfile);
    x_df0 = getdunefoot(x0, z0);
    x_mkl0 = getmkl(x0, z0);

    % read profile after 1, 5 and 10 years of simulation
    fname = fullfile(D(i).pathname, D(i).name);
    [x, z, t] = readdata(fname);
    
    if ~issorted(t)
        idx(i) = false;
        continue
    end

    % calculate indicators
    x_df = cellfun(@(zt) getdunefoot(x, zt), num2cell(z, 2));
    x_mkl = cellfun(@(zt) getmkl(x, zt), num2cell(z, 2));
    

    
%     if any(isnan(x_df))
%         dbstopcurrent
%     end
    
    % calculate nourishment volume
    V = jarkus_getVolume(x00, z00,...
        'x2', x0,...
        'z2', z0);
    % find the maximum level of the nourishment
    znmax = max(z00(z00 ~= z0));
    
    for j = 1:length(x_df)
        if ~any(isnan([x_df(j) x_mkl(j)])) &&...
                max(abs(z(end,:) - z(1,:)))<6 && ...% exclude cases with instable very high peaks in profile
                all(z(j,:) <= max(z0)+1)
            data(end+1) = struct(...
                'Time_horizon', t(j)/365.24,...
                'Nourishment_Type', znmax>0,...
                'Nourishment_Volume', V,...
                'Dune_foot_change', x_df00-x_df(j),...
                'MKL_change', x_mkl00-x_mkl(j),...
                'irun', str2num(irun{1}(4:end)) ...
                );
        end
    end
    fprintf('%s\t%.2f\t%.2f\t%i\n', irun{1}, V, znmax, znmax>0);
end

sum(idx)/length(idx)

function [x, z, t] = readdata(fname)
fi = qpfopen(fname);
DataProps = qpread(fi);
% read x grid
Grid = qpread(fi, DataProps(1), 'griddata');
% read bottom profile after 1, 5 and 10 years of simulation
data = qpread(fi, DataProps(3), 'data', [122, 609, 1218]);
x = Grid.X;
z = data.Val;
t = data.Time;

function x_df = getdunefoot(x, z)
xcr = findCrossings(x, z(:), minmax(x'), [3 3]);
if isempty(xcr)
    xcr = NaN;
end
x_df = min(xcr);

function x_mkl = getmkl(x, z)
if all(z<3)
    x_mkl = NaN;
else
    x_mkl = jarkus_getMKL(flipud(-x), flipud(z(:)), 3, -4.4);
    x_mkl = -x_mkl;
end