% PFRAC
%
%   Function PFRAC returns constants of partial fractions.
%
%   If the roots of input denominator DEN have no multiple roots, then
%   CONSTANTS represents following partial fractions:
%
%   NUM(s)       CONST_1                     CONST_N
%   ------ = ----------------- + ... +  -----------------
%   DEN(s)   (s - ROOTS_DEN_1)          (s - ROOTS_DEN_N)
%
%   where
%   N     - is number of roots
%
%   If the roots of input denominator DEN have multiple roots, then
%   CONSTANTS represents following partial fractions:
%
%   NUM(s)       CONST_1                     CONST_M
%   ------ = ----------------- + ... +  -----------------
%   DEN(s)   (s - ROOTS_DEN)^1          (s - ROOTS_DEN)^M
%
%   where
%   M     - is number of multiple roots
%
%   If the roots of input denominator DEN have roots with non-zero
%   imaginary part, then CONSTANTS represents following partial
%   fractions: 
%
%   NUM(s)   CONST_1 + CONST_2*s
%   ------ = -----------------------
%   DEN(s)   s^2 + COEF_1*s + COEF_2
%
%   Type:
%
%   [const, roots_den, const_matrix] = pfrac(num, den, tol_diff, tol_imag))
%
%   where
%
%   const           - is output vector of required constants of partial
%                     fractions
%   roots_den       - is output vector of roots of denominator DEN
%                     associated with output vector CONST
%   const_matrix    - is output matrix for calculating the coeficients:
%                     CONST = INV(CONST_MATRIX)*NUM_WITH_CORRECT_DEGREE'
%   num             - is input vector of numerator of input fraction
%   den             - is input vector of denominator of input fraction
%   tol_diff        - is obtinal tolarance of "different" roots (DEFAULT = 1e-4)
%   tol_imag        - is obtinal tolarance of "imaginary" roots (DEFAULT = 1e-4)
%
%   juraj.oravec@stuba.sk
%
%   since   2011.04.07.
%   rev.    2013.02.25.
%
function [const, rts, coef_mtx] = pfrac(num, den,tol_diff,tol_imag)


% ToDo:
%
% 1) Multiple Imaginary Roots:
%    [const, rts, coef_mtx] = pfrac(2,poly([1-i,1+i,1-i,1+i]))
%
% 2) Roots with only Imaginary part:
%    [const, rts, coef_mtx] = pfrac([1 3 10 27],[1,0,13,0,36]))



% ------------------------------------------ %
%
% Input parse
%
% ------------------------------------------ %


% Check type of inputs
%
% Temporary variables of input class
%
class_num = class(num);
class_den = class(den);
%
if (nargin < 2)|((isempty(num)==1 | isempty(den)==1)) | (class_num(1) ~='d') | (class_den(1) ~='d')
    error(' PFRAC: No empty two input vectors of numbers required!')
end

% Check obtional inputs
%
if((nargin > 2)&(isempty(tol_diff) == 1)) | (nargin == 2)

    % Tolerance of "fake" different roots
    %
    tol_diff = 1e-4;
        
end % if

if(nargin < 4)
    
    % Tolerance of "fake" zero imaginary part
    %
    tol_imag = 1e-4;

end % if
tol_diff; % Verbose
tol_imag; % Verbose


% Ensure column-vector sizes of inputs
%
[num,den] = force_column_vector(num,den);
den = den';


% ------------------------------------------ %
%
% Preliminary Stage of the roots processing
%
% ------------------------------------------ %


% Get roots
%
rts_raw = roots(den);


% Remove "fake" imaginary part
%
for k = 1 : length(rts_raw)
    if(abs(imag(rts_raw(k))) <= tol_imag)
        rts_raw_noimag(k,1) = real(rts_raw(k));
    else
        rts_raw_noimag(k,1) = rts_raw(k);
    end % if
end % for k
rts_raw_noimag; % Verbose


% Remove "fake" different roots
%
for k = 1 : length(rts_raw_noimag )
    rts_raw_diff(k,1) = round(rts_raw_noimag(k)*(1/tol_diff))*tol_diff;
end % for k
rts_raw_diff; % Verbose


% Sort "clean" roots
%
rts = sort(rts_raw_diff); % Verbose

% Count unique roots
%
rts_unique = sort(unique(rts));
for k = 1 : length(rts_unique)
    rts_cnt(k,1) = sum(rts_unique(k) == rts);
end % for k
rts_cnt; % Verbose


% Count multiple roots
%
cnt = 0;
for k = 1 : length(rts_unique)
    for j = 1 : rts_cnt(k)
        cnt = cnt + 1;
        rts_idx(cnt,1) = j;
    end % for j
end % for k
rts_idx; % Verbose


% ------------------------------------------ %
%
% Generate polynomials of partial fractions
%
% ------------------------------------------ %

rts; % Verbose
cnt = 0;
rts_idx = [];
for k = 1 : length(rts_unique)
    
    for j = 1 : rts_cnt(k)
        
        % Find first J indexes of current root to be removed
        %
        rts_idx = find(rts_unique(k) == rts);
        rts2remove = rts_idx(1:j);

        if(imag(rts_unique(k)) < 0)
            
            % Find first J indexes of conjugated current root to be removed
            %
            rts_unique_conj = conj(rts_unique(k));
            rts_idx_conj = find(rts_unique_conj == rts);
            rts2remove_conj = rts_idx_conj(1:j);
            rts2remove = [rts2remove; rts_idx_conj(1:j)];
            
        end % if
        
        % Find indexes of roots that generate the coefficient polynomial
        %
        rts_pol_idx = index_not(rts2remove, [1:length(rts)]);
        rts_pol = rts(rts_pol_idx);
        
        if(imag(rts_unique(k)) == 0)

            % Generate polynomial for current root with no imaginary part
            %
            cnt = cnt + 1;
            coef_pol{cnt,1} = poly(rts_pol);
            
        end
        
        if(imag(rts_unique(k)) < 0)

            % Generate polynomials for current root with imaginary part
            %
            cnt = cnt + 1;
            coef_pol{cnt,1} = poly([rts_pol]);
            
            cnt = cnt + 1;
            coef_pol{cnt,1} = poly([rts_pol; 0]);
        
        end % if
        
    end % for j
    
end % for k


% Get maximal polynomial degree
%
for k = 1 : length(coef_pol)
    coef_deg(k,1) = length(coef_pol{k});
end % for k
coef_deg_max = max(coef_deg);


% Ensure same degree
%
for k = 1 : length(coef_pol)
    coef_pol_deg{k,1} = [zeros(1,coef_deg_max - coef_deg(k)),coef_pol{k}];
end % for k
num_deg = [zeros(1,coef_deg_max - length(num)),num];
num_deg = force_column_vector(num_deg);


% ------------------------------------------ %
%
% Compute the coefficients of partial fractions
%
% ------------------------------------------ %


% Coefficients' matrix
%
for k = 1 : length(coef_pol)
    coef_mtx(:,k) = coef_pol_deg{k,1}';
end % for k
coef_mtx; % Verbose
num; % Verbose

const = inv(coef_mtx)*num_deg;

[rts, const]; % Verbose


end % function


% INDEX_NOT
%
%   Function INDEX_NOT remove required elements from input array of type
%   DOUBLE or CELL.
%
%   idx_not = index_not(idx,idx_all)
%
%   where:
%
%   idx_not     - is output array of same type, as input IDX_ALL with removed
%                 elements listed in input IDX
%   idx         - is input list of elements to be removed from array IDX_ALL
%   idx_all     - is input arary of elements to be filtered
%
%   Example #1:
%
%   Type:
%   idx_not = index_not([10 20 30 40],[20 40 60])
%
%   Results:
%   idx_not =
% 
%        60
%
%   Example #2:
%
%   Type:
%   idx_all{1} = 10; idx_all{2} = 20; idx_all{3} = 30;
%   idx_not = index_not({20},idx_all)
%
%   Results:
%   idx_not = 
% 
%       [10]
%       [30]
%
%   juraj.oravec@stuba.sk
%
%   2012.07.17.
%
function idx_not = index_not(idx,idx_all)

if(nargin < 2)
   error(sprintf('INDEX_NOT:ERROR: Too inputs required!\n\n idx_not = index_not(idx,idx_all)'))
end

if(nargin > 2)
   error(sprintf('INDEX_NOT:ERROR: Too many inputs!\n\n idx_not = index_not(idx,idx_all)'))
end

n_idx = length(idx);
n_idx_all = length(idx_all);

% if(n_idx_all == 0)
%    error('INDEX_NEG:ERROR: Empty input "IDX_ALL"!')
% end %if

% if(n_idx_all < n_idx )
%     error('INDEX_NEG:ERROR: Too much elements in input "IDX" according to input "IDX_ALL"!')
% end % if

idx_not = []; % Initial value of output
cnt = 0;

if(iscell(idx) == 1)

    % For arrays class CELL
    %
    for k = 1 : n_idx_all

        chk = 0;

        for n = 1 : n_idx

            if(isequal(idx_all{k}, idx{n}))
                chk = 1;
            end

        end % for n

        if(chk == 0)
            cnt = cnt + 1;
            idx_not{cnt,1} = idx_all{k};
        end % if

    end % for k

    
else
    
    % For arrays class DOUBLE
    %
    for k = 1 : n_idx_all

        chk = 0;

        for n = 1 : n_idx

            if(idx_all(k) == idx(n))
                chk = 1;
            end

        end % for n

        if(chk == 0)
            cnt = cnt + 1;
            idx_not(cnt,1) = idx_all(k);
        end % if

    end % for k

end % if

end % function


% FORCE_COLUMN_VECTOR
%
%   Function FORCE_COLUMN_VECTOR returnes the column vectors.
%
%   [v_out_1,...,v_out_N] = force_column_vector(v_in_1,...,v_in_N)
%
%   v_in_1,...,v_in_N       - input row or column vectors
%   v_out_1,...,v_out_N     - output column vectors
%
%   Example:
%
%   [v1,v2,v3] = force_column_vector([1,1,1],[2;2],[1,2,3,4,5])
%
%   juraj.oravec@stuba.sk
%
%   2013.01.15.
%
function varargout = force_column_vector(varargin)

% Check if there is output required
%
if (nargout == 0)
    return
else
    %
    % Number fo required outputs
    %
    n_out = nargout;
end % if

% Check if there is input inserted
%
if (nargin == 0)
    %
    % Ensure emtpy outputs
    %
    for k = 1 : n_out
        varargout{k} = [];
    end % for k
    return
else
    %
    % Number of input vectors
    %
    n_in = nargin;
end % if

% Check if number of required outputs is equal to the number of inserted
% inputs
%
if (nargin ~= nargout)
    error('FORCE_COLUMN_VECTOR: Number of input vectors must be the same as the number of output vectors!')
end % if

% Check if there is no matrix on the input
%
for k = 1 : n_in
    %
    % Matrix
    %
    if(size(varargin{k},1) > 1) & (size(varargin{k},2) > 1)
        error('FORCE_COLUMN_VECTOR: Matrix found! Required inputs are only the row/column vector!')
    end % if
    
end % for k


% Ensure column vector
%
for k = 1 : n_in
    
    if(size(varargin{k},2) > 1)
        varargout{k} = varargin{k}';
    else
        varargout{k} = varargin{k};
    end % if
    
end % for k

end % function
