function lyapunov(filename)
%This function will plot a set of data, measure it's 
%Fast Fourier Transfor, and find the Lyapunov and the 
%Correlation Dimension for the set.

% load in data file
    data = load(filename);
% calculate number of data points
    N = length(data);
% create another array for return map
    nextdata = zeros(N,1);
    for (i=1:N-1)
     nextdata(i)=data(i+1);
    end
     nextdata(N)=nextdata(N-1);

% now plot it (except last point)
    plot(data(1:N-1),nextdata(1:N-1),'.');
% give it a title
    title(['Return map; Data = ', filename ]);
%save jpeg
    name = [ filename,'_return.jpg'];
    print ('-djpeg(100)', name)

% Pause and wait for user input ...
    display('press any key to continue');
    pause;

% calculate the fast fourier transform
    fftcoeff = fft(data,N);
    t = 2:N;
% plot the absolute value fourier coefficients
    plot(t,abs(fftcoeff(2:N)));
% give it a title
    title(['Fourier Coefficients; Data = ', filename ]);
%save jpeg
    name = [ filename,'_fast_fourier_coeff.jpg'];
    print ('-djpeg(100)', name)

% Pause and wait for user input
    display('press any key to continue');
    pause;

    N2 = floor(N/2);
    N4 = floor(N/4);
% find mid point of orbit sequence
    k=N2;
% create space for exponents
    exponent = zeros(N4,1);
%   look at 1/4 of the points
    for (j=1:N4)
        % set distance initially
        d = abs(data(k+1)-data(k));
        index = k+1;
        for (i=2:N-1)
            % see if there is a closer point
            if (i ~= k) && (abs(data(i)-data(k)))<d
                % if so, store index and distance
                d = abs(data(i)-data(k));
                index = i;
            end
        end
        % write log of quotient as difference of logs to get better accuracy!
        if (data(k) ~= data(index)) && (data(k+1) ~= data(index+1))
            exponent(j) = log( abs(data(k+1)-data(index+1)))-log(abs(data(k)-data(index)));
        end
    % repeat with the next point
    k = k+1;
    end

% now plot the lyapunov exponents
    t = 1:N4;
    lyapunov = exponent(1:N4);
    exp_avg = 0.0;

% find the average value for lyapunov exponent
    for (i=1:N4)
        exp_avg = exp_avg + exponent(i);
    end
% plot the exponents, the average, and the baseline
    exp_avg = exp_avg/N4;
    plot(t,lyapunov,t,0,t,exp_avg);
    TITLE(['Lyapunov Exponents; Data = ',filename, '; exponent avg = ', num2str(exp_avg)]);
%save jpeg
    name = [ filename,'_lyp_exps.jpg'];
    print ('-djpeg(100)', name)

% Pause and wait for user input
    display('press any key to continue');
    pause;

display 'average value for lyapunov exponent is:'
    exp_avg

%Call separate function for an approximation of the Correlation Dimension 
    corrdim(filename);

end

function corrdim(filename)
%This function uses the given data set to find the average 
%correlation dimension as the average of the slope of the plot

% load in data file
data = load(filename);

%measure length of data sets and number of pairs of points.
N = length(data);
p = N*(N-1)/2;

%initialize matrix to store values of r and C(r), the Correlation function
C_r = zeros(8,1);
R = zeros (8,1);

for (k=1:8)
    % set distance 
    r = 2^(-k);
    %initialize c=0 for total number of pairs {w_1,w_2} with |w_1-w_2|<r
    c = 0;
    %Compare distances between each pair of points
    for (i=1:N-1)
        for (j=i+1:N)
          % find pairs of points with distance less than d
            if (abs(data(j)-data(i)))<r
             % if so, increment c
             c=c+1;
            end
        end   
    end
    %Calculate and store log2 of radius and Correlation Function C(r)=C_r 
    x(k)=log2(r);
    y(k)=log2(c/p);  
    
    % repeat for next power of 2
    k = k+1;    
end

for (i=1:8)
dim(i)=(y(i)/x(i));
end

C_dim=mean(dim);
plot (x,y,'.')
TITLE(['Correlation Dimension; Data = ',filename '; cordim(S) avg =', num2str(C_dim)]);


%save jpeg
name = [ filename,'_cordim.jpg'];
print ('-djpeg(100)', name);

display 'Correlation Dimension= ',C_dim

end



