function T = compute_T_Rayleigh_force(m, delta, F, unitType)
% compute_T_Rayleigh_force
%
% Rayleigh force–displacement method for fundamental period
%
%   ω² = g * Σ(F_i * δ_i) / Σ(m_i * δ_i²)
%   T  = 2π / ω
%
% UNIT REQUIREMENTS:
%   English:
%       m      = lbm
%       delta  = inches
%       F      = lbf
%       g      = 386.4 in/s^2
%
%   Metric:
%       m      = kg
%       delta  = meters
%       F      = N
%       g      = 9.80665 m/s^2
%
% OUTPUT:
%   T = fundamental period (seconds)

    arguments
        m (:,1) double {mustBePositive}
        delta (:,1) double
        F (:,1) double
        unitType (1,:) char {mustBeMember(unitType,{'english','metric'})}
    end

    n = length(m);

    if length(delta) ~= n || length(F) ~= n
        error('m, delta, and F must be vectors of equal length.');
    end

    % Remove any rigid-body or zero displacement cases
    if all(abs(delta) < eps)
        error('All displacements are zero — cannot compute Rayleigh period.');
    end

    % Gravity (consistent with displacement units)
    if strcmpi(unitType,'english')
        g = 386.4;        % in/s^2
    else
        g = 9.80665;      % m/s^2
    end

    % Rayleigh numerator and denominator
    forceTerm = sum(F .* delta);        % Σ(F_i δ_i)
    massTerm  = sum(m .* delta.^2);     % Σ(m_i δ_i²)

    if massTerm <= 0 || forceTerm <= 0
        error('Invalid Rayleigh terms. Check stiffness or load direction.');
    end

    % Angular frequency squared
    omega2 = g * forceTerm / massTerm;

    if omega2 <= 0
        error('Computed negative or zero ω². Check input units.');
    end

    % Period
    T = 2*pi / sqrt(omega2);

end