function [f_vec, S_aa, S_ff_unnotched, S_ff_limit, S_aa_notched, Mass_eff] = ...
    calc_force_limit_psd(psd_spec, f_am, M_A_mag, df, C, M, Phi, iu, mode_flag)

% CALC_FORCE_LIMIT_PSD  Force-limited acceleration & force PSD per NASA-HDBK-7004C.
%
% iu        : 1 = English (lbf, in, sec, g = 386.089 in/s^2)
%             2 = Metric  (N,   m,  sec, g = 9.80665 m/s^2)
%
% mode_flag : 1 -> Simple Mc^2 method, flat envelope
%                  S_FF(f) = (C * Mtot)^2 * max(S_AA) * g^2
%             2 -> Mc^2 with frequency-varying input
%                  S_FF(f) = (C * Mtot)^2 * S_AA(f) * g^2
%             3 -> Apparent-mass method (preferred)
%                  S_FF(f) = C^2 * |M_A(f)|^2 * S_AA(f) * g^2
%
% Mtot in modes 1 and 2 is the payload effective mass (excludes the shaker
% mass at DOF 1), computed from the modal participation factors.
%
% Input/output PSD are in g^2/Hz.  Force PSD is in lbf^2/Hz or N^2/Hz.

%% 0. Unit system
if iu == 1
    g = 386.089;
elseif iu == 2
    g = 9.80665;
else
    error('iu must be 1 (English) or 2 (Metric)');
end

if nargin < 9 || isempty(mode_flag)
    mode_flag = 3;
end

%% 1. Frequency vector (uniform linear grid)
f_min = psd_spec(1,1);
f_max = psd_spec(end,1);
f_vec = (f_min:df:f_max)';

%% 2. Log-log interpolation of input PSD onto f_vec
log_f_spec  = log10(psd_spec(:,1));
log_S_spec  = log10(psd_spec(:,2));
log_f_query = log10(f_vec);
S_aa = 10.^(interp1(log_f_spec, log_S_spec, log_f_query, 'linear', 'extrap'));

%% 3. Apparent mass interpolation
M_A_interp = interp1(f_am, M_A_mag, f_vec, 'linear', 'extrap');
M_A_interp(M_A_interp <= 0) = eps;

%% 4. Convert input PSD to physical acceleration units
S_aa_phys = (g^2) * S_aa;

%% 5. Effective payload mass for Mc-style methods.
%    Works for any modal normalization (mass-normalized or otherwise).
switch mode_flag
    case 2
        r = ones(size(M,1), 1);            % base translation
    otherwise
        r = ones(size(M,1), 1);
        r(1) = 0;                          % exclude shaker DOF
end

Gamma   = (Phi' * M) * r;
Mi      = diag(Phi' * M * Phi);
Mass_eff = sum( Gamma.^2 ./ max(Mi, eps) );

%% 6. Unnotched force PSD  (apparent-mass prediction)
S_ff_unnotched = (M_A_interp.^2) .* S_aa_phys;

%% 7. Force-limit PSD
switch mode_flag
    case 1
        S_aa_ref_phys = (g^2) * max(S_aa);
        S_ff_limit = (C * Mass_eff)^2 * S_aa_ref_phys * ones(size(f_vec));
    case 2
        S_ff_limit = (C * Mass_eff)^2 * S_aa_phys;
    case 3
        S_ff_limit = (C^2) * (M_A_interp.^2) .* S_aa_phys;
    otherwise
        error('mode_flag must be 1, 2, or 3');
end

%% 8. Notched acceleration PSD
S_aa_limit_phys = S_ff_limit ./ max(M_A_interp.^2, eps);
S_aa_limit      = S_aa_limit_phys / (g^2);

S_aa_notched = min(S_aa, S_aa_limit);

% Gentle log-domain smoothing
S_aa_notched = log_smooth(S_aa_notched, 9);

% Re-enforce ceilings after smoothing
S_aa_notched = min(S_aa_notched, S_aa);
S_aa_notched = min(S_aa_notched, S_aa_limit);

end


function y = log_smooth(x, window)
x = max(x, eps);
logx = log10(x);
N    = numel(logx);
half = floor(window/2);
logy = logx;
for k = 1:N
    a = max(1, k - half);
    b = min(N, k + half);
    logy(k) = mean(logx(a:b));
end
y = 10.^logy;
end
