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

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2014 Deltares
%       Joost den Bieman
%
%       joost.denbieman@deltares.nl
%
%       P.O. Box 177
%       2600 MH 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: 19 Mar 2014
% Created with Matlab version: 8.2.0.701 (R2013b)

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

%% Settings
OPT = struct(...
    'JarkusID',     [],             ...
    'Years',        [],             ...
    'ExtrPeriod',   6,              ...
    'ExtrRFID',     [],             ...
    'ReturnFreqs',  [],             ...
    'LGW',          [],             ...
    'plotBKL',      false,           ...
    'plotMKL',      true,           ...
    'Extrapolate',  true,           ...         
    'NormShift',    0,              ...
    'plotTRDA2006', false,          ...
    'TRDAmarker',   'm^',     ...
    'FigName',      'R-t_diagram',  ...
    'MarkerSize',   50,             ...
    'PrintDir',     '',             ...
    'PrintEN',      false,          ...
    'Print',        false);

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

%% Initializing
Years       = OPT.Years;
ReturnFreqs = OPT.ReturnFreqs;
R           = NaN(numel(Years),numel(ReturnFreqs));
Converged   = false(numel(Years),numel(ReturnFreqs));
MKL         = NaN(numel(Years),1);
FreqColors  = ['c';'g';'m'];

%% Plotting
hFig    = figure;
hAx     = gca;
hold on

idxPeriod = 1:(numel(Years)-OPT.ExtrPeriod);

for iFreq = 1:numel(ReturnFreqs)
    for iYear = 1:numel(Years)%(numel(Years)-OPT.ExtrPeriod)
        for iCountDown = numel(DuneErosionResult(iYear,iFreq,:)):-1:1
            if ~isempty(DuneErosionResult{iYear,iFreq,iCountDown})
                R(iYear,iFreq)  = DuneErosionResult{iYear,iFreq,iCountDown}.settings.x2zVariables{12};
                Converged(iYear,iFreq)  = abs(DuneErosionResult{iYear,iFreq,iCountDown}.Output.P_f-ReturnFreqs(iFreq))<0.1*ReturnFreqs(iFreq);
                break
            end
        end
    end
    
    % Plot R point for every return frequency
    if OPT.PrintEN
        RPointName = ['R point 1/' num2str(1/ReturnFreqs(iFreq))];
    else
        RPointName = ['R punt 1/' num2str(1/ReturnFreqs(iFreq))];
    end
    
    if iFreq == OPT.ExtrRFID
        scatter(hAx,Years(Converged(:,iFreq)),R(Converged(:,iFreq),iFreq),OPT.MarkerSize,FreqColors(1+mod(iFreq,numel(FreqColors))),'Filled', ...
            'MarkerEdgeColor','k','DisplayName',RPointName)
        if sum(~Converged(:,iFreq))>0
            scatter(hAx,Years(~Converged(:,iFreq)),R(~Converged(:,iFreq),iFreq),OPT.MarkerSize,['h' FreqColors(1+mod(iFreq,numel(FreqColors)))],'Filled', ...
                'MarkerEdgeColor','k','DisplayName',[RPointName ' not converged'])
        end
    else
        scatter(hAx,Years(idxPeriod(Converged(idxPeriod,iFreq))),R(idxPeriod(Converged(idxPeriod,iFreq)),iFreq),OPT.MarkerSize,FreqColors(1+mod(iFreq,numel(FreqColors))),'Filled', ...
            'MarkerEdgeColor','k','DisplayName',RPointName)
        if sum(~Converged(:,iFreq))>0
            scatter(hAx,Years(idxPeriod(~Converged(idxPeriod,iFreq))),R(idxPeriod(~Converged(idxPeriod,iFreq)),iFreq),OPT.MarkerSize,['h' FreqColors(1+mod(iFreq,numel(FreqColors)))],'Filled', ...
                'MarkerEdgeColor','k','DisplayName',[RPointName ' not converged'])
        end
    end
end

PolyFit         = polyfitn(Years(idxPeriod(Converged(idxPeriod,iFreq))),R(idxPeriod(Converged(idxPeriod,iFreq)),OPT.ExtrRFID),1);
LinearFit       = polyvaln(PolyFit,Years);

% Include MKL point per year
if OPT.plotMKL
    scatter(hAx,Years,MKL,OPT.MarkerSize,'rs','Filled','MarkerEdgeColor','k','DisplayName','MKL')
end

xlim(hAx,[floor(min(get(hAx,'XLim'))-1) ceil(max(get(hAx,'XLim'))+1)])
if OPT.Extrapolate
    SDFit           = std(R(idxPeriod(Converged(idxPeriod,OPT.ExtrRFID)),OPT.ExtrRFID)-LinearFit(idxPeriod(Converged(idxPeriod,OPT.ExtrRFID))));
    plot(hAx,Years,LinearFit,'k--','LineWidth',1.5,'DisplayName',['Linear fit and extrapolation, \sigma = ' num2str(SDFit,'%.1f') ' m'])
    % scatter(hAx, Years(end), LinearFit(end), 'ko', 'filled', 'DisplayName', 'Prediction')
    hEB =   errorbar(hAx, Years(end), LinearFit(end), SDFit, 'ko', 'MarkerFaceColor', 'k', 'DisplayName', ['Prediction, ' num2str(LinearFit(end),'%.1f') ' m']);
    errorbar_tick(hEB, 50)
    % plot(hAx, [Years(end) Years(end)], [LinearFit(end) LinearFit(end)-SDFit], '-ko', 'MarkerFaceColor', 'k', 'DisplayName', 'Prediction & limit')
else
    Rselection  = R(idxPeriod(Converged(idxPeriod,OPT.ExtrRFID)),OPT.ExtrRFID);
    Rranked     = sort(Rselection, 1, 'ascend');
    line(get(hAx,'XLim'),[Rranked(3)+OPT.NormShift Rranked(3)+OPT.NormShift],'Color','k','LineStyle','--','LineWidth',1.5,'DisplayName',['Normative point, ' num2str(Rranked(3)+OPT.NormShift,'%.1f') ' m'])
end


% Plot critical retreat line
if ~isempty(OPT.LGW)
    if OPT.PrintEN
        LGWName = 'Landward boundary flood defense';
    else
        LGWName = 'Landwaartse grens waterkering';
    end
    line(get(hAx,'XLim'),[OPT.LGW OPT.LGW],'Color','b','LineStyle','--','LineWidth',2,'DisplayName',LGWName)
end

% Include BKL line
if OPT.plotBKL
    BKL = jarkus_getBKL(OPT.JarkusID);
    line(get(hAx,'XLim'),[BKL BKL],'Color','b','LineStyle','-','LineWidth',2,'DisplayName','BKL')
end

% Include TRDA2006 results
if OPT.plotTRDA2006
    [R06, y06] = calc_Rt_trda2006('id',OPT.JarkusID);
    if OPT.PrintEN
        TRDAName = 'TRDA 2006 results';
        TRDAresultname = 'TRDA 2006 normative point';
    else
        TRDAName = 'TRDA 2006 resultaten';
        TRDAresultname = 'TRDA 2006 maatgevend afslagpunt';
    end

    yearFilter = false(size(y06));
    for iYr = 1:numel(Years)
        yearFilter = (y06==Years(iYr)) + yearFilter;
    end
    yearFilter = logical(yearFilter);
    
%     idStart         = find(y06 == Years(idxPeriod(1)));
    
    y06selection = y06(yearFilter);
    R06selection = R06(yearFilter);
    filter      = ~isnan(R06selection) ;
    y06filtered = y06selection(filter);
    R06filtered = R06selection(filter);
    
    idStop          = find(y06filtered == Years(idxPeriod(end)));
    
    scatter(hAx, y06filtered, R06filtered, OPT.TRDAmarker, 'filled', 'MarkerEdgeColor', 'k','DisplayName',TRDAName)
    
    PolyFit06       = polyfitn(y06filtered(1:idStop),R06filtered(1:idStop),1);
    LinearFit06     = polyvaln(PolyFit06,Years);
    SDFit06         = std(R06filtered(1:idStop)'-LinearFit06(1:(end-OPT.ExtrPeriod)));
    
    if OPT.Extrapolate
        plot(hAx,Years,LinearFit06,'k:','LineWidth',1.5,'DisplayName',['Fit ' TRDAName ', \sigma = ' num2str(SDFit06,'%.1f') ' m'])
        
        %     scatter(hAx, Years(end), LinearFit06(end), 'k^', 'filled', 'DisplayName', 'Prediction (TRDA 2006)')
        %     plot(hAx, [Years(end) Years(end)], [LinearFit06(end) LinearFit06(end)-SDFit06], '-k^', 'MarkerFaceColor', 'k', 'DisplayName', 'Prediction & limit (TRDA 2006)')
        hEB06           = errorbar(hAx, Years(end), LinearFit06(end), SDFit06, 'r^', 'MarkerFaceColor', 'r', 'DisplayName', ['Prediction (TRDA 2006), ' num2str(LinearFit06(end),'%.1f') ' m']);
        errorbar_tick(hEB06, 50)
        
        %     plot(hAx, [Years((end-OPT.ExtrPeriod):end)' Years((end-OPT.ExtrPeriod):end)'], [LinearFit06((end-OPT.ExtrPeriod):end)-2*SDFit06 LinearFit06((end-OPT.ExtrPeriod):end)+2*SDFit06], 'r', 'LineWidth', 1, 'DisplayName', ['Prediction (TRDA 2006), ' num2str(LinearFit06(end),'%.1f') ' m']);
        
        % search for 3rd worst (op-2-na-slechtste) punt to compare with
        % TRDA2006
    end
    
    R06ranked   = sort(R06filtered(1:idStop), 2, 'ascend');
    line(get(hAx,'XLim'),[R06ranked(3) R06ranked(3)],'Color','k','LineStyle','-.','LineWidth',1.5,'DisplayName',[TRDAresultname ', ' num2str(R06ranked(3),'%.1f') ' m'])
end

legend(hAx,'show','location','SouthOutside')

ylims           = get(hAx,'YLim');
Ymin            = min(ylims);
Ymax            = max(ylims);
Xmin            = Years(numel(Years)-OPT.ExtrPeriod)+0.5;
Xmax            = ceil(max(get(hAx,'XLim'))+1);
patch([Xmin Xmax Xmax Xmin],[Ymax Ymax Ymin Ymin],[0.2 0.2 0.2],'FaceAlpha',0.2,'EdgeColor','none','DisplayName','none');

hold off

if OPT.PrintEN
    xlabel('Year')
    ylabel('Horizontal distance w.r.t. RSP [m]')
    title(['R-t diagram for transect ' num2str(OPT.JarkusID)])
else
    xlabel('Jaar')
    ylabel('Horizontale afstand t.o.v. RSP [m]')
    title(['R-t diagram voor raai ' num2str(OPT.JarkusID)])
end



if OPT.Print
    print('-dpng','-r600',fullfile(OPT.PrintDir,OPT.FigName))
    print('-dpdf',fullfile(OPT.PrintDir,OPT.FigName))
end