function [result] = xb_analyse_basic_output(pointtime, point_zs, pointx,fcutoff,advanced)
    fres = 0.01;
    for ii=1:length(pointx)
        % --- total waves
        int_time                            = linspace(pointtime(1),pointtime(end),length(pointtime)); % not completly correct; wl should be interpolated as well!
        [result(ii).fx, result(ii).vy]      = spectrumsimple(int_time, point_zs(:,ii), 'fresolution', fres,'detrend',1, 'figures', 0);

        m0                  = trapz(result(ii).fx,result(ii).vy);
        result(ii).Hm0_tot  = 4.*sqrt(m0);
        result(ii).Hrms     = sqrt(8*nanmean((point_zs(:,ii)-nanmean(point_zs(:,ii))).^2));
        result(ii).eta      = mean(point_zs(:,ii),'omitnan');
        
        % --- LF
        f_l                 = result(ii).fx(result(ii).fx>fcutoff/100 & result(ii).fx<fcutoff);
        vy_l                = result(ii).vy(result(ii).fx>fcutoff/100 & result(ii).fx<fcutoff);
        result(ii).Hm0_LF    = 4.*sqrt(trapz(f_l, vy_l));
		
        % moments
        [m0 ,fminlim,fmaxlim]   = computeMoment(  0, result(ii).fx', result(ii).vy, 0.03, 2, 'no' );
        [m_1,fminlim,fmaxlim]   = computeMoment( -1, result(ii).fx', result(ii).vy, 0.03, 2, 'no' );
        [m1,fminlim,fmaxlim]    = computeMoment( 1, result(ii).fx', result(ii).vy, 0.03, 2, 'no' );
        [m2,fminlim,fmaxlim]    = computeMoment( 2, result(ii).fx', result(ii).vy, 0.03, 2, 'no' );
        result(ii).Tm01                  = m0  / m1;
        result(ii).Tm02                  = sqrt( m0 / m2 );
        result(ii).Tm_10                 = m_1 / m0;
        
        
        if advanced
            wl  =point_zs(:,ii);
            result(ii).skewness             = nanmean((wl-nanmean(wl)).^3)/nanmean((wl-nanmean(wl)).^2).^(1.5);
            result(ii).A                    = abs(hilbert(wl)-nanmean(wl));
            result(ii).H                    = imag(hilbert(wl)-nanmean(wl));
            result(ii).groupiness           = sqrt(std(result(ii).H))/mean(result(ii).A);
            % assymetry: Elgar and Guza (1985)
            result(ii).asymmetry            = nanmean(result(ii).H.^3)/nanmean( ( (wl-nanmean(wl) ).^2 ).^(1.5) );
        end
        
    end
    
    
end

function [moment,fminlim,fmaxlim] = computeMoment( n,f,P,fmin,fmax, ltail );
% --------------------------------------------------------------
% This function computes the moment, which is equal to the integral
% of function P*f^n over frequency range [fmin,fmax]. The 
% trapezoidal rule is employed. An analytical tail is added, provided
% ltail = 'yes'
%
% This function is tested to give the correct results. It is
% tested for linear functions, for which it gives, as it should,
% the exact solution. 
%
% See also document FrequencyIntegration.pdf
% --------------------------------------------------------------

% === Some useful quantities
imin  = find( f <= fmin, 1, 'last' );
imax  = find( f >= fmax, 1, 'first');
if ( isempty(imin) )
    imin=1;
end
if ( isempty(imax) )
    imax=length(f);
end
f0    = f(imin);
f1    = f(imin+1);
fN_1  = f(imax-1);
fN    = f(imax);

% === Take possibility of fmin<f0 and fmax>fN into account
fminlim = max( f0, fmin );
fmaxlim = min( fN, fmax );

% === Spectrum weighted with frequency
Pf = P .* f.^n;

% === Integration over part I 
df = f1 - f0;
moment_pt1 = 0.5 * ( f1   - fminlim ) * ( Pf(imin)*(f1 - fminlim) + Pf(imin+1)*(fminlim - f0 + df) ) / df;

% === Integration over part III
df = fN - fN_1;
moment_pt3 = 0.5 * ( fmaxlim - fN_1 ) * ( Pf(imax-1)*(df + fN - fmaxlim) + Pf(imax)*(fmaxlim - fN_1) ) / df;

% === Integration over part II
moment_pt2 = 0;
for i = (imin+1):(imax-2)
    moment_pt2 = moment_pt2 ...
               + 0.5 * ( f(i+1) - f(i) ) * ( Pf(i) + Pf(i+1) );
end

% === Add diagnostic tail
if ( strcmp(ltail,'yes') & (fmax>fN) )
    moment_tail = P(imax) * ( fmax^(n-3) - fN^(n-3) ) / (n - 3);
else
    moment_tail = 0;
end

% === Add all contributions
moment = moment_pt1 + moment_pt2 + moment_pt3 + moment_tail;
end