function varargout = dunebreach_probability(varargin)
%DUNEBREACH_PROBABILITY  derive the breaching probability of dunes
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = dunebreach_probability(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   dunebreach_probability
%
%   See also 

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2010 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 Lesser General Public
%   License as published by the Free Software Foundation, either
%   version 2.1 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
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library. If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://OpenEarth.nl">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: 26 Oct 2010
% Created with Matlab version: 7.10.0.499 (R2010a)

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

%%
OPT = struct(...
    'jarkus_url', 'http://opendap.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc',...
    'HBC_url', 'HBC_VTV2006prob.nc',...
    'grainsize_url', 'http://opendap.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/grainsize_vtv/grainsize.nc',...
    'landboundary_url', 'http://opendap.deltares.nl:8080/thredds/dodsC/opendap/deltares/landboundaries/holland.nc',...
    'ncpath', fileparts(which(mfilename)),...
    'verticalextend', [-5 5],... % minimum requirement for inclusion of transect
    'errortablefile', [mfilename '.sqlite'],...
    'runbasedir', fullfile(fileparts(which(mfilename)), '..', 'runs'),...
    'NrSamples', 1000,...
    'year', 2008);

try
    OPT = setproperty(OPT, varargin{:});
end

ncfile = fullfile(OPT.ncpath, [mfilename '_' num2str(OPT.year) '.nc']);

%%
% retrieve netcdf info
info = nc_info(OPT.jarkus_url);
% filter variable names
varnames = {info.Dataset.Name};

% retrieve transect id's
id = nc_varget(OPT.jarkus_url, 'id');

% find variable applicable variable name for 'time' (depends on datafile
% (old/new))
time = varnames{ismember(varnames, {'year' 'time'})};
altitude = varnames{ismember(varnames, {'height' 'altitude'})};

% define start and count for dimension time
timeid = round(nc_varget(OPT.jarkus_url, time)/365.24)+1970 == OPT.year;
[start1 count1] = deal(find(timeid)-1, 1);

% get altitude
z = nc_varget(OPT.jarkus_url, altitude, [start1 0 0], [count1 -1 -1]);
fillvalue = nc_attget(OPT.jarkus_url, altitude, '_FillValue');
z(z == fillvalue) = NaN;

% get cross-shore coordinates
cross_shore = varnames{ismember(varnames, {'seaward_distance' 'cross_shore'})};
x = nc_varget(OPT.jarkus_url, cross_shore);

if ~exist(ncfile, 'file')
    
    % get kustvak codes
    areacode = nc_varget(OPT.jarkus_url, 'areacode');

    % get alongshore coordinates
    y = nc_varget(OPT.jarkus_url, varnames{ismember(varnames, {'coastward_distance' 'alongshore'})});

    % get angle
    Angle = nc_varget(OPT.jarkus_url, 'angle');

    xRD = nc_varget(OPT.jarkus_url, 'x');
    xRD0 = xRD(:,1);
    yRD = nc_varget(OPT.jarkus_url, 'y');
    yRD0 = yRD(:,1);
    xRD(isnan(z)) = NaN;
    yRD(isnan(z)) = NaN;

    %% find crossings of extended transects with connection lines between
    %% reference points (steunpunten)
    xstp = nc_varget(OPT.HBC_url, 'x');
    ystp = nc_varget(OPT.HBC_url, 'y');
    xr = [min(xstp) max(xstp)];
    [a b xcr ycr] = deal(nan(size(xRD0)));
    for i = 1:length(xRD0)
        [a(i) b(i)] = xydegN2ab(xRD0(i), yRD0(i), Angle(i));
        try
            if ismember(Angle(i), [0 180])
                % north or south directed transect, no y = ax+b approach possible
                xcr(i) = xRD0(i);
                ycr(i) = interp1(xstp, ystp, xcr(i));
            else
                if any(isnan(polyval([a(i) b(i)], xr)))
                    dbstopcurrent
                end
                [xcr(i) ycr(i)] = findCrossings(xr, polyval([a(i) b(i)], xr), xstp, ystp);
            end
        catch
            [xcr(i) ycr(i)] = deal(NaN);
        end
    end

    %% select only transects that face seaward
    % derive angles of lines between steunpunten
    % positive clockwise 0 north
    angle_stp = xy2degN(xstp(1:end-1), ystp(1:end-1), xstp(2:end), ystp(2:end));
    stp1id = NaN(size(xcr));
    stp_angle = NaN(size(xcr));
    for i = 1:length(xcr)
        if ~isnan(xcr(i))
            try
                stp1id(i) = find(xcr(i) - xstp > 0, 1, 'last');
                stp_angle(i) = angle_stp(stp1id(i));
            catch
                stp1id(i) = NaN;
            end
        end
    end
    stp2id = stp1id + 1; % other steunpunt is the neighbouring one

    angle_tr = xy2degN(xRD0, yRD0, xcr', ycr');

    angle_range = 21;
    angleid = Angle > 270-angle_range + stp_angle' & Angle < 270+angle_range + stp_angle'  | ...
        Angle < -90+angle_range + stp_angle' & Angle < -90-angle_range + stp_angle';

    %% exclude areas with dikes/hard structures
    % Pettemer Sea Defence
    selid = id < 7e6 + 2041 | id > 7e6 + 2582;
    % Dike near head of Den Helder
    selid = selid & (id < 7e6 | id >= 7e6 + 110);

    selid = (selid & angleid)';
    
    %% check for vertical extend criteria
    [minz maxz] = deal(NaN(size(selid)));
    for transect = find(selid)
        nonnanid = ~isnan(z(transect,:));
        if ~any(nonnanid)
            selid(transect) = any(nonnanid);
        else
            maxz(transect) = max(z(transect,nonnanid));
            minz(transect) = min(z(transect,nonnanid));
            selid(transect) = ...
                maxz(transect) > max(OPT.verticalextend) &&...
                minz(transect) < min(OPT.verticalextend);
        end
        
        if selid(transect)
            lwb(transect) = jarkus_identify_first_dunerow(...
                x(nonnanid),...
                z(transect,nonnanid)');
        end
    end
    
    %% plot selected transects
    %{
    plot(nc_varget(OPT.landboundary_url, 'x'), nc_varget(OPT.landboundary_url, 'y'), 'k',...
        xstp, ystp, 'o')
    hold on
    rgb = repmat({'r' 'g' 'b'}, 1, 6);
    for i = find(P_f > 1e-5)'%find(selid)%1:size(xRD,1)
    %     plot([xRD0(i) xcr(i)], [yRD0(i) ycr(i)], 'color', rgb{areacode(i)});
        nonnanid = ~isnan(z(i,:));
        plot(xRD(i,nonnanid), yRD(i,nonnanid), 'color', rgb{areacode(i)});
        plot([xRD(i,find(nonnanid, 1, 'last')) xcr(i)], [yRD(i,find(nonnanid, 1, 'last')) ycr(i)],...
            'LineStyle', ':',...
            'color', rgb{areacode(i)});
    end
    %}

    %% derive lambda
    [Angle stp1id stp2id xcr ycr selectedid lwb] = deal(Angle(selid), stp1id(selid), stp2id(selid), xcr(selid), ycr(selid), id(selid), lwb(selid));
    Lambda = NaN(size(Angle));
    for i = 1:length(Lambda)
        idinterp = stp1id(i)+[0 1];
        if stp2id(i) == 5
            rangelambda = [1 .57];
            stp2id(i) = stp2id(i) + 1;
        elseif stp1id(i) == 5
            rangelambda = [.57 0];
            stp1id(i) = stp1id(i) - 1;
        else
            rangelambda = [1 0];
        end
        Lambda(i) = interp1(xstp(idinterp), rangelambda, xcr(i));
    end

    %% find grainsizes for selected transects
    [meanD50 sigmaD50] = rws_grainsize2jarkus(...
        'output', {'meanD50' 'sigmaD50'},...
        'jarkus_url', OPT.jarkus_url,...
        'grainsize_url', OPT.grainsize_url);
    [meanD50 sigmaD50] = deal(meanD50(selid), sigmaD50(selid));

    %% convert x y coordinates to lat lon
    [lon lat] = convertCoordinates(xcr, ycr, ...
        'CS1.name', 'Amersfoort / RD New',...
        'CS2.code', 4326);

    %% write netCDF file with boundary locations and lambdas
    nc_create_empty(ncfile)

    % add global atributes
    nc_attput(ncfile, nc_global, 'title', 'VTV2006 probabilistic Hydraulic Boundary conditions interpolation information');
    % nc_attput(ncfile, nc_global, 'institution', 'Rijkswaterstaat');
    nc_attput(ncfile, nc_global, 'history', ['converted to netCDF on ' date ' with script with $Id$']);
    % nc_attput(ncfile, nc_global, 'references', 'Walstra, D. J. R., Diermanse, F. L. M., van Geer, P. F. C., van Thiel de Vries, J., Groeneweg, J., and van de Graaff, J. (2008). SBWDuinen2-Ontwikkeltraject. Technical report, Deltares.');
    nc_attput(ncfile, nc_global, 'year', OPT.year')
    % add dimensions
    nc_add_dimension(ncfile, 'alongshore', length(xcr))

    s = struct(...
        'Name', {
            'x'
            'y'
            'lat'
            'lon'
            'referencepoint1'
            'referencepoint2'
            'lambda1'
            'muD50'
            'sigmaD50'
            'landwardboundary'
            'selectedid'
            },...
        'Nctype', {
            nc_double
            nc_double
            nc_double
            nc_double
            nc_double
            nc_double
            nc_double
            nc_double
            nc_double
            nc_double
            nc_int
            },...
        'Dimension', {
            {'alongshore'}
            {'alongshore'}
            {'alongshore'}
            {'alongshore'}
            {'alongshore'}
            {'alongshore'}
            {'alongshore'}
            {'alongshore'}
            {'alongshore'}
            {'alongshore'}
            {'alongshore'}
            });

    %% add variables
    for i = 1:length(s)
        nc_addvar(ncfile, s(i));
    end

    nc_varput(ncfile, 'x', xcr);
    nc_varput(ncfile, 'y', ycr);
    nc_varput(ncfile, 'lat', lat);
    nc_varput(ncfile, 'lon', lon);
    nc_varput(ncfile, 'referencepoint1', stp1id);
    nc_varput(ncfile, 'referencepoint2', stp2id);
    nc_varput(ncfile, 'lambda1', Lambda);
    nc_varput(ncfile, 'muD50', meanD50);
    nc_varput(ncfile, 'sigmaD50', sigmaD50);
    nc_varput(ncfile, 'landwardboundary', lwb);
    nc_varput(ncfile, 'selectedid', selectedid);
else
    selectedid = nc_varget(ncfile, 'selectedid');
    Lambda = nc_varget(ncfile, 'lambda1');
    stp1id = nc_varget(ncfile, 'referencepoint1');
    stp2id = nc_varget(ncfile, 'referencepoint2');
    meanD50 = nc_varget(ncfile, 'muD50');
    sigmaD50 = nc_varget(ncfile, 'sigmaD50');
    lwb = nc_varget(ncfile, 'landwardboundary');
end

%% check of last run is complete
% if so: create a new one
% else: continue run
basedir = OPT.runbasedir;
runfolders = dir(fullfile(basedir, 'run*'));
if ~isempty(runfolders)
    runfolder = fullfile(basedir, runfolders(end).name);
else
    runfolder = '';
end
ncfile = fullfile(runfolder, ['year_' num2str(OPT.year) '.nc']);
continuerun = false;
if ~isempty(runfolders) && exist(ncfile, 'file')
    run_id = nc_varget(ncfile, 'run_id');
    continuerun = run_id(end) < selectedid(end);
    if continuerun
        % remove transects already done from id, in order to continue process
        selectedid = selectedid(selectedid ~= run_id(1:end-1));
        % set run_id of last run to 0, to prevent repeated run_id's
        run_id_cut = run_id;
        run_id_cut(run_id_cut > run_id(end-1)) = 0;
        nc_varput(ncfile, 'run_id', run_id_cut)
    end
end
if ~continuerun
    runfolder = fullfile(basedir, sprintf('run%03i', length(runfolders)));
    mkdir(runfolder)
    ncfile = fullfile(runfolder, ['year_' num2str(OPT.year) '.nc']);
end

%% prepare sqlite database to store all DUROS messages
database = fullfile(runfolder, OPT.errortablefile);
table = 'errortable';

if exist(database, 'file')
    mksqlite('open', database);
    mksqlite(['DELETE FROM ' table ' WHERE profile > ' num2str(run_id(end-1))]);
    mksqlite('close', database);
else
    % Open database
    mksqlite('open', database);
    mksqlite(['create table ' table ' (profile int, sample int, errorcode int, message char(255))']);
    mksqlite('close', database);
end

%% Monte Carlo
Locations = {
    'Hoek van Holland'...
    'IJmuiden'...
    'Den Helder'...
    'Eierlandse Gat'...
    'Steunpunt Waddenzee'...
    'Borkum'
    };

P_f = zeros(size(Lambda));

% loop over all selected transects
for i = 1:length(selectedid)
    transectid = selectedid(i);
    % define stochastic variables
    stochast = createStochastVar_2stp(...
        'station', [Locations(stp1id(selectedid == transectid)) Locations(stp2id(selectedid == transectid))],...
        'lambda', Lambda(selectedid == transectid),...
        'D50', [meanD50(selectedid == transectid) sigmaD50(selectedid == transectid)]);
    zi = z(id == transectid,:);
    % id of non-NaN and seaward of landwardboundary
    nonnanid = ~isnan(zi) & x' >= lwb(selectedid == transectid);
    if any(nonnanid)
        [xi zi] = deal(x(nonnanid), zi(nonnanid)');
        
        % run Monte Carlo
        MCresult = MC(stochast,...
            'x2zFunction', @x2zBreach,...
            'NrSamples', OPT.NrSamples,...
            'ISvariable', 'WL_t',...
            'f1', 1,...
            'f2', 1e-10,...
            'variables', {'xInitial' xi 'zInitial' zi 'database' database 'transectid' transectid 'table' table});
        P_f(i) = MCresult.Output.P_f;
        % write results to netcdf file
        prob_MCresult2netcdf(MCresult,...
            'ncfile', ncfile,...
            'run_id', transectid);
    end
end