% coded by Yue Hua, 5/21/14
%               Stanford University
%               yuehua@stanford.edu
% 
% NGA-West2 Model for Estimating Average Horizontal Values of Pseudo-Absolute 
% Spectral Accelerations Generated by Crustal Earthquakes
% Idriss, 2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT
%
%   M               = moment magnitude
%   Vs30(m/sec)     = average shear wave velocity over the top 30m below
%   the ground surface
%   T               = period of vibration
%                 Use 1000 for output the array of Sa with period
%   R_RUP           = closest distance to fault rupture (rupture distance)(km)
%   F               = 0 strike slip & other non-reverse faulting
%                   = 1 reverse
%
% OUTPUT   
%
%   Sa              = median spectral acceleration prediction
%   sigma           = logarithmic standard deviation of spectral acceleration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Sa sigma period1] = I_2014_nga(M, Vs30, T, Rrup, F)
%
% for the given period T, get the index for the constants
period   = [ 0.01	0.02	0.03	0.05	0.075	0.1	0.15	0.2	0.25	0.3	0.4	0.5	0.75	1	1.5	2	3	4	5	7.5	10	0];

if Vs30 > 1200
    Vs30 = 1200;
end

if length (T) == 1 && T == 1000; % Compute Sa and sigma with pre-defined period
    Sa=zeros(1,length(period)-1);
    sigma=zeros(1,length(period)-1);
    period1=period(1:end-1);
    for ip=1:length(period)-1
        [Sa(ip),sigma(ip)] = I_2014_sub(M, ip, Rrup, Vs30,F);
    end
else                            % Compute Sa and sigma with user-defined period
    Sa=zeros(1, length(T));
    sigma=zeros(1, length(T));
    period1=T;
    for i=1:length(T)
        Ti = T(i);
        if ( isempty(find(abs(period-Ti) < 0.0001))) % The user defined period requires interpolation
            T_low = max(period(period < Ti));
            T_high = min(period(period > Ti));
            ip_low  = find(period==T_low);
            ip_high = find(period==T_high);
            
            [Sa_low, sigma_low] = I_2014_sub(M, ip_low, Rrup, Vs30,F);
            [Sa_high, sigma_high] = I_2014_sub(M, ip_high, Rrup, Vs30,F);
            x = [log(T_low) log(T_high)];
            Y_sa = [log(Sa_low) log(Sa_high)];
            Y_sigma = [sigma_low sigma_high];
            Sa(i) = exp(interp1(x, Y_sa, log(Ti)));
            sigma(i) = interp1(x, Y_sigma, log(Ti));
        else
            ip_T = find(abs((period- Ti)) < 0.0001);
            [Sa(i), sigma(i)] = I_2014_sub(M, ip_T, Rrup, Vs30,F);
        end
    end
end

function [Sa, sigma]=I_2014_sub(M, ip, R_RUP, Vs30,F)
% Constant
period = [ 0.01	0.02	0.03	0.05	0.075	0.1	0.15	0.2	0.25	0.3	0.4	0.5	0.75	1	1.5	2	3	4	5	7.5	10	0];

if M <= 6.75  % VS30 > 450
a1=[7.0887	7.1157	7.2087	6.2638	5.9051	7.5791	8.0190	9.2812	9.5804	9.8912	9.5342	9.2142	8.3517	7.0453	5.1307	3.3610	0.1784	-2.4301	-4.3570	-7.8275	-9.2857	7.0887];
a2=[0.2058	0.2058	0.2058	0.0625	0.1128	0.0848	0.1713	0.1041	0.0875	0.0003	0.0027	0.0399	0.0689	0.1600	0.2429	0.3966	0.7560	0.9283	1.1209	1.4016	1.5574	0.2058];
b1=[2.9935	2.9935	2.9935	2.8664	2.9406	3.0190	2.7871	2.8611	2.8289	2.8423	2.8300	2.8560	2.7544	2.7339	2.6800	2.6837	2.6907	2.5782	2.5468	2.4478	2.3922	2.9935];
b2=[-0.2287	-0.2287	-0.2287	-0.2418	-0.2513	-0.2516	-0.2236	-0.2229	-0.2200	-0.2284	-0.2318	-0.2337	-0.2392	-0.2398	-0.2417	-0.2450	-0.2389	-0.2514	-0.2541	-0.2593	-0.2586	-0.2287];
elseif M > 6.75
a1	=[ 9.0138	9.0408	9.1338	7.9837	7.7560	9.4252	9.6242	11.1300	11.3629	11.7818	11.6097	11.4484	10.9065	9.8565	8.3363	6.8656	4.1178	1.8102	0.0977	-3.0563	-4.4387	9.0138];
a2	=[ -0.0794	-0.0794	-0.0794	-0.1923	-0.1614	-0.1887	-0.0665	-0.1698	-0.1766	-0.2798	-0.3048	-0.2911	-0.3097	-0.2565	-0.2320	-0.1226	0.1724	0.3001	0.4609	0.6948	0.8393	-0.0794];
b1	=[ 2.9935	2.9935	2.9935	2.7995	2.8143	2.8131	2.4091	2.4938	2.3773	2.3772	2.3413	2.3477	2.2042	2.1493	2.0408	2.0013	1.9408	1.7763	1.7030	1.5212	1.4195	2.9935];
b2	=[ -0.2287	-0.2287	-0.2287	-0.2319	-0.2326	-0.2211	-0.1676	-0.1685	-0.1531	-0.1595	-0.1594	-0.1584	-0.1577	-0.1532	-0.1470	-0.1439	-0.1278	-0.1326	-0.1291	-0.1220	-0.1145	-0.2287];
end

a3	=[ 0.0589	0.0589	0.0589	0.0417	0.0527	0.0442	0.0329	0.0188	0.0095	-0.0039	-0.0133	-0.0224	-0.0267	-0.0198	-0.0367	-0.0291	-0.0214	-0.0240	-0.0202	-0.0219	-0.0035	0.0589];
zeta	=[ -0.854	-0.854	-0.854	-0.631	-0.591	-0.757	-0.911	-0.998	-1.042	-1.030	-1.019	-1.023	-1.056	-1.009	-0.898	-0.851	-0.761	-0.675	-0.629	-0.531	-0.586	-0.854];
gamma	=[ -0.0027	-0.0027	-0.0027	-0.0061	-0.0056	-0.0042	-0.0046	-0.0030	-0.0028	-0.0029	-0.0028	-0.0021	-0.0029	-0.0032	-0.0033	-0.0032	-0.0031	-0.0051	-0.0059	-0.0057	-0.0061	-0.0027];
phi	=[ 0.08	0.08	0.08	0.08	0.08	0.08	0.08	0.08	0.08	0.08	0.08	0.08	0.08	0.06	0.04	0.02	0.02	0	0	0	0	0.08];

ln_PSA = a1(ip) + a2(ip) * M + a3(ip)* (8.5 - M)^2 - (b1(ip) + b2(ip) * M) * log(R_RUP + 10) + zeta(ip) * log(Vs30) + gamma(ip) * R_RUP + phi(ip)*F;
Sa = exp(ln_PSA);

M_SE = max(min(M, 7.5),5);

if period(ip) < 0.05
sigma = 1.18 + 0.035 * log(0.05) - 0.06 * M_SE;
elseif period(ip) > 3
sigma = 1.18 + 0.035 * log(3) - 0.06 * M_SE;
else
sigma = 1.18 + 0.035 * log(period(ip)) - 0.06 * M_SE;
end
    





