function [K,Kinv] = referenceImpedanceArrays(data,z0)
%referenceImpedanceArrays Make 3d reference impedance arrays.
% [K,Kinv] = referenceImpedanceArrays(z0,dataSize) makes reference
% impedance arrays of size dataSize, handling either a scalar z0, a
% per-port vector z0 of length dataSize(1) (or dataSize(2)), and for
% backward compatibility, a per-frequency vector z0 of length dataSize(3).
% The use of per-frequency z0 is NOT recommended because of it could result
% in non-causality.
%
% In the unlikely case that the number of ports (dataSize(1) or
% dataSize(2)) equals the number of frequencies (dataSize(3)), per-port z0
% is assumed unless z0 is specified precisely as a 3-dimensional vector.

%   Copyright 2022 The MathWorks, Inc.

arguments
    data {mustBeNumeric,mustBeNonempty,mustBeNonNan, ...
        rf.internal.mustBeNxNxM}
    z0 {mustBeNumeric,mustBeNonempty,mustBeFinite, ...
        rf.internal.mustBeRealValued, ...
        rf.internal.mustBeDataSized(z0,data)} = 50
end

if isscalar(z0) || ...
        (size(data,1) == size(data,3) && ...
        ndims(z0) == 3 && numel(z0) == size(data,3))
    % Note rare per-frequency exception if data is NxNxN and z0 is 1x1xN
    sqz = sqrt(z0);
elseif numel(z0) == size(data,1) % z0 per-port
    sqz = sqrt(z0(:));
else % z0 per-frequency
    sqz = sqrt(reshape(z0,[1,1,size(data,3)]));
end

Z = zeros(size(data)); % preallocate storage
I = eye(size(data,1));
% Carefully use implicit expansion for all shape combinations
K = sqz .* I + Z;
Kinv = (1./sqz) .* I + Z;
