% Created by: Mohsen Kohrangi, 23Sep2014
% GMPM of average spectral acceleration
% The code will get as input the selected periods, 
% Magnitude, distance and all other parameters of the 
% selected GMPM (e.g. Boore & Atkinson, 2008 used here) and 
% will return the median and logarithmic Spectral 
% acceleration of the product of the Spectral
% accelerations at selected periods; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Variables
% T = Array of interested Periods (sec); 
% preferred GMPM
% 1: BA_2008_nga
% 2: CY_2008_nga
%
% Output Variables
% Sa: Median average spectral acceleration prediction
% sigma: logarithmic standard deviation of average 
% spectral acceleration prediction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Sa sigma] = Sa_avg (M, T, Rjb, Fault_Type, Vs30)

n = length(T);
SPa_Med = zeros(n,1);
SPa_STD = zeros(n,1);

for i = 1:n
[SPa_Med(i) SPa_STD(i)] = BA_2008_nga (M, T(i), Rjb, Fault_Type, Vs30);
    % convert to sigma_arb
    % One should uncomment this line if the arbitary component is used for
    % record selection.
%   ro_xy=0.79-0.023*log(T(i)); 
    ro_xy = 1;
    SPa_STD(i) = (SPa_STD(i).^2*(2./(1+ro_xy))).^0.5;
end
MoC = zeros(n);
for i = 1:n
       for j = 1:n;
       [rho] = baker_jayaram_correlation(T(i), T(j));  
       MoC (i,j) = rho;
       end
end

SPa_avg_meanLn = (1/n) *sum(log(SPa_Med)); % logarithmic mean of  Sa,avg

SPa_avg_STD = 0;
for i = 1:n
       for j = 1:n;
       SPa_avg_STD     = SPa_avg_STD  + (MoC(i,j) *SPa_STD(i)*SPa_STD(j)); % logarithmic Var of the Sa,avg
       end
end
SPa_avg_STD = (1/n).^2*SPa_avg_STD;
% compute median of average spectral acceleration
% and logarithmic standard deviation of 
% spectral acceleration prediction
Sa       = exp(SPa_avg_meanLn);
sigma = sqrt(SPa_avg_STD);
end