classdef (Hidden = true) ltidata 
   % Class definition for @ltidata

   %   Author(s): P. Gahinet
   %   Copyright 1986-2022 The MathWorks, Inc.
   properties
      Delay  % delay structure
      Ts     % sampling time
   end

   % Methods
   methods
      
      %% Characteristics
      function boo = issparse(~)
         boo = false;
      end
      
      function pf = hasSectorResp(D,W1,W2,Focus)
         % Checks if system has an R-index plot
         D = D.ss();
         checkCompatibility(D,'sector');
         [a,b,c,d,~,e] = getABCDE(D);
         Ts = D.Ts; %#ok<*PROPLC>
         rD = ltipack.util.gsvmax(W1'*d,W2'*d,0);
         if isfinite(rD)
            if ~D.Scaled
               % Scale before calling TZERO
               [a,b,c,e] = ltipack.xscale(a,b,c,d,e,Ts,'Focus',Focus,'Warn',false);
            end
            z = ltipack.tzero(a,b,W2'*c,W2'*d,e);
            pf = (Ts==0 && all(real(z)<0)) || (Ts~=0 && all(abs(z)<1));
         else
            pf = false;
         end
      end

      function pf = hasIFPOFPResp(D,type,Focus)
         % Checks if system has an IFP or OFP plot
         D = D.ss();
         checkCompatibility(D,'sector');
         [a,b,c,d,~,e] = getABCDE(D);
         Ts = D.Ts;
         if ~D.Scaled
            % Scale before calling TZERO
            [a,b,c,e] = ltipack.xscale(a,b,c,d,e,Ts,'Focus',Focus,'Warn',false);
         end
         nu = size(b,2);
         if type==1
            % IFP
            % Check I+delta*G has no RHP zeros for small delta
            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
            pf = (Ts==0 && all(real(z)<0)) || (Ts~=0 && all(abs(z)<1));
         else
            % OFP
            tD = ltipack.getOFP(d);
            if tD==-Inf
               pf = false;
            else
               % 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);
               pf = (Ts==0 && all(real(z)<0)) || (Ts~=0 && all(abs(z)<1));
            end
         end
      end

      %% Conversions
      function GData = getPIDPlantData(GData)
         % Format plant data for PID tuning
         % Default implementation: convert to ZPKDATA
         GData = zpk(GData);
      end
      
      function [D,SFLAG] = sparss(D)
         % Conversion to SPARSS
         [D,SFLAG] = D.ss();
         D = sparss(D);
      end
      
      function [D,SFLAG] = mechss(D)
         % Conversion to MECHSS
         [D,SFLAG] = D.ss();
         D = mechss(D);
      end
      
      function Dlft = genss(D)
         % Conversion to ltipack.lftdataSS
         Dlft = ltipack.lftdataSS(D.ss(),ltipack.LFTBlockWrapper.emptyBlockList());
      end
      
      function Dlft = genfrd(D,varargin)
         % Converts to ltipack.lftdataFRD
         IC = frd(D,varargin{:});
         Dlft = ltipack.lftdataFRD(elimDelay(IC,[],[],IC.Delay.IO),...
            ltipack.LFTBlockWrapper.emptyBlockList());
      end
      
      function Ds = idss(D, varargin)
         % Conversion to idpack.ssdata. Default implementation.
         %
         % Split: If TRUE, the last Ny input channels are treated as originating
         % from noise sources resulting in non-zero K matrix (default: FALSE). For
         % using this option, the model must have Nu+Ny input channels, where Nu =
         % number of measured inputs.
         Ds = idss(D.ss(), varargin{:});
      end
      
      function Dp = idpoly(D, varargin)
         % Conversion to polydata. Default implementation.
         Warn = ctrlMsgUtils.SuspendWarnings('Control:ltiobject:UseSSforInternalDelay'); %#ok<NASGU>
         Dp = idpoly(tf(D),varargin{:});
      end
      
      function Dp = idtf(D, varargin)
         % Conversion to idpack.tfdata. Default implementation.
         Warn = ctrlMsgUtils.SuspendWarnings('Control:ltiobject:UseSSforInternalDelay'); %#ok<NASGU>
         Dp = idtf(D.tf(),varargin{:});
      end
      
      function Df = idfrd(D, freqs)
         % Conversion to idpack.frddata. 
         % Default implementation. Overridden for identified models.
         if ischar(freqs)
            [~,~,freqs,FocusInfo] = freqresp(D,3,[],false);
            freqs = ltipack.util.roundFreqFocus(NaN(1,2),FocusInfo.Focus,freqs);
         end
         Df = D.frd(freqs);
         % protect against ill-conditioned D leading to complex response values at 0/Nyquist
         if isreal(D)
            Ts = Df.Ts;
            w = Df.Frequency;
            if Ts==0
               i0 = w==0;
            else
               i0 = w==0 | abs(w-pi/Ts)<1e4*eps | abs(w+pi/Ts)<1e4*eps;
            end
            if any(i0)
               Df.Response(:,:,i0) = real(Df.Response(:,:,i0));
            end
         end
         Df = idfrd(Df);
      end
      
      %% Analysis
      function [z,nrk] = tzero(D)
         % Default implementation
         [z,nrk] = tzero(D.ss());
      end
      
      function [gpeak,fpeak] = norminf(D,varargin)
         [gpeak,fpeak] = norminf(D.ss(),varargin{:});
      end
      
      function [Index,fIndex] = sectorbnd(D,varargin)
         [Index,fIndex] = sectorbnd(D.ss(),varargin{:});
      end
      
      function [Index,fIndex] = ifpofp(D,IndexType,tol,fBand)
         [Index,fIndex] = ifpofp(D.ss(),IndexType,tol,fBand);
      end
      
      function n = normh2(D)
         n = normh2(D.ss());
      end
      
      function wc = getGainCrossover(D,g)
         wc = getGainCrossover(D.ss(),g);
      end
      
      function wc = getSectorCrossover(D,varargin)
         wc = getSectorCrossover(D.ss(),varargin{:});
      end
      
      function fb = bandwidth(D,varargin)
         fb = bandwidth(D.ss(),varargin{:});
      end

      function D = scaleTime(D,sf)
         % Time vector scaling tnew = sf * told
         Ts = D.Ts;
         if Ts>0
            D.Ts = sf * Ts;
         elseif Ts==0
            % Update delay values
            Delay = D.Delay;
            Delay.Input = sf * Delay.Input;
            Delay.Output = sf * Delay.Output;
            if isfield(Delay,'IO')
               Delay.IO = sf * Delay.IO;
            else
               Delay.Internal = sf * Delay.Internal;
            end
            D.Delay = Delay;
         end
      end
      
      function [D,xinit,ComputeX] = checkCompatibility(D,PlotType,xinit,varargin)
         % Check compatibility of model with a given plot type
         if any(strcmp(PlotType,{'step','impulse','lsim'}))
            if ~isproper(D)
               error(message('Control:analysis:NotSupportedSimulationImproperSys'))
            end
            % Always assume simulation starts from steady state for non-SS
            % models
            xinit = [];  ComputeX = false;
         elseif strcmp(PlotType,'sector')
            % System must be proper and delay free in CT
            if ~isproper(D)
               error(message('Control:analysis:getSectorIndex10'))
            elseif D.Ts==0 && hasdelay(D)
               error(message('Control:analysis:getSectorIndex9'))
            end
         end
      end

      %% Transformations
      function D = exp(D) %#ok<*MANU>
         % Default: not implemented
         error(message('Control:transformation:exp1'))
      end
      

      %% Delay management
      % Generic access to delay values (used by user-facing objects)
      function Value = getInputDelay(D)
         Value = D.Delay.Input;
      end
      
      function Value = getOutputDelay(D)
         Value = D.Delay.Output;
      end
      
      function iod = getIODelay(D,~)
         %   IOD = getIODelay(D) returns the I/O delay matrix.
         %   For state-space models, getIODelay tests if internal
         %   delays are equivalent to I/O delays and returns the
         %   equivalent I/O delay or NaN for each I/O pair.
         %
         %   IOD = getIODelay(D,'total') returns the total I/O delay
         %   including the input and output delays.
         iod = D.Delay.IO;  % default implementation
         if nargin>1
            [ny,nu] = size(iod);
            iod = iod + D.Delay.Input(:,ones(1,ny)).' + D.Delay.Output(:,ones(1,nu));
         end
      end
      
      function D = setInputDelay(D,Value,CheckFlag)
         % Sets input delay value
         if CheckFlag
            [~,nu] = iosize(D);
            Value = ltipack.util.checkInputDelay(Value,nu,D.Ts);
         end
         if any(Value<0)
            error(message('Control:ltiobject:setLTI1','InputDelay'))
         end
         D.Delay.Input = Value;
      end

      function D = setOutputDelay(D,Value,CheckFlag)
         % Sets output delay value
         if CheckFlag
            [ny,~] = iosize(D);
            Value = ltipack.util.checkOutputDelay(Value,ny,D.Ts);
         end
         D.Delay.Output = Value;
      end
      
      function D = setIODelay(D,Value,CheckFlag)
         % Sets I/O delay value
         if CheckFlag
            [ny,nu] = iosize(D);
            Value = ltipack.util.checkIODelay(Value,ny,nu,D.Ts);
         end
         D.Delay.IO = Value;
      end
      
      function boo = hasdelay(D)
         % Returns TRUE if model has delays.
         Delay = D.Delay;
         boo = any(Delay.Input) || any(Delay.Output) || any(Delay.IO(:));
      end
      
      function boo = hasDelayDynamics(varargin)
         % Returns T if the dynamics (poles) depend on the delays.
         boo = false;
      end
      
      function boo = hasInternalDelay(D)
         % Returns T if model has internal delays.
         % Treal I/O delays as internal for consistency with state space
         boo = any(D.Delay.IO(:));
      end

      % Delay manipulations
      [D,Tdio] = elimDelay(D,id,od,iod)
      D = checkDelay(D)
      D = utFoldDelay(D,InputDelay,OutputDelay)
      
   end

   %% RCTB methods
   methods
            
      function [isStable,omega] = getMUGrid(D,PadeN,FreqRange)
         % Checks stability and generates appropriate frequency grid
         % for mu analysis or synthesis.
         if hasInternalDelay(D)
            Dnd = pade(D,0,0,PadeN);
         else
            Dnd = D;
         end
         isStable = isstable(Dnd);
         if isStable
            omega = sigmaGrid(Dnd,FreqRange);
         else
            omega = [];
         end
      end
      
   end
   
   
   % Protected methods (delay management)
   methods(Access=protected)

      [Ni,No,Nf] = checkPadeOrders(D,Ni,No,Nf)
      [D1,D2] = catDelay(dim,D1,D2,Zero1,Zero2)
      [D1,D2,ElimFlag] = plusDelay(D1,D2,Zero1,Zero2)
      [D1,D2,ElimFlag] = mtimesDelay(D1,D2,ScalarFlag)
      boo = utIsIOScaling(D)

      function [w,Focus] = parseFreqSpec(D,wspec,REALPLOT)
         % Interprets W argument of BODE, NYQUIST,...
         Focus = ltipack.getFreqFocus(wspec,D.Ts,'log'); % log-scale focus, may contain NaNs
         if isempty(wspec) || iscell(wspec)
            w = [];  % no frequency grid
         else
            % User-defined grid
            w = wspec(:); % e.g., for LTI Viewer when freq. grid defined in Prefs
            if REALPLOT
               % Map w<0 to w>0 for plots of systems with real coefficients
               w = sort(abs(w));
               w(diff(log2(w))<1e3*eps,:) = [];
            end
         end
      end
      
      function MaxDelay = getMaxDelay(D)
         % Estimates max I/O delay.
         MaxDelay = 0;
         for ct=1:length(D)
            Delay = D(ct).Delay; %#ok<*PROP>
            iod = bsxfun(@plus,bsxfun(@plus,Delay.IO,Delay.Output),Delay.Input.');
            MaxDelay = max([MaxDelay;iod(:)]);
         end
      end

      function D1 = timesDelay(D1,D2)
         % Delay management for D=D1.*D2
         Delay1 = D1.Delay;  Delay2 = D2.Delay;
         Delay1.Input = Delay1.Input + Delay2.Input;
         Delay1.Output = Delay1.Output + Delay2.Output;
         Delay1.IO = Delay1.IO + Delay2.IO;
         D1.Delay = Delay1;
      end

      function D = elimZeroDelay(D)
         % Nothing to do in general
      end

   end
   
   %% Static methods
   methods(Static, Access=protected)
      
      [Delay,fiod] = discretizeDelay(Delay,Ts)
      
      [wout,magout,phout] = fShowPhaseShift(Grade,w,mag,ph,Ts,Td,Focus)
            
   end

end