% This function calculates the stress distributed around the borehole in anisotropic or isotropic formations.
% The inputs includes:
%     i) Radius rbh and orientation of the borehole traj;
%     ii) Magitude and orientation of the far-field stress: sigmas and orientS;
%     iii) The borehole fluid pressure pw;
%     iv) Magnitude and orientation of the elastic stiffness matrix: Cm and orientM;
%     v) The size and mesh of the medium geoM
% The outputs are the stress in cylindrical co-ordinates including:
%     sigmarr- ¦Òrr;
%     sigmass- ¦Òss;
%     sigmazz- ¦Òzz;
%     sigmasz- ¦Òsz;
%     sigmarz- ¦Òrz;
%     sigmars- ¦Òrs.
function [sigmarr, sigmass, sigmazz, sigmasz, sigmarz, sigmars] = StressDistribution(rbh, traj, sigmas, orientS, pw, Cm, orientM, geom)
%% Convert the input angles from degrees into radians
traj=traj/180*pi;
orientM=orientM/180*pi;
orientS=orientS/180*pi;
%% Convert the geology angles into Euler rotation angles
% Please look at Fig. 3 in sectioin 2.1. See also equations (4) and (8).
alphazm1= orientM(1) - pi/2;% The rotation angle of the medium symmetry frame around D-axis
alphaym= orientM(2);% The rotation angle of the medium symmetry frame around y'-axis
alphazm2= pi/2 - orientM(3);% The rotation angle of the medium symmetry frame around zm-axis
alphazs1= orientS(1) - pi/2;% The rotation angle of the far-field stress frame around D-axis
alphays= orientS(2);% The rotation angle of the far-field stress frame around y''-axis
alphazs2= pi/2 - orientS(3);% The rotation angle of the far-field stress frame around zs-axis
%% Rotate the elastic stiffnesses from the medium symmetric coordinate frame to borehole coordinate frame
C=StiffnessRotation3D(Cm, alphazm1, alphaym, alphazm2, traj);
%% Rotate the stress tensor from the principle coordinate frame to borehole coordinate frame
is=StressRotation3D(sigmas, alphazs1, alphays, alphazs2, traj);
istr=[is(1,1) is(2,2) is(3,3) is(2,3) is(1,3) is(1,2)];
%% General solution
rmedia=geom(1);% radius of the media (unit, m)
n1=geom(2); % the mesh of radius
n2=geom(3); % the mesh of azimuth angle as defined in Fig.2(b).
A=C^(-1); % compliance matrix, Eqn.(15)
B=A-A(:,3)*A(:,3)'/A(3,3); % the matrix of reduced strain coefficients, Eqn.(14).
weili4=[B(1,1), -2*B(1,6), 2*B(1,2)+B(6,6), -2*B(2,6), B(2,2)]; % Eqn.(C.2)
weili3=[B(1,5), -(B(1,4)+B(5,6)), B(2,5)+B(4,6), -B(2,4)]; % Eqn.(C.2)
weili2=[B(5,5), -2*B(4,5), B(4,4)]; % Eqn.(C.2)
% for general cases
sixorder=conv(weili4,weili2)-conv(weili3, weili3); % Eqn.(C.1)
rootsmiu=roots(sixorder);
mius=rootsmiu(imag(rootsmiu)>=0); % solutions of the sixorder equation
l4=polyval(weili4,mius);% Eqn.(C.2)
l3=polyval(weili3,mius);% Eqn.(C.2)
l2=polyval(weili2,mius);% Eqn.(C.2)
if ( all(abs(l3)<10^(-6)) ) % for special cases (see section 2.5)
    rootsmiu1=roots(weili4); % Eqn.(C.2)
    roots1=rootsmiu1(imag(rootsmiu1)>=0);
    rootsmiu2=roots(weili2); % Eqn.(C.2)
    roots2=rootsmiu2(imag(rootsmiu2)>=0);
    mius=zeros(3,1);
    mius(1)=roots1(1);
    mius(2)=roots1(2);
    mius(3)=roots2(1);
end
lamda=[-l3(1)/l2(1); -l3(2)/l2(2); -l3(3)/l4(3)];% lamdas, see Eqn.(18)
laplace=mius(2)-mius(1)+lamda(3)*lamda(2)*(mius(1)-mius(3))+lamda(1)*lamda(3)*(mius(3)-mius(2)); % Eqn.(22)
theta=-pi:2*pi/(n2-1):pi;
radius=0:rmedia/(n1-1):rmedia;
dphi=zeros(3,1);
sigmazz=zeros(n1,n2);
sigmarr=zeros(n1,n2);
sigmass=zeros(n1,n2);
sigmasz=zeros(n1,n2);
sigmarz=zeros(n1,n2);
sigmars=zeros(n1,n2);
for ii=1:n1
    for jj=1:n2
        if (radius(ii)<rbh) % let all stress components inside the borehole being NaN.
            sigmazz(ii,jj)=NaN;
            sigmarr(ii,jj)=NaN;
            sigmass(ii,jj)=NaN;
            sigmasz(ii,jj)=NaN;
            sigmarz(ii,jj)=NaN;
            sigmars(ii,jj)=NaN;
        else
            xx=cos(theta(jj))*radius(ii);
            yy=sin(theta(jj))*radius(ii);
            zz=xx+mius*yy; % Eqn.(17)
            % The relationship between eta and zk in the theory is not
            % unique. There are two possible sign (positive or negative).
            % We choose the sign according to which one offers a larger magnitude of eta.
            % This is referred to the discussion section.
            eta01=((zz/rbh).^2-1-mius.*mius).^0.5; % the part in the square root of Eqn.(D.6)
            eta02=-((zz/rbh).^2-1-mius.*mius).^0.5;
            eta1=( zz/rbh+ eta01 )./(1-mius*1i); % Eqn.(D.6)
            eta2=( zz/rbh+ eta02 )./(1-mius*1i);
            X1=2*laplace*eta1.*eta01; % Eqn.(21) which is the dominator of Eqn.(20)
            X2=2*laplace*eta2.*eta02;
            eta0=eta01;
            eta=eta1;
            X=X1;
            for kkk= 1:3
                if ( abs(eta1(kkk))<abs(eta2(kkk)) )
                    eta0(kkk)=eta02(kkk);
                    eta(kkk)=eta2(kkk);
                    X(kkk)=X2(kkk);
                end
            end 
            phi1=1i*istr(6)-istr(2)+pw;
            phi2=istr(6)-1i*istr(1)+1i*pw;
            phi3=istr(4)-1i*istr(5);
            % Eqn.(20)
            dphi(1)=-( phi1*(mius(2)-lamda(2)*lamda(3)*mius(3))+...
                phi2*(lamda(2)*lamda(3)-1)+...
                phi3*lamda(3)*(mius(3)-mius(2)) )/...
                X(1);
            dphi(2)=-( -phi1*(mius(1)-lamda(1)*lamda(3)*mius(3))-...
                phi2*(lamda(1)*lamda(3)-1)+...
                phi3*lamda(3)*(mius(1)-mius(3)) )/...
                X(2);
            dphi(3)=-( phi1*(mius(1)*lamda(2)-mius(2)*lamda(1))+...
                phi2*(lamda(1)-lamda(2))+...
                phi3*(mius(2)-mius(1)) )/...
                X(3);
            % Eqn.(23)
            istrbore=[ 2*real( mius(1)^2*dphi(1)+ mius(2)^2*dphi(2)+ lamda(3)*mius(3)^2*dphi(3) ) ...
                2*real( dphi(1)+ dphi(2)+ lamda(3)*dphi(3) ) ...
                0 ...
                -2*real( lamda(1)*dphi(1)+ lamda(2)*dphi(2)+ dphi(3) ) ...
                2*real( lamda(1)*mius(1)*dphi(1)+ lamda(2)*mius(2)*dphi(2)+ mius(3)*dphi(3) ) ...
                -2*real( mius(1)*dphi(1)+ mius(2)*dphi(2)+ lamda(3)*mius(3)*dphi(3) )];
            tt=-A(3,:)*istrbore.'/A(3,3); % Eqn.(B6)
            istrbore(3)=tt;
            % Eqn.(23)
            Carstr=istr+istrbore;
            % Eqn.(24) used to convert the stress into cylindrical co-ordinate
            Cylinstr=Cartocylin(Carstr,theta(jj));
            Cylinstr=roundn(Cylinstr,-3);% round the results
            sigmarr(ii,jj)=Cylinstr(1);
            sigmass(ii,jj)=Cylinstr(2);
            sigmazz(ii,jj)=Cylinstr(3);
            sigmasz(ii,jj)=Cylinstr(4);
            sigmarz(ii,jj)=Cylinstr(5);
            sigmars(ii,jj)=Cylinstr(6);
        end
    end
end
%% Plots
figure(1)
hold on
nnn=int16(rbh/rmedia*(n1-1)+1);
plot(theta/pi*180, sigmarr(nnn,:), 'black');
plot(theta/pi*180, sigmass(nnn,:), 'blue');
plot(theta/pi*180, sigmazz(nnn,:), 'r');
plot(theta/pi*180, sigmasz(nnn,:), 'green');
plot(theta/pi*180, sigmarz(nnn,:), '--green');
plot(theta/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
% figure(2)
nthe=4;
nrr=8;
PolarPcolorWeili(theta,radius,sigmarr,rmedia,n2,nthe,nrr);
title('sigma rr','position',[-rmedia,rmedia],'FontSize',12);
% caxis([5,30]); % These colorbar ranges are used to plot Fig.7. We leave them here for the convinience in plotting.
% figure(3)
PolarPcolorWeili(theta,radius,sigmazz,rmedia,n2,nthe,nrr);
title('sigma zz','position',[-rmedia,rmedia],'FontSize',12);
% caxis([15,30]);
% figure(4)
PolarPcolorWeili(theta,radius,sigmass,rmedia,n2,nthe,nrr);
title('sigma ss','position',[-rmedia,rmedia],'FontSize',12);
% caxis([-15,75]);
% figure(5)
PolarPcolorWeili(theta,radius,sigmasz,rmedia,n2,nthe,nrr);
title('sigma sz','position',[-rmedia,rmedia],'FontSize',12);
% caxis([-10,10]);
% figure(6)
PolarPcolorWeili(theta,radius,sigmarz,rmedia,n2,nthe,nrr);
title('sigma rz','position',[-rmedia,rmedia],'FontSize',12);
caxis([-5,5]);
% figure(7)
PolarPcolorWeili(theta,radius,sigmars,rmedia,n2,nthe,nrr);
title('sigma rs','position',[-rmedia,rmedia],'FontSize',12);
% caxis([-15,15]);
toc