function [Sa sigma period1] = I_2008_nga(M, Vs30, T, Rrup, F)
%
% by Yoshifumi Yamamoto,  5/8/10
%       Stanford University, yama4423@stanford.edu
%
% An NGA Empirical Model for Estimating the Horizontal Spectral Values
% Generated By Shallow Crustal Earthquakes
% I. M. Idriss
% Earthquake Spectra, Volume 24, No.1, pages 217-242, February 2008
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT
%
%   M               = moment magnitude
%   Vs30
%   T               = period of vibration
%                 Use 1000 for output the array of Sa with period
%   Rrup            = closest distance to fault rupture (rupture distance)(km)
%   F               = Flag for normal faulting earthquakes
%                   = 0 strike slip
%                   = 1 reverse
%
% OUTPUT   
%
%   Sa              = median spectral acceleration prediction
%   sigma           = logarithmic standard deviation of spectral acceleration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% for the given period T, get the index for the constants
period   = [   0.01,    0.02,    0.03,    0.04,    0.05,    0.06,    0.08,     0.1,    0.15,     0.2,    0.25,     0.3,    0.35,     0.4,    0.45,     0.5,     0.6,     0.7,     0.8,     0.9,       1,     1.5,       2,       3,       4,       5,       6,       7,       8,       9,      10];
    
nT=length(T);
iflg=0;
if(nT==1);
    if(T==1000);
        iflg=1;
        nperi=length(period);
        Sa=zeros(1,nperi);
        sigma=zeros(1,nperi);
        period1=period(1:end);
        
        for index=1:1:nperi;
            % get constants for the given index value
            V = get_idriss_constants(index,M,Vs30);
            Sa(index) = exp(calc_val(M, Rrup, F, V));
            [sigma(index)] = idriss_sigma(M, V);
        end;
    end;
end;

if(iflg==0);
    Sa=zeros(1,nT);
    sigma=zeros(1,nT);
    period1=T;
    for it=1:1:nT;
        if T(it)==0; T(it)=0.01; end;
        Teach=T(it);
        if Teach>period(end); Teach = period(end); end;
    % interpolate between periods if neccesary  
        if (length(find(period == Teach)) == 0);
            T_low = max(period(find(period<Teach)));
            T_hi = min(period(find(period>Teach)));

            [sa_low, sigma_low] = I_2008_nga(M, Vs30, T_low, Rrup, F);
            [sa_hi, sigma_hi] = I_2008_nga(M, Vs30, T_hi, Rrup, F);

            x = [log(T_low) log(T_hi)];
            Y_sa = [log(sa_low) log(sa_hi)];
            Y_sigma = [sigma_low sigma_hi];
            Sa(it) = exp(interp1(x,Y_sa,log(Teach)));
            sigma(it) = interp1(x,Y_sigma,log(Teach));
            
        else
            index = find(abs((period - Teach)) < 0.0001); % Identify the period

            % get constants for the given index value
            V = get_idriss_constants(index,M,Vs30);
            Sa(it) = exp(calc_val(M, Rrup, F, V));
            [sigma(it)] = idriss_sigma(M, V);

        end;
    end;
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    local functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
function [X] = calc_val(M, Rrup, F, constants)
% calculate predicted value

	X = constants.a1 + constants.a2*M - (constants.b1+constants.b2*M)*log(Rrup+10) ...
        + constants.gamma*Rrup + constants.phi*F;
  
function [constants] = get_idriss_constants(index,M,Vs30)
% get relevant constants

% arrays with values by index
period   = [   0.01,    0.02,    0.03,    0.04,    0.05,    0.06,    0.08,     0.1,    0.15,     0.2,    0.25,     0.3,    0.35,     0.4,    0.45,     0.5,     0.6,     0.7,     0.8,     0.9,       1,     1.5,       2,       3,       4,       5,       6,       7,       8,       9,      10];
a1_m1    = [ 3.7066,  3.7066,  3.7566,  3.8066,  4.1248,  4.4681,  4.4853,  4.4592,  3.4793,  3.2354,  2.7628,  2.3813,  2.0302,  1.7037,   1.394,  1.0893,  0.5308,   0.024, -0.4141, -0.8184, -1.2290, -2.9168, -4.2783, -6.2431, -7.6967, -8.8110, -9.7232,-10.4706,-11.0814,-11.5896,-12.0149];
a2_m1    = [-0.1252, -0.1252, -0.1252, -0.1252, -0.1781, -0.2228, -0.1949, -0.1624, -0.0188,  0.0346,  0.0791,  0.1187,  0.1545,  0.1873,  0.2177,  0.2461,  0.2979,  0.3443,  0.3866,  0.4255,  0.4615,  0.6103,  0.7246,  0.8935,  1.0137,  1.1027,  1.1696,  1.2197,  1.2566,  1.2826,  1.2995];
b1_m1    = [ 2.9832,  2.9832,  2.9832,  2.9832,  3.0156,  3.0708,  3.1071,  3.1212,  2.8609,  2.8739,  2.8203,  2.8126,  2.8056,  2.7992,  2.7932,  2.7876,  2.7772,  2.7677,   2.759,   2.751,  2.7434,  2.7112,  2.6851,  2.6437,   2.611,  2.5839,  2.5607,  2.5406,  2.5228,   2.507,  2.4928];
b2_m1    = [-0.2339, -0.2339, -0.2339, -0.2339, -0.2445, -0.2536, -0.2576, -0.2570, -0.2267, -0.2282, -0.2292, -0.2301, -0.2309, -0.2317, -0.2324, -0.2330, -0.2342, -0.2353, -0.2363, -0.2373, -0.2381, -0.2418, -0.2447, -0.2493, -0.2529, -0.2558, -0.2582, -0.2603, -0.2621, -0.2636, -0.2650];
gamma_m1 = [0.00047, 0.00047, 0.00047, 0.00047,       0,       0,       0,       0,       0,       0,-0.00049, 0.00052, 0.00099, 0.00112, 0.00114, 0.00132, 0.00154,  0.0017, 0.00152, 0.00157, 0.00188,  0.0025, 0.00268,  0.0005,-0.00248,-0.00453,-0.00566,-0.00633,-0.00671,-0.00689,-0.00709];
phi_m1   = [   0.12,    0.12,   0.12,   0.12,      0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.11,     0.1,    0.06,    0.04,       0,       0,       0,       0,       0,       0,       0,       0];
a1_m2    = [ 5.6315,  5.6315,  5.6815,  5.7315,  5.8447,  6.0362,  6.4307,  6.3053,  5.0845,  5.0842,  4.5453,  4.2719,  4.0174,  3.7792,  3.5519,  3.3235,  2.9047,  2.5222,  2.1972,  1.8971,  1.5822,  0.2888, -0.7737, -2.3037, -3.4564, -4.3563, -5.1145, -5.7538, -6.2921, -6.7588, -7.1679];
a2_m2    = [-0.4104, -0.4104, -0.4104, -0.4104, -0.4329, -0.4551, -0.4831, -0.4359, -0.2566, -0.2393, -0.1850, -0.1614, -0.1399, -0.1202, -0.1020, -0.0849, -0.0538, -0.0258, -0.0003,  0.0232,   0.045,  0.1354,  0.2054,  0.3099,  0.3855,  0.4427,  0.4868,  0.5209,  0.5471,  0.5669,  0.5814];
b1_m2    = [ 2.9832,  2.9832,  2.9832,  2.9832,  2.9487,  2.9494,  2.9788,  2.9153,  2.4829,  2.5066,  2.3687,  2.3475,  2.3284,  2.3105,  2.2937,  2.2793,  2.2507,   2.225,  2.2014,  2.1786,  2.1588,   2.072,  2.0027,  1.8938,  1.8091,  1.7401,  1.6825,  1.6327,    1.59,  1.5532,  1.5201];
b2_m2    = [-0.2339, -0.2339, -0.2339, -0.2339, -0.2346, -0.2356, -0.2386, -0.2265, -0.1707, -0.1738, -0.1623, -0.1612, -0.1602, -0.1593, -0.1584, -0.1577, -0.1562, -0.1549, -0.1537, -0.1525, -0.1515, -0.1471, -0.1436, -0.1382, -0.1341, -0.1308, -0.1281, -0.1258, -0.1239, -0.1223, -0.1209];
gamma_m2 = [0.00047, 0.00047, 0.00047, 0.00047,       0,       0,       0,       0,       0,       0,-0.00049, 0.00052, 0.00099, 0.00112, 0.00114, 0.00132, 0.00154,  0.0017, 0.00152, 0.00157, 0.00188,  0.0025, 0.00268,  0.0005,-0.00248,-0.00453,-0.00566,-0.00633,-0.00671,-0.00689,-0.00709];
phi_m2   = [   0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.12,    0.11,     0.1,    0.06,    0.04,       0,       0,       0,       0,       0,       0,       0,       0];
a1_v1m1  = [ 3.7066,  3.7066,  3.7566,  3.8066,  4.1248,  4.4681,  4.4853,  4.4592,  3.4793,  3.2354,  2.7628,  2.3813,  2.0302,  1.7037,   1.394,  1.0893,  0.5308,   0.024, -0.4141, -0.8184, -1.2290, -2.9168, -4.2783, -6.2431, -7.6967, -8.8110, -9.7232,-10.4706,-11.0814,-11.5896,-12.0149];
a1_v1m2  = [ 5.6315,  5.6315,  5.6815,  5.7315,  5.8447,  6.0362,  6.4307,  6.3053,  5.0845,  5.0842,  4.5453,  4.2719,  4.0174,  3.7792,  3.5519,  3.3235,  2.9047,  2.5222,  2.1972,  1.8971,  1.5822,  0.2888, -0.7737, -2.3037, -3.4564, -4.3563, -5.1145, -5.7538, -6.2921, -6.7588, -7.1679];
da1      = [-0.1492, -0.1492, -0.1492, -0.1492, -0.1492, -0.1657, -0.1922, -0.2126, -0.2475, -0.2693, -0.2842, -0.2948, -0.3028, -0.3088, -0.3136, -0.3173, -0.3227, -0.3262, -0.3285, -0.3299, -0.3308, -0.3297, -0.3248, -0.3110, -0.2956, -0.2797, -0.2638, -0.2481, -0.2327, -0.2175, -0.2026];
a1_v2m1  = [ 3.5574,  3.5574,  3.6074,  3.6574,  3.9756,  4.3024,  4.2931,  4.2466,  3.2318,  2.9661,  2.4786,  2.0865,  1.7274,  1.3949,  1.0804,   0.772,  0.2081, -0.3022, -0.7426, -1.1483, -1.5598, -3.2465, -4.6031, -6.5541, -7.9923, -9.0907, -9.9870,-10.7187,-11.3141,-11.8071,-12.2175];
a1_v2m2  = [ 5.4823,  5.4823,  5.5323,  5.5823,  5.6955,  5.8705,  6.2385,  6.0927,   4.837,  4.8149,  4.2611,  3.9771,  3.7146,  3.4704,  3.2383,  3.0062,   2.582,   2.196,  1.8687,  1.5672,  1.2514, -0.0409, -1.0985, -2.6147, -3.7520, -4.6360, -5.3783, -6.0019, -6.5248, -6.9763, -7.3705];


constants.period = period(index);
if Vs30>=450 && Vs30<=900 && M<=6.75;
    constants.a1    = a1_m1(index);
    constants.a2    = a2_m1(index);
    constants.b1    = b1_m1(index);
    constants.b2    = b2_m1(index);
    constants.gamma = gamma_m1(index);
    constants.phi   = phi_m1(index);
elseif Vs30>=450 && Vs30<=900 && M<=8.5;
    constants.a1    = a1_m2(index);
    constants.a2    = a2_m2(index);
    constants.b1    = b1_m2(index);
    constants.b2    = b2_m2(index);
    constants.gamma = gamma_m2(index);
    constants.phi   = phi_m2(index);
end;



function [sigma] = idriss_sigma(M, V)
% calculate the sigma

Mt=max(5,M);
Mt=min(7.5,Mt);
peri=max(0.05,V.period);
peri=min(3,peri);
sigma = 1.28+0.05*log(peri) - 0.08*Mt;

