function   varargout = ltpda_parfit_derivative(varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LTPDA_PARFIT_DERIVATIVE performs the numeric time derivative using the method
% of 5 points parabolic fit [1]. The function can perform the second
% derivative, the first derivative and the point estimation. 
% 
%
%
% CALL:                 Deriv = ltpda_parfit_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:
%               - 'ZERO' perform the point estimation using the
%               coefficients vector d0 = [-3 12 17 12 -3]./35
%               - 'FIRST' perform the first derivative using the
%               couefficients vector d1 = [-2 -1 0 1 2]./(10*T)
%               - 'SECOND' perform the second derivative using the
%               coefficients vector d2 = [2 -1 -2 -1 2]./(7*T^2)
%           
%           NOTE1: T is the sampling period
%           NOTE2: The default option for 'ORDER' is 'SECOND'
%              
%
%
% OUTPUTS:     
%           - Deriv is the resulting Ao after the operation defined by the
%           input parameters
%
% REFERENCES:
%           [1] L. Carbone et al., Physical Review D 75, 042001 (2007)
%
% VERSION:     $Id: ltpda_parfit_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_parfit_derivative('Params')
%
% The following call returns a string that contains the routine CVS version:
%
% >> version = ltpda_parfit_derivative('Version')
%
% The following call returns a string that contains the routine category:
%
% >> category = ltpda_parfit_derivative('Category')
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% Standard history variables

ALGONAME = mfilename;
VERSION  = '$Id: ltpda_parfit_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 'ZERO'
        Coeffs = [-3 12 17 12 -3]./35;
    case 'FIRST'
        Coeffs = [-2 -1 0 1 2]./10;
    case 'SECOND'
        Coeffs = [2 -1 -2 -1 2]./7;
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 'ZERO'
            Deriv = (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 '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);
    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','Parabolic Fit Derivative Method');


% END