% Author: Wei Li
% Please report any bug to weili.rockphysics@gmail.com or liwei2014@cugb.edu.cn.
% Referred manuscript title: "Stress state around an inclined borehole in the anisotropic formation".
% The ASCIB3D algorithm is short for the arbitrary stress concentrated around inclined borehole in 3-dementions.
% The algorithm calculates the stress distributed around the borehole in anisotropic or isotropic formations.
% The outputs are the stress in cylindrical co-ordinates including:
%     sigmarr- ¦Òrr;
%     sigmass- ¦Òss;
%     sigmazz- ¦Òzz;
%     sigmasz- ¦Òsz;
%     sigmarz- ¦Òrz;
%     sigmars- ¦Òrs.
clf;
clear all;
close all;
tic
%% 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=[45,45];%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
sigmas=[20 0 0;% The applied far-field stress tensor in units of MPa (Eqn. 8)
        0 10 0;
        0 0 30];
% Orientation of far-field stress frame defined referred to the NED (North(x)-East(y)-Downward(z)) coordinate system
orientS=[90,0,90];%orientS=[strikeS,dipS,rakeS]=(¦×s,¦Äs,¦Ãs):a 1 X 3 vector containing the strike, dip, and rake of the stress tensor in degrees
% strikeS- strike of ¦Òxx-¦Òyy plane; dipS- dip angle of ¦Òxx-¦Òyy plane;
% rakeS- the angle of line rake indicating the direction of ¦Òxx within ¦Òxx-¦Òyy plane
%% iii) The borehole fluid pressure
pw=5.0; % magnitude of wellbore pressure in units of MPa
%% iv) The media elastic property
% Magitude of the components
% example 0: Isotropic case referred to Fig. 5(a), this matrix is calculated with Young's modulus and Poisson's ratio equaling 20 and 0.079, respectively.
% Cm=[20.275,1.739,1.739,0,0,0;%Cm is the elastic constant matrix in the material coordinate frame in GPa
%     1.739,20.275,1.739,0,0,0;
%     1.739,1.739,20.275,0,0,0;
%     0,0,0,9.268,0,0;
%     0,0,0,0,9.268,0;
%     0,0,0,0,0,9.268];
% example 1: Transversely isotropic sample(Gaede et al., 2012) referred to Fig. 5(b) and (c)
Cm=[45.2,16.4,19.67,0,0,0;%Cm is the elastic constant matrix in the material coordinate frame in GPa
    16.4,45.2,19.67,0,0,0;
    19.67,19.67,28,0,0,0;
    0,0,0,7.05,0,0;
    0,0,0,0,7.05,0;
    0,0,0,0,0,14.4];
% example 2: Transversely isotropic mylonite(Li et al., 2019) referred to the Fig. 6
% Cm=[62.885,9.672,-0.053,0,0,0;%Cm is the elastic constant matrix in the material coordinate frame in GPa
%     9.672,62.885,-0.053,0,0,0;
%     -0.053,-0.053,33.427,0,0,0;
%     0,0,0,19.239,0,0;
%     0,0,0,0,19.239,0;
%     0,0,0,0,0,26.607];
% example 3: orthorhombic sample referred to Fig. 5(d) and Fig. 7
% Cm=[35.6,8.2,7.9,0,0,0;%Cm is the elastic constant matrix in the material coordinate frame in GPa
%     8.2,22.15,3.85,0,0,0;
%     7.9,3.85,16.97,0,0,0;
%     0,0,0,6,0,0;
%     0,0,0,0,7,0;
%     0,0,0,0,0,8];
% Orientation of medium symmetry frame defined referred to the NED (North(x)-East(y)-Downward(z)) coordinate system
orientM=[120,30,90];%orientM=[strikeM,dipM,rakeM]=(¦×m,¦Äm,¦Ãm): a 1 X 3 vector containing the strike, dip, and rake of the stiffness matrix in degrees. 
%strikeM -strike of the foliation; dipM -dip angle of the foliation; rakeM -the angle of line rake indicating the direction of lineation
%% v) The size and mesh of the medium
geom=[4*rbh,181,361];%geom = [rmedia, n1, n2]: a 1 X 3 vector that describes the geometry of the set of calculations to be performed
% rmedia -the maximu radial distance into the medium that the calculations are to extend;
% n1 -the number of points calculated along the radial distance;
% n2 -the number of points calculated along the azimuth.
%% General solution
% Function GeneralSolution solves the stress distributed on the contour of and near the borehole.
[sigmarr, sigmass, sigmazz, sigmasz, sigmarz, sigmars] = StressDistribution(rbh, traj, sigmas, orientS, pw, Cm, orientM, geom);
toc