function [fit,err] = rationalfit(varargin)
%RATIONALFIT Perform rational fitting to complex frequency-dependent data.
%   RATIONALFIT uses vector fitting with complex frequencies
%   S = j*2*pi*FREQ and DATA to construct a rational function.
%
%            C(1)     C(2)           C(n)
%   F(S) = (------ + ------ + ... + ------ + D) * EXP(-S*DELAY*DELAYFACTOR)
%           S-A(1)   S-A(2)         S-A(n)
%
%   FREQ is a vector of nonnegative frequencies in Hz. DATA is a single
%   vector, a 2-D array (column vectors), or a 3-D array of network
%   parameter data.  The length of the DATA vectors must be equal to the
%   length of FREQ. The RATIONALFIT function returns an RFMODEL.RATIONAL
%   object. For more information, see rfmodel.rational documentation page. 
%
%   FIT = RATIONALFIT constructs a non-reflective one-port.
%
%   FIT = RATIONALFIT(FREQ,DATA) returns a rational fit for the given FREQ
%   and DATA.
%
%   FIT = RATIONALFIT(S) for N-port S-parameters object S fits all NxN
%   elements of S and returns FIT as an NxN array of RFMODEL.RATIONAL
%   objects.  This is equivalent to
%   FIT = RATIONALFIT(S.Frequencies,S.Parameters).
%
%   FIT = RATIONALFIT(FILENAME) constructs a fit for network data
%   supplied in FILENAME (*.snp,*.ynp,etc ...).
%
%   FIT = RATIONALFIT(...,TOL) constructs a fit that attempts to satisfy a
%   relative error tolerance specified in dB by the optional argument TOL.
%   The default value of TOL is -40 dB, i.e. 20*log10(1e-2).
%
%   FIT = RATIONALFIT(S,I,J) for S-parameters object S fits Sij using
%   FREQ = S.Frequencies and DATA = rfparam(S,I,J), returning a scalar FIT.
%
%   FIT = RATIONALFIT(...,NAME=VALUE) constructs a rational fit using
%   name-value arguments.
%
%   [FIT,ERR] = RATIONALFIT(...) returns the error of the fit in dB as the
%   second output argument.
%
%
%   Name-value pair arguments:
%   --------------------------
%
%   'Tolerance'         Relative error tolerance in decibels, specified as
%                       a scalar less than or equal to zero.
%
%                       Default: -40
%
%   'TendsToZero'       A logical variable specifying the behavior of the
%                       rational fit F(S) for large S. When true, the D
%                       term of the fit is set to zero so that the rational
%                       fit F(S) tends to zero as S approaches infinity.
%                       When false, a nonzero D will be allowed.
%
%                       Default: true (set D to zero so that F(S) -> 0)
%
%   'NPoles'            Either a scalar nonnegative integer specifying the
%                       number of poles, or a two-element vector of
%                       nonnegative integers specifying a search range. The
%                       number of poles might be less than specified if any
%                       computed residues are zero.
%
%                       Default: [0 48]
%
%   'Weight'            A nonnegative frequency-weighting array equal in
%                       size to the DATA array. WEIGHT values less than one
%                       de-emphasize the DATA at corresponding FREQ values.
%
%                       Default: ones(size(DATA)) (use equal weighting)
%
%   'DelayFactor'       A real scalar between 0 and 1 inclusive, specifying
%                       the fraction of estimated delay to extract from the
%                       DATA before rational fitting,
%                       NEW_DATA = DATA * EXP(S*DELAY*DELAYFACTOR).
%
%                       Default: 0 (do not remove any delay)
%
%   'IterationLimit'    The maximum number of rationalfit iterations,
%                       specified as a positive integer. Provide a
%                       two-element vector to specify the minimum and
%                       maximum.
%
%                       Default: [4 24]
%
%   'WaitBar'           A logical variable specifying whether to display a
%                       waitbar.
%
%                       Default: false (run without a waitbar)
%
%   'PoleSharing'       'All'       : Poles of each S-parameter term are
%                                     identical
%                       'Column'    : Poles of S-parameter column terms
%                                     are identical
%                       'Individual': Poles of each S-parameter term are
%                                     individual
%
%                       Default: 'All'
%
%   EXAMPLE:
%
%   S = sparameters('passive.s2p');
%   fit = rationalfit(S,1,1,'TendsToZero',false);
%
%   See also S2TF, SNP2SMP, SPARAMETERS/RATIONALFIT, RFMODEL.RATIONAL

%   B. Gustavsen and A. Semlyen, "Rational Approximation of Frequency
%   Domain Responses by Vector Fitting", IEEE Transactions on Power
%   Delivery, Vol. 14, No. 3, pp. 1052-1061, July 1999.
%
%   B. Gustavsen, "Improving the Pole Relocating Properties of Vector
%   Fitting", IEEE Transactions on Power Delivery, Vol. 21, No. 3, pp.
%   1587-1592, July 2006.
%
%   D. Deschrijver, M. Mrozowski, T. Dhaene, D. De Zutter, "Macromodeling
%   of Multiport Systems Using a Fast Implementation of the Vector Fitting
%   Method", IEEE Microwave and Wireless Components Letters, Vol. 18, No.
%   6, pp. 383-385, June 2008.

%   Copyright 2006-2024 The MathWorks, Inc.

    [varargin{:}] = convertStringsToChars(varargin{:});

if isempty(varargin)
    % The default rational function has one data point, 0, at
    % frequency 0, which fits perfectly with no poles and a direct
    % term of 0.
    freqs = 0;
    data = 0;
    varargin = {'TendsToZero',1};
elseif isa(varargin{1},'char') || isa(varargin{1},'string')
    % Assumes a valid s-parameter file else error
    sObj = sparameters(varargin{1});
    freqs = sObj.Frequencies;
    data  = sObj.Parameters;
    if (nargin>1)
        if isnumeric(varargin{2})
            varargin = {'Tolerance' varargin{2:end}};
        else
            varargin = varargin(2:end);
        end
    else
        varargin = cell.empty;
    end
elseif isa(varargin{1},'numeric')  &&  nargin>=2    &&                  ...
       isa(varargin{2},'numeric')
    % Cases covered by this case:
    % 1) Sparameter Objects enters via "rationalfit(sObject)",
    %    See toolbox/rf/rf/@sparameters/rationalfit
    % 2) Syntax: rationalfit(freq,data,varargin)
    if ~iscolumn(varargin{1})  &&  ~isrow(varargin{1})
        error(message('rf:rationalfit:NotRowOrColumn'))
    end
    freqs = varargin{1}(:);
    numfreq = length(freqs);
    data = varargin{2};
    datasize = size(data);
    if length(datasize) == 3
        % S-parameter data 3d matrix
        ndatafreq = datasize(3);
    elseif length(datasize) == 2
        % S-parameter data 2d matrix
        if numfreq == 1                       %#ok<ISCL>
            numfreq = 1;
            ndatafreq = 1;
        else
            if isrow(data) % silently accept a single row vector
                data = data(:);
                datasize = size(data);
            end
            ndatafreq = datasize(1);
        end
    elseif isscalar(datasize)
        % S-parameter data 1
        ndatafreq = datasize;
    else
        error(message('rf:rationalfit:WrongFreqOrDataInput'))
    end
    if numfreq ~= ndatafreq
        error(message('rf:rationalfit:WrongFreqOrDataInput'))
    end
    uniqueFreqs = unique(freqs);
    if numel(uniqueFreqs) < numfreq
        error(message('rf:rationalfit:FrequenciesNotUnique'))
    end
    if (nargin>2)
        if isnumeric(varargin{3})
            varargin = {'Tolerance' varargin{3:end}};
        else
            varargin = varargin(3:end);
        end
    else
        varargin = cell.empty;
    end
else
    error(message('rf:rationalfit:InvalidArgumentPattern'))
end

if isempty(varargin)
    poleShareVal = 'All';
else
    poleShareLogical = strcmpi('PoleSharing',varargin);
    if poleShareLogical(end)==true
        error(message('rf:rationalfit:PoleSharingNotSpecified'))
    elseif any(poleShareLogical)
        shareOpt = varargin(circshift(poleShareLogical,1));
        poleShareVal = shareOpt{end};
        if ~any(strcmpi(poleShareVal,{'All','Column','Individual'}))
                error(message('rf:rationalfit:PoleSharingInvalid'))
        end
    else
        poleShareVal = 'All';
    end
end

dataDimensions = size(data);
numDimensions = size(dataDimensions,2);
switch lower(poleShareVal)
    case 'column'
        numCols = dataDimensions(2);
        ratTerms = rfmodel.rational;
        err = NaN(numCols,1);
        if numDimensions==3
            numRows = dataDimensions(1);
            fit = ratTerms(ones(numRows,numCols));
            for colIdx = 1:numCols
                [partFit,err(colIdx)] =                                 ...
                    performFit(freqs,data(:,colIdx,:),varargin{:});
                fit(numRows*(colIdx-1)+(1:numRows)) = partFit(1:numRows);
            end
        else
            fit = ratTerms(ones(1,numCols));
            for colIdx = 1:numCols
                [fit(colIdx),err(colIdx)] =                             ...
                    performFit(freqs,data(:,colIdx),varargin{:});
            end
        end
    case 'all'
        [fit,err] = performFit(freqs,data,varargin{:});
    case 'individual'
        numCols = dataDimensions(2);
        ratTerms = rfmodel.rational;
        if numDimensions==3
            numRows = dataDimensions(1);
            fit = ratTerms(ones(numRows,numCols));
            err = NaN(numRows,numCols);
            for colIdx = 1:numCols
                for rowIdx = 1:numRows
                    [fit(rowIdx,colIdx),err(rowIdx,colIdx)] =         ...
                        performFit(freqs,data(rowIdx,colIdx,:),varargin{:});
                end
            end
        else
            fit = ratTerms(ones(1,numCols));
            err = NaN(numCols,1);
            for colIdx = 1:numCols
                [fit(colIdx),err(colIdx)] =                             ...
                    performFit(freqs,data(:,colIdx),varargin{:});
            end
        end
end
end  % End of rationalfit

function [fit,err] = performFit(freqs,data,varargin)
[a,c,d,delay,outsize,err] = rationalfitcore(freqs,data,varargin{:});
temp = rfmodel.rational;
fit = temp(ones(outsize)); % array of identical handles
for k = 1:numel(fit)
    fit(k) = rfmodel.rational('a',a,'c',c(:,k),'d',d(k),'delay',delay(k));
end
end
