% Sabetta & Pugliese ground motion prediction model, 1996
% Marko Brozovi�
% 2.8.2011
%
% Useful for 4.6<=M<=6.8 and R<=100 km
% Only for larger horizontal component
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input variables
% T = period (sec)
% M = magnitude
% R = fault distance (km)
% Site Geology
% S1: 1-shallow alluvium, 0-otherwise; 1-plitva naplavina H<20m, 0-skala(400  < stri�na hitrost < 800m/s), 0-skala(stri�na hitrost > 800m/s)
% S2: 1-deep alluvium, 0-otherwise 1-globoka naplavina H>20m, 0-skala(400  < stri�na hitrost < 800m/s), 0-skala(stri�na hitrost > 800m/s)
% coef: 1-fault distance, raw; 2-fault distance, smooth; 3-epicentral
%       distance, raw; 4-epicentral distance, smooth;
%
% Output variables
% Sa: median spectral acceleration prediction
% sigma: logarithmic standard deviation of spectral acceleration
%          prediction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Sa sigma period]=SP_1996(T,M,R,S1,S2,coef)


% Period and frequency
if T<0.04
    T=0.04;
end
period=T; % (s)
f=1/T; % (Hz)



% Regression coefficients
if coef==1
    % fault distance, raw coefficients (Table 1)
    fSP=[0.25 0.33 0.5 0.67 1 1.33 2 2.5 3.33 5 6.67 10 15 25]; % Frequency (Hz)
    aSP=[-2.32 -2.017 -1.935 -1.498 -1.112 -0.753 -0.496 -0.077 0.467 0.489 0.499 0.29 -0.033 -0.505]; % Constant term a
    bSP=[0.653 0.637 0.667 0.614 0.56 0.503 0.458 0.397 0.301 0.283 0.255 0.244 0.249 0.275]; % Magnitude coefficient b
    cSP=-1*ones(1,14); % Distance coefficient c
    e1SP=[0 0 0 0 0 0.12 0.24 0.196 0.122 0.103 0.135 0.168 0.112 0.147]; % Site coefficient e1
    e2SP=[0.134 0.149 0.194 0.2 0.242 0.232 0.156 0 0 0 0 0 0 0]; % Site coefficient e2
    hSP=[2.3 2.3 3.1 4 4.2 4.9 4.6 4.9 6.1 5.7 6.5 7.3 7.2 5.8]; % h
    sigmaSP=log(10.^([0.275 0.317 0.3 0.242 0.246 0.261 0.279 0.268 0.245 0.215 0.207 0.2 0.181 0.182])); % sigma
elseif coef==2
    % fault distance, smooth coefficients (Table 1)
    fSP=[0.25 0.33 0.5 0.67 1 1.33 2 2.5 3.33 5 6.67 10 15 25]; % Frequency (Hz)
    aSP=[-2.4 -2.17 -1.8 -1.51 -1.12 -0.85 -0.4 -0.077 0.4 0.55 0.5 0.29 -0.035 -0.505]; % Constant term a
    bSP=[0.685 0.675 0.65 0.62 0.57 0.53 0.455 0.4 0.315 0.273 0.255 0.245 0.255 0.273]; % Magnitude coefficient b
    cSP=-1*ones(1,14); % Distance coefficient c
    e1SP=[0 0 0 0.01 0.05 0.12 0.22 0.21 0.165 0.13 0.13 0.13 0.13 0.13]; % Site coefficient e1
    e2SP=[0.13 0.151 0.184 0.21 0.242 0.232 0.156 0.115 0.065 0 0 0 0 0]; % Site coefficient e2
    hSP=[2.1 2.5 3.1 3.5 4.0 4.4 5.0 5.3 5.7 6.3 6.7 7.3 7.2 5.8]; % h
    sigmaSP=log(10.^([0.297 0.297 0.297 0.295 0.29 0.283 0.271 0.264 0.245 0.215 0.2 0.19 0.182 0.182])); % sigma
elseif coef==3
    % epicentral distance, raw coefficients (Table 2)
    fSP=[0.25 0.33 0.5 0.67 1 1.33 2 2.5 3.33 5 6.67 10 15 25]; % Frequency (Hz)
    aSP=[-2.409 -2.146 -2.082 -1.647 -1.297 -1.006 -0.65 -0.281 0.202 0.296 0.222 -0.019 -0.312 -0.817]; % Constant term a
    bSP=[0.684 0.675 0.706 0.654 0.605 0.557 0.498 0.442 0.354 0.326 0.311 0.304 0.304 0.336]; % Magnitude coefficient b
    cSP=-1*ones(1,14); % Distance coefficient c
    e1SP=[0 0 0 0 0 0.12 0.249 0.219 0.147 0.126 0.161 0.194 0.137 0.174]; % Site coefficient e1
    e2SP=[0.111 0.109 0.157 0.165 0.209 0.208 0.124 0 0 0 0 0 0 0]; % Site coefficient e2
    hSP=[2.7 2.7 3.4 4.4 4.2 4.4 5.1 4.8 5.5 5.8 5.8 6.3 6.5 4.7]; % h
    sigmaSP=log(10.^([0.303 0.344 0.32 0.27 0.271 0.278 0.301 0.285 0.26 0.238 0.22 0.214 0.194 0.195])); % sigma
else
    % epicentral distance, smooth coefficients (Table 2)
    fSP=[0.25 0.33 0.5 0.67 1 1.33 2 2.5 3.33 5 6.67 10 15 25]; % Frequency (Hz)
    aSP=[-2.5 -2.25 -1.9 -1.647 -1.28 -1 -0.595 -0.281 0.1 0.296 0.222 -0.019 -0.312 -0.817]; % Constant term a
    bSP=[0.725 0.715 0.687 0.66 0.612 0.57 0.5 0.445 0.377 0.323 0.31 0.304 0.304 0.33]; % Magnitude coefficient b
    cSP=-1*ones(1,14); % Distance coefficient c
    e1SP=[0 0 0 0.01 0.05 0.12 0.23 0.222 0.185 0.161 0.161 0.161 0.161 0.161]; % Site coefficient e1
    e2SP=[0.1 0.108 0.15 0.175 0.208 0.19 0.124 0.078 0.02 0 0 0 0 0]; % Site coefficient e2
    hSP=[2.6 3 3.6 4 4.4 4.7 5 5.2 5.4 5.7 5.9 6.2 6.3 4.7]; % h
    sigmaSP=log(10.^([0.319 0.319 0.319 0.315 0.308 0.303 0.29 0.28 0.26 0.234 0.22 0.208 0.2 0.195])); % sigma
end
% Interpolation of coefficients

a=interp1(fSP,aSP,f);
b=interp1(fSP,bSP,f);
c=-1;
e1=interp1(fSP,e1SP,f);
e2=interp1(fSP,e2SP,f);
h=interp1(fSP,hSP,f);
sigma=interp1(fSP,sigmaSP,f);

% Calculation of spectral acceleration
g=9.81; % gravitational acceleration (m/s2)
Sa=2*pi/T/100/g*10^(a+b*M+c*log10(sqrt(R^2+h^2))+e1*S1+e2*S2); % (g)