function Divergence=div_calc(state,ws_sec,fs,period_sec,progress)
% calculate divergence curve (needed for max lyapunov exponent). 
% Input: 
%   state: state space
%   ws_sec: window size over which divergence will be calculated(in seconds)
%   fs: sample frequency 
%   period_sec: indicates period of time-wise near samples to exclude in
%   nearest neighbour search
%   progress: show progress or not
% 
% Output: 
%   Divergence: the divergence curve
%
% History:
%   August 2011: v1, Sietse Rispens based on version KvS en SMB
%   1 November 2012, Sietse Rispens: Use div_calc_shorttimeseries for short
%   time series

if size(state,1) <= 10000
    Divergence=div_calc_shorttimeseries(state,ws_sec,fs,period_sec,progress);
    return;
end

[N,D]=size(state);
ws=round(ws_sec*fs);
Period=round(period_sec*fs);
if N<ws, error('ws shall not be larger than N'); end
NCompleteWindow = N-ws+1;

Divergence_sum=zeros(1,ws);
Divergence_count=zeros(1,ws);

SumStd = sqrt(sum(std(state,1).^2));
DistLimStart = SumStd*nthroot(1/NCompleteWindow,D);
DistLim = DistLimStart;

TreeRoot=kdtree_build(state(1:NCompleteWindow,:));
k0 = 2; % The initial number of nearest neighbors to look for
kmax = min(2*Period + 2, NCompleteWindow); % The maximum number of nearest neighbors to look for 
for i = 1:NCompleteWindow 
    StateI = state(i,:);
    if ~isnan(StateI)
        DistLim = DistLim*nthroot(1/5,D);
        DistLim = min(DistLim,DistLimStart);
        Index = [];
        k=k0;
        kmaxreached = 0;
        while isempty(Index) && ~kmaxreached
            [Idx, Dist] = kdtree_k_nearest_neighbors(TreeRoot,StateI,k);
            Dist(abs(i-Idx) < Period)=nan;
            Index = Idx(find(Dist==min(Dist),1,'first'));
            if k >= kmax
                kmaxreached = 1;
            elseif isempty(Index)
                k = min(kmax,k*2);
            end
        end
        if ~isempty(Index)
            if i+ws>N || Index+ws>N
                % do not use these data
            else
                if ~isempty(Index)
                    DistCurve = log(sqrt(sum((state(i:i+ws-1,:)-state(Index:Index+ws-1,:)).^2,2)));
                    NotNan = ~isnan(DistCurve);
                    Divergence_sum(NotNan) = Divergence_sum(NotNan) + DistCurve(NotNan)';
                    Divergence_count(NotNan) = Divergence_count(NotNan) + 1;
                end
            end
        end
    end
    if progress > 0
        if mod(i,round(progress))==0
            if i>round(progress)
                fprintf('\b\b\b\b\b\b\b\b\b\b');
            end
            fprintf('i=%8d', i);
        end
    end
end
kdtree_delete(TreeRoot); % Free the pointer to k-d-tree
if progress > 0
    fprintf('\n');
end
Divergence= (Divergence_sum./Divergence_count);