function [Zp, Yp, k, D, mu, fc_crit] = true_point_impedance(fc, E, nu, h, rho, a, eta, r0, p, c0)
% TRUE_POINT_IMPEDANCE
% True mechanical point impedance of a thin plate, smoothly blended
% (by frequency ratio) to the Cremer–Heckl high-frequency constant.
%
%   [Zp, Yp, k, D, mu, fc_crit] = true_point_impedance( ...
%       fc, E, nu, h, rho, a, eta, r0, p, c0)
%
% Inputs
%   fc    - frequency array in Hz (vector)
%   E     - Young's modulus (Pa)
%   nu    - Poisson's ratio (-)
%   h     - plate thickness (m)
%   rho   - mass density (kg/m^3)
%   a     - effective contact radius for the logarithmic cutoff (m)
%   eta   - OPTIONAL hysteretic loss factor for bending stiffness (default 0)
%   r0    - OPTIONAL freq-ratio where blend centers (default 1.2 ~ near fc)
%   p     - OPTIONAL blend steepness exponent (default 6; higher = sharper)
%   c0    - OPTIONAL sound speed in surrounding fluid (m/s), default 343
%
% Outputs
%   Zp      - blended point impedance (N·s/m), complex
%   Yp      - driving-point mobility (m/N·s) = 1./Zp
%   k       - bending wavenumber (1/m)
%   D       - bending stiffness (N·m)
%   mu      - mass per unit area (kg/m^2)
%   fc_crit - coincidence/critical frequency in the surrounding fluid (Hz)
%
% Formulas
%   D  = E*h^3 / [12*(1 - nu^2)],   mu = rho*h
%   k  = (( (2*pi*fc).^2 * mu / D ))^(1/4)
%   Z_true = (8*pi*D*k.^2) ./ ( j*2*pi*fc .* ( 1 + j*(2/pi)*log(k*a) ) )
%   Z_CH   = 8*sqrt(D*mu) * (1 + j)     (Cremer–Heckl high-frequency constant)
%   fc_crit = (c0^2 / (2*pi)) * sqrt(mu / D)
%
% Blend (in frequency ratio r = f/fc_crit):
%   r  = max(fc(:),eps) / fc_crit;
%   r0 = 1.2; p = 6;
%   w  = 1 ./ (1 + (r./r0).^p);           % low f: w≈1 (true), high f: w≈0 (C&H)
%   Zp = w .* Z_true + (1 - w) .* Z_CH;

    % ---- Defaults for optional inputs ----
    if nargin < 7 || isempty(eta),  eta  = 0;   end
    if nargin < 8 || isempty(r0),   r0   = 1.2; end
    if nargin < 9 || isempty(p),    p    = 6;   end
    if nargin < 10 || isempty(c0),  c0   = 343; end   % m/s (air, ~20°C)

    % Preserve input shape
    fc_shape = size(fc);
    fc = fc(:);                                 % column vector

    % ---- Section properties ----
    D  = E*h^3 / (12*(1 - nu^2));               % bending stiffness (N·m)
    if eta ~= 0
        D = D * (1 + 1j*eta);                   % hysteretic loss in D
    end
    mu = rho * h;                               % mass per unit area (kg/m^2)

    % ---- Critical (coincidence) frequency (Hz) in surrounding fluid ----
    fc_crit = (c0^2 / (2*pi)) * sqrt(mu / D);

    % ---- Bending wavenumber ----
    omega = 2*pi*fc;                            % rad/s
    k = (omega.^2 .* mu ./ D).^(1/4);           % 1/m (principal branch)

    % ---- True (logarithmic) point impedance (avoid log(0)) ----
    ka = k .* a;
    ka(ka == 0) = eps;
    bracket = 1 + 1j*(2/pi)*log(ka);            % 1 + j*(2/pi)*ln(ka)
    Z_true = (8*pi*D .* k.^2) ./ (1j*omega .* bracket);

    % ---- Cremer–Heckl constant (semi-infinite/high-f asymptote) ----
    Z_CH = 8 * sqrt(D*mu) * (1 + 1j);

    % ---- Smooth frequency-ratio blend around fc_crit ----
    % After computing fc_crit and Z_true, Z_CH:
    r  = max(fc(:), eps) / fc_crit;             % frequency ratio
    % You can also vectorize without max(:):
    % r = max(fc, eps) ./ fc_crit; r = r(:);
    % Use the given/ default r0 and p:
    w  = 1 ./ (1 + (r./r0).^p);                 % low f: w≈1 (true), high f: w≈0 (C&H)
    Zp = w .* Z_true + (1 - w) .* Z_CH;

    % ---- Mobility & DC limits ----
    Yp = 1 ./ Zp;
    Zp(fc == 0) = Inf;
    Yp(fc == 0) = 0;

    % ---- Restore original shapes ----
    Zp = reshape(Zp, fc_shape);
    Yp = reshape(Yp, fc_shape);
    k  = reshape(k,  fc_shape);
end
