function [DataOut,HeaderOut,DSpecOut,HSpecOut] = scalarstats(w,scalar,DataIn,descriptor,DSpecIn,sf,Nf,index,flag)
%function [DataOut,HeaderOut,DSpecOut,HSpecOut] = fscalar(w,scalar,DataIn,descriptor,DSpecIn,sf,Nf,index)
% calculate scalar flux statistics and spectra.
% accepts an input Nf that can be the number of frequency bins to trim to,
% or, alternatively, can be a reference frequency vector and the spectra
% will be interpolated (linearly in log space) to these frequencies

% flag can be set to 'detrend' for detrending timeseries
if nargin<9
  flag=[];
end

w=w(:);
scalar=scalar(:);

if nargin<5
  DSpecOut=[];
  HSpecOut=[];
  if nargin==5
    sf=DSpecIn;
    DSpecIn=[];
    HSpecIn=[];
  end
end

if isempty(w)
  Data=NaN*ones(5,1);
else

  % statistics
  means=mean([w scalar]);
  if strcmp(flag,'detrend')
    covmatrix=cov([detrend(w) detrend(scalar)]);
  else
    covmatrix=cov([w scalar]);
  end

  Data(1) = means(1);
  Data(2) = means(2);
  Data(3) = covmatrix(1,1);
  Data(4) = covmatrix(2,2);
  Data(5) = covmatrix(1,2);

  Header = [descriptor ': w mean'];
  Header = strvcat(Header,[descriptor ': scalar mean']);
  Header = strvcat(Header,[descriptor ': w variance']);
  Header = strvcat(Header,[descriptor ': scalar variance']);
  Header = strvcat(Header,[descriptor ': w-scalar covariance']);

  %%% process the spectra
  if nargout > 2
    if strcmp(flag,'detrend')
      Pwc=hspec_ss(detrend(w),detrend(scalar),sf);
    else
      Pwc=hspec_ss(w,scalar,sf);
    end

    [dum,dum2]=ogive(Pwc(:,1),Pwc(:,2)); Pwc(:,6)=dum;
    [dum,dum2]=ogive(Pwc(:,1),Pwc(:,3)); Pwc(:,7)=dum;
    [dum,dum2]=ogive(Pwc(:,1),Pwc(:,4)); Pwc(:,8)=dum; %high to low freq
    [dum,dum2]=ogive(Pwc(:,1),Pwc(:,4),'reverse'); Pwc(:,9)=dum; %low to high freq

    if nargin<7|isempty(Nf)
      Nf=size(Pwc,1);
    elseif length(Nf)>1

      % need to integrate these spectra to this frequency axis
      % do a linear interpolation in log spaced frequencies
      Faxis=Nf(:); Nf=length(Faxis);

      Pwc_=Pwc;clear Pwc
      Pwc(:,1)=Faxis; Pwc(:,9)=flipud(Faxis);
      Pwc(:,2)=interp1(log(Pwc_(:,1)),Pwc_(:,2),log(Faxis),'linear');
      Pwc(:,3)=interp1(log(Pwc_(:,1)),Pwc_(:,3),log(Faxis),'linear');
      Pwc(:,4)=interp1(log(Pwc_(:,1)),Pwc_(:,4),log(Faxis),'linear');
      Pwc(:,5)=interp1(log(Pwc_(:,1)),Pwc_(:,5),log(Faxis),'linear');
      Pwc(:,6)=interp1(log(Pwc_(:,1)),Pwc_(:,6),log(Faxis),'linear','extrap');
      Pwc(:,7)=interp1(log(Pwc_(:,1)),Pwc_(:,7),log(Faxis),'linear','extrap');
      Pwc(:,8)=interp1(log(Pwc_(:,1)),Pwc_(:,8),log(Faxis),'linear','extrap');
      Pwc(:,9)=interp1(log(Pwc_(:,1)),Pwc_(:,9),log(Faxis),'linear','extrap');
    end
    Spec=NaN*ones(Nf*9,1);
    Spec=[ Pwc(1:Nf,1) ; Pwc(1:Nf,2) ; Pwc(1:Nf,3) ; Pwc(1:Nf,4) ; Pwc(1:Nf,5) ; Pwc(1:Nf,6) ; Pwc(1:Nf,7) ; Pwc(1:Nf,8) ; Pwc(1:Nf,9)];
  end
end

if isempty(DataIn)
  DataOut   = [ Data(:) ];
elseif ~isempty(index)
  DataOut=DataIn;
  DataOut(:,index) = Data(:);
else
  DataOut   = [ DataIn  Data(:) ];
end
HeaderOut = Header;

if nargout > 2
  SpecHeader = [descriptor ': frequency, Nbins = ' int2str(Nf) ];
  SpecHeader = strvcat(SpecHeader,[descriptor ': w power spectrum']);
  SpecHeader = strvcat(SpecHeader,[descriptor ': scalar power spectrum']);
  SpecHeader = strvcat(SpecHeader,[descriptor ': w-scalar cospectrum']);
  SpecHeader = strvcat(SpecHeader,[descriptor ': w-scalar quadspectrum']);
  SpecHeader = strvcat(SpecHeader,[descriptor ': w ogive']);
  SpecHeader = strvcat(SpecHeader,[descriptor ': scalar ogive']);
  SpecHeader = strvcat(SpecHeader,[descriptor ': w-scalar ogive high to low']);
  SpecHeader = strvcat(SpecHeader,[descriptor ': w-scalar ogive low to high']);

  if isempty(DSpecIn)
    DSpecOut   = [ Spec(:) ];
  elseif ~isempty(index)
    DSpecOut=DSpecIn;
    DSpecOut(:,index) = Spec(:);
  else
    DSpecOut = [DSpecIn  Spec(:)];
  end
  HSpecOut = SpecHeader;
end

return