function T = compute_T_Rayleigh(m, k, u, unitType)
% compute_T_Rayleigh
% Computes fundamental period using Rayleigh method
%
% INPUTS:
%   m        = vector of lumped story masses
%              (lbm if 'english', kg if 'metric')
%   k        = vector of story stiffnesses
%              (lbf/in if 'english', N/m if 'metric')
%   u        = assumed mode shape (inches if 'english', meters if 'metric')
%   unitType = 'english' or 'metric'
%
% OUTPUT:
%   T = fundamental period (seconds)

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

    n = length(m);

    if length(k) ~= n || length(u) ~= n
        error('m, k, and u must have same length.');
    end

    % Unit correction factor
    switch lower(unitType)
        case 'english'
            % Using lbm, lbf, inches
            gc = 386.088;   % lbm*in/(lbf*s^2)
        case 'metric'
            gc = 1.0;       % no correction needed
    end

    % Rayleigh numerator
    num = sum(m .* u.^2);

    % Rayleigh denominator
    den = 0;

    for i = 1:n
        if i == 1
            delta = u(i);  % base displacement = 0
        else
            delta = u(i) - u(i-1);
        end
        den = den + k(i) * delta^2;
    end

    % Natural circular frequency
    omega = sqrt(gc * den / num);

    % Period
    T = 2*pi / omega;

end
