function analyze_this(info,dirs)

%% Part 1: measurements
%%
% Getting right information out of the data
day         = 13;
id          = 16;   % hour interval (6 corresponds to 7th hour)
iy          = 41;

% Read measurements
hrms_meas   =   load(fullfile(dirs.data,'del_hrms.mat'));
                load(fullfile(dirs.data,['sp13_' num2str(id) '_' num2str(id+1) ]));
                load(fullfile(dirs.data,'vmean.mat'));

% Adapt H_hf
Hrms_m      = interp2(hrms_meas.xmeter,hrms_meas.thm,hrms_meas.Hrms,hrms_meas.xmeter,day+(id+0.5)/24);

% Adapth H_VLF
fi          = linspace(0, 0.004, 1000); 
Sl          = interp1(fm,Sn,fi);
Hrms_VLFm   = sqrt(8*trapz(fi,Sl));   % assuming the transfer function equals one for infragravity waves

% Adapt H_lf
fi          = linspace(0.004, (0.0989/2), 1000);
Sl          = interp1(fm,Sn,fi);
Hrms_lo     = sqrt(8*trapz(fi,Sl));   % assuming the transfer function equals one for infragravity waves

% Change vmean
jm          = find(x2>day+id/24 & x2<day+(id+1)/24);
jm          = jm(1);
vmean3i(738:1016,3) = 0.5*(vmean3i(738:1016,2) + vmean3i(738:1016,4));
vmean3i(360:500,7)  = 0.5*(vmean3i(360:500,6)  + vmean3i(360:500,8));
vmean3i(600:670,7)  = 0.5*(vmean3i(600:670,5)  + vmean3i(600:670,8));
vmean3i(500:1016,6) = 0.5*(vmean3i(500:1016,5) + vmean3i(500:1016,7));
vm = 0.5*(vmean3i(jm,:)+vmean3i(jm+1,:))/100;

% Change grid
xmeasurement   = 950 - xmeter;

% Remove station 60, which contains no data
Hrms_lo(6)  = [];
Hrms_m(6)   = [];
Hrms_VLFm(6) = [];
xmeasurement(6) = [];
vm(6) = [];

%% Part 2: XB simulation

globalx = nc_varget([pwd filesep 'xboutput.nc'],'globalx');
%ve_mean = nc_varget([pwd filesep 'xboutput.nc'],'ve');

% read dimensions
d           = xb_read_dims;

% read output
xbo1        = xb_read_output(pwd, 'vars', 'point_*');
xbo2        = xb_read_output(pwd, 'vars', {'H' 'zs' 've'});

% change waves in structure
zs              = xs_get(xbo2,'zs');
H              = xs_get(xbo2,'H');
x               = xs_get(xbo2, 'DIMS.globalx_DATA');       
y               = xs_get(xbo2, 'DIMS.globaly_DATA');   
time            = xbo2.data(1).value.data(19).value;
[nx ny]         = size(x);

%  Split waves
for i = 1:nx;
    for j = 1:ny
        
        Testz                       = zs(:,i,j);
        Testzf                      = detrend(Testz);
        rem                         = Testz-Testzf;

        [wav_high1 , wav_low1]      = filterhjb(Testzf , 0.004 , 1/3568);
        [wav_high2 , wav_low2]      = filterhjb(wav_high1 , 0.0989/2 , 1/3568);
       
        var_wave                    = var(wav_low2);
        Hrms_lf(:, i,j)             = squeeze(sqrt(8*abs(var_wave)));
        var_wave                    = var(wav_high2);
        Hrms_hf1(i,j)               = squeeze(sqrt(8*abs(var_wave)));
      
    end
end

Hrms_hf           = (Hrms_hf1.^2 + squeeze(mean(H)).^2).^0.5;
Hrms_t            = (Hrms_hf.^2 + squeeze(Hrms_lf).^2).^0.5;
vmean             = squeeze(mean(xs_get(xbo2, 've')));

% Make structure of hydro
xbw = xs_empty();
xbw = xs_set(xbw, 'DIMS', xs_get(xbo2, 'DIMS'));

if ~isscalar(Hrms_hf);  xbw = xs_set(xbw, 'Hrms_hf',    squeeze( mean (Hrms_hf([20:40],:)))' ); end 
if ~isscalar(Hrms_lf);  xbw = xs_set(xbw, 'Hrms_lf',    squeeze( mean (Hrms_lf(1,[20:40],:))) ); end 
if ~isscalar(Hrms_t);   xbw = xs_set(xbw, 'Hrms_t',     squeeze( mean (Hrms_t([20:40],:)))' ); end 
if ~isscalar(vmean);    xbw = xs_set(xbw, 'vmean',      squeeze( mean (vmean([20:40],:)))' ); end 



% Remove station 60, which contains no data
Hrms_lo(6)  = NaN;
Hrms_m(6)   = NaN;

Hrms_t_xb = squeeze( mean (Hrms_t([20:40],:)))';
Hrms_lf_xb = squeeze( mean (Hrms_lf(1,[20:40],:)));
vmean_xb = squeeze( mean (vmean([20:40],:)));
figure;
subplot('position',[0.1 0.6 0.6 0.3]);
plot(x(1,:), Hrms_t_xb,'r-','linewidth',2)
hold on
plot(x(1,:), Hrms_lf_xb,'b-','linewidth',2)
plot(xmeasurement, Hrms_m(:),'ko')
plot(xmeasurement, Hrms_lo(:),'k<')
grid on
xlabel('Distance [m]');ylabel('Height [m]'); title('Wave heights and water levels')
legend({'Computed H_{m0,HF}','Computed H_{m0,LF}','Observed \eta','Observed H_{m0,HF}','Observed H_{m0,LF}','Observed \eta', },'Location',[.75 .6 .2 .3])
subplot('position',[0.1 0.1 0.6 0.3]);
plot(x(1,:), vmean_xb,'r-','linewidth',2)
hold on
plot(xmeasurement, vm(:),'ko')

grid on
xlabel('z [m]');ylabel('Distance [m]');
legend({'v mean xb','v mean obs'},'Location',[.75 .1 .2 .3])
xb_write_plot(gcf, dirs.output, 'hydro');

point       = cell(1,sum(strfilter({xbo1.data.name}, 'point_*')));
[point{:}]  = xs_get(xbo1, 'point_*');
point       = cell2mat(point);
point       = point(:,[1 3 6 9]);

xbs         = xb_get_spectrum(...
    point(1:1200,:),...
    'sfreq',1,...
    'fsplit',0.05,...
    'fcutoff',.25,...
    'detrend',true);

xb_plot_spectrum(xbs, ...
    'measured', [fm(:) Sn]);

xb_write_plot(gcf, dirs.output, 'spectra');

% write skill table
xb_write_skilltable(                                        ...
    xbw,                                                    ...
    {[xmeasurement(:) Hrms_lo(:)] [xmeasurement(:) Hrms_m(:)] [xmeasurement(:) vm(:)]},      ...
    'vars', {'Hrms_lf' 'Hrms_hf' 'vmean'},                  ...
    'file', fullfile(dirs.output, 'skills.tex')                 );

%% Private function
    function [wav_high , wav_low] = filterhjb(wav , f_cutoff , df , av_yes);

    % function [wav_tot_high , wav_tot_low] = filterhjb(wav , f_cutoff , df , av_yes);
    %
    % This function divides a wave signal in a high frequency part and a low frequency part
    % separated by a cutoff frequency.
    % Input variables:
    %   - wav       : wave signal, series in columns, must be even number of points;
    %   - f_cutoff  : cut-off frequency;
    %   - df        : frequency step;
    %   - av_yes    : include average value yes = 1, no = 0 (default is no);
    %
    % Author: Henk Jan Bakkenes

    if nargin <4
        av_yes = 0;
    end

    ii_cutoff = round(f_cutoff/df); % calculates integer for cut off frequency

    A_wav = fft(wav);

    A_high = zeros(size(wav));
    A_low  = zeros(size(wav));

    %size(wav);
    %length(wav);
    A_high(ii_cutoff+2:round(length(wav)/2),:) = A_wav(ii_cutoff+2:round(length(wav)/2),:);

    if av_yes
        A_low(1,:) = 0.5*A_wav(1,:);
        A_low(2:ii_cutoff+1,:) = A_wav(2:ii_cutoff+1,:);
    else
        A_low(2:ii_cutoff+1,:) = A_wav(2:ii_cutoff+1,:);
    end

    wav_high = 2*real(ifft(A_high)); % 2 compensates for the missing 'mirror' fourier coefficients
    wav_low = 2*real(ifft(A_low));
