function [filter,C]=power_law_cov(time,sigma,spectral_index)

%    [filter,C]=power_law_cov(time,sigma,spectral_index)
%        or
%    [filter]=power_law_cov(time,sigma,spectral_index)
%
%    Compute the covariance matrix and filter matrix (from white)
%    for any power-law process using fractional differencing.
%    Correlated data is then equal to filter*randn(n,1).
%    If you only want the filter, use the second version, it is MUCH faster.
%    For large number of data (N=length(time)) this can use a lot of
%    memory quickly: size(C)=n^2*8 bytes.
%
%    filter         filter matrix to change white noise to correlated
%    C              (optional) covariance matrix
%
%    time           time vector of data points (years) not necessarily evenly
%                   spaced.
%    sigma          standard deviation of noise (meters)
%    spectral_index exponent on frequency for power-law process,
%                   all values work, not just integers.
%                   (e.g., -1 flicker, -2 random walk)

% Hadley Johnson, hjohnson@ucsd.edu, 12-Oct-1998 version

time=time(:);
num_data=length(time);
delta_time=[0;diff(time)];
 constant=sigma*delta_time'.^(-spectral_index/4);
% constant=sigma*realpow(delta_time',(-spectral_index/4));
if spectral_index>0,
  constant(1)=0;
end
C=constant(ones(length(constant),1),:);
filter=fractional_diff(num_data,spectral_index/2);
filter=filter.*C;

%For some reason using the sparce matrix routines limits the size of memory
%useage for the C matrix by a LOT!  This is a MATLAB problem.  Don't use them.

% reuse the C variable to save memory space, this is just a
% temporary variable.

%   C=sparse(filter);
%   C=full(C*C');
% -or-
%   C=filter*filter';

if nargout>=2,
  C=filter*filter';
end
