function   varargout = ltpda_series_derivative(varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LTPDA_SERIES_DERIVATIVE performs the numeric time derivative
% using the method of five points stencil Taylor series expansion [1].
% The function can perform the first, second, third and fourth derivetive
% of a series of input data.
% 
%
%
% CALL:                 Deriv = ltpda_series_derivative(a, ... , pl)
%                      
%
% INPUTS:      
%           - Aos containing the data to be differentiated.
%           - pl: is a parameter list containing the options for the
%           calculation. Parameters are:
%           - 'ORDER' whose accepeted options are:
%               - 'FIRST' perform the first derivative using the
%               couefficients vector d1 = [1 -8 0 8 -1]./(12*T)
%               - 'SECOND' perform the second derivative using the
%               coefficients vector d2 = [-1 16 -30 16 -1]./(12*T^2)
%               - 'THIRD' perform the third derivative using the
%               coefficients vector d3 = [-1 2 0 -2 1]./(2*T^3)
%               - 'FOURTH' perform the third derivative using the
%               coefficients vector d4 = [1 -4 6 -4 1]./(T^4)
%                
%           NOTE1: T is the sampling period
%           NOTE2: The default option for 'ORDER' is 'SECOND'
%              
%
%
% OUTPUTS:     
%           - Deriv is avector containing the resulting Aos after
%           differentiation
% 
% 
% REFERENCES:
%           [1] S. E. Koonin and D. C. Meredith, COMPUTATIONAL PHYSICS -
%           Fortran Version, 1990, Westview Press
%
%
% VERSION:     $Id: ltpda_series_derivative.m,v 1.1 2008/04/24 16:13:12 luigi Exp $
% 
% HISTORY:     18-03-2008 L Ferraioli
%                 Creation
%
% The following call returns a parameter list object that contains the
% default parameter values:
%
% >> pl = ltpda_series_derivative('Params')
%
% The following call returns a string that contains the routine CVS version:
%
% >> version = ltpda_series_derivative('Version')
%
% The following call returns a string that contains the routine category:
%
% >> category = ltpda_series_derivative('Category')
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% Standard history variables

ALGONAME = mfilename;
VERSION  = '$Id: ltpda_series_derivative.m,v 1.1 2008/04/24 16:13:12 luigi Exp $';
CATEGORY = 'Signal Processing';

%% Check if this is a call for parameters, the CVS version string
% or the method category
if nargin == 1 && ischar(varargin{1})
  in = char(varargin{1});
  if strcmpi(in, 'Params')
    varargout{1} = getDefaultPL();
    return
  elseif strcmpi(in, 'Version')
    varargout{1} = VERSION;
    return
  elseif strcmpi(in, 'Category')
    varargout{1} = CATEGORY;
    return
  end
end

%% Collect input ao's, plist's and ao variable names

in_names = {};

for ii = 1:nargin
    in_names{end+1} = inputname(ii); % Collect the inputs names
end

[as, pl, invars] = collect_inputs(varargin, in_names);

% produce one parameter list
if isa(pl, 'plist')
  pl = combine(pl, getDefaultPL());
else
  pl = getDefaultPL();
end

%% Initialize outputs
bs = [];

%% Go through analysis objects

% Assigning coefficients values based on the input options
switch find(pl, 'ORDER')
    case 'FIRST'
        Coeffs = [1 -8 0 8 -1]./12;
    case 'SECOND'
        Coeffs = [-1 16 -30 16 -1]./12;
    case 'THIRD'
        Coeffs = [-1 2 0 -2 1]./2;
    case 'FOURTH'
        Coeffs = [1 -4 6 -4 1];     
end

for jj = 1:numel(as)
    
    a = as(jj);
    phi = a.data.y;
    fs = get(a, 'fs');
    T = 1/fs;
    
    %% Building vectors for calculation
    phi_temp = [phi(1);phi(1);phi(1);phi(1);phi;phi(end);phi(end);phi(end);phi(end)];

    % Switching between the input options
    switch find(pl, 'ORDER')
        case 'FIRST'
            Deriv = (1/T).*(Coeffs(1)*phi_temp(1:end-4) + Coeffs(2)*phi_temp(2:end-3) + Coeffs(3)*phi_temp(3:end-2) + Coeffs(4)*phi_temp(4:end-1) + Coeffs(5)*phi_temp(5:end));
            Deriv = Deriv(3:end-2);
        case 'SECOND'
            Deriv = (1/T^2).*(Coeffs(1)*phi_temp(1:end-4) + Coeffs(2)*phi_temp(2:end-3) + Coeffs(3)*phi_temp(3:end-2) + Coeffs(4)*phi_temp(4:end-1) + Coeffs(5)*phi_temp(5:end));
            Deriv = Deriv(3:end-2);
        case 'THIRD'
            Deriv = (1/T^3).*(Coeffs(1)*phi_temp(1:end-4) + Coeffs(2)*phi_temp(2:end-3) + Coeffs(3)*phi_temp(3:end-2) + Coeffs(4)*phi_temp(4:end-1) + Coeffs(5)*phi_temp(5:end));
            Deriv = Deriv(3:end-2);
        case 'FOURTH'
            Deriv = (1/T^4).*(Coeffs(1)*phi_temp(1:end-4) + Coeffs(2)*phi_temp(2:end-3) + Coeffs(3)*phi_temp(3:end-2) + Coeffs(4)*phi_temp(4:end-1) + Coeffs(5)*phi_temp(5:end));
            Deriv = Deriv(3:end-2);
    end
    
    %% Building output Aos
    
    % Creating new time series data
    Deriv = tsdata(Deriv,fs);

    Deriv = set(Deriv, 'name', ...
    sprintf('ltpda_parfit_derivative(%s)', a.data.name), ...
    'xunits', a.data.xunits, 'yunits', 'm/s^2');

    % make new history objects
    h = history(ALGONAME, VERSION, pl, a.hist);
    h = set(h, 'invars', invars);

    % make output analysis objects
    b = ao(Deriv, h);
    clear Deriv

    % name, mfilename, description for these objects
    b = setnh(b, 'name', sprintf('ltpda_parfit_derivative(%s)', invars{jj}),...
       'description', find(pl,'description'));

    %% Output the data

    bs = [bs b];
    clear a b

end

%% Reshape the ouput to the same size of the input
varargout{1} = reshape(bs, size(as));

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION:    getDefaultPL
%
% DESCRIPTION: Get default params
%
% HISTORY:     01-02-2008 M Hueller
%                 Creation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plo = getDefaultPL()

plo = plist('NAME', 'ltpda_parfit_derivative', 'ORDER', 'SECOND', ...
    'DESCRIPTION','Series Expansion Derivative Method');


% END