function out=lyapunov(ptime,ppnts)
%
%  Program to calculate Lyanpunov Exponents
%

if isempty(ptime) % nargin < 2 || isempty(ptime)
    ptime=[1:length(ppnts(:,1))]';   
end
ptime=ptime-ptime(1);

%tdiff=ptime(2:end)-ptime(1:end-1);

eps=1e-8; % pad number. 
n=length(ppnts(:,1));
idx=nearpoint(ppnts(1:end-1,:));
idx1=[idx;n];
nppnts=ppnts(idx1,:);
sep=0;
for i=1:n-1
    sep=sep+log((eps+norm(ppnts(i+1,:)-ppnts(idx1(i)+1,:)))/(eps+norm(ppnts(i,:)-nppnts(i,:))));
end
out=sep/(ptime(end-1)-ptime(1))/log(2); % Lyapunov exponent, log base 2.

function out=nearpoint(point)
%
m=length(point(:,1));
out=zeros(m,1);

for j=1:m
    diff=repmat(point(j,:),m,1)-point;
    if length(point(1,:))>1
        %dist=sqrt(sum([diff.*diff]'))';
               dist=sqrt(sum(diff.*diff,2));
    else
        dist=abs(diff);
    end
    [~,bb]=sort(dist);
    if bb(1)==j
        out(j)=bb(2);
    else
        out(j)=bb(1);
    end
        
end
end
end