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;
