function [T, M, K, Phi_mn, Gamma, modalMassRatio, omega] = ...
    compute_shear_building_eig(m, kStory, unitType)

% compute_period_eig
%
% Modal analysis of an n-story shear building.
%
% INPUTS:
%   m        = story masses
%              English: lbm
%              Metric : kg
%
%   kStory   = story stiffnesses
%              English: lbf/in
%              Metric : N/m
%
%   unitType = 'english' or 'metric'
%
% OUTPUTS:
%   T               = fundamental period (sec)
%   M               = assembled mass matrix
%   K               = assembled stiffness matrix
%   Phi_mn          = mass-normalized eigenvectors
%   Gamma           = modal participation factors
%   modalMassRatio  = modal mass ratios
%   omega           = natural circular frequencies (rad/sec)

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

n = length(m);

%% ==========================================================
% MASS MATRIX
%% ==========================================================

if strcmpi(unitType,'english')
    
    gc = 386.0886;  % lbm*in/(lbf*s^2)
    M = diag(m / gc);   % dynamic mass (lbf*s^2/in)

else  % metric
    
    M = diag(m);        % kg
    
end

%% ==========================================================
% STIFFNESS MATRIX
%% ==========================================================

K = zeros(n);

for i = 1:n
    ki   = kStory(i);
    kip1 = 0;
    if i < n
        kip1 = kStory(i+1);
    end

    K(i,i) = ki + kip1;

    if i > 1
        K(i,i-1) = -ki;
        K(i-1,i) = -ki;
    end

    if i < n
        K(i,i+1) = -kip1;
        K(i+1,i) = -kip1;
    end
end

%% ==========================================================
% GENERALIZED EIGENVALUE SOLUTION
%% ==========================================================

[Phi,Omega2] = eig(K,M);

w2 = real(diag(Omega2));
[w2,sortIdx] = sort(w2);   % sort ascending
Phi = Phi(:,sortIdx);

omega = sqrt(w2);          % natural circular frequencies

%% ==========================================================
% MASS NORMALIZATION
% Enforce: Phi_mn' * M * Phi_mn = I
%% ==========================================================

Phi_mn = zeros(size(Phi));

for i = 1:n
    scale = sqrt(Phi(:,i)' * M * Phi(:,i));
    Phi_mn(:,i) = Phi(:,i) / scale;
end

%% ==========================================================
% PARTICIPATION FACTORS (base excitation)
% Gamma_i = (phi_i^T M 1)
% since mass-normalized → denominator = 1
%% ==========================================================

r = ones(n,1);   % influence vector for base excitation

Gamma = zeros(n,1);
modalMass = zeros(n,1);

for i = 1:n
    Gamma(i) = Phi_mn(:,i)' * M * r;
    modalMass(i) = Gamma(i)^2;  % because mass-normalized
end

%% ==========================================================
% MODAL MASS RATIOS
%% ==========================================================

totalMass = sum(diag(M));
modalMassRatio = modalMass / totalMass;

%% ==========================================================
% FUNDAMENTAL PERIOD
%% ==========================================================

omega1 = omega(1);

T = 2*pi ./ omega;   % vector of periods (sec)

end