function [cdf] = InvCharacteristic(poles,thr,eps)

% Description: InvCharacteristic(poles,thr,eps) inverts a general 
%              charateristic function which consists of both simple 
%              and multiple poles.
%
% Inputs:  poles - A vector containing all poles of the characteristic 
%                  function including both left plane and right plane 
%                  poles. If the cdf of the decision variable is required, 
%                  the pole at zero should be included in the vector.
%          thr   - The threshold such that the cdf of the decision variable,
%                  f, is given by P(f<thr).
%          eps   - This gives the accuracy of an exponential which is represented
%                  using the Maclaurin series expansion. If eps = [], the default
%                  value of 1e-4 is used. Note that when thr = 0, eps is not used
%                  in the calculations, and the solution is exact.
% Outputs: cdf   - The is the numerical result for P(f<thr)

pl=poles(find(poles<=0));     % Left plane (LP) poles
pr=poles(find(poles>0));      % Right plane (RP) poles
A = prod(poles(find(poles))); % The multiplication factor.

% If thr>0, the LP and RP poles are swapped since the residues will be
% calculated over the LP poles.
if thr>0
   t=pl;
   pl=pr;
   pr=t;
end
pr=sort(pr);

% If eps = [], use the default value of 1e-4
if isempty(eps) == 1
   eps=1e-4;
end

nRPpoles = length(pr);  % The number of RP poles
nLPpoles = length(pl);  % The number of LP poles

% The RP poles are represented as pr(i) = porg + pdiff(i)
porg = min(pr)/2;
pdiff = pr - porg;

% P1 represents polynomial #1;
P1=1;
for c1=1:nLPpoles
   P1 = conv(P1,[porg-pl(c1),1]);
end
LP1=length(P1); % Length of polynomial #1

% The Maclaurin series expansion is used to create P2, polynomial #2.
% Note that exp(+j) converges faster than exp(-j). Therefore, if
% thr < 0, q is used to shift the exponential so that the
% Maclaurin series expansion is used to represent exp(thr*pdiff(i)+q)
% instead of exp(thr*pdiff(i)). This shift is reversed when P2 is later 
% multiplied by exp(-q).
if thr < 0
   q=abs(max(abs(pdiff))*thr);
else
   q = 0;
end
w=1;
P2=1;
LP2=1;
Ptmp=1;
[th_cmp,thi]=max(q+thr.*pdiff);
Pdtmat = tril(ones(LP2,LP2),0)+triu(pdiff(thi).*ones(LP2,LP2),1);
if thr~=0
   while ((abs(1 - sum(prod([P2.*Pdtmat(1,1:LP2);Pdtmat(2:LP2,1:LP2)],1))/ ...
         exp(th_cmp)) > eps) | (w<LP1) | (w<nRPpoles))
      Ptmp = conv(Ptmp,[q/thr,1]);
      P2 = [P2,0] + prod(repmat(thr,1,w)./[1:w]).*Ptmp;
      LP2 = length(P2);
      Pdtmat = tril(ones(LP2,LP2),0)+triu(pdiff(thi).*ones(LP2,LP2),1);
      w=w+1;
   end
end

% Normalize polynomial #2 to ensure numerical stability and reverse
% the shift of exp(q);
P2 = A*prod(1./(pr(nRPpoles)-pl))*exp(porg*thr-q).*P2;

LP2=length(P2); % Length of polynomial #2

% LP is the maximum length of polynomials #1 or #2. The length,
% LP, must be at least nRPpoles since the procedure requires nRPpoles-1
% iterations to invert a characteristic function with nRPpoles RP poles.
LP = max(max(LP1,LP2),nRPpoles);

for cnt=2:length(pdiff)
   % Normalize polynomial #1 to ensure numerical stability.
   P1t=P1./prod(pr(cnt-1)-pl);
   
   % Extend polynomials #1 and #2 by the required amount.
   P1t=[P1t,zeros(1,LP-LP1)];
   P2t=[P2,zeros(1,LP-LP2)];

   % The multiplicative factors for the upper and lower vectors.
   Pdtmat = tril(ones(LP,LP),0)+triu(pdiff(cnt-1).*ones(LP,LP),1);
   
   % The vector, P3, contains the coefficients for P2 in the following
   % iteration. This uses the upper and lower vectors of P1t and P2t.
   P3=zeros(1,LP-1);
   for c1=LP-1:-1:1
      P3(c1) = sum(sum((((prod([P1t(1:c1).*Pdtmat(1,1:c1);Pdtmat(2:LP,1:c1)],1)).')* ...
          (prod([P2t(c1+1:LP).*Pdtmat(1,1:LP-c1);Pdtmat(2:LP,1:LP-c1)],1))- ...
          ((prod([P2t(1:c1).*Pdtmat(1,1:c1);Pdtmat(2:LP,1:c1)],1)).')* ...
          (prod([P1t(c1+1:LP).*Pdtmat(1,1:LP-c1);Pdtmat(2:LP,1:LP-c1)],1)))));
   end
   
   P2=P3;
   LP2=length(P2);
   LP = max(LP1,LP2);
end

% The numerical result for P(f<thr)
Pdtmat = tril(ones(LP2,LP2),0)+triu(pdiff(length(pdiff)).*ones(LP2,LP2),1);
cdf = abs(sum(prod([P2.*Pdtmat(1,1:LP2);Pdtmat(2:LP2,1:LP2)],1)));

