function [strandbreedte,duingrens,zeegrens] = ...
    bepaalStrandbreedte(maindir,datadir,raai,jaar,strandtype)

%% Script voor bepalen van strandbreedtes op basis van Jarkus-data
%
% Programmer: Vincent Vuik
% Company:    HKV lijn in water
% Project:    PR2063.20 - Kustlijnzorg
% Date:       19 mei 2011
%
% Stuur dit script aan vanuit runStrandbreedte.m
% Dit script bepaalt op basis van JARKUS-data de strandbreedte per raai en
% jaar. Er worden snijpunten bepaald met het duinvoetniveau en het GLW- of
% GHW-niveau. De snijpunten die de begrenzingen van het strand aan zee- en
% duinzijde representeren worden geselecteerd. Het verschil tussen deze
% begrenzingen is gelijk aan de strandbreedte.
%
% Wijzig eventueel:
% - locaties en instellingen
% - definities duinvoet, duinhoogte, zeediepte etc

%% Locaties

figdir     = [maindir,'Strandbreedte\figuren\strandbreedte_raaien\'];
locGHWGLW  = 'file';      % 'file' of 'OET'

%% Instellingen

plotfigures    = 'nee';   % ja/nee (bij nee alleen strandbreedtes wegschrijven, gaat snelst)
savefigures    = 'nee';   % ja/nee
closefigures   = 'nee';   % ja/nee
figure_visible = 'nee';   % ja/nee (nee gaat sneller, dan wel savefigures='ja', anders plotfigures zinloos)
printscreen    = 'ja';    % ja/nee (strandbreedtes per raai & jaar naar command window wegschrijven)
zoom           = 'str';   % 'all' (compleet profiel) of 'str' (alleen strand) --> keuze xlim, ylim

fig_xlim_all   = [-3500 3500];                                      % m vanaf nulpunt raai
fig_ylim_all   = [-15 20];                                          % m+NAP
fig_xtick_all  = 500;
fig_ytick_all  = 4;
fig_xlim_str   = [-500 500];                                        % m vanaf nulpunt raai
fig_ylim_str   = [-10 10];                                          % m+NAP
fig_xtick_str  = 100;
fig_ytick_str  = 2;

fig_xlim       = eval(['fig_xlim_',zoom]);
fig_ylim       = eval(['fig_ylim_',zoom]);
fig_xtick      = eval(['fig_xtick_',zoom]);
fig_ytick      = eval(['fig_ytick_',zoom]);

duinvoet       = 3;                                             % m+NAP, voor definitie duinvoet

% duinhoogte: opgeven minimale duinhoogte. wat beschouw je daadwerkelijk
% als duin (h>duinhoogte) en wat als slechts een hobbel cq oneffenheid
% in het strand (h<duinhoogte)?
duinhoogte     = duinvoet+1;                                    % m+NAP, default: duinhoogte = duinvoet + 1 m

if strcmp(locGHWGLW,'file')==1
    [GHW,GLW] = bepaalGHWGLW(maindir,raai);                     % aanroepen functie bepaalGHWGLW.m
    hw_lijn   = GHW;
    lw_lijn   = GLW;
elseif strcmp(locGHWGLW,'OET')==1
    check_oet;
    url       = jarkus_url;
    vakIDs    = nc_varget(url,'id');
    GHW       = nc_varget(url,'mean_high_water');
    GLW       = nc_varget(url,'mean_low_water');
    vakIDnr   = find(vakIDs==raai);
    hw_lijn   = GHW(vakIDnr);                                  % m+NAP, voor definitie gem. hoogwaterlijn (GHW)
    lw_lijn   = GLW(vakIDnr);                                  % m+NAP, voor definitie gem. laagwaterlijn (GLW)
end

% zeediepte: opgeven minimale zeediepte. wat beschouw je daadwerkelijk als
% zee (h<zeediepte) en wat slechts als een kuil cq oneffenheid in het
% strand (h>zeediepte)? definitie afhankelijk van strandtype droog/totaal
if strcmp(strandtype,'droog')==1
    zeediepte = (hw_lijn+lw_lijn)/2;                            % m+NAP, default: zeediepte = (ghw+glw)/2
elseif strcmp(strandtype,'totaal')==1
    zeediepte = lw_lijn-0.5;                                    % m+NAP, default: zeediepte = glw-0.5
end

%% Inlezen data

data          = load([datadir,num2str(jaar),'.txt']);
data_raai     = data(find(data(:,1)==raai),:);
if isempty(data_raai)
    strandbreedte = NaN;
    duingrens     = NaN;
    zeegrens      = NaN;
    return;
end
distance_hor  = flipdim(-data_raai(:,3),1);
distance_vert = flipdim( data_raai(:,4),1);

%% A) Bepalen alle snijpunten met duinvoetniveau

i_duinvoet = []; % maak lege i_duinvoet aan

% als er minimaal 1 snijpunt met duinvoet-niveau is:
if ~isempty(find(distance_vert>duinvoet))==1 || ~isempty(find(distance_vert<duinvoet))==1
    a=1; % teller
    for i=2:length(distance_vert)
        if distance_vert(i-1)>duinvoet && distance_vert(i)<=duinvoet
            i_duinvoet(a) = i; % locaties waar profiel duinvoethoogte snijdt
            a=a+1;
        elseif distance_vert(i-1)>=duinvoet && distance_vert(i)<duinvoet
            i_duinvoet(a) = i; % locaties waar profiel duinvoethoogte snijdt
            a=a+1;
        end
    end
end

%% B) Bepalen alle snijpunten met waterlijn
% Droog strand:  snijpunten met GHW
% Totaal strand: snijpunten met GLW

if strcmp(strandtype,'droog')==1
    waterlijn = hw_lijn;
elseif strcmp(strandtype,'totaal')==1
    waterlijn = lw_lijn;
end

i_waterlijn = []; % maak lege i_waterlijn aan

% als er minimaal 1 snijpunt met geldende waterlijn-niveau is:
if ~isempty(find(distance_vert>=waterlijn))==1 || ~isempty(find(distance_vert<=waterlijn))==1
    b=1; % teller
    for i=2:length(distance_vert)
        if distance_vert(i-1)>waterlijn && distance_vert(i)<=waterlijn
            i_waterlijn(b) = i; % locaties waar profiel duinvoethoogte snijdt
            b=b+1;
        elseif distance_vert(i-1)>=waterlijn && distance_vert(i)<waterlijn
            i_waterlijn(b) = i; % locaties waar profiel duinvoethoogte snijdt
            b=b+1;
        end
    end
end

%% C) Bepalen duingrens van het strand

if isempty(i_duinvoet)==1 % lege i_duinvoet: geen enkel snijpunt...
    if max(distance_vert)<=duinvoet % gehele profiel onder duinvoetniveau
        def_duinvoet = 1; % dan duingrens aan begin van profiel
        duingrens    = distance_hor(def_duinvoet);
    elseif min(distance_vert)>=duinvoet % gehele profiel boven duinvoetniveau
        def_duinvoet = length(distance_hor); % dan duingrens aan eind van profiel
        duingrens    = distance_hor(def_duinvoet);
    end
else % wel snijpunt(en)...
    a=1;
    def_duinvoet = i_duinvoet(a); % definitie duinvoet
    while ~isempty(find(distance_vert(def_duinvoet:end)>duinhoogte)) % als er nog bulten>duinhoogte zeewaarts van def_duinvoet liggen...
        a=a+1;
        if length(i_duinvoet)-a>=0 % als a=length(i_duinvoet) of kleiner...
            def_duinvoet = i_duinvoet(a); % def_duinvoet=i_duinvoet(1...length(i_duinvoet))
        else
            break;
        end
    end
    duingrens = (distance_vert(def_duinvoet-1)-duinvoet)/(distance_vert(def_duinvoet-1)-distance_vert(def_duinvoet))...
        *(distance_hor(def_duinvoet)-distance_hor(def_duinvoet-1))+distance_hor(def_duinvoet-1); % lineaire interpolatie
end

%% D) Bepalen zeegrens van het strand

if isempty(i_waterlijn)==1  % lege i_waterlijn: geen enkel snijpunt...
    if max(distance_vert)<=waterlijn % gehele profiel ligt onder waterlijn
        def_waterlijn = 1; % dan zeegrens aan begin van profiel
        zeegrens      = distance_hor(def_waterlijn);
    elseif min(distance_vert)>=waterlijn % gehele profiel ligt boven waterlijn
        def_waterlijn = length(distance_hor); % dan zeegrens aan eind van profiel
        zeegrens      = distance_hor(def_waterlijn);
    else
        strandbreedte = NaN;
        zeegrens      = NaN;
        return;
    end
else % wel snijpunt(en)...
    select_waterlijn = find(i_waterlijn>=def_duinvoet);
    i_waterlijn      = i_waterlijn(select_waterlijn);     % alleen die elementen die zeewaarts van duinvoet liggen. Voorkomt neg. strandbreedtes
    if isempty(i_waterlijn)==1  % na deze bewerking geen enkel snijpunt meer over...
        def_waterlijn = length(distance_hor); % dan zeegrens aan eind van profiel
        zeegrens      = distance_hor(def_waterlijn);
    else
        b=0;
        def_waterlijn    = i_waterlijn(end);
        if def_waterlijn==def_duinvoet
            def_waterlijn    = i_waterlijn(end);
        else
            while ~isempty(find(distance_vert(def_duinvoet:def_waterlijn)<zeediepte)) % als er nog kuilen<zeediepte landwaarts van def_waterlijn liggen...
                b=b+1;
                if length(i_waterlijn)-b>0 % als b<length(i_waterlijn)...
                    def_waterlijn = i_waterlijn(end-b); % end-b>0
                else
                    break;
                end
            end
        end
        zeegrens = (distance_vert(def_waterlijn-1)-waterlijn)/(distance_vert(def_waterlijn-1)-distance_vert(def_waterlijn))...
            *(distance_hor(def_waterlijn)-distance_hor(def_waterlijn-1))+distance_hor(def_waterlijn-1); % lineaire interpolatie
    end
end

strandbreedte = zeegrens-duingrens;
if strcmp(printscreen,'ja')==1
    disp(['Strandbreedte raai ',num2str(raai),', ',num2str(jaar),' = ',num2str(strandbreedte),' m']);
end

%% Plotten

if strcmp(plotfigures,'ja')==1

    if strcmp(figure_visible,'ja')==1
        figure(volgnr);
    elseif strcmp(figure_visible,'nee')==1
        figure('Visible','off');
    end

    % plot profiel en grenzen
    plot(distance_hor,distance_vert,'.-','Color','k','LineWidth',0.5);
    hold on;
    plot(duingrens,duinvoet,'o','Color',[0 0.42 0.24]);
    plot(zeegrens,waterlijn,'bo');
    hold off;

    % lijn voor weergeven waarde voor strandbreedte in figuur
    height_line = fig_ylim(2)-0.25*(fig_ylim(2)-fig_ylim(1));

    % voeg lijnen en labels in
    line(fig_xlim               ,[duinvoet duinvoet]      ,'Color',[0 0.42 0.24],'Linestyle','-.');
    text(fig_xlim(2)-fig_xtick/5,duinvoet+fig_ytick/5     ,'Duinvoet','HorizontalAlignment','left');
    line(fig_xlim               ,[hw_lijn hw_lijn]        ,'Color','b','Linestyle','-.');
    text(fig_xlim(2)-fig_xtick/5,hw_lijn+fig_ytick/5      ,'GHW','HorizontalAlignment','left');
    line(fig_xlim               ,[lw_lijn lw_lijn]        ,'Color','b','Linestyle','-.');
    text(fig_xlim(2)-fig_xtick/5,lw_lijn+fig_ytick/5      ,'GLW','HorizontalAlignment','left');
    text(fig_xlim(2)-fig_xtick/5,fig_ylim(2)-fig_ytick/2  ,['GHW: ',num2str(hw_lijn),' m+NAP'],'HorizontalAlignment','left');
    text(fig_xlim(2)-fig_xtick/5,fig_ylim(2)-fig_ytick    ,['GLW: ',num2str(lw_lijn),' m+NAP'],'HorizontalAlignment','left');
    line([zeegrens duingrens]   ,[height_line height_line],'Color','k');
    text((zeegrens)+fig_xtick/10,height_line+fig_ytick/5,...
        ['Strandbreedte = ',num2str(round(strandbreedte)),' m'],'HorizontalAlignment','right');
    line([zeegrens zeegrens]    ,fig_ylim                 ,'Color','r','LineStyle','-.');
    line([duingrens duingrens]  ,fig_ylim                 ,'Color','r','LineStyle','-.');

    % set axes properties
    xlim(fig_xlim);
    ylim(fig_ylim);
    set(gca,'box','on','XDir','reverse');
    set(gca,'XTick',[fig_xlim(1):fig_xtick:fig_xlim(2)],'XtickLabel',{-1*[fig_xlim(1):fig_xtick:fig_xlim(2)]});
    set(gca,'YTick',[fig_ylim(1):fig_ytick:fig_ylim(2)]);
    xlabel('Afstand vanaf nulpunt raai (m) \rightarrow ');
    ylabel('Bodemhoogte (m+NAP) \rightarrow');
    title(['Dwarsdoorsnede o.b.v. Jarkus-metingen en AHN - raai ',num2str(raai),...
        ', ',num2str(jaar)]);

    % save figures
    printname = [figdir,strandtype,'\strandbreedte_zoom_',zoom,'_raai_',n2svol(raai,8),'_',num2str(jaar),'.png'];
    if strcmp(savefigures,'ja')==1
        print(printname,'-dpng');
    end

    % close figures
    if strcmp(closefigures,'ja')==1
        close all;
    end

end