%
%  Run getDesignSpectrumASCE7_22_LatLon.m before running this script.
%
function [T, fn, Sa] = Sa_multiPeriod_fromSpecData_VAFB(specData)
% Sa_multiPeriod_fromSpecData_LatLon
%
% Compute and plot the Multi-Period horizontal design spectrum
% per §11.4.5.1 of ASCE/SEI 7-22 using LATITUDE/LONGITUDE metadata.
%
% Inputs:
%   specData = struct returned by getDesignSpectrumASCE7_22_LatLon
%              Required fields:
%                 .periods
%                 .ordinates
%                 .TL
%                 .siteClass
%                 .riskCategory
%                 .latitude
%                 .longitude
%
% Outputs:
%   T   = period vector (sec)
%   fn  = frequency vector (Hz)
%   Sa  = spectral acceleration (g)

%% ---------------- Extract Data ----------------

Tvals     = specData.periods(:);
Saved     = specData.ordinates(:);
TL        = specData.TL;
siteClass = specData.siteClass;
riskCat   = specData.riskCategory;
lat       = specData.latitude;
lon       = specData.longitude;

if isempty(Tvals) || isempty(Saved)
    error('specData does not contain multi-period spectrum data.');
end

%% ---------------- Sort Data ----------------

[Tvals_sorted, sortIdx] = sort(Tvals);
Saved_sorted = Saved(sortIdx);

T = Tvals_sorted;

%% ---------------- Build Spectrum ----------------

Sa = nan(size(T));

% 1) Interpolate for T ≤ max(periods)
maskInterp = (T <= Tvals_sorted(end));
Sa(maskInterp) = interp1(Tvals_sorted, Saved_sorted, ...
                         T(maskInterp), 'linear');

% 2) Tail: T > max(periods) & ≤ TL
maskTail1 = (T > Tvals_sorted(end)) & (T <= TL);
Sa(maskTail1) = Saved_sorted(end) .* (10 ./ T(maskTail1));

% 3) Tail: T > TL
maskTail2 = (T > TL);
Sa(maskTail2) = Saved_sorted(end) .* ...
                (10 .* TL ./ (T(maskTail2).^2));

% Remove T=0 if present
if T(1) == 0
    T(1)  = [];
    Sa(1) = [];
end

%% ---------------- Generate 1/3-Octave Frequencies ----------------

f_min = 0.1;
f_max = 100;
ratio = 2^(1/6);     % 1/3 octave

frequencies = f_min;
f = f_min;

while f * ratio <= f_max
    f = f * ratio;
    frequencies(end+1,1) = f; %#ok<SAGROW>
end

if frequencies(end) < f_max
    frequencies(end+1,1) = f_max;
end

T_fine = 1 ./ frequencies;

%% ---------------- Resample to 1/3-Octave Grid ----------------

[T, Sa] = interpCubicSpline2D(T, Sa, T_fine);

fn = frequencies;

%% ---------------- Plotting ----------------

fig_num = 1;

% Reverse for period plots
ppp = flipud([T Sa]);

x_label = 'Period T (sec)';
y_label = 'Spectral Acceleration Sa (g)';

t_string = sprintf(['Multi-Period Design Spectrum  5%% Damping\n' ...
                    'Vandenberg AFB  |  Site %s  |  Risk %s'], ...
                    siteClass, riskCat);

% --- Log-Log vs Period
[fig_num,~] = plot_loglog_function_h2(fig_num, ...
                x_label, y_label, t_string, ...
                ppp, min(T), max(T));

% --- Log-Lin vs Period
[fig_num,~] = plot_loglin_function_h2(fig_num, ...
                x_label, y_label, t_string, ...
                ppp, min(T), max(T));

% --- Frequency Plots
ppp = [fn Sa];
x_label = 'Natural Frequency fn (Hz)';

% Log-Log vs Frequency
[fig_num,~] = plot_loglog_function_h2(fig_num, ...
                x_label, y_label, t_string, ...
                ppp, min(fn), max(fn));

% Log-Lin vs Frequency
[fig_num,~] = plot_loglin_function_h2(fig_num, ...
                x_label, y_label, t_string, ...
                ppp, min(fn), max(fn));

end
