function out = seismic_elf_shear(m, h, kStory, SDS, SD1, R, Ie, Cd, unitType, periodMethod, riskCategory)
% seismic_elf_shear
%
% ASCE 7-22 Equivalent Lateral Force (ELF) solver for shear buildings.
%
% Units:
%   English: m = lbm, h = ft, kStory = lbf/in
%   Metric : m = kg,  h = m,  kStory = N/m

arguments
    m (:,1) double {mustBePositive}
    h (:,1) double {mustBePositive}
    kStory (:,1) double {mustBePositive}
    SDS (1,1) double {mustBePositive}
    SD1 (1,1) double {mustBePositive}
    R   (1,1) double {mustBePositive}
    Ie  (1,1) double {mustBePositive}
    Cd  (1,1) double {mustBePositive}
    unitType (1,:) char {mustBeMember(unitType,{'english','metric'})}
    periodMethod (1,:) char {mustBeMember(periodMethod, ...
        {'trialRayleigh','iterateMode','eig'})}
    riskCategory (1,:) char {mustBeMember(riskCategory,{'I','II','III','IV'})} = 'II'
end

n = length(m);

%% ==========================================================
% UNIT HANDLING
%% ==========================================================

if strcmpi(unitType,'english')

    % Convert geometry
    hU = 12*h;                    % ft → in

    % Constants
    g  = 386.0886;                % in/s^2
    gc = 386.0886;                % lbm*in/(lbf*s^2)

    % Seismic weight (lbf)
    W = m * g / gc;

    % Dynamic mass matrix (lbf*s^2/in)
    Mdyn = diag(m / gc);

else  % metric

    hU = h;                      % m
    g  = 9.80665;                % m/s^2

    W = m * g;                   % N
    Mdyn = diag(m);              % kg

end

Wtotal = sum(W);

%% ==========================================================
% ASSEMBLE SHEAR BUILDING 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

%% ==========================================================
% PERIOD CALCULATION
%% ==========================================================

switch periodMethod

    case 'trialRayleigh'

        phi = hU(:);
        phi = phi / sqrt(phi.'*Mdyn*phi);

        omega2 = (phi.'*K*phi) / (phi.'*Mdyn*phi);
        T = 2*pi / sqrt(omega2);

    case 'iterateMode'

        phi = hU(:);
        phi = phi / sqrt(phi.'*Mdyn*phi);

        for it = 1:200
            y = K \ (Mdyn*phi);
            phi_new = y / sqrt(y.'*Mdyn*y);

            if norm(phi_new - phi) < 1e-8
                phi = phi_new;
                break
            end

            phi = phi_new;
        end

        omega2 = (phi.'*K*phi) / (phi.'*Mdyn*phi);
        T = 2*pi / sqrt(omega2);

    case 'eig'

        [Phi,Omega2] = eig(K,Mdyn);
        w2 = real(diag(Omega2));

        [w2min,idx] = min(w2);
        phi = real(Phi(:,idx));
        phi = phi / sqrt(phi.'*Mdyn*phi);

        omega2 = w2min;
        T = 2*pi / sqrt(omega2);
end

%% ==========================================================
% ASCE Cs CALCULATION
%% ==========================================================

Tshort = SD1 / SDS;

if T <= Tshort
    Cs_raw = SDS * Ie / R;
else
    Cs_raw = SD1 * Ie / (T * R);
end

Cs_min1 = 0.044 * SDS * Ie;
Cs_min2 = 0.01  * Ie;

Cs = max([Cs_raw, Cs_min1, Cs_min2]);

V = Cs * Wtotal;

%% ==========================================================
% DESIGN RESPONSE SPECTRUM (g LEVELS)
%% ==========================================================

% Transition period
Ts = SD1 / SDS;

% Spectral acceleration at building period (in g)
if T <= Ts
    Sa_T = SDS;          % constant acceleration region
else
    Sa_T = SD1 / T;      % constant velocity region
end

% Peak spectral acceleration of two-period design spectrum
Sa_peak = SDS;           % always the maximum ordinate (in g)


%% ==========================================================
% VERTICAL FORCE DISTRIBUTION
%% ==========================================================

if T <= 0.5
    kexp = 1;
elseif T >= 2.5
    kexp = 2;
else
    kexp = 1 + (T - 0.5)/2;
end

Fx = V * (W .* hU.^kexp) / sum(W .* hU.^kexp);

%% ==========================================================
% DISPLACEMENTS
%% ==========================================================

deltaElastic = K \ Fx;

deltaDesign = (Cd/Ie) * deltaElastic;

storyHeights = [hU(1); diff(hU)];

driftElastic = [deltaElastic(1); diff(deltaElastic)];
driftDesign  = [deltaDesign(1);  diff(deltaDesign)];

driftRatioElastic = driftElastic ./ storyHeights;
driftRatioDesign  = driftDesign  ./ storyHeights;

%% ==========================================================
% DRIFT LIMIT
%% ==========================================================

switch riskCategory
    case {'I','II'}
        driftLimit = 0.020;
    case 'III'
        driftLimit = 0.015;
    case 'IV'
        driftLimit = 0.010;
end

maxDriftRatio = max(driftRatioDesign);
maxDriftOK = maxDriftRatio <= driftLimit;

%% ==========================================================
% P-DELTA STABILITY COEFFICIENT
%% ==========================================================

theta = zeros(n,1);

for i = 1:n
    Vstory = sum(Fx(i:end));
    Pstory = sum(W(i:end));

    theta(i) = (Pstory * abs(driftDesign(i))) / ...
               (Vstory * storyHeights(i));
end

%% ==========================================================
% OUTPUT STRUCT
%% ==========================================================

out = struct();

out.T = T;
out.Cs = Cs;
out.V = V;
out.kexp = kexp;

out.Fx = Fx;

out.deltaElastic = deltaElastic;
out.deltaDesign  = deltaDesign;

out.driftRatioElastic = driftRatioElastic;
out.driftRatioDesign  = driftRatioDesign;

out.driftLimit = driftLimit;
out.maxDriftRatio = maxDriftRatio;
out.maxDriftOK = maxDriftOK;

out.theta = theta;

out.phi = phi;
out.periodMethod = periodMethod;
out.riskCategory = riskCategory;

out.Sa_T = Sa_T;         % spectral acceleration at T (g)
out.Sa_peak = Sa_peak;   % peak spectrum acceleration (g)


%% ==========================================================
% PLOT TWO-PERIOD DESIGN RESPONSE SPECTRUM
%% ==========================================================

% Generate smooth period vector
Tvec = logspace(-2,1,1000)';   % 0.01 to 10 sec

% Characteristic periods
T0 = 0.2 * SD1 / SDS;
TS = SD1 / SDS;
TL = 8;   % You may replace with input if desired

% Piecewise spectrum
Sa_curve = zeros(size(Tvec));

idx1 = (Tvec < T0);
Sa_curve(idx1) = SDS .* (0.4 + 0.6 .* (Tvec(idx1) ./ T0));

idx2 = (Tvec >= T0) & (Tvec <= TS);
Sa_curve(idx2) = SDS;

idx3 = (Tvec > TS) & (Tvec <= TL);
Sa_curve(idx3) = SD1 ./ Tvec(idx3);

idx4 = (Tvec > TL);
Sa_curve(idx4) = (SD1 / TL) .* (TL ./ Tvec(idx4)).^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% ==========================================================
% PLOT TWO-PERIOD DESIGN RESPONSE SPECTRUM
%% ==========================================================

%% ==========================================================
% PLOT TWO-PERIOD DESIGN RESPONSE SPECTRUM
%% ==========================================================

% Generate smooth period vector
Tvec = logspace(-2,1,1000)';   % 0.01 to 10 sec

% Characteristic periods
T0 = 0.2 * SD1 / SDS;
TS = SD1 / SDS;
TL = 8;   % Replace with input if desired

% Piecewise spectrum
Sa_curve = zeros(size(Tvec));

idx1 = (Tvec < T0);
Sa_curve(idx1) = SDS .* (0.4 + 0.6 .* (Tvec(idx1) ./ T0));

idx2 = (Tvec >= T0) & (Tvec <= TS);
Sa_curve(idx2) = SDS;

idx3 = (Tvec > TS) & (Tvec <= TL);
Sa_curve(idx3) = SD1 ./ Tvec(idx3);

idx4 = (Tvec > TL);
Sa_curve(idx4) = (SD1 / TL) .* (TL ./ Tvec(idx4)).^2;

% ----------------------------------------------------------
% Create Figure
% ----------------------------------------------------------

figure
semilogx(Tvec, Sa_curve, 'b-', 'LineWidth', 1.6)
hold on
grid on

% Minor gridlines dotted
set(gca,'XMinorGrid','on','YMinorGrid','on')
set(gca,'MinorGridLineStyle',':')
set(gca,'GridLineStyle','-')

xlabel('Period T (sec)')
ylabel('Spectral acceleration, S_a(T) (g)')

title(sprintf(['ASCE 7-22 Two-Period DESIGN Spectrum (5%% Damping)\n' ...
               'Risk Category %s'], riskCategory))

xlim([1e-2 10])

% Slightly above highest point (10% margin)
ymax = max(Sa_curve);
ylim([0 ymax*1.10])

% ----------------------------------------------------------
% Marker at Building Period
% ----------------------------------------------------------

semilogx(T, Sa_T, 'bo', ...
    'MarkerSize', 9, ...
    'LineWidth', 2, ...
    'MarkerFaceColor', 'none')

% ----------------------------------------------------------
% Fixed Annotation (Top-Right Corner)
% ----------------------------------------------------------

txt = sprintf('T = %.3f s\nS_a = %.3f g', T, Sa_T);

annotation('textbox', ...
    [0.72 0.75 0.28 0.12], ...  % normalized [x y width height]
    'String', txt, ...
    'FontSize', 10, ...
    'Color', 'b', ...
    'BackgroundColor', 'white', ...
    'EdgeColor', 'b', ...
    'FitBoxToText', 'on');

hold off


end