function s_params_new = s2s(s_params,z0,z0_new)
%S2S Convert S-parameters to S-parameters with a different reference impedance
%   S_PARAMS_NEW = S2S(S_PARAMS, Z0, Z0_NEW) converts the scattering
%   parameters S_PARAMS with the reference impedance Z0 into the scattering
%   parameters S_PARAMS_NEW with the reference impedance Z0_NEW.
%
%   S_PARAMS is a complex NxNxM array of M N-port S-parameters. Z0 is the
%   reference impedance of the input S-parameters. Z0_NEW is the reference
%   impedance of the output S-parameters, the default is 50 ohms. The
%   output S_PARAMS_NEW is a complex NxNxM array of M N-port S-parameters.
%
%   Z0 and Z0_NEW are each specified as one of the following:
%   - Scalar used for all ports and frequencies
%   - Vector of length size(S_PARAMS,1), for per-port impedances
%   - Vector of length size(S_PARAMS,3), for per-frequency impedances
%
%   If size(S_PARAMS,1) == size(S_PARAMS,3) == numel(Z0), then Z0 specifies
%   per-port impedances, unless Z0 is a 3-D array.  Similarly for Z0_NEW.
%
%   Note that Z0 may be complex, but Z0_NEW must be real-valued, i.e. have
%   zero imaginary part.
%
%   See also S2Y, S2Z, S2H, S2ABCD, Y2S, Z2S, H2S, ABCD2S

%   Copyright 2003-2023 The MathWorks, Inc.

arguments
    s_params {mustBeNumeric,mustBeNonempty,mustBeNonNan, ...
        rf.internal.mustBeVectorOrNxNxM}
    z0 {mustBeNumeric,mustBeNonempty,mustBeFinite, ...
        rf.internal.mustBeDataSized(z0,s_params)}
    z0_new {mustBeNumeric,mustBeNonempty,mustBeFinite, ...
        rf.internal.mustBeRealValued, ...
        rf.internal.mustBeDataSized(z0_new,s_params)} = 50
end
if isvector(s_params)
    s_params = reshape(s_params,1,1,[]);
end

if all(z0 == z0_new)
    s_params_new = s_params;
    return
end

% Calculate the new S-parameters if any change in reference impedance
[K1,K1inv] = referenceImpedanceArrays(s_params,z0);
[K2,K2inv] = referenceImpedanceArrays(s_params,z0_new);
alpha = pagemtimes(K2inv,K1);
alphaInv = pagemtimes(K1inv,K2);
P = alpha + alphaInv;
M = alpha - alphaInv;
sM = pagemtimes(s_params,M);
sP = pagemtimes(s_params,P);
s_params_new = pagemldivide(P + sM,M + sP);
