function [Jerk1, Jerk2, Jerk3,MaxAx1,MaxAx2,MaxAx3,MaxJerk1,MaxJerk2,MaxJerk3] = accelAnalysis(Ax1,Ax2,Ax3,Time,TimeS,PlotTime,YAxisRange);

% use as [Jerk1, Jerk2,Jerk3,MaxAx1,MaxAx2,MaxAx3,MaxJerk1,MaxJerk2,MaxJerk3] = accelAnalysis(Ax1,Ax2,Ax3,Time,TimeS,PlotTime,YAxisRange);
% Inputs:   Ax1 : acceleration values on axis 1 in g's (array of integer)
%           Ax2 : acceleration values on axis 2 in g's(array of integer)
%           Ax3 : acceleration values on axis 3 in g's(array of integer)
%           Time : time of each measurement in matlab time(array of integer)
%           TimeS : time of each measurement in seconds, from time of
%           first measurement(array of integer)
%           PlotTime : if set to 0, time in seconds in used for plots, if
%           set to 1, matlab time is used for plots.
%           XAxisRange : The range between the upper and the lower Y axis
%                   limits on the plots.
%
% Outputs:  Jerk1 : jerk values on axis 1 in g's (array of integer)
%           Jerk2 : jerk values on axis 2 in g's(array of integer)
%           Jerk3 : jerk values on axis 3 in g's(array of integer)
%           MaxAx1, MAxAx2,MaxAx3: These n by 4 arrays of integers
%           contains n largest acceleration values for each axis on 
%           column 1. Columns 2, 3 and 4 are the asscociated jerk,
%           time and matlab time. Default n is 20 but can easily be
%           modified in the code.
%           MaxJerk1, MaxJerk2,MaxJerk3: These n by 4 arrays of integers
%           contains n largest jerk values for each axis on column 1. 
%           Columns 2, 3 and 4 are the asscociated acceleration,
%           time and matlab time. Default n is 20 but can easily be
%           modified in the code.
%           
% This function processes data generated by the function 'accel'.
%
% Francois Cazenave, august 2006
%
%This function is to be used to process data generated by the function
%'accel'. The matlab workspace should contain the variable above.
%This program will generate the variables containing an approximation of
%the rate of change of acceleration in G/s.
%
% Example:
% [Acc1, Acc2, Acc3,TotalTime,TotalTimeS] = accel('44AB7180','44AB73D8',5,1,4);
% In the line above, 44AB7180 is the name on the first file and 44AB73D8 the
% name of the second file. 5 files will be processed, including 44AB7180.DAT
% and 44AB73D8.DAT
% X axis on plots will be time in matlab time.
% The Y axis range will be 4 g's.


% Calculate the sampling Frequency
Frequency = 1 / (TimeS(3)-TimeS(2));

% Create vector containing the absolute value of an approximation of the jerk (derivative of acceleration)
% The jerk vector are 1 unit shorter that acceleration vector
Jerk1=abs(diff(Ax1))*Frequency;
Jerk2=abs(diff(Ax2))*Frequency;
Jerk3=abs(diff(Ax3))*Frequency;

% Remove the last value so that time, acceleration and jerk vectors have the same size 
Ax1=Ax1(1:(length(Ax1)-1));
Ax2=Ax2(1:(length(Ax2)-1));
Ax3=Ax3(1:(length(Ax3)-1));
Time=Time(1:(length(Time)-1));
TimeS=TimeS(1:(length(TimeS)-1));

figure;

if PlotTime==0; % use time is seconds from first measurement for x axis
    % Plot the Jerk on 3 subplots
    % Plot jerk along axis 1
    subplot(3,1,1)
    plot(TimeS,Jerk1);
    title(['Jerk along 3 axis from ' datestr(Time(1)) ' to ' datestr(Time(end))]);
    ylim([0 YAxisRange]);
    ylabel('axis 1 (m/s3)');

    % Plot jerk along axis 2
    subplot(3,1,2)
    plot(TimeS,Jerk2);
    ylim([0 YAxisRange]);
    ylabel('axis 2 (m/s3)');

    % Plot jerk along axis 3
    subplot(3,1,3)
    plot(TimeS,Jerk3);
    ylim([0 YAxisRange]);
    xlabel('time (s)')
    ylabel('axis 3 (m/s3)');


elseif PlotTime==1; % use matlab time for x axis
    % Plot the Jerk on 3 subplots
    % Plot Jerk along axis 1
    subplot(3,1,1)
    plot(Time,Jerk1);
    datetick;
    title(['Jerk along 3 axis from ' datestr(Time(1)) ' to ' datestr(Time(end))]);
    ylim([0 YAxisRange]);
    ylabel('axis 1 (m/s3)');

    % Plot Jerk along axis 2
    subplot(3,1,2)
    plot(Time,Jerk2);
    datetick;
    ylim([0 YAxisRange]);
    ylabel('axis 2 (m/s3)');

    % Plot Jerk along axis 3
    subplot(3,1,3)
    plot(Time,Jerk3);
    datetick;
    ylim([0 YAxisRange]);
    xlabel('time (s)')
    ylabel('axis 3 (m/s3)');
end

% Sort Jerk and Acceleration in descending order
% First for acceleration
[SortedAccel1,IndexAccel1]=sort(abs(Ax1),'descend'); % use the absolute value for sorting
SortedAccel1=Ax1(IndexAccel1);  % use the real value for returning to user
[SortedAccel2,IndexAccel2]=sort(abs(Ax2),'descend');
SortedAccel2=Ax2(IndexAccel2);
[SortedAccel3,IndexAccel3]=sort(abs(Ax3),'descend');
SortedAccel3=Ax3(IndexAccel3);

% Then for jerk
[SortedJerk1,IndexJerk1]=sort(Jerk1,'descend');
[SortedJerk2,IndexJerk2]=sort(Jerk2,'descend');
[SortedJerk3,IndexJerk3]=sort(Jerk3,'descend');


%Create matrices which contain max values
v=20; % v is the number of maximum values that will be placed in the result arrays
% Max acceleration on axis 1
MaxAx1(:,1)=SortedAccel1;
MaxAx1(:,2)=Jerk1(IndexAccel1);
MaxAx1(:,3)=TimeS(IndexAccel1);
MaxAx1(:,4)=Time(IndexAccel1);
MaxAx1=MaxAx1(1:v,1:4);
% Max acceleration on axis 2
MaxAx2(:,1)=SortedAccel2;
MaxAx2(:,2)=Jerk2(IndexAccel2);
MaxAx2(:,3)=TimeS(IndexAccel2);
MaxAx2(:,4)=Time(IndexAccel2);
MaxAx2=MaxAx2(1:v,1:4);
% Max acceleration on axis 3
MaxAx3(:,1)=SortedAccel3;
MaxAx3(:,2)=Jerk3(IndexAccel3);
MaxAx3(:,3)=TimeS(IndexAccel3);
MaxAx3(:,4)=Time(IndexAccel3);
MaxAx3=MaxAx3(1:v,1:4);

% Max Jerk on axis 1
MaxJerk1(:,1)=SortedJerk1;
MaxJerk1(:,2)=Ax1(IndexJerk1);
MaxJerk1(:,3)=TimeS(IndexJerk1);
MaxJerk1(:,4)=Time(IndexJerk1);
MaxJerk1=MaxJerk1(1:v,1:4);
% Max Jerk on axis 2
MaxJerk2(:,1)=SortedJerk2;
MaxJerk2(:,2)=Ax2(IndexJerk2);
MaxJerk2(:,3)=TimeS(IndexJerk2);
MaxJerk2(:,4)=Time(IndexJerk2);
MaxJerk2=MaxJerk2(1:v,1:4);
% Max Jerk on axis 3
MaxJerk3(:,1)=SortedJerk3;
MaxJerk3(:,2)=Ax3(IndexJerk3);
MaxJerk3(:,3)=TimeS(IndexJerk3);
MaxJerk3(:,4)=Time(IndexJerk3);
MaxJerk3=MaxJerk3(1:v,1:4);

