function [ccov] = lagcov(data,refer,lags,level,degf);
% function [ccov] = lagcov(data,refer,[lags],[level],[degf]);
%
% calculate lagged covariances
%
% input  :	data		- data array
%		refer		- reference series
%		[lags]	[0]	- vector of lags (reference leads for positive)
%
% output :	ccov		- covariance
%
% uses :	nsum.m	nmean.m
%
% version 0.1.0		last change 27.01.1999

% G.Krahmann, LDEO Palisades, Jan 1999


% history:
%


%
% check arguments
%
if nargin<2
  error('Not enough input arguments.')
end
refer = squeeze(refer);
lr = length(refer);
if ~(size(refer,1)==1 | size(refer,2)==1)
  error('Reference series must be a vector.')
end
if sum( lr==size(data) )==0
  error('Data must contain a dimension with the length of the reference.')
end
if nargin<3
  lags = 0;
end
if sum(isnan(refer))>1
  error('Reference series must not contain NaN''s')
end


%
% check for simple case
%
simple = 0;
if length(lags)==1
  if lags==0
    simple = 1;
  end
end


%
% prepare output arrays
%
ds = size(data);
di = find(ds==lr);
di = di(1);
ns = ds;
ns(di) = length(lags);
ccoef = repmat(0,ns);


% 
% subtract mean and enlarge reference
%
dr = ones(1,length(ns));
dr(di) = size(data,di);
nr = ds;
nr(di) = 1;
data = data - repmat(nmean(data,di),dr);
refer = refer - nmean(refer);
ref = repmat(0,size(data));
ref = shiftdim(ref,di-1);		% I was not able to get repamt
for n=1:length(refer)			% to work properly here
  ref(n,:) = ref(n,:) + refer(n);
end
ref = shiftdim(ref,length(ds)-di+1);


%
% do covariances
%
if simple==1

%  ref2 = refer.^2;
%  mref2 = mean(ref2);
  ccov = nsum(data.*ref,di);
  nccov = nsum( (data*0+1).*(ref*0+1),di);
  ccov = ccov./nccov;

else

  % shift dimensions to get the series into the first
  ref = shiftdim(ref,di-1);
  data = shiftdim(data,di-1);
  ccoef = shiftdim(ccoef,di-1);

  % reshape to get only two dimensions
  sr = size(ref);
  ref = reshape(ref,[sr(1),prod(sr(2:length(sr)))]);
  data = reshape(data,[sr(1),prod(sr(2:length(sr)))]);
  sc = size(ccoef);
  ccov = reshape(ccoef,[sc(1),prod(sc(2:length(sc)))]);

%  data2 = data.^2;
  for n=1:length(lags)
    if lags(n)<0
%      ref2 = refer(1-lags(n):length(refer)).^2;
%      mref2 = mean(ref2);
      d = data(1:size(data,1)+lags(n),:);
      r = ref(1-lags(n):size(ref,1),:);
      ccov(n,:) = nsum(d.*r,1);
      nccov(n,:) = nsum( (d*0+1).*(r*0+1),di);
      ccov(n,:) = ccov(n,:)./nccov(n,:);
    elseif lags(n)>0
%      ref2 = refer(1:length(refer)-lags(n)).^2;
%      mref2 = mean(ref2);
      d = data(1+lags(n):size(data,1),:);
      r = ref(1:size(ref,1)-lags(n),:);
      ccov(n,:) = nsum(d.*r,1);
      nccov(n,:) = nsum( (d*0+1).*(r*0+1),di);
      ccov(n,:) = ccov(n,:)./nccov(n,:);
    elseif lags(n)==0
%      ref2 = refer.^2;
%      mref2 = mean(ref2);
      d = data;
      r = ref;
      ccov(n,:) = nsum(d.*r,1);
      nccov(n,:) = nsum( (d*0+1).*(r*0+1),di);
      ccov(n,:) = ccov(n,:)./nccov(n,:);
    end
  end 

  % reshape back
  ref = reshape(ref,sr);
  data = reshape(data,sr);
  ccov = reshape(ccov,sc);

  % shift dimensions back
  ref = shiftdim(ref,length(ds)-di+1);
  data = shiftdim(data,length(ds)-di+1);
  ccov = shiftdim(ccov,length(ds)-di+1);

end
