% Last edited by L. Shen @ Feb 27, 2019. Contact luyi@ualberta.ca for bugs 
% and glitches.
% This is part of the supplementary materisl for Shen et al., 2019, 
% Quantitative constraints to the complete state of stress from combined 
% borehole and focal mechanism inversions: Fox Creek, Alberta 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% before running please make sure the following files are in the same
% folder: read_in_AZ.m read_in_PP.m read_in_Sh.m read_in_SV.m
% FormationPP.txt MinimumSh.txt SH_azimuth.txt shape_ratio_statistics.mat
% specify input coordinates
% here is an example input for 11 earthquake's locations 
Lat = [
    54.3845
54.4618
54.4972
54.5076
54.414057
54.438115
54.3536
54.3543
54.3261
54.1572
54.3096
];
Lon = [-117.4
-117.266
-117.382
-117.318
-117.309173
-117.340393
-117.32
-117.222
-116.706
-116.886
-117.62
];
Depth = [ 3.5 3.3 3.4 3.3 3.4 3.4 3.5 3.4 3.1 3.3 3.9]; % depth in km
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Depth =  Depth(:);
Lat = Lat(:);
Lon = Lon(:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s_ratio = 1.72; 
% 1.72 is the rough ratio between great circle distances of 1 degree 
% latitude and 1 degree longitude at the latitude of 52 degree. This ratio is
% used for 2D interpolation
%
n = 10000; %% number of random simulations to for to invert for SH
%% input gradient used to extrapolate stress values at formations near the Duvernay
%% see data and methods section in the main text for details
SV_gradient = 25.5; % KPa/m
SV_gradient_uncertainty = 2; %KPa/m
Sh_gradient = 32.1;%KPa/m
Sh_gradient_uncertainty = 3.1;%KPa/m
PP_gradient = 29.1;%KPa/m
PP_gradient_uncertainty = 7.2;%KPa/m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output stress tensor T
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Do not modify below this line
%%
T_SH = [];
T_SH_uncertainty = [];
T_SV = [];
T_SV_uncertainty = [];
T_Sh = [];
T_Sh_uncertainty = [];
T_AZ = [];
T_Az_confidence = [];
T_PP = [];
T_PP_uncertainty = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% load shape ratio statistics from stress inversion
load shape_ratio_statistics
% read in SV
read_in_SV;
Lat_SV = VerticalSV.Lat;
Lon_SV = VerticalSV.Lon;
Dep_SV = VerticalSV.Dep;
SV = VerticalSV.SV;
SV_uncertainty = VerticalSV.Uncertainty;
Dep_SV(isnan(SV))=[];
Lat_SV(isnan(SV))=[];
Lon_SV(isnan(SV))=[];
SV_uncertainty(isnan(SV))=[];
SV(isnan(SV))=[];
%
Lat_SV_scale = Lat_SV*s_ratio; 
%
%% creat scatter interpolant
%
F_SV = scatteredInterpolant(Lat_SV_scale,Lon_SV,SV);
F_SV_uncertainty = scatteredInterpolant(Lat_SV_scale,Lon_SV,SV_uncertainty);
F_dep = scatteredInterpolant(Lat_SV_scale,Lon_SV,Dep_SV);
%
SV_eq = F_SV(Lat*s_ratio, Lon);
SV_eq_uncertainty = F_SV_uncertainty(Lat*s_ratio, Lon);
SV_eq_dep = F_dep(Lat*s_ratio, Lon);
Dep_diff = SV_eq_dep(:) - Depth(:)*1000;
SV_dep_correction = Dep_diff*SV_gradient/1000;
SV_dep_correction_uncertainty = Dep_diff*SV_gradient_uncertainty/1e6;
T_SV = SV_eq - SV_dep_correction;
T_SV_uncertainty = SV_eq_uncertainty + SV_dep_correction_uncertainty;
% read in Sh map
read_in_Sh;
%%%%%%%%%%%%%%%%%%%
Lat_Sh = MinimumSh.Lat;
Lon_Sh = MinimumSh.Lon;
Dep_Sh = MinimumSh.Dep;
Sh = MinimumSh.Sh;
Sh_uncertainty = MinimumSh.Sh_uncertainty;

Dep_Sh(isnan(Sh))=[];
Lat_Sh(isnan(Sh))=[];
Lon_Sh(isnan(Sh))=[];
Sh_uncertainty(isnan(Sh))=[];
Sh(isnan(Sh))=[];
%
Lat_Sh_scale = Lat_Sh*s_ratio; 
%
F_Sh = scatteredInterpolant(Lat_Sh_scale,Lon_Sh,Sh);
F_Sh_uncertainty = scatteredInterpolant(Lat_Sh_scale,Lon_Sh,Sh_uncertainty);
F_Sh_dep = scatteredInterpolant(Lat_Sh_scale,Lon_Sh,Dep_Sh);
%
Sh_eq = F_Sh(Lat*s_ratio, Lon);
Sh_eq_uncertainty = F_Sh_uncertainty(Lat*s_ratio, Lon);
Sh_eq_dep = F_dep(Lat*s_ratio, Lon);
Dep_diff = Sh_eq_dep(:) - Depth(:)*1000;
Sh_dep_correction = Dep_diff*Sh_gradient/1e3;
Sh_dep_correction_uncertainty = Dep_diff*Sh_gradient_uncertainty/1e6;
T_Sh = Sh_eq - Sh_dep_correction;
T_Sh_uncertainty = Sh_eq_uncertainty + Sh_dep_correction_uncertainty;
%
%
% compute SH
for i = 1:length(T_SV)
    for j = 1:n
        SV_input = (rand - 0.5)*2 * T_SV_uncertainty(i) + T_SV(i);
        Sh_input = (rand - 0.5)*2 * T_Sh_uncertainty(i) + T_Sh(i);
        SR_input = shape_ratio_statistics(round(rand*(length(shape_ratio_statistics)-1)+1));
        T_SH_stat(i,j) = (SV_input - SR_input*Sh_input)/(1-SR_input);
    end
    T_SH(i) = median(T_SH_stat(i,:));
    pd = fitdist(T_SH_stat(i,:)','kernel');
    Beta_fit = icdf(pd,[0.05 0.95]);
    T_SHmax(i) = Beta_fit(2);
    T_SHmin(i) = Beta_fit(1);
end
% read in Az map
read_in_AZ;
Lat_AZ = SHAzimuth.Lat;
Lon_AZ = SHAzimuth.Lon;
AZ = SHAzimuth.Az;
AZ_Confidence = SHAzimuth.Confidence;
Lat_AZ(isnan(AZ))=[];
Lon_AZ(isnan(AZ))=[];
AZ_Confidence(isnan(AZ))=[];
AZ(isnan(AZ))=[];
Lat_AZ_scale = Lat_AZ*s_ratio; 
F_AZ = scatteredInterpolant(Lat_AZ_scale,Lon_AZ,AZ);
F_AZ_Confidence = scatteredInterpolant(Lat_AZ_scale,Lon_AZ,AZ_Confidence);
T_AZ = F_AZ(Lat*s_ratio, Lon);
T_AZ_Confidence = F_AZ_Confidence(Lat*s_ratio, Lon);
% read Pp
read_in_PP;
%%%%%%%%%%%%
Lat_PP = FormationPP.Lat;
Lon_PP = FormationPP.Lon;
Dep_PP = FormationPP.Depth;
PP = FormationPP.Pp;
PP_uncertainty = FormationPP.Uncertainty;
Dep_PP(isnan(PP))=[];
Lat_PP(isnan(PP))=[];
Lon_PP(isnan(PP))=[];
PP_uncertainty(isnan(PP))=[];
PP(isnan(PP))=[];
%
Lat_PP_scale = Lat_PP*s_ratio; 
F_PP = scatteredInterpolant(Lat_PP_scale,Lon_PP,PP);
F_PP_uncertainty = scatteredInterpolant(Lat_PP_scale,Lon_PP,PP_uncertainty);
F_PP_dep = scatteredInterpolant(Lat_PP_scale,Lon_PP,Dep_PP);
PP_eq = F_PP(Lat*s_ratio, Lon);
PP_eq_uncertainty = F_Sh_uncertainty(Lat*s_ratio, Lon);
PP_eq_dep = F_dep(Lat*s_ratio, Lon);
Dep_diff = PP_eq_dep(:) - Depth(:)*1000;
PP_dep_correction = Dep_diff*PP_gradient/1e3;
PP_dep_correction_uncertainty = Dep_diff*PP_gradient_uncertainty/1e6;
T_PP = PP_eq - PP_dep_correction;
T_PP_uncertainty = PP_eq_uncertainty + PP_dep_correction_uncertainty;
%%
%%%%%%%%
%%%%%%%%
%%
%% start printout
%%
fprintf(' Lat     Lon        Depth  SH     SHmax    SHmin  Sh     Uncertainty     SV   Uncertainty  PP  Uncetainty   SH_Azimuth     Azimuh_confidence \n')
format short
datasave=[];

for i = 1:length(Depth)
    fprintf('% 2.4f % 3.4f % 5.3f % 5.2f % 5.2f % 5.2f % 5.2f %5.2f           %5.2f  %5.2f       %5.2f %5.2f         %5.0f        %5.1f\n', Lat(i), Lon(i), Depth(i), T_SH(i), T_SHmax(i), T_SHmin(i), T_Sh(i), T_Sh_uncertainty(i), T_SV(i), T_SV_uncertainty(i), T_PP(i), T_PP_uncertainty(i), T_AZ(i), T_AZ_Confidence(i));
end
% figure
% 
% boxplot(T_SH_stat')
% title('Stress Values')
% ylabel('Stress (MPa)')
% xlabel('Enquiry point')
% 
% hold on
% errorbar(1:length(T_SV),T_SV,T_SV_uncertainty,'O')
% errorbar(1:length(T_Sh),T_Sh,T_Sh_uncertainty,'O')
% errorbar(1:length(T_PP),T_PP,T_PP_uncertainty,'O')

