%
%  Run getDesignSpectrumForLocation.m before running the following script.
%
function [T,fn,Sa] = Sa_multiPeriod_fromSpecData(specData)
% Sa_multiPeriod_fromSpecData  Compute and plot the Multi-Period horizontal design spectrum
%   using specData (from getDesignSpectrumForLocation) per §11.4.5.1 of ASCE/SEI 7-22.
%
%   Sa = Sa_multiPeriod_fromSpecData(T, specData)
%
%   Inputs:
%     T        = vector of periods (s) where you want Sa
%     specData = struct returned by getDesignSpectrumForLocation, with fields:
%                   .periods    – vector of discrete periods (s)
%                   .ordinates  – vector of spectral accelerations (g) at those periods
%                   .TL         – long­period transition period (s)
%                   .siteClass  – e.g., 'A','B',…
%                   .city       – location name
%                   (plus other metadata)
%
%   Output:
%     Sa       = vector of spectral acceleration (g) at each period in T
%
%   Notes:
%     - The discrete values specData.periods / specData.ordinates are used
%       for exact points for T ≤ max(specData.periods).  
%     - For T > max(specData.periods) & ≤ specData.TL: tail = ordinates(end) * (10/T).  
%     - For T > specData.TL: tail = ordinates(end) * (10 * specData.TL / T^2).  
%     - Four standard plotting formats are included.
%
%   Reference: §11.4.5.1 of ASCE/SEI 7-22. :contentReference[oaicite:0]{index=0}

% Extract discrete data
Tvals  = specData.periods(:);
Saved  = specData.ordinates(:);
TL     = specData.TL;
siteClass = specData.siteClass;

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


% generate frequencies from 0.1 to 100 Hz in 1/12-octave steps

T=Tvals_sorted;


% Pre-allocate Sa
Sa = nan(size(T));

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

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

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

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


f_min = 0.1;       % minimum frequency (Hz)
f_max = 100;       % maximum frequency (Hz)
ratio  = 2^(1/6);  % 1/3-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

% Display result
disp('1/3-octave frequencies from 0.1 Hz to 100 Hz:');
disp(frequencies);

T_fine=1./frequencies;


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


% --- Plotting four formats ---

fig_num=1;

% 1. Log-Log vs Period

ppp=[T Sa];

ppp = flipud(ppp);

x_label='Period T (sec)';
y_label='Spectral acceleration Sa (g)';
t_string=sprintf('Multi-Period Spectrum  5%% Damping\n (Site %s, Class %s, Risk Category %s)', specData.city, siteClass,specData.riskCategory);
[fig_num,~]=plot_loglog_function_h2(fig_num,x_label,y_label,t_string,ppp,min(T),max(T));

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




% 3. Log-Log vs Frequency


ppp=[frequencies Sa];


x_label='Natural Frequency fn (Hz)';

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

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

fn=frequencies;

end
