% Sa_vertical  Compute and plot the vertical design response spectrum
%   per ASCE/SEI 7-22 §11.9 (Eqs. 11.9-1..3).
%
%   Sa_v = Sa_vertical(Cv, SMS)
%
%   Inputs:
%     Tv  = vector of vertical periods (s), all > 0
%     Cv  = vertical coefficient, either numeric OR letter 'A'..'E' (site class)
%     SMS = numeric MCE_R short-period spectral acceleration (g), e.g., 1.5
%
%   Output:
%     Sa_v = spectral acceleration (g)
%
%  The vertical spectrum in §11.9 uses fixed breakpoints at 
%  0.025 s, 0.05 s, 0.15 s, and 2.0 s, and then holds constant beyond 2.0 s
%
function [T,fn,Sa_v] = Sa_vertical(Cv, SMS)
%
    t_string=sprintf('Vertical Design Response Spectrum 5%% Damping \n Site Class %s (§11.9)',Cv);

    f_min = 0.1;       % minimum frequency (Hz)
    f_max = 100;       % maximum frequency (Hz)
    ratio  = 2^(1/12);  % 1/12-octave frequency ratio
    
    frequencies = f_min;
    f = f_min;
    while f * ratio <= f_max
        f = f * ratio;
        frequencies(end+1,1) = f;   %#ok<SAGROW>
    end
    
    % Ensure the last value reaches or exceeds f_max (optional)
    if frequencies(end) < f_max
        frequencies(end+1,1) = f_max;
    end

    Tv=1./frequencies;

    Tv=sortrows(Tv);


    % --- Validate/resolve inputs ---
    validateattributes(Tv, {'numeric'}, {'vector','real','finite','positive'});
    Cv  = resolveCv(Cv);
    validateattributes(SMS, {'numeric'}, {'scalar','real','finite','positive'});

    % --- Compute Sa_v (piecewise per §11.9) ---
    Sa_v = zeros(size(Tv));

    idxA = (Tv < 0.025);
    Sa_v(idxA) = 0.3 * Cv * SMS;

    idxB = (Tv >= 0.025) & (Tv < 0.05);
    Sa_v(idxB) = 20 * Cv * SMS .* (Tv(idxB) - 0.025) + 0.3 * Cv * SMS;

    idxC = (Tv >= 0.05) & (Tv < 0.15);
    Sa_v(idxC) = 0.8 * Cv * SMS;

    idxD = (Tv >= 0.15) & (Tv < 2.0);
    Sa_v(idxD) = 0.8 * Cv * SMS .* (0.15 ./ Tv(idxD)).^0.75;

    idxE = (Tv >= 2.0);
    Sa_v(idxE) = 0.8 * Cv * SMS .* (0.15 / 2.0).^0.75;

    % === Plotting four formats ===

    fig_num=1;

    y_label='Spectral acceleration, Sa_v(T_v) (g)';
   
    f_v = 1 ./ Tv;

    ppp=[Tv Sa_v];

    T=Tv;
    fn=f_v;


    x_label='Period T (sec)';

    [fig_num,~]=plot_loglog_function_h2(fig_num,x_label,y_label,t_string,ppp,min(Tv),max(Tv));
    [fig_num,~]=plot_loglin_function_h2(fig_num,x_label,y_label,t_string,ppp,min(Tv),max(Tv));

    ppp=flipud([f_v Sa_v]);

    x_label='Natural Frequency fn (Hz)';
    
    [fig_num,~]=plot_loglog_function_h2(fig_num,x_label,y_label,t_string,ppp,min(f_v),max(f_v));
    [fig_num,~]=plot_loglin_function_h2(fig_num,x_label,y_label,t_string,ppp,min(f_v),max(f_v)); 

end

% --------- Helper ---------
function Cv_num = resolveCv(Cv_in)
    if isnumeric(Cv_in)
        Cv_num = Cv_in;
        validateattributes(Cv_num, {'numeric'}, {'scalar','positive','finite'});
        return;
    end
    if ~(ischar(Cv_in) || isstring(Cv_in))
        error('Cv must be numeric or a site-class letter A..E.');
    end
    key = upper(string(strtrim(Cv_in)));
    switch key
        case "A", Cv_num = 0.9;
        case "B", Cv_num = 1.0;
        case "C", Cv_num = 1.1;
        case "D", Cv_num = 1.2;
        case "E", Cv_num = 1.3;
        otherwise
            error('Unrecognized Cv letter "%s". Use A, B, C, D, or E.', key);
    end
end
