function Analyse_Buckley(sim,name2,output_folder, data_folder)

addpath 'p:\xbeach\msc_studies\Menno_deRidder\Msc\Scripts\Matlab\Menno_Tools\'

%% input
name        = 'sim';
idtime      = [500    2499];
fcutoff     = 0.2;           
fres        = 0.025;
rough       = 1; %1 smooth; 2: rough
%% load data
ppath = data_folder;
load([ppath '/data']);

for i=1:20
    time(i,:) = M(rough).R(sim).g(i).ts(:,1);
    id        = time(i,:) > 100 & time(i,:) < 100000000;
    wl(i,:)   = M(rough).R(sim).g(i).ts(:,2);
    [fx1, data.vy(i,:)]       = spectrumsimple(0.025, wl(i,id), 'fresolution', fres,'detrend',1, 'figures', 0);
    %% Hs
    f                         = fx1(fx1>fcutoff);
    var                       = data.vy(i,fx1>fcutoff);
    data.Hs(i)                = swan_hs(f,var);
    %% spectra properties

    [m0 ,fminlim,fmaxlim]  = computeMoment(  0, fx1(2:end), data.vy(i,2:end), 0.03, 1.2, 'no' );
    [m_1,fminlim,fmaxlim]  = computeMoment( -1, fx1(2:end), data.vy(i,2:end), 0.03, 1.2, 'no' );
    [m1,fminlim,fmaxlim]  = computeMoment( 1, fx1(2:end), data.vy(i,2:end), 0.03, 1.2, 'no' );
    [m2,fminlim,fmaxlim]  = computeMoment( 2, fx1(2:end), data.vy(i,2:end), 0.03, 1.2, 'no' );
    
    data.Tm01(i)  = m0  / m1;
    data.Tm02(i)  = sqrt( m0 / m2 );
    data.Tm_10(i) = m_1 / m0;
    
    q = waveProperties(time(i,:),wl(i,id));
    data.skew(i)      = q.skewness;
    data.group(i)     = q.groupiness; 
    data.asymmetry(i) = q.asymmetry;
    %% LF waves
    f_l                         = fx1(fx1>fcutoff/10&fx1<fcutoff);
    var_l                       = data.vy(i,fx1>fcutoff/10&fx1<fcutoff);
    data.Hs_low(i)              = swan_hs(f_l, var_l);

    data.eta(i) = M(1).R(sim).g(i).SetUpZero;
    data.time(i,:) = time(i,:);
    data.wl(i,:) = wl(i,:);
    data.fx = fx1;
    [data.wl_l(i,:), HH] = low_filter(data.time(i,:) ,data.wl(i,:),fcutoff/10,fcutoff);
end
EMS = [12,15,17,19];
for i=1:4
     udat   = M(rough).R(sim).g(EMS(i)).ts(:,3);
     vdat   = M(rough).R(sim).g(EMS(i)).ts(:,4);
     bed         = M(rough).R(sim).gz00;
     [data.zsin(EMS(i),:), data.zsout(EMS(i),:), uin, uout, reflc] = Guza_split_waves_menno(data.time(EMS(i),:),data.wl(EMS(i),:),udat,-1*bed(EMS(i)),'boundin',fcutoff,0.5, 1);
     %%
     
     [fx1_in, vy_in]       = spectrumsimple(0.025, data.zsin(EMS(i),:), 'fresolution', 0.05,'detrend',1, 'figures', 0);
     [fx1_out, vy_out]       = spectrumsimple(0.025, data.zsout(EMS(i),:), 'fresolution', 0.05,'detrend',1, 'figures', 0);
     
     f_l                         = fx1(fx1>fcutoff/10&fx1<fcutoff);
     var_l                       = vy_in(fx1>fcutoff/10&fx1<fcutoff);
     data.Hs_low_in(EMS(i))              = swan_hs(f_l, var_l);
     var_l                       = vy_out(fx1>fcutoff/10&fx1<fcutoff);
     data.Hs_low_out(EMS(i))              = swan_hs(f_l, var_l);
end
data.Hs_low_in(20) = 0;
data.Hs_low_out(20) = 0;

figure
plot(time(6,:),wl(6,:))
ylabel('wl [m]')

save_folder = [output_folder];

saveas(gcf, [save_folder, '\wl.png']);

%% read xboutput
% points
fname            = [output_folder '\xboutput.nc'];
point_zs         = nc_varget(fname, 'point_zs'); 
point_u          = nc_varget(fname, 'point_u'); 
point_zb         = nc_varget(fname, 'point_zb'); point_zb = squeeze(point_zb(1,:));
globalx          = nc_varget(fname, 'globalx'); 
t                = nc_varget(fname, 'globaltime');
% global
zb               = nc_varget(fname, 'zb'); zb = squeeze(zb(1,:,:));
zs               = nc_varget(fname, 'zs');
u                = nc_varget(fname, 'u');
hh               = nc_varget(fname, 'hh');
point_time       = nc_varget(fname, 'pointtime');
pointx           = nc_varget(fname, 'pointx');

new_time = linspace(0,point_time(end-1),length(point_time));
point_time = interp1(point_time(1:end-1),point_time(1:end-1),new_time);
idwanted    = point_time > idtime(1) & point_time < idtime(2);




%% create spectrum

for i = 1:length(pointx)
    % Guza filter
    result.time(i,:)    = point_time(idwanted);
    result.wl(i,:)      = squeeze(point_zs(idwanted,i));
    result.eta(i)       = mean(result.wl(i,:));
    [result.zsin(i,:), result.zsout(i,:), uin, uout, result.reflc(i)] = Guza_split_waves_menno(point_time(idwanted),squeeze(point_zs(idwanted,i)),squeeze(point_u(idwanted,i)),squeeze(point_zb(i)),'boundin',fcutoff,0.5, 1);

    
    % compute Hrms from serie
    result.Hrms(i)      = sqrt(8*mean(result.zsin(i,:).^2));
    % spectra
    [fx_in, result.vy_in(i,:)]     = spectrumsimple(point_time(idwanted), result.zsin(i,:), 'fresolution', fres,'detrend',1, 'figures', 0);
    [fx, result.vy(i,:)]        = spectrumsimple(point_time(idwanted), squeeze(point_zs(idwanted,i)), 'fresolution', fres,'detrend',1, 'figures', 0);
    [fx_out, result.vy_out(i,:)]    = spectrumsimple(point_time(idwanted), result.zsout(i,:), 'fresolution', fres,'detrend',1, 'figures', 0);
    result.fx = fx;
    
    [m0 ,fminlim,fmaxlim]  = computeMoment(  0, fx(2:end), result.vy(i,2:end), 0.03, 1.2, 'no' );
    [m_1,fminlim,fmaxlim]  = computeMoment( -1, fx(2:end), result.vy(i,2:end), 0.03, 1.2, 'no' );
    [m1,fminlim,fmaxlim]  = computeMoment( 1, fx(2:end), result.vy(i,2:end), 0.03, 1.2, 'no' );
    [m2,fminlim,fmaxlim]  = computeMoment( 2, fx(2:end), result.vy(i,2:end), 0.03, 1.2, 'no' );

    result.Tm_10(i) =  m_1 / m0;
    result.Tm01(i) = m0  / m1;
    result.Tm02(i) = sqrt( m0 / m2 );
    
    q = waveProperties(point_time(idwanted),result.wl(i,:))
    
    result.skew(i)      = q.skewness;
    result.group(i)     = q.groupiness;   
    result.asymmetry(i) = q.asymmetry;
    
    % Hs
    f                         = fx(fx>fcutoff);
    var                       = result.vy(i,fx>fcutoff);
    result.Hs(i)                = swan_hs(f,var);
    result.Hs_in(i)             = swan_hs(fx_in,result.vy_in(i,:));
    result.Hs_out(i)            = swan_hs(fx_in,result.vy_out(i,:));
    
    
    % low freq
    f_l                         = fx(fx>fcutoff/10&fx<fcutoff);
    var_l                       = result.vy(i,fx>fcutoff/10&fx<fcutoff);
    result.Hs_low(i)            = swan_hs(f_l, var_l);
    
    f_l                         = fx_in(fx_in>fcutoff/10&fx_in<fcutoff);
    var_l                       = result.vy_out(i,fx_in>fcutoff/10&fx_in<fcutoff);
    result.Hs_low_out(i)        = swan_hs(f_l, var_l);
    var_l                       = result.vy_in(i,fx_in>fcutoff/10&fx_in<fcutoff);
    result.Hs_low_in(i)         = swan_hs(f_l, var_l);
    
    [result.wl_l(i,:), HH]      = low_filter(result.time(i,:),result.wl(i,:),fcutoff/10,fcutoff);
end


if rough==2 & canopy==1 
    %% analyse
    ucan = nc_varget(fname, 'ucan');
    point_ucan = nc_varget(fname, 'point_ucan'); 
    Fvegu = nc_varget(fname, 'Fvegu');
    point_Fvegu = nc_varget(fname, 'point_Fvegu'); 
    dzsdx = nc_varget(fname, 'point_dzsdx'); 
    u = nc_varget(fname, 'point_u'); 
    h = nc_varget(fname, 'point_hh'); 
    
    result.avg_F = mean(point_Fvegu(idwanted,:));
    result.rms_ucan = sqrt(mean(point_ucan(idwanted,:).^2));
    result.rms_u = sqrt(mean(point_u(idwanted,:).^2));
    result.rms_F = sqrt(mean(point_Fvegu(idwanted,:).^2));
    result.avg_dzsdx = mean(dzsdx(idwanted,:));
    result.avg_ucan = mean(point_ucan(idwanted,:));
    result.avg_u = mean(point_u(idwanted,:));
    result.inertia = mean(diff( (point_ucan(idwanted,:))/(point_time(2)-point_time(1))));
    result.drag = mean(point_ucan(idwanted,:).*abs(point_ucan(idwanted,:)));
    result.shear = mean((point_u(idwanted,:)- point_ucan(idwanted,:)).*abs(point_u(idwanted,:)- point_ucan(idwanted,:)));
    
    result.inertia_rms = sqrt(mean(( diff( (point_ucan(idwanted,:))/(point_time(2)-point_time(1)))).^2));
    result.drag_rms = sqrt( mean((point_ucan(idwanted,:).*abs(point_ucan(idwanted,:))).^2));
    result.shear_rms = sqrt(mean(((point_u(idwanted,:)- point_ucan(idwanted,:)).*abs(point_u(idwanted,:)- point_ucan(idwanted,:))).^2));
    result.rms_dzsdx = sqrt(mean(dzsdx(idwanted,:).^2));
    n=0.05;
    Cf = sqrt(9.81*n^2./(h.^(1.0/12)));
    fric = 1000* Cf .* abs(u) .* u;
    result.fric_rms = sqrt(mean(fric(idwanted,:).^2));
    result.fric = mean(fric(idwanted,:));
end

%% Spectra
loc = [1,10];

figure('pos',[100,100,1200,600])
for i=1:length(loc)
    subplot(2,1,i)
    plot(fx,result.vy(loc(i),:),'Linewidth',2);
    hold on
    plot(fx_in,result.vy_in(loc(i),:),'Linewidth',2);
    hold on
    plot(fx_in,result.vy_out(loc(i),:),'Linewidth',2);
    hold on
    plot(fx1, data.vy(i,:),'Linewidth',2);
    xlim([0,3])
    %ylim([0, 0.03])
    title(pointx(loc(i)))
    xlabel('Freqentie [Hz]'); ylabel('Variance density [m^2]');
    legend('Total','Incoming','Outgoing','data')
    saveas(gcf, [save_folder, '\freq_variance.png']);
end
%% plot Hrms
figure('pos',[100,100,1200,600])
subplot(3,1,1)
plot(globalx(:,1), zb(:,1),'linewidth',2)
xlim([5 40])

subplot(3,1,2)

plot(pointx, result.Hs(:)/sqrt(2),'linewidth',2)
hold on
plot(pointx, result.Hs_in(:)/sqrt(2),'linewidth',2)
hold on
plot(pointx, result.Hs_out(:)/sqrt(2),'linewidth',2)
hold on
plot(pointx, data.Hs(:)/sqrt(2),'linewidth',2)
ylabel('H_{rms}[m]')
xlabel('Distance [m]')
legend('Total','Incoming','Outgoing','Data')

subplot(3,1,3)
plot(pointx, result.Hs_low(:),'linewidth',2)
hold on
plot(pointx, result.Hs_low_in(:),'linewidth',2)
hold on
plot(pointx, result.Hs_low_out(:),'linewidth',2)
hold on
plot(pointx, data.Hs_low(:),'linewidth',2)
hold on; plot(pointx, data.Hs_low_in(:),'<')
hold on; plot(pointx, data.Hs_low_out(:),'>')


ylabel('Hs[m]')
xlabel('Distance [m]')
legend('Total','Incoming','Outgoing','Data')
title('low frequency')
saveas(gcf, [save_folder, '\Hs_analyse.png']);


figure;
plot(pointx, data.eta,'linewidth',2)
hold on
plot(pointx, result.eta,'linewidth',2)
legend('data','XBeach')
saveas(gcf, [save_folder, '\eta_vali.png']);




save([save_folder  '\data.mat'], 'data')
save([save_folder  '\result.mat'], 'result')
model.zs = zs;
model.zb = zb;
model.pointx = pointx;
model.globalx = globalx;
model.hh = hh;
save([save_folder '\model.mat'], 'model')
clear all
close all
