function [R_dof, rho, Ri_modes] = cqc_from_spectrum(Phi, fn, zeta, S, varargin)
% CQC_FROM_SPECTRUM
% Combine modal peaks into per-DOF peaks using CQC (equal damping).
%
% Inputs
%   Phi   : [nDOF x nModes] mass-normalized eigenvectors (columns = modes)
%   fn    : [nModes x 1] natural frequencies in Hz
%   zeta  : scalar damping ratio (e.g., 0.05)
%   S     : [nModes x 1] response spectrum value at each mode's fn
%           (must correspond to the SAME response type you want to combine
%            in the DOFs, e.g., modal relative displacement, velocity, accel)
%
% Name-Value (optional)
%   'UseAbsPhi' : logical (default true). If true, uses abs(Phi) when
%                 mapping modal peaks to DOFs (sign is irrelevant for
%                 envelope-type peaks).
%   'Gamma'     : [nModes x 1] participation factors to scale S before
%                 mapping to DOFs. Default = ones(nModes,1).
%                 (Useful if S is a ground-accel spectrum for base excitation
%                  and you want modal-coordinate peaks scaled by Gamma.)
%
% Outputs
%   R_dof   : [1 x nDOF] CQC-combined peak for each DOF
%   rho     : [nModes x nModes] modal correlation matrix
%   Ri_modes: [nModes x nDOF] modal peak at each DOF used in CQC
%
% Notes
%   • CQC correlation (equal damping; Der Kiureghian/Chopra form):
%       r_ij = ω_i/ω_j
%       ρ_ij = (8 ξ^2 r (1+r)^2)/((1 - r^2)^2 + 4 ξ^2 r (1+r)^2)
%   • If you need absolute responses from base excitation, ensure S (and
%     optionally Gamma) represent the correct modal coordinate peak quantity
%     you intend to map to the physical DOF response via Φ.
%
% Example
%   [R_dof, rho, Ri_modes] = cqc_from_spectrum(Phi, fn, 0.05, S, 'Gamma', Gamma);
%

    % ---- input checks ----
    if ~isscalar(zeta) || zeta < 0
        error('zeta must be a nonnegative scalar (e.g., 0.05).');
    end
    [nDOF, nModes] = size(Phi);
    fn = fn(:);
    S  = S(:);

    if numel(fn) ~= nModes
        error('Length(fn) must equal number of columns in Phi (nModes).');
    end
    if numel(S) ~= nModes
        error('Length(S) must equal number of modes.');
    end
    if any(fn <= 0)
        error('All natural frequencies (fn) must be positive.');
    end

    % ---- parse options ----
    p = inputParser;
    p.addParameter('UseAbsPhi', true, @(x)islogical(x) && isscalar(x));
    p.addParameter('Gamma', ones(nModes,1), @(x)isvector(x) && numel(x)==nModes);
    p.parse(varargin{:});
    useAbs = p.Results.UseAbsPhi;
    Gamma  = p.Results.Gamma(:);

    % ---- frequencies (rad/s) and correlation matrix rho ----
    omega = 2*pi*fn;
    r = omega ./ omega.';  % nModes-by-nModes
    xi  = max(zeta, 1e-8); % guard tiny damping
    num = 8*(xi^2) .* r .* (1 + r).^2;
    den = (1 - r.^2).^2 + 4*(xi^2) .* r .* (1 + r).^2;
    rho = num ./ den;
    rho(1:nModes+1:end) = 1.0;
    rho = 0.5*(rho + rho.');  % symmetrize

    % ---- modal peaks mapped to each DOF ----
    % Scale the modal spectrum by participation if provided
    S_eff = S .* Gamma;      % [nModes x 1]

    % Map modal coordinate peaks to physical DOFs via eigenvector components.
    % Row j of Ri_modes corresponds to mode j; column k to DOF k.
    if useAbs
        Ri_modes = (abs(Phi)).';  % nModes x nDOF
    else
        Ri_modes = (Phi).';
    end
    Ri_modes = Ri_modes .* S_eff; % multiply each row j by S_eff(j)

    % ---- CQC per DOF (column-wise combination) ----
    R_dof = zeros(1, nDOF);
    for k = 1:nDOF
        Rik = Ri_modes(:,k);      % [nModes x 1]
        Sij = Rik * Rik.';        % outer product
        R_dof(k) = sqrt( sum(rho .* Sij, 'all') );
    end
end
