function AnalyseXBeachOutput_Hydro(info,dirs)
%% Notes
% Hydraulische analyse output

%% ------------------------------------------------------------------ %%
%  ---------------------------- Settings ---------------------------- %
%  ------------------------------------------------------------------ %

% XBeach model settings
runnameXB       = info.run;
refdate         = datetime(2015,01,10,07,00,00); % Startmoment XBeach som (na inspeeltijd)
tstart          = 1200;                          % [s] Start simulatie (/output) na inspeeltijd
tstop           = 133200;  % 4899600 voor run09  % [s] Eind simulatie 
dtXB_pnt        = 0.1;                           % [s] Puntuitvoerinterval 
dtXB_mean       = 1800;                          % [s] Mean output interval
NrOut           = round((tstop-tstart)./dtXB_pnt)+1; % Aantal tijdstappen uitvoer (points)
tpeakstart      = datetime(2015,01,11,00,30,00); % Tijd start stormpiek voor in dwarsdoorsnede
tpeakend        = datetime(2015,01,11,02,00,00); % Tijd eind stormpiek voor in dwarsdoorsnede

% Tijdinstellingen meetdata
tMeasStart      = datetime(2014,11,01,00,00,00);
tMeasInt        = 15*60;                        % [s] (interval van 15 min)
tMeasEnd        = datetime(2015,01,31,17,00,00);

HFreq           = [0.05 1];
LFreq           = [0.005 0.05];                 % Dus zonder very low frequency waves, net als in gemeten data

%% ------------------------------------------------------------------ %%
%  ------------------------- Preprocessing -------------------------- %
%  ------------------------------------------------------------------ %

% Time
endref          = refdate + seconds(tstop);                 % Eind simulatie
ts              = 1:tMeasInt/dtXB_pnt:NrOut;                % Indices in output van start van bin gelijk aan meetbins voor simulatie (zonder inspeeltijd) t.o.v. start wegschrijven output.
ts2             = [ts NrOut];                               % Idem, maar nu incl. index eind laatste bin
tsReal          = (ts-1).*dtXB_pnt;                         % [s] Tijdserie van modeltijd per bin

% Paths     
fname           = [dirs.output filesep 'xboutput.nc'];

%% ------------------------------------------------------------------ %%
%  --------------------------- Load data ---------------------------- %
%  ------------------------------------------------------------------ %
% Load XBeach grid
x = load([dirs.output filesep 'x.grd']);

% Load XBeach output
pointtime   = nc_varget(fname,'pointtime');
point_zs    = nc_varget(fname,'point_zs');
point_zb    = nc_varget(fname,'point_zb');
point_H     = nc_varget(fname,'point_H');
pointx      = nc_varget(fname,'pointx');
pointy      = nc_varget(fname,'pointy');
meantime    = nc_varget(fname,'meantime');
H_mean      = squeeze(nc_varget(fname,'H_mean'));
H_min       = squeeze(nc_varget(fname,'H_min'));
H_max       = squeeze(nc_varget(fname,'H_max'));
zs1_var     = squeeze(nc_varget(fname,'zs1_var'));
zs_min      = squeeze(nc_varget(fname,'zs_min'));
zs_max      = squeeze(nc_varget(fname,'zs_max'));
zs_mean     = squeeze(nc_varget(fname,'zs_mean'));
zb          = squeeze(nc_varget(fname,'zb_mean'));  zb = zb(1,:);   % Bed level is alle tijdstappen hetzelfde, dus maar 1 tijdstap nodig

% Load measurement data
load([dirs.data filesep 'obs_hydro.mat']);

%% ------------------------------------------------------------------ %%
%  ---------------- Processing/Analyse model output ----------------- %
%  ------------------------------------------------------------------ %

%% Timeseries analysis measurement locations
for ii = 1:size(point_zs,2)   % Loop over observation points
    
    %% Determine wave characteristics
    for jj = 1:length(ts)    % Loop over timebins (1 minder dan timesteps)
        mean_zb(ii,jj)  = mean(point_zb(ts2(jj):ts2(jj+1),ii));
        mean_wl(ii,jj)  = mean(point_zs(ts2(jj):ts2(jj+1),ii));
        
        % Wave height high frequency waves
        Hrms_hf(ii,jj)  = rms(point_H(ts2(jj):ts2(jj+1),ii));   
        Hm0_hf(ii,jj)   = sqrt(2).*Hrms_hf(ii,jj);                           % Alleen geldig in diep water?

        % Wave height low frequency/infragravity waves
        zstmp           = point_zs(ts2(jj):ts2(jj+1),ii);
        %    Methode 2: O.b.v. energiespectrum
        zs_IG           = lfhf2(detrend(zstmp), LFreq(1), LFreq(2), dtXB_pnt);     % Split surface elevation time series into a high/low frequency signal after fast fourier transform of the detrended surface elevation signal
        [freqS,PwrS]    = spectrumsimple(1:length(zs_IG), zs_IG, 'figures',0); % Generates simple spectrum of y(x) with hann filter and welch method
        Hm0_IG(ii,jj)   = spec2hs(freqS,PwrS);                                 % Calculate significant wave height from wave spectrum without high frequency tail
    end
    
    %% Determine goodness of fit measured-modeled for storm period
    % Bepaal indices in meetdata binnen storm
    idStart     = find(time.(obs_points{ii}) <= refdate, 1, 'last');        % index start storm in meetdata
        if isempty(idStart); idStart = 1; end
    idEnd       = find(time.(obs_points{ii}) >= endref, 1, 'first');        % index eind storm in meetdata
        if isempty(idEnd); idEnd = length(time.(obs_points{ii})); end
    id{ii}      = idStart:idEnd;                                            % alle indices tijdens storm in meetdata
    id{ii}      = id{ii}(~isnan(data.(obs_points{ii})(idStart:idEnd,7)));   % verwijder indices van tijdstippen met NaNs in data (o.b.v. Hm0 LF)  
        
    % Te vergelijken data (d = measured data, m = model data)
    d_time          = datenum(time.(obs_points{ii})(id{ii})); 
    d_datahh{ii}    = data.(obs_points{ii})(id{ii},1);      % water depth
    d_dataHhf{ii}   = data.(obs_points{ii})(id{ii},3);      % Hm0_hf
    d_dataHig{ii}   = data.(obs_points{ii})(id{ii},7);      % Hm0_lf (zonder very low frequency waves)
    
    m_timetmp       = datenum(refdate+seconds(tsReal));
    [~, id3{ii}]    = ismember(d_time, m_timetmp);
    m_time          = m_timetmp(id3{ii})';                       
    m_datahh{ii}    = (mean_wl(ii,id3{ii}) - mean_zb(ii,id3{ii}))'; 
    m_dataHhf{ii}   = Hm0_hf(ii,id3{ii})';
    m_dataHig{ii}   = Hm0_IG(ii,id3{ii})';

    % Bepaal goodness of fit voor waterdepth, Hm0 hf en Hm0 ig
    [GoF_WL.r2(ii),  GoF_WL.sci(ii),  GoF_WL.relbias(ii),  GoF_WL.bss(ii),  GoF_WL.rms(ii),  GoF_WL.bias(ii)]   = xb_skill_BOI([d_time,d_datahh{ii} ], [m_time,m_datahh{ii} ]);
    [GoF_Hhf.r2(ii), GoF_Hhf.sci(ii), GoF_Hhf.relbias(ii), GoF_Hhf.bss(ii), GoF_Hhf.rms(ii), GoF_Hhf.bias(ii)]  = xb_skill_BOI([d_time,d_dataHhf{ii}], [m_time,m_dataHhf{ii}]);
    [GoF_Hig.r2(ii), GoF_Hig.sci(ii), GoF_Hig.relbias(ii), GoF_Hig.bss(ii), GoF_Hig.rms(ii), GoF_Hig.bias(ii)]  = xb_skill_BOI([d_time,d_dataHig{ii}], [m_time,m_dataHig{ii}]);

    clear idStart idEnd d_time m_timetmp  
end

% GoF all data
obs_hh_all  = cell2mat(d_datahh');
obs_Hhf_all = cell2mat(d_dataHhf');
obs_Hig_all = cell2mat(d_dataHig');
mod_hh_all  = cell2mat(m_datahh');
mod_Hhf_all = cell2mat(m_dataHhf');
mod_Hig_all = cell2mat(m_dataHig');
[GoF_WL.r2(size(point_zs,2)+1),  GoF_WL.sci(size(point_zs,2)+1),  GoF_WL.relbias(size(point_zs,2)+1),  GoF_WL.bss(size(point_zs,2)+1),  GoF_WL.rms(size(point_zs,2)+1),  GoF_WL.bias(size(point_zs,2)+1)]  = xb_skill_BOI([[1:length(obs_hh_all)]',obs_hh_all],   [[1:length(mod_hh_all)]',mod_hh_all]);
[GoF_Hhf.r2(size(point_zs,2)+1), GoF_Hhf.sci(size(point_zs,2)+1), GoF_Hhf.relbias(size(point_zs,2)+1), GoF_Hhf.bss(size(point_zs,2)+1), GoF_Hhf.rms(size(point_zs,2)+1), GoF_Hhf.bias(size(point_zs,2)+1)] = xb_skill_BOI([[1:length(obs_Hhf_all)]',obs_Hhf_all], [[1:length(mod_Hhf_all)]',mod_Hhf_all]);
[GoF_Hig.r2(size(point_zs,2)+1), GoF_Hig.sci(size(point_zs,2)+1), GoF_Hig.relbias(size(point_zs,2)+1), GoF_Hig.bss(size(point_zs,2)+1), GoF_Hig.rms(size(point_zs,2)+1), GoF_Hig.bias(size(point_zs,2)+1)] = xb_skill_BOI([[1:length(obs_Hig_all)]',obs_Hig_all], [[1:length(mod_Hig_all)]',mod_Hig_all]);

% GoF Table
GoFtable = [GoF_WL.rms' GoF_WL.sci' GoF_WL.bias' GoF_WL.relbias' GoF_Hhf.rms' GoF_Hhf.sci' GoF_Hhf.bias' GoF_Hhf.relbias' GoF_Hig.rms' GoF_Hig.sci' GoF_Hig.bias' GoF_Hig.relbias'];
save([dirs.output filesep 'skill.mat'], 'GoFtable');

dlm = {'$WL$\newline$_{RMSE}$','$WL$\newline$_{SCI}$','$WL$\newline$_{bias}$','$WL$\newline$_{relbias}$','$HF$\newline$_{RMSE}$','$HF$\newline$_{SCI}$',...
    '$HF$\newline$_{bias}$','$HF$\newline$_{relbias}$','$IG$\newline$_{RMSE}$','$IG$\newline$_{SCI}$','$IG$\newline$_{bias}$','$IG$\newline$_{relbias}$'};
matrix2latex(GoFtable, 'filename', 'skill.tex', ...
'rowjustification','p{0.75cm}',...
'rowlabel', {'P1','P2','P3','P4','P5','P6','P7','P8','P9','P10', 'Mean'}, ...
'collabel', dlm, ...
'format',	{'%s', '%4.2f', '%4.2f','%4.2f','%4.2f','%4.2f', '%4.2f','%4.2f','%4.2f','%4.2f', '%4.2f','%4.2f','%4.2f'});

% Linear fit
linfit_WL_p     = polyfit(mod_hh_all,obs_hh_all,1);
linfit_WL_f     = polyval(linfit_WL_p,mod_hh_all); 
linfit_Hhf_p    = polyfit(mod_Hhf_all,obs_Hhf_all,1);
linfit_Hhf_f    = polyval(linfit_Hhf_p,mod_Hhf_all); 
linfit_Hig_p    = polyfit(mod_Hig_all,obs_Hig_all,1);
linfit_Hig_f    = polyval(linfit_Hig_p,mod_Hig_all); 

%% Cross-section during storm peak
tpeakmean           = [find(refdate+seconds(meantime-tstart) == tpeakstart+seconds(dtXB_mean)) : find(refdate+seconds(meantime-tstart) == tpeakend)];  % Tijdstappen voor piek van storm in time-averaged output 
tpeakpntmeas        = [find(time.(obs_points{1}) == tpeakstart) : find(time.(obs_points{1}) == tpeakend)-1];  % Tijdsbins voor piek van storm in measured wave characteristics per point

% Golfdata langs cross-section
Hm0_hf_StormPiek    = rms(H_mean(tpeakmean,:)) .* sqrt(2);
Hm0_hf_StormPiek_min= min(H_min(tpeakmean,:) .* sqrt(2));
Hm0_hf_StormPiek_max= max(H_max(tpeakmean,:) .* sqrt(2));
Hm0_lf_StormPiek    = rms(4.*sqrt(zs1_var(tpeakmean,:)));

% Gemeten golfdata op meetlocaties tijdens piek storm
for i = 1:10
    Hm0_hf_pntmeas_StormPiek(i)    = rms(data.(obs_points{i})(tpeakpntmeas,3));
    Hm0_IG_pntmeas_StormPiek(i)    = rms(data.(obs_points{i})(tpeakpntmeas,7));
    WL_pntmeas_StormPiek(i)        = mean(data.(obs_points{i})(tpeakpntmeas,1))+InstrLocProj.z(i);
end
% Hm0_hf_pnt_StormPiek    = rms(Hm0_hf(:,tpeakpntbin)');
% Hm0_IG_pnt_StormPiek    = rms(Hm0_IG(:,tpeakpntbin)');
% WL_pnt_StormPiek        = mean(point_zs(tpeakpnt,:));

%% ------------------------------------------------------------------ %%
%  -------------------------- Plot figures -------------------------- %
%  ------------------------------------------------------------------ %

%% Compare model results with measurement data

j = 1;  % Alleen voor P1
figure(101)
    set(gcf, 'Position', [200 200 900 500])
    subplot(4,1,1)
        hold on
        a = plot(time.(obs_points{j}), data.(obs_points{j})(:,1), '.k', 'markersize', 3);
        b = plot(refdate+seconds(tsReal), mean_wl(j,:)-mean_zb(j,:), 'linewidth', 1);
        plot([refdate refdate], [-10 10], 'k--')
        plot([endref endref], [-10 10], 'k--')
        ylabel('Water depth [m]')
        set(gca,'XTicklabel',[]);
        legend([a,b], {'data', 'model'}, 'location', 'northeast')
        xlim([min(time.(obs_points{j})), max(time.(obs_points{j}))])
        ylim([0 2])

    subplot(4,1,2)
        hold on
        plot(time.(obs_points{j}), data.(obs_points{j})(:,3), '.k', 'markersize', 3)
        plot(refdate+seconds(tsReal), Hm0_hf(j,:), 'linewidth', 1)
        plot([refdate refdate], [0 7], 'k--')
        plot([endref endref], [0 7], 'k--')
        ylabel('H_{m0,G} [m]')
        set(gca,'XTicklabel',[]);
        xlim([min(time.(obs_points{j})), max(time.(obs_points{j}))])
        ylim([0 1])

    subplot(4,1,3)
        hold on
        plot(time.(obs_points{j}), data.(obs_points{j})(:,7), '.k', 'markersize', 3)
        plot(refdate+seconds(tsReal), Hm0_IG(j,:), 'linewidth', 1)
        plot([refdate refdate], [0 1.5], 'k--')
        plot([endref endref], [0 1.5], 'k--')
        ylabel('H_{m0,IG} [m]')
        set(gca,'XTicklabel',[]);     
        xlim([min(time.(obs_points{j})), max(time.(obs_points{j}))])
        ylim([0 0.5])        

    subplot(4,1,4)
        hold on
        plot(time.(obs_points{j}), data.(obs_points{j})(:,8), '.k', 'markersize', 3)
        plot([refdate refdate], [0 25], 'k--')
        plot([endref endref], [0 25], 'k--')
        ylabel('Wave period T_{m1,0} [s]')
        xlim([min(time.(obs_points{j})), max(time.(obs_points{j}))])
        ylim([0 15])
        xlabel('Date')    

     xb_write_plot(gcf, dirs.output, ['Comparison_XBeach-measured_EntireMeasurementPeriod_' obs_points{j}])



%%  
    InstrLocProj.Name = {'P1','P2','P3','P4','P5','P6','P7','P8','P9','P10'};
    % Zoom transect VOOR IN MEMO
    figure(); 
    set(gcf, 'Position', [200 200 900 500])
    h1 = subplot(2,1,1); hold on; box on; grid on;
        a = fill([x'; flipud(x')] ,[Hm0_hf_StormPiek_min'; flipud(Hm0_hf_StormPiek_max')], [0.6 0.6 0.6], 'FaceAlpha', 0.5, 'EdgeColor', 'none', 'displayname', 'Min. & max. H_{m0,hf} - Model');
        b = plot(x, Hm0_hf_StormPiek, '-b',                             'displayname', 'Mean H_{m0,hf} - Model', 'linewidth', 1.5);
        b2 = plot(InstrLocProj.R(1:10),Hm0_hf_pntmeas_StormPiek, '.b',  'displayname', 'Mean H_{m0,hf} - Measured', 'markersize', 10);
        c = plot(x, Hm0_lf_StormPiek, '-k',                             'displayname', 'Mean H_{m0,lf} - Model', 'linewidth', 1.5);
        c2 = plot(InstrLocProj.R(1:10),Hm0_IG_pntmeas_StormPiek, '.k',  'displayname', 'Mean H_{m0,ig} - Measured', 'markersize', 10);
        legend([c c2 b b2 a],'location', 'northeast')
        xlim([-1800 -100]);
        yticks([0:0.5:3])
        ylim([0 3]);
        ylabel('H_{m0} [m]')
        pos1 = get( h1, 'Position' );
        h1.Position = [pos1(1) pos1(2)-0.06 pos1(3) pos1(4)+0.04]; % [x y breedte hoogte]   
        h1.XTickLabel = [];
        set(gca, 'FontSize', 11)
    h2 = subplot(2,1,2); hold on; box on; grid on;
        fill([-8000; x'; 200] ,[-30; zb'; -30], [255 200 0]./255, 'FaceAlpha', 0.4, 'EdgeColor', 'none', 'HandleVisibility','off')
        fill([x'; flipud(x')] ,[mean(zs_mean(tpeakmean,:))'; flipud(zb')], [59 171 227]./255, 'FaceAlpha', 0.5, 'EdgeColor', 'none', 'HandleVisibility','off')
        plot(x, zb, 'displayname', 'Bed level', 'linewidth', 1.5)
        plot(x, min(zs_min(tpeakmean,:)), '--b', 'displayname', 'Min. & max. water level - Model', 'linewidth', 1.5)
        plot(x, max(zs_max(tpeakmean,:)), '--b', 'HandleVisibility','off', 'linewidth', 1.5)
        plot(x, mean(zs_mean(tpeakmean,:)), '-b', 'displayname', 'Mean water level - Model', 'linewidth', 1.5)
        plot(InstrLocProj.R(1:10), WL_pntmeas_StormPiek, '.b', 'displayname', 'Mean water level - Measured', 'markersize', 10)
        plot(InstrLocProj.R(1:10), InstrLocProj.z(1:10), '.k', 'markersize', 15, 'HandleVisibility','off')
        text(InstrLocProj.R(1:10)-20, InstrLocProj.z(1:10)+0.4, InstrLocProj.Name(1:10));
        xlim([-1800 -100]);
        ylim([-2 4])
        ylabel('Elevation [m +NAP]')
        xlabel('Distance along transect [m]')
        legend('location', 'southeast')
        pos1 = get( h2, 'Position' );
        h2.Position = [pos1(1) pos1(2) pos1(3) pos1(4)+0.04]; % [x y breedte hoogte]   
        set(gca, 'FontSize', 11)


    xb_write_plot(gcf, dirs.output, ['Cross-section_XBeach-results_PeakStorm_Zoom'])

    
    %% Figuren per meetlocatie
    for j=1:numel(obs_points)
        idStart = find(time.(obs_points{j}) <= refdate, 1, 'last');
        if isempty(idStart); idStart = 1; end
        
        idEnd = find(time.(obs_points{j}) >= endref, 1, 'first');
        if isempty(idEnd); idEnd = length(time.(obs_points{j})); end
        
        % Figuur met per meetlocatie de meetdata en modeldata over tijd (WL, Hm0 HF en Hm0 LF) VOOR IN MEMO
        figure(j)
        set(gcf, 'Position', [200 200 900 500])
        h1 = subplot(3,1,1);
            hold on; box on; grid on;
            plot(time.(obs_points{j})(idStart:idEnd), data.(obs_points{j})(idStart:idEnd,1), 'ko', 'markersize', 5)
            plot(refdate+seconds(tsReal), mean_wl(j,:)-mean_zb(j,:), 'color', h1.ColorOrder(1,:), 'linewidth', 1.5)
            ylims = [0 2.5];    txty = 0.9.*(ylims(2)-ylims(1));     txtdy = 0.12.*(ylims(2)-ylims(1));         
            ylim(ylims);
            xlim([time.(obs_points{j})(idStart) time.(obs_points{j})(idEnd)]);
            text(0.32,txty,['RMSE = ' num2str(GoF_WL.rms(j), '%.2f')]);           
            text(0.32,txty-txtdy,['sci = ' num2str(GoF_WL.sci(j), '%.2f')]); 
            text(0.32,txty-2*txtdy,['bias = ' num2str(GoF_WL.bias(j), '%.2f')]); 
            text(0.32,txty-3*txtdy,['relbias = ' num2str(GoF_WL.relbias(j), '%.2f')]); 
            ylabel('Water depth [m]')
            set(gca,'XTicklabel',[]);
            pos1 = get( h1, 'Position' );
            h1.Position = [pos1(1) pos1(2) pos1(3) pos1(4)+0.06]; % [x y breedte hoogte]    
        h2 = subplot(3,1,2);
            hold on; box on; grid on;
            plot(time.(obs_points{j})(idStart:idEnd), data.(obs_points{j})(idStart:idEnd,3), 'ko', 'markersize', 5)
            plot(refdate+seconds(tsReal), Hm0_hf(j,:), 'color', h2.ColorOrder(1,:), 'linewidth', 1.5)
            ylims = [0 1.2];    txty = 0.9.*(ylims(2)-ylims(1));     txtdy = 0.12.*(ylims(2)-ylims(1));         
            ylim(ylims);
            xlim([time.(obs_points{j})(idStart) time.(obs_points{j})(idEnd)]);
            text(0.32,txty,['RMSE = ' num2str(GoF_Hhf.rms(j), '%.2f')]);           
            text(0.32,txty-txtdy,['sci = ' num2str(GoF_Hhf.sci(j), '%.2f')]); 
            text(0.32,txty-2*txtdy,['bias = ' num2str(GoF_Hhf.bias(j), '%.2f')]); 
            text(0.32,txty-3*txtdy,['relbias = ' num2str(GoF_Hhf.relbias(j), '%.2f')]); 
            ylabel('H_{m0,G} [m]')
            set(gca,'XTicklabel',[]);
            pos1 = get( h2, 'Position' );
            h2.Position = [pos1(1) pos1(2) pos1(3) pos1(4)+0.06]; % [x y breedte hoogte]    
        h3 = subplot(3,1,3);
            hold on; box on; grid on;
            a = plot(time.(obs_points{j})(idStart:idEnd), data.(obs_points{j})(idStart:idEnd,7), 'ko', 'markersize', 5);
            b = plot(refdate+seconds(tsReal), Hm0_IG(j,:), 'color', h3.ColorOrder(1,:), 'linewidth', 1.5);
            ylims = [0 0.6];    txty = 0.9.*(ylims(2)-ylims(1));     txtdy = 0.12.*(ylims(2)-ylims(1));         
            ylim(ylims);
            xlim([time.(obs_points{j})(idStart) time.(obs_points{j})(idEnd)]);
            text(0.32,txty,['RMSE = ' num2str(GoF_Hig.rms(j), '%.2f')]);           
            text(0.32,txty-txtdy,['sci = ' num2str(GoF_Hig.sci(j), '%.2f')]); 
            text(0.32,txty-2*txtdy,['bias = ' num2str(GoF_Hig.bias(j), '%.2f')]); 
            text(0.32,txty-3*txtdy,['relbias = ' num2str(GoF_Hig.relbias(j), '%.2f')]); 
            ylabel('H_{m0,IG} [m]')
            xlabel('Date')
            leg = legend([a,b], {'data', 'model'}, 'location', 'northeast');
            leg.ItemTokenSize = [15 20];
            pos1 = get( h3, 'Position' );
            h3.Position = [pos1(1) pos1(2) pos1(3) pos1(4)+0.06]; % [x y breedte hoogte]    

        xb_write_plot(gcf, dirs.output, ['Comparison_XBeach-measured_ModelPeriod_Timeseries_' obs_points{j} ])
        clear idStart idEnd
    end 
    

    figure(201)
    set(gcf, 'Position', [200 200 1000 300])
    subplot(1,4,1); hold on; axis equal; box on; grid on;
        cmap= jet(numel(obs_points));
        for j = 1:numel(obs_points)
            scatter(data.(obs_points{j})(id{j},1),mean_wl(j,id3{j})-mean_zb(j,id3{j}), 5, cmap(j,:), 'filled')
        end
        plot([0 10], [0 10], ':k','linewidth',1.5)
        plot(linfit_WL_f,mod_hh_all,'k', 'linewidth',1.5);
        xlabel('Measured water depth [m]')
        ylabel('Modeled water depth [m]')
        xlim([0 2.5]); ylim([0 2.5]);
        set(gca,'xtick',0:0.5:2.5); set(gca,'ytick',0:0.5:2.5); 
%                 set(gca,'fontsize',6)
    subplot(1,4,2); hold on; axis equal; box on; grid on;
        for j = 1:numel(obs_points)
            scatter(data.(obs_points{j})(id{j},3),Hm0_hf(j,id3{j}), 5, cmap(j,:), 'filled')
        end
        plot([0 10], [0 10], ':k','linewidth',1.5)
        plot(linfit_Hhf_f,mod_Hhf_all,'k','linewidth',1.5);
        xlabel('Measured H_{m0,G} [m]')
        ylabel('Modeled H_{m0,G} [m]')
        xlim([0 1.2]); ylim([0 1.2]);
        set(gca,'xtick',0:0.3:1.2); set(gca,'ytick',0:0.3:1.2);
%                 set(gca,'fontsize',6)
    subplot(1,4,3); hold on; axis equal; box on; grid on;
        for j = 1:numel(obs_points)
            scatter(data.(obs_points{j})(id{j},7),Hm0_IG(j,id3{j}), 5, cmap(j,:), 'filled')
        end
        plot([0 10], [0 10], ':k','linewidth',1.5)
        plot(linfit_Hig_f,mod_Hig_all,'k','linewidth',1.5);
        xlabel('Measured H_{m0,IG} [m]')
        ylabel('Modeled H_{m0,IG} [m]')
        xlim([0 0.6]); ylim([0 0.6]);
        set(gca,'xtick',0:0.2:0.6); set(gca,'ytick',0:0.2:0.6);
%                 set(gca,'fontsize',6)
    subplot(1,4,4); hold on; axis off  % omslachtige manier om legenda toe te voegen
        for j = 1:numel(obs_points)
            scatter(data.(obs_points{j})(id{j},7),Hm0_IG(j,id3{j}), 5, cmap(j,:), 'filled')
        end
        plot([0 10], [0 10], ':k','linewidth',1.5)
        plot(linfit_Hig_f,mod_Hig_all,'k','linewidth',1.5);
        ylim([100 101]); xlim([100 101])
        legend([obs_points, 'best fit', 'linear fit'], 'location', 'west')
            
            
    xb_write_plot(gcf, dirs.output, ['Comparison_XBeach-measured_ModelPeriod_Scatterplot_AllObservationPoints'])
    
end