function [rho] = markhvida_et_al_spatial_correlation(T1, T2, h, var_flag)

% Created by Maryia Markhvida, 08/12/2018
% Minor edits by Jack Baker, 8/20/2019
% Minor edits by Maryia Markhvida - hard coded model coefficients, 10/10/2019
%
% Compute spatial cross-correlation of epsilons for NGA-West2 ground motion 
% models, using principal components 
%
% The function is strictly empirical, fitted over the range the range 0.01s <= T1, T2 <= 5s
%
% Documentation is provided in the following document:
%   Markhvida, M., Ceferino, L., & Baker, J. W. (2018).
%   Modeling spatially correlated spectral accelerations at multiple periods
%   using principal component analysis and geostatistics. Earthquake
%   Engineering & Structural Dynamics, 47(5), 1107-1123.
%
% INPUT
%
%   T1, T2      = The two periods of interest. The periods may be equal,
%                 and there is no restriction on which one is larger.
%
%   h           = [1 x n] vector of separation distances between two sites (units of km)
%
%  (optional)
%   var_flag    = 1  assumes that variance (covariance) of the random field
%                    is the sill of the nested exponential semi-variogram.
%                    This is in line with the assumption made in the PCA
%                    spatial cross-correlation model (Markhvida et al, 2018)
%               = 0  assumes that variance of the random field is 1 if T1 == T2.
%                    If T1 /= T2, then the covariance is assumed to be the
%                    cross-variance at 125km. This is the median distance at which
%                    semi-variance across different periods reaches 1.
%
% OUTPUT
%
%   rho         = [1 x n] vector with the predicted correlation coefficients

% Use all 19 principal components by default
nPCs = 19;

% Define the period vector, T, in accordance with the coefficients
T = [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];

PCAcoefs = get_PCA_coefs();
modelVario = get_model_varios();

% Check the validity of input arguments
if min(T1,T2)<0.01
    error('The periods must be greater or equal to 0.01s')
end
if max(T1,T2)>5
    error('The periods must be less than or equal to 5s')
end
if any(h<0)
    error('The separation distance must be positive')
end


% Determine PCA coefficients for periods of interest
PCAcoefs_T1 = interp1(T,PCAcoefs,T1);
PCAcoefs_T2 = interp1(T,PCAcoefs,T2);

% Get the theoretical variance of different principal components
for i = 1:nPCs
    [~,varVector(i)] =  getSemiVariance(0,modelVario{i});
end

% Get the theoretical variance for normalized residuals of different
% periods
var_T1 = (PCAcoefs_T1.*PCAcoefs_T1)* varVector';
var_T2 = (PCAcoefs_T2.*PCAcoefs_T2)* varVector';


% Get the auto-covariance and auto-correlation for all distances in vector
% h and periods T1 and T2
for j = 1:length(h)
    % Get semi-variances and variances for difference principal components
    for i = 1:nPCs
        [semivarVector(i),var_PCA(i)] =  getSemiVariance(h(j),modelVario{i});
    end
    
    % Convert the semi-variance and variance of principal components into
    % semi-variance (or cross-variance) and variance (or covariance) in the
    % original basis
    semivar_o(j) = (PCAcoefs_T1.*PCAcoefs_T2)*semivarVector';
    var_o(j)     = (PCAcoefs_T1.*PCAcoefs_T2)*var_PCA';
    
    
    % calculate the auto-covariance and auto-correlation in the original basis
    if (nargin < 4 || var_flag == 1)
        cov_o(j) = var_o(j) - semivar_o(j);
        rho(j) = cov_o(j)/(var_T1*var_T2)^0.5;
        
    else
        if T1 == T2
            cov_T1_T2 = 1;
        else
            h_cutoff = 125;
            for i_pca = 1:nPCs
                [crossvarVector(i_pca),~] =  getSemiVariance(h_cutoff,modelVario{i_pca});
            end
            cov_T1_T2 = (PCAcoefs_T1.*PCAcoefs_T2)*crossvarVector';
        end
        cov_o(j) = max(cov_T1_T2 - semivar_o(j),0); % calculate covariance > 0
        rho(j) = cov_o(j)/(1*1)^0.5;       
    end
       
end

end

% Additional helper functions to calculate the semi-variance/cross-variance
function [SEMIVAR,var] = getSemiVariance(DIST,modelVario)
switch modelVario.Type
    case 'iso nest'
        [SEMIVAR,var] = getIsoNestedVario(modelVario,DIST);
    case 'nug'
        [SEMIVAR,var] = getNugVario(modelVario,DIST);
    otherwise
        disp('NOT RIGHT MODEL')
end
end

function [SEMIVAR,var] = getIsoNestedVario(varModel,DIST)
var = varModel.Cn+varModel.C1+varModel.C2;
SEMIVAR = varModel.Cn+varModel.C1*(1-exp(-3*DIST/varModel.a1))+varModel.C2*(1-exp(-3*DIST/varModel.a2));
SEMIVAR(DIST == 0) = 0;
end

function [SEMIVAR,var] = getNugVario(varModel,DIST)
var = varModel.Cn;
SEMIVAR = var;
SEMIVAR(DIST == 0) = 0;
end

% Additional function to extract PCA coefficients
function PCAcoefs = get_PCA_coefs()
PCAcoefs = [0.270963956240109,-0.139418157111539,0.0690420060759825,-0.106094866059058,-0.0922880747536406,-0.113489976128860,-0.188935371416588,0.153956801827048,-0.160082932478587,-0.0485878661500656,0.106169114128661,0.0545367125260002,-0.0842347288819985,0.00206507178166717,0.233666515544382,-0.0444106080542835,-0.298766213268524,-0.527588859528322,-0.580349072958488;...
    0.270185457409379,-0.141734438837191,0.0770156687411159,-0.116393534099990,-0.103464378298138,-0.124082463392091,-0.199840301009291,0.155452551301531,-0.157024101304166,-0.0511781532022337,0.102685985245256,0.0534091781987030,-0.0785807732589693,0.00538047460862058,0.220317828901249,-0.0394525931272570,-0.257172668625987,-0.150994889177740,0.781868928002583;...
    0.266716484131893,-0.150918021372557,0.101241750461614,-0.144620230365225,-0.128327845057672,-0.150413273486784,-0.217519114614170,0.154533128422060,-0.144555133487480,-0.0493784913365888,0.0865380808702005,0.0369034473261574,-0.0551975904000634,0.00787249481774084,0.149651509850332,-0.0232087969820605,-0.0284597879541604,0.808901723567414,-0.226437334732122;...
    0.251688240452975,-0.184642998841202,0.178879968100826,-0.221328310597117,-0.175557525753130,-0.176668887066346,-0.188651351366600,0.0424749058524636,-0.0455090101766787,-0.0291885726697761,-0.0315764520880226,-0.0605411454951765,0.0935508235649957,0.0224423933374744,-0.299181352126819,0.0599168858817312,0.754350402692246,-0.206472972801382,0.0231093445072980;...
    0.236434540660266,-0.218922079202474,0.237254183941850,-0.234559034122274,-0.133267087790244,-0.0431828093963935,0.119447151365187,-0.272310554909179,0.238192698302355,0.100676333288767,-0.263315034080162,-0.121207177353811,0.202769694434657,0.00661936093850581,-0.493306767118467,0.116246172684757,-0.475923822815555,0.0367733765077595,-0.00643955732220411;...
    0.232994643271854,-0.228087987254185,0.230554572919478,-0.160443024111330,0.0400564218872398,0.181657484726543,0.427112684239685,-0.324579279763868,0.263780433255153,0.142634796459082,-0.0813780347714383,0.0465305509298986,-0.151801546733413,-0.0833183140197916,0.534198178434453,-0.184596749983802,0.210357916588565,-0.00284225563610220,0.00331108666921881;...
    0.238919759244457,-0.211905954003063,0.132646222385762,0.0820453502922968,0.327946972937890,0.393273105011282,0.325836316409324,0.162029624546621,-0.182164846060428,-0.138319895254022,0.470111475270552,0.177876087511963,-0.111256500094100,0.0883177907230135,-0.291143253148395,0.262494244929844,-0.00152291509197000,0.0154770015927024,0.00150080927064067;...
    0.247247200513419,-0.174053609784519,-0.00825743327819653,0.277382297057791,0.403271333843203,0.220437620135438,-0.0837312940200531,0.224796020087722,-0.171941472633753,-0.0292747130158189,-0.381524121322450,-0.237244495814365,0.356271008578838,-0.0850494996785505,-0.0125434811410974,-0.442559354727306,0.0151125187962444,0.0110964548458275,0.00194187520085148;...
    0.253677096925723,-0.122375885147353,-0.148595585559558,0.365271223153003,0.253186443805678,-0.0612442510569459,-0.283389735157485,-0.0811437931924326,0.212210668293107,0.143362426718511,-0.275727618896954,-0.0411307164923357,-0.202014525375018,0.0228459038515480,0.155257213551213,0.632145331896681,0.0455130187869974,0.000690596761983894,0.000469634928915885;...
    0.254921191707502,-0.0713194463533730,-0.237030888421264,0.359073100317565,0.0401080106567685,-0.248766601944991,-0.141859037791190,-0.286692239119504,0.300971237941779,0.0579993716202019,0.328411991836789,0.208361703130758,-0.194768507716927,0.0324295008946133,-0.258822160901248,-0.477244327079575,0.00191094743796522,0.00662185258908696,0.000188266912551749;...
    0.252458254214951,0.0125091293591534,-0.327121080864820,0.226053196913636,-0.261297620204730,-0.216236975254179,0.344080559560279,-0.121230620609225,-0.0602714322805403,-0.219189670580381,0.211470671417798,-0.128634841134357,0.576234521196323,-0.0550760430415841,0.197333043823494,0.205766455381630,0.0236621449982310,0.00370370109158968,0.000185197342559174;...
    0.245944240653730,0.0799604139803053,-0.358449872812475,0.0640998104960830,-0.341792253779399,0.0224967545657750,0.388717982330275,0.177122203733985,-0.255990758059885,-0.00644303562267801,-0.375390122964730,-0.0762061466572891,-0.502002035655987,0.0183662355165662,-0.176351451500561,-0.0686345486566678,0.0154233464320777,0.00541822807684994,0.00127868850390311;...
    0.225758567264058,0.191381035473567,-0.335176303224685,-0.216152633771253,-0.165178954634192,0.423011571619087,-0.144255461671014,0.187567292296388,0.149360081833295,0.530105300516035,0.0417962320677310,0.326784764375550,0.274609836157270,0.0558724755215897,0.00442737635788848,0.0111352503084675,0.0244045338837114,0.000682889258092345,0.000198665379271490;...
    0.211097168683674,0.259405649803992,-0.243643584807687,-0.325745718814885,0.0763484285808871,0.330279571319951,-0.220001696329049,-0.117395307490011,0.271296590718256,-0.438774328341411,0.148165305851337,-0.484545679436120,-0.143691433970383,-0.0389147466608019,0.00893381613114769,-0.0205549231060488,-0.00611584426449596,0.00380246484862132,-0.000389229000503269;...
    0.188387436860863,0.329799740851314,-0.0946692611641819,-0.273646378894757,0.356651426359913,-0.153161129681719,-0.000682050969765767,-0.329897663448327,-0.267361749795915,-0.279382155927600,-0.263740500608820,0.528987613257166,0.0703281456138742,-0.0835258438506085,-0.0254953444060441,0.0271139318663258,0.0126194670232453,0.00385357399236438,-0.000809502038431328;...
    0.176395533106338,0.357332294420881,0.0554387508851632,-0.155161957798969,0.354513035089274,-0.343041979311133,0.161895880026785,-0.0275355960813715,-0.207859409120980,0.506833313170992,0.205450556183935,-0.413916086467618,-0.0407497250821277,0.168472840930614,-0.00235465742383666,-0.00588317616785561,-0.00334222371106789,0.00197868686362638,0.00170392027445259;...
    0.165469018345669,0.360040619637271,0.260392170105291,0.0669803672081155,0.0572701975680616,-0.220913199055638,0.181062550106522,0.519913777311151,0.462086625105155,-0.104489884552583,-0.0219632037127792,0.119011798891668,-0.00429988165275597,-0.417545575417374,-0.0401081045672218,0.0200607203872710,-0.00526025219733044,-0.00580167364426092,0.000458230034086340;...
    0.159580892256856,0.347927159738324,0.348346295207843,0.239394140984691,-0.157211928082521,0.0928779108728929,-0.00501602103568442,0.0169759687957175,0.109978222336614,-0.182603153808833,-0.121233489669211,0.0711306930722663,0.0620582733857731,0.750450107378881,0.0785420660886401,-0.0521102089376762,0.00751936646888490,-0.00188268149870413,-0.00185190445021904;...
    0.148832920730500,0.332847695392180,0.365098051686060,0.331227694894700,-0.281334582456685,0.283334016381439,-0.182848344502475,-0.325817997238157,-0.310946153889112,0.128556113560638,0.0837173663270596,-0.0703828760736262,-0.0461924125750960,-0.441499963638099,-0.0405906261145499,0.0337161618200980,0.00116291918248088,0.00401605033019270,0.000505570344725481];

end

% Additional function to extract semi-variogram model coefficients
function modelVario = get_model_varios()
Cn = [2.5, 0.5, 0.15, 0.15, 0.314321867136085,...
    0.190749535511365, 0.137846758971697, 0.111283843493347, 0.096499281204429,...
    0.071736796680044, 0.064816215163269, 0.054076635653157, 0.051188751201166,...
    0.043316419835114, 0.041398046055658, 0.034663671578289, 0.018796994070527,...
    0.002856941187582, 3.606453961200486e-04];
C1 = [4.52, 1.4, 0.42, 0.225];
C2 = [6.78, 2.6, 0.63, 0.225];
a1 = [15, 10, 15, 10];
a2 = [250, 160, 160, 120];
Type = {'iso nest','iso nest', 'iso nest', 'iso nest','nug','nug','nug','nug',...
    'nug','nug','nug','nug','nug','nug','nug','nug','nug','nug','nug'};
modelVario = {};
for i = 1:19
    modelVario{i,1}.Cn = Cn(i);
    modelVario{i,1}.Type = Type{i};
    if  strcmp(Type{i},'iso nest')
        modelVario{i,1}.C1 = C1(i);
        modelVario{i,1}.C2 = C2(i);
        modelVario{i,1}.a1 = a1(i);
        modelVario{i,1}.a2 = a2(i);
    end
end
end