 function MCresult = dutch_dune_coast_MC(varargin)
%DUTCH_DUNE_COAST_MC  Monte Carlo dune erosion simulation
%
%   More detailed description goes here.
%
%   Syntax:
%   MCresult = dutch_dune_coast_MC(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   MCresult =
%
%   Example
%   dutch_dune_coast_MC
%
%   See also 

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2011 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: 30 Mar 2011
% Created with Matlab version: 7.7.0.471 (R2008b)

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

%% default input

NC_ID = double(nc_varget('http://opendap.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc','id'));

OPT = struct(...
    'jarkus_url', 'http://opendap.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc',...
    'HBC_url', fullfile(fileparts(which(mfilename)), 'HBC_VTV2006prob.nc'),...
    'grainsize_url', 'http://opendap.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/grainsize_vtv/grainsize.nc',...
    'runbasedir', fullfile(fileparts(which(mfilename)), '..', 'runs'),...
    'errortablefile', 'messages.sqlite',...
    'NrSamples', 100,...
    'id', NC_ID(849:1441,1),...
    'year', 2010,...
    'zRef', 5,...
    'zRef1',3,...
    'zRef2',-4,...
    'ncfile', 'MCresult.nc');
    
OPT = setproperty(OPT, varargin{:});

global zVOL 
global zERO 

%% obtain transects
tr = jarkus_transects(...
    'id', OPT.id,...
    'year', OPT.year);
tr.altitude = squeeze(tr.altitude);

for i = 1:length(OPT.id)
    QQ(i) = isequal(sum(isnan(tr.altitude(i,:))),1925);
end
aantal_transects_without_measurement_for_the_given_year = sum(QQ)

% find id's of nan-transects
ZZ = [tr.id' QQ']; 
nan_transects = ZZ(find(ZZ(:,2)==1),1);

if ~isequal(OPT.id(:), tr.id(:))
    warning('Not all requested transects have been found.')
end

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

%% boundary conditions
station_names = cellstr(nc_varget(OPT.HBC_url, 'name'));
x_station = nc_varget(OPT.HBC_url, 'x');
y_station = nc_varget(OPT.HBC_url, 'y');
[stationids lambda id angle_within_range transects_facing_seaward transects_in_range] = jarkus_derive_HB_stations(OPT.id, x_station, y_station,...
    'virtual', [false false false false true false]);

%%
selected_transects = angle_within_range & transects_facing_seaward & transects_in_range;
[stationids lambda id] = deal(stationids(selected_transects,:), lambda(selected_transects), id(selected_transects));

%% create runfolder
basedir = OPT.runbasedir;
runfolders = dir(fullfile(basedir, 'run*'));
runfolder = fullfile(basedir, sprintf('run%03i', length(runfolders)));
mkdir(runfolder)

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

% Open database
mksqlite('open', database);
mksqlite(['create table ' table ' (profile int, sample int, errorcode int, message char(255))']);
mksqlite('close', database);

%% Monte Carlo
P_f = NaN(size(id));
for i = 1:length(P_f)
    
    transectid = id(i);
    % jarkusid is the index of the jarkus file that corresponds to i of
    % this loop
    jarkusid = id(i) == tr.id;
    
    % define stochastic variables
    stochast = createStochastVar_2stp(...
        'station', station_names(stationids(i,:))',...
        'lambda', lambda(i),...
        'D50', [meanD50(jarkusid) sigmaD50(jarkusid)]);
    zi = tr.altitude(jarkusid,:);
            
    % id of non-NaN
    nonnanid = ~isnan(zi);
    if any(nonnanid)
        [xi zi] = deal(tr.cross_shore(nonnanid), zi(nonnanid));
        if OPT.zRef  > max(zi)
            continue
        end
        
        % derive cross-shore reference location
        xcr = findCrossings(xi, zi, [min(xi) max(xi)], ones(1,2)*OPT.zRef);
        xRef = max(xcr);
        
        %characterize profile: find volume and slope between both z-ref's
        xcr1 = findCrossings(xi, zi, [min(xi) max(xi)], ones(1,2)*OPT.zRef1);
        xRef1 = max(xcr1);      
        
        xcr2 = findCrossings(xi, zi, [min(xi) max(xi)], ones(1,2)*OPT.zRef2);
        xRef2 = min(xcr2);
       
        VolA = getVolume(...
            'x',xi,...
            'z',zi,...
            'UpperBoundary',OPT.zRef1,...
            'LowerBoundary',OPT.zRef2,...
            'LandwardBoundary',xRef1,...
            'SeawardBoundary',xRef2);
        
        % Compute the slope of the 'foreshore' as defined above.
        if transectid == 8005625
            Slope = 0;
        else
            if transectid == 8005650
                Slope = 0;
            else
                Slope = (OPT.zRef1 - OPT.zRef2) / (xRef2 - xRef1);
            end
        end
        
          
        % run Monte Carlo
        MCresult = MC(...
            'stochast', stochast,...
            'x2zFunction', @x2z_duros,...
            'NrSamples', OPT.NrSamples,...
            'variables',    {'xInitial' xi,...
                             'zInitial' zi,...
                             'xRef' xRef,...
                             'transectid', transectid,...
                             'database' database,...
                             'table' 'errortable'},...
            'ISvariable', 'WL_t',...
            'f1', 1,...
            'f2', 1e-7);
        P_f(i) = MCresult.Output.P_f;
         
        % importance sampling settings
        %             'isvariable', 'wl_t',...
        %             'f1', 1,...
        %             'f2', 1e-7,...
        
        % write results to netcdf file
        
        prob_MCresult2netcdf(MCresult,...
            'ncfile', fullfile(runfolder, OPT.ncfile),...
            'run_id', transectid,...
            'volume', VolA,...
            'slope', Slope,...
            'zVOL', zVOL,...
            'zERO', zERO);
    end
end