function [Index,fIndex] = ifpofp(D,IndexType,tol,fBand)
% Compute IFP or OFP index for state-space model G(s)

%  Copyright 1986-2017 The MathWorks, Inc.

% Default values
Index = -Inf;
fIndex = NaN;
Ts = abs(D.Ts);

% Limitations
if hasdelay(D)
   if Ts==0
      error(message('Control:analysis:getSectorIndex9'))
   else
      D = elimDelay(D,D.Delay.Input,D.Delay.Output,D.Delay.Internal);
   end
end

% Simplify and scale realization
D = sminreal(D);
[a,b,c,d,~,e] = getABCDE(D);
if ~D.Scaled
   % Rescale for freq. response accuracy
   [a,b,c,e] = ltipack.xscale(a,b,c,d,e,Ts,'Warn',false);
end

% Handle descriptor case
if ~isempty(e)
   Ds = ltipack.ssdata(a,b,c,d,e,0);  % D may contain residual internal delays
   Ds.Scaled = true;
   [isP,Ds] = isproper(Ds);
   if isP
      % System is proper. Use reduced model with nonsingular E
      [a,b,c,d,~,e] = getABCDE(Ds);
   else
      error(message('Control:analysis:getSectorIndex10'))
   end
end

% Compute poles and test frequencies
if isempty(e)
   p = eig(a);
else
   p = eig(a,e);
end
if Ts>0
   % Compute equivalent continuous-time poles
   p = log(p(p~=0))/Ts;
end
fTest = unique(abs(p));
fTest = [-fTest ; fTest];  % for complex data

% Set up directional index computation
[nx,nu] = size(b);
switch lower(IndexType(1))
   case 'i'
      % IFP
      % Check I+delta*G has no RHP zeros for small delta
      % Note: Pick delta to correspond to t=-1e8 for scaled G
      Gscale = norm(d,1) + norm(b,1)^2/(1+norm(a,1));
      if Gscale==0
         z = ltipack.sspole(a,e);
      else
         delta = sqrt(eps)/Gscale;   aux = sqrt(delta);
         z = ltipack.tzero(a,aux*b,aux*c,eye(nu)+delta*d,e);
      end
      if (Ts==0 && any(real(z)>=0)) || (Ts>0 && any(abs(z)>=1))
         return
      end
      % Compute index from FDI
      % Note: Custom code needed to correctly handle pure integrators
      [Index,fIndex] = ltipack.ifpConicIndex(a,b,c,d,e,Ts,tol,-Inf,fTest,fBand);

   case 'o'
      % OFP
      % D must be nonsingular for the FDI to hold strictly at s=Inf
      % Note: tFD>=tD always holds since W'*H has a causal, stable inverse
      %       and ||X||oo >= ||X(oo)||
      tD = ltipack.getTmax(-(d+d'),d'*d)/2;
      if tD==-Inf
         return
      end
      % Check G+delta*I has no RHP zeros for small delta
      delta = sqrt(eps)*norm(d,1);
      z = ltipack.tzero(a,b,c,d+delta*eye(nu),e);
      if (Ts==0 && any(real(z)>=0)) || (Ts>0 && any(abs(z)>=1))
         return
      end
      % Compute index from FDI
      Q = [zeros(nu),-eye(nu)/2;-eye(nu)/2,zeros(nu)];
      W = [eye(nu);zeros(nu)];
      c = [c ; zeros(nu,nx)];
      d = [d ; eye(nu)];
      [Index,fIndex] = ltipack.dirConicIndex(a,b,c,d,e,Ts,Q,W,tol,-Inf,fTest,fBand);
end
