function [COR,CORp]=lag_cor(D1,D2,lag,varargin)
% [COR,CORp]=lag_cor(D1,D2,lag,varargin)
% berechnet die correlatation mit Verschiebungen zwischen den
% Reihen D1 und D2.
% D1,D2 muessen den gleichen Anfangszeitpunkt haben
% und den gleichen Increment, brauchen aber nicht gleich
% lang zu sein.
% D1 und D2 duerfen NaNs haben; D1 und D2 bisher nur 1-dimensional
%
% lag: Werte fuer lag Eingabe (Verschiebung um soundsoviel Punkte)
%      gibt die Verschiebung von D2 an (z.B. bei monatlicher
%      Aufloesung gibt ein Wert von +12 an, dass bei Cor=1 ein
%      Signal bei D1 genau ein Jahr nachher bei D2 auftritt)
% Als variable Argumente können angegeben werden:
%  "plot"       : dann wird ein plot erstellt
%  "run",width  : die Korrelation wird jeweils über Teilstücke der
%                 Länge width durchgeführt. Ausgabe ist dann z.B.
%                 die zeitlich veränderliche Korrelation zweier
%                 Zeitreihen
%  "runstep",i  : hiermit kann festgelegt werden, dass die running
%                 Korrelation nicht bei jeden Zyklus/Zeitschritt
%                 sondern nur jeden i'ten durchgeführt wird
%  "minlength",p: Normalerweise werden NaN gesetzt, wenn D2 weniger
%                 als 70% Daten als D1 besitzt, bzw. in D1 weniger
%                 als 70% von der running width Daten hat. Hiermit
%                 kann man einen andern Prozentteil angeben. Bei 0
%                 werden NaN nur gesetzt falls keine Daten da, bei
%                 100 muss D2 mindestens gleich lang wie D1 sein
%                 (bei schon entfernten NaNs). bzw. immer running
%                 width Daten in D1 vorhanden sein.
%  "startshift",i: gibt an dass die Zeitreihe D2 erst bei den i'ten
%                 Werte des Zeitreihe D1 anfangt (bei negativen
%                 entsprechend andersrum)
%
%  Ausgabe:
%  COR  : Korrelationskoeffizent als Funktion der Zeit
%         (1-dimension) und der Verschiebung. Wenn keine running
%         correlation berechnet wird ist die erste Dimension=1
%  CORp : wie COR, nur zeiter Ausgabeparameter von corrcoef
%
% COPYRIGHT J.Holfort; under GNU public license

PLOTTEN=0; irun=0; minlength=0.7;runstep=1;

iD1=size(D1);
if length(iD1)>2; D1=squeeze(D1); iD1=size(D1); end
if length(iD1)>2; error('just one-D data allowed'); end
if iD1(1)>1 & iD1(2)>1; error('just one-D data allowed'); end 
iD2=size(D2);
if length(iD2)>2; D2=squeeze(D2); iD2=size(D2); end
if length(iD2)>2; error('just one-D data allowed'); end
if iD2(1)>1 & iD2(2)>1; error('just one-D data allowed'); end 
if iD2(1)==1 & iD1(1)~=1; D2=D2'; end

ii=1;
while ii<=nargin-3
    TT=varargin{ii}; ii=ii+1;
    if ~ischar(TT); error('wrong variable argument list'); end
    switch lower(TT)
        case {'plot'}
            PLOTTEN=1;
        case {'run'}
            irun=1;
            if ii>nargin-3; error('you must give a running width');end
            nrun=varargin{ii}; ii=ii+1;
            if nrun<2; 
                warning(['length of running correlation width' ...
                        ' to small, set to 20']);
                nrun=20
            end
        case {'runstep'}
            if ii>nargin-3; error('with runstep you have to give a value');end
            runstep=round(varargin{ii}); ii=ii+1;
            if runstep<1; 
                warning(['step for running correlation to' ...
                        ' small, set to 1']);
                runstep=1
            end	    
        case {'minlength'}
            if ii>nargin-3; error('with minlength you have to give a value');end
            minlength=varargin{ii}/100.; ii=ii+1;
            if minlength<0 or minlength>1;
                warning('minimal data length outside 0 to 100 Percent,set to 70%');
                minlength=.7
            end
        case {'startshift'}
            if ii>nargin-3; error('with startshift you have to give a value');end
            idum=varargin{ii}-1; ii=ii+1;
            if idum>0
                if iD1(1)==1; D2=[zeros(1,idum)+NaN D2]; else
                    D2=[zeros(idum,1)+NaN;D2]; end
            elseif idum<0
                if iD1(1)==1; D1=[zeros(1,-idum)+NaN D1]; else
                    D1=[zeros(-idum,1)+NaN;D1]; end
            end
        otherwise
            fprintf([TT '\n']);
            error('unknown variable argument')
        end
end

iD1=length(D1);
iD2=length(D2);

if exist('lag')~=1 | isempty(lag)
    lag=max(iD1,iD2)/1.5;
    lag=[floor(-lag):ceil(lag)];
end

ir2=1;
if irun==1; 
  ir2=length(D1); 
else
  ir2=1; nrun=length(D1); 
end
COR =zeros(length(1:runstep:ir2),length(lag))+NaN;
if nargout>1 | PLOTTEN==1; 
    CORp=zeros(length(1:runstep:ir2),length(lag))+NaN;
end

irx=0;
for ir=1:runstep:ir2
  irx=irx+1;
  ii=find(~isnan(D1));
  % ii(ii<ir-nrun/2 | ii>ir+nrun/2)=[]; %Aenderung 8/2004 nur falls irun==1
  if irun==1 
     ii(ii<ir-nrun/2 | ii>ir+nrun/2)=[];
     if length(ii)<nrun*minlength; continue; end
  end
     
  for iversatz=1:length(lag)
      jj1=ii; 
      jj2=ii+lag(iversatz);
      kk=jj2>0 & jj2<iD2;
      jj1=jj1(kk); jj2=jj2(kk);
      kk=~isnan(D2(jj2)); jj1=jj1(kk); jj2=jj2(kk);
      if length(jj2)> length(ii)*minlength  % nur bei genuegender Datenabdeckung
          if std(D1(jj1))==0 | std(D2(jj2))==0
              if std(D1(jj1))==0 & std(D2(jj2))==0
                  a=ones(2,2); p=a*0;
              else
                  jj2=[]; continue
              end
          end
          if nargout==1 & ~PLOTTEN==1
             [a]=corrcoef(D1(jj1),D2(jj2)); 
             COR (irx,iversatz)=a(1,2);
          else
             [a,p]=corrcoef(D1(jj1),D2(jj2)); 
             COR (irx,iversatz)=a(1,2);
             CORp(irx,iversatz)=p(1,2);
          end
      end    
  end
end

if PLOTTEN==1
  if irun~=1
    figure
    plot(lag,COR(1,:))
    hold on
    h=CORp(1,:)<0.05;
    plot(lag(h),COR(1,h),'ro')
    xlabel('lag');ylabel('correlation')
  else
    F=jet(length(lag));
    figure; hold on
    for ii=1;length(lag)
        plot(COR(:,ii),'color',F(ii,:))
    end
    xlabel('cycle');ylabel('correlation')
  end  
end