function varargout = desb_xbeach_artificial_jarkus(varargin)
%DESB_XBEACH_ARTIFICIAL_JARKUS  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = desb_xbeach_artificial_jarkus(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   desb_xbeach_artificial_jarkus
%
%   See also 

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2012 Delft University of Technology
%       Kees den Heijer
%
%       C.denHeijer@TUDelft.nl	
%
%       Faculty of Civil Engineering and Geosciences
%       P.O. Box 5048
%       2600 GA Delft
%       The Netherlands
%
%   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: 02 Mar 2012
% Created with Matlab version: 7.7.0.471 (R2008b)

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

%%
OPT = struct(...
    'id1', [4001740 7001503 7003775],...
    'id2', [4001740 7001503 7003775],...
    'year', 2010,...
    'contour', -4,... % contour to connect 2 profile parts
    'zmin', -20,... % deep water boundary
    'zmin_extend', 'first',... % put boundary at most nearshore ('first') or most offshore ('last') zmin crossing
    'irun', 2,...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations'),...
    'WL', 5,...
    'Hsig', 9,...
    'Tp', 12,...
    'D50', 225e-6,...
    'duration', 5*3600,...
    'morfac', 5);

OPT = setproperty(OPT, varargin);

%%
id = unique([OPT.id1(:); OPT.id2(:)]);
[n1 n2] = deal(length(OPT.id1), length(OPT.id2));
n = n1*n2;

[xr zr] = deal(repmat({[]}, length(id), 1));
for i = 1:length(id)
    if id ~= 0
        tr = jarkus_extend_with_vaklodingen(id(i), OPT.year, 'jarkus_extend', true, 'debug', true);
        xr{i} = tr.cross_shore;
        zr{i} = tr.altitude;
    else
        [xr{i} zr{i}] = SimpleProfile('zmin', OPT.zmin);
    end
end

% create combinations
id1 = repmat(OPT.id1(:)', n/n1, 1);
id2 = repmat(OPT.id2(:), 1, n/n2);
[id1 id2] = deal(id1(:), id2(:));

[xtmp ztmp] = deal(repmat({[]}, n, 1));
xc = zeros(size(xtmp));
[x z] = deal([]);
for i = 1:n
    [xtmp{i} ztmp{i}] = jarkus_combine_profiles(...
        xr{id==id1(i)}, zr{id==id1(i)},...
        xr{id==id2(i)}, zr{id==id2(i)});
    idcontour = find(ztmp{i}<OPT.contour, 1, 'first') + [-1 0];
    xc(i) = interp1(ztmp{i}(idcontour), xtmp{i}(idcontour), OPT.contour);
    xtmp{i} = xtmp{i} - xc(i);
    x = unique([x(:); xtmp{i}(:)]);
end

%%
for i = 1:n
    z(:,i) = interp1(xtmp{i}, ztmp{i}, x);
    if strcmp(OPT.zmin_extend, 'first')
        xtoe(i) = x(find(z(:,i)<OPT.zmin, 1, 'first')-1);
    elseif strcmp(OPT.zmin_extend, 'last')
        xtoe(i) = x(find(z(:,i)>OPT.zmin, 1, 'last'));
    else
        error
    end
    % extend the seaward side of the profile using the zmin level
    z(x>xtoe(i),i) = OPT.zmin;
    % extend the landward side of the profile using the most landward
    % available level
    z(isnan(z(:,i)),i) = z(find(~isnan(z(:,i)), 1, 'first'),i);
end
z(x>max(xtoe)+150,:) = [];
x(x>max(xtoe)+150) = [];

%%
[x ix] = sort(-x);
z = z(ix,:);
zm = nanmax(z,2);

xb_setpref('grid_finalise', {'zmin', OPT.zmin})

xbi = xb_generate_model(...
    'bathy', {'x', x, 'z', zm, 'xgrid', {'dxmin', 2}, 'finalise', {'seaward_extend'}, 'superfast', true},...
    'tide', {'front' OPT.WL},...
    'waves', {...
    'type', 'jonswap',...
    'Hm0', OPT.Hsig,...
    'Tp', OPT.Tp,...
    'duration', 1800*OPT.morfac},...
    'wavegrid', {'nbins' 1},...
    'settings', {...
    'morfac', OPT.morfac,...
    'posdwn', -1,...
    'tintg', 100,...
    'tstart', 4000,...
    'morstart', 4000,...
    'tstop', 4000+ceil(OPT.duration/OPT.morfac)*OPT.morfac,...
    'D50', OPT.D50,...
    'D90', OPT.D50*1.5,...
    'random' 0,...
    'outputformat', 'netcdf',...
    'meanvars', {'H' 'zs' 'zb'},...
    });
% remove swtable from structure
xbi = xs_del(xbi, 'swtable');
xbi = xs_set(xbi, 'bcfile.duration', xb_get(xbi, 'tstop'));

xori = xs_get(xbi, 'xori');
xfile = xs_get(xbi, 'xfile.xfile');
depfile = xs_get(xbi, 'depfile.depfile');

%%
xb = repmat(xbi, size(z,2), 1);
for islope = 1:n
    depfilea = interp1(x-xori, z(:,islope), xfile);
    depfilea(isnan(depfilea)) = min(depfilea);
    
    xb(islope) = xb_set(xb(islope), 'depfile.depfile', depfilea);
    
    basepath = fullfile(OPT.maindir, sprintf('%08i_%08i', id1(islope), id2(islope)), 'xbeach', sprintf('run%02i', OPT.irun));
    casename = sprintf('%s%.0f', OPT.zmin_extend, OPT.zmin);
    calcdir = fullfile(basepath, casename);
    mkpath(calcdir)
    xb_write_input(calcdir, xb(islope));
end