clear all
close all



%% input
sim         = 12;
point_id = 15;
idtime      = [500    2499];

name        = 'sim';

rough       = 1; %1 smooth; 2: rough
%% load data
ppath = 'data';

load([ppath '/data'])

time = M(rough).R(sim).g(point_id).ts(:,1);
id        = time > 100 & time < 100000000;
wl        = M(rough).R(sim).g(point_id).ts(:,2);

%%
folder = 'sims_default\';
name2 = '_nh+2DV' 
fname            = [folder name num2str(sim) name2 '\xboutput.nc'];
point_zs         = nc_varget(fname, 'point_zs'); 
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);
wl2         = squeeze(point_zs(idwanted,point_id));
%%



figure;
plot(time,wl-mean(wl))
hold on
plot(point_time(idwanted), wl2-mean(wl2))
title(pointx(point_id))
legend('data','xb')

[fx, vy] = spectrumsimple(0.025,wl,'fresolution', 0.05, 'figures', 0)
[fx2, vy2] = spectrumsimple(point_time(idwanted),wl2,'fresolution', 0.05, 'figures', 0)

figure;
plot(fx,vy)
hold on
plot(fx2,vy2)
title(pointx(point_id))


q1 = waveProperties(time,wl);
q2 = waveProperties(point_time(idwanted),wl2);

[m0 ,fminlim,fmaxlim]  = computeMoment(  0, fx(2:end)', vy(2:end), 0.03, 1.2, 'no' );
[m_1,fminlim,fmaxlim]  = computeMoment( -1, fx(2:end)', vy(2:end), 0.03, 1.2, 'no' );
[m1,fminlim,fmaxlim]  = computeMoment( 1, fx(2:end)', vy(2:end), 0.03, 1.2, 'no' );
[m2,fminlim,fmaxlim]  = computeMoment( 2, fx(2:end)', vy(2:end), 0.03, 1.2, 'no' );

Tm_10 =  m_1 / m0;
Tm01 = m0  / m1;
Tm02 = sqrt( m0 / m2 );

