% This script calculates the stress concentration around the borehole wall in anisotropic or isotropic formations.
clf;
clear all;
close all;
tic
method= 2;% =1 Kirsch's solution (2D); =2 Hiramatsu's solution (3D)
%% i) Radius and trajectory of the borehole defined referred to the NED (North(x)-East(y)-Downward(z)) coordinate system
rbh=0.1;% radius of borehole (unit, m)
traj=[90,0];%traj = [boreholeI boreholeA]:  A 1 X 2 vector describing inclination and azimuth of the borehole in degrees
% boreholeI -inclination of the borehole; boreholeA -azimuth of the borehole
%% ii) The far-field stress state applied to the borehole
% Magitude of the components
sigmaV=30.0; % magnitude of vertical stress, MPa
sigmaH=20.0; % magnitude of maximum horizontal stress, MPa
sigmah=10.0; % magnitude of minimum horizontal stress, MPa
% The direction of stress coordinates. We followed the original paper to
% use the rotation angles instead of the geology angles.
principle1=0;% azimuth angle of principle stress tensor coordinates
principle2=0;% tile angle of principle stress tensor coordinates
%% iii) The media elastic property
Poisson=0.079;% Poisson's ratio
%% Convert the input angles from degrees into radians
borehole1=traj(2);% azimuth angle of borehole coordinates
borehole2=traj(1);% tilt angle of borehole coordinates
principle1=principle1/180*pi;% azimuth angle of principle stress tensor coordinates
principle2=principle2/180*pi;% tile angle of principle stress tensor coordinates
borehole1=borehole1/180*pi;% azimuth angle of borehole coordinates
borehole2=borehole2/180*pi;% tilt angle of borehole coordinates
%% Main function
n1=101; % the mesh of radius
rmedia=4*rbh;% radius of the media (unit, m)
n2=361; % the mesh of azimuth angle
radius=0:rmedia/(n1-1):rmedia;
sita=-pi:2*pi/(n2-1):pi;
dphi=zeros(3,1);
sigmax=zeros(n1,n2);
sigmay=zeros(n1,n2);
sigmaz=zeros(n1,n2);
taoyz=zeros(n1,n2);
taoxz=zeros(n1,n2);
taoxy=zeros(n1,n2);
sigmarr=zeros(n1,n2);
sigmass=zeros(n1,n2);
sigmazz=zeros(n1,n2);
sigmars=zeros(n1,n2);
sigmasz=zeros(n1,n2);
sigmarz=zeros(n1,n2);
dcs=StrRMatrixY(principle2)*StrRMatrixZ(principle1); % direction cosines
for ii=1:n1
    for jj=1:n2
        if (radius(ii)<rbh) % let all stress components inside the borehole to be NaN.
            sigmax(ii,jj)=NaN;
            sigmay(ii,jj)=NaN;
            sigmaz(ii,jj)=NaN;
            taoyz(ii,jj)=NaN;
            taoxz(ii,jj)=NaN;
            taoxy(ii,jj)=NaN;
            sigmarr(ii,jj)=NaN;
            sigmass(ii,jj)=NaN;
            sigmasz(ii,jj)=NaN;
        else
        xx=cos(sita(jj))*radius(ii);
        yy=sin(sita(jj))*radius(ii);
        if (method == 1)
            sigmass(ii,jj)= sigmaH*(1+rbh^2/radius(ii)^2)/2- sigmaH*(1+3*rbh^4/radius(ii)^4)*cos(2*sita(jj))/2;
            sigmarr(ii,jj)= sigmaH*(1-rbh^2/radius(ii)^2)/2+ sigmaH*(1+3*rbh^4/radius(ii)^4-4*rbh^2/radius(ii)^2)*cos(2*sita(jj))/2;
            sigmars(ii,jj)= -sigmaH*(1-3*rbh^4/radius(ii)^4+2*rbh^2/radius(ii)^2)*sin(2*sita(jj))/2;
        end
        if (method == 2)
            coalpha1=0.5*( (dcs(1,1)^2*sin(borehole2)^2+ dcs(1,2)^2+ dcs(1,3)^2*cos(borehole2)^2- 2*dcs(1,3)*dcs(1,1)*sin(borehole2)*cos(borehole2))*sigmaH+...
                           (dcs(2,1)^2*sin(borehole2)^2+ dcs(2,2)^2+ dcs(2,3)^2*cos(borehole2)^2- 2*dcs(2,3)*dcs(2,1)*sin(borehole2)*cos(borehole2))*sigmah+...
                           (dcs(3,1)^2*sin(borehole2)^2+ dcs(3,2)^2+ dcs(3,3)^2*cos(borehole2)^2- 2*dcs(3,3)*dcs(3,1)*sin(borehole2)*cos(borehole2))*sigmaV );
            coalpha2=0.5*( (-dcs(1,1)^2*sin(borehole2)^2+ dcs(1,2)^2- dcs(1,3)^2*cos(borehole2)+ 2*dcs(1,3)*dcs(1,1)*sin(borehole2)*cos(borehole2))*sigmaH+...
                           (-dcs(2,1)^2*sin(borehole2)^2+ dcs(2,2)^2- dcs(2,3)^2*cos(borehole2)^2+ 2*dcs(2,3)*dcs(2,1)*sin(borehole2)*cos(borehole2))*sigmah+...
                           (-dcs(3,1)^2*sin(borehole2)^2+ dcs(3,2)^2- dcs(3,3)^2*cos(borehole2)^2+ 2*dcs(3,3)*dcs(3,1)*sin(borehole2)*cos(borehole2))*sigmaV );
            coalpha3=( dcs(1,2)*dcs(1,3)*cos(borehole2)-dcs(1,1)*dcs(1,2)*sin(borehole2) )*sigmaH+...
                     ( dcs(2,2)*dcs(2,3)*cos(borehole2)-dcs(2,1)*dcs(2,2)*sin(borehole2) )*sigmah+...
                     ( dcs(3,2)*dcs(3,3)*cos(borehole2)-dcs(3,1)*dcs(3,2)*sin(borehole2) )*sigmaV;
            cobeta1=(dcs(1,1)^2*cos(borehole2)^2+ dcs(1,3)^2*sin(borehole2)^2+ 2*dcs(1,3)*dcs(1,1)*sin(borehole2)*cos(borehole2))*sigmaH+...
                    (dcs(2,1)^2*cos(borehole2)^2+ dcs(2,3)^2*sin(borehole2)^2+ 2*dcs(2,3)*dcs(2,1)*sin(borehole2)*cos(borehole2))*sigmah+...
                    (dcs(3,1)^2*cos(borehole2)^2+ dcs(3,3)^2*sin(borehole2)^2+ 2*dcs(3,3)*dcs(3,1)*sin(borehole2)*cos(borehole2))*sigmaV;
            cogamma1=(-dcs(1,1)^2*sin(borehole2)*cos(borehole2)+ dcs(1,3)^2*cos(borehole2)*sin(borehole2)+ dcs(1,3)*dcs(1,1)*(cos(borehole2)^2-sin(borehole2)^2))*sigmaH+...
                     (-dcs(2,1)^2*sin(borehole2)*cos(borehole2)+ dcs(2,3)^2*cos(borehole2)*sin(borehole2)+ dcs(2,3)*dcs(2,1)*(cos(borehole2)^2-sin(borehole2)^2))*sigmah+...
                     (-dcs(3,1)^2*sin(borehole2)*cos(borehole2)+ dcs(3,3)^2*cos(borehole2)*sin(borehole2)+ dcs(3,3)*dcs(3,1)*(cos(borehole2)^2-sin(borehole2)^2))*sigmaV;
            cogamma2=(-dcs(1,2)*dcs(1,3)*sin(borehole2)-dcs(1,1)*dcs(1,2)*cos(borehole2))*sigmaH+...
                     (-dcs(2,2)*dcs(2,3)*sin(borehole2)-dcs(2,1)*dcs(2,2)*cos(borehole2))*sigmah+...
                     (-dcs(3,2)*dcs(3,3)*sin(borehole2)-dcs(3,1)*dcs(3,2)*cos(borehole2))*sigmaV;
            sigmarr(ii,jj)=coalpha1*(1-rbh^2/radius(ii)^2)+...
                           coalpha2*(1-4*rbh^2/radius(ii)^2+ 3*rbh^4/radius(ii)^4)*cos(2*sita(jj))+...
                           coalpha3*(1-4*rbh^2/radius(ii)^2+ 3*rbh^4/radius(ii)^4)*sin(2*sita(jj));
            sigmass(ii,jj)=coalpha1*(1+rbh^2/radius(ii)^2)+...
                           coalpha2*(-1- 3*rbh^4/radius(ii)^4)*cos(2*sita(jj))+...
                           coalpha3*(-1- 3*rbh^4/radius(ii)^4)*sin(2*sita(jj));
            sigmazz(ii,jj)=cobeta1- 4*Poisson*(coalpha2*rbh^2/radius(ii)^2*cos(2*sita(jj))+ coalpha3*rbh^2/radius(ii)^2*sin(2*sita(jj)));
            sigmasz(ii,jj)=cogamma1*(1+rbh^2/radius(ii)^2)*cos(sita(jj))+ cogamma2*(1+rbh^2/radius(ii)^2)*sin(sita(jj));
            sigmarz(ii,jj)=cogamma1*(1-rbh^2/radius(ii)^2)*sin(sita(jj))- cogamma2*(1-rbh^2/radius(ii)^2)*cos(sita(jj));
            sigmars(ii,jj)=coalpha2*(-1-2*rbh^2/radius(ii)^2+3*rbh^4/radius(ii)^4)*sin(2*sita(jj))+...
                           coalpha3*(1+2*rbh^2/radius(ii)^2-3*rbh^4/radius(ii)^4)*cos(2*sita(jj));
        end    
        end
    end
end
%% Plots
ppx=radius.'*cos(sita);
ppy=radius.'*sin(sita);
figure(1)
hold on
nnn=int16(rbh/rmedia*(n1-1)+1);
plot(sita/pi*180, sigmarr(nnn,:), 'black');
plot(sita/pi*180, sigmass(nnn,:), 'blue');
plot(sita/pi*180, sigmazz(nnn,:), 'r');
plot(sita/pi*180, sigmasz(nnn,:), 'green');
plot(sita/pi*180, sigmarz(nnn,:), '--green');
plot(sita/pi*180, sigmars(nnn,:), '--black');
legend({'sigmarr', 'sigmass', 'sigmazz', 'sigmasz', 'sigmarz', 'sigmars'},'Location','EastOutside'); 
xlabel('Azimuth angle','FontSize',12,'FontName','Times new roman');
ylabel('Stress, MPa','FontSize',12,'FontName','Times new roman');
set(gca,'xlim',[-180 180]);
set(gca,'Xtick',-180:90:180);
hold off