% This function calculates the calibration matrix entries for two cracks of
% length cl1, cl2 and angles theta1, theta2 respectively, placed in a wedge
% of internal angle 2*alpha, all angles measured in radians. This is done 
% using the number of points on the free surfaces, n, the multiplicative
% constant facl and the number of points on the crack faces N. The outputs
% of the function are the entries of the calibration matrix set out in 
% equation 25.
% 
% This script was written in MATLAB R2017b.
%
% The results produced are not displayed
% 
% Daniel Riddoch, 2021-03-02
% University of Oxford, Department of Engineering Science
function [a,b,c,d,e,f,g,h]=TwoCrackCalib(alpha,theta1,theta2,cl1,cl2,n,facl,N)
%% Eigensolution
% We begin by calculating the eigensolution for both cracks
[EigenVal1,EigenVec1]=WilliamsEigenSol(alpha,theta1);  % Williams eigensolution
[~,EigenVec2]=WilliamsEigenSol(alpha,theta2);  % Williams eigensolution
%% Creating parameter structure
% This section folds the variables needed by other functions into a single
% structure for ease.
parameter.N=N;
parameter.theta1=theta1;
parameter.theta2=theta2;
parameter.n=n;
parameter.facl=facl;
parameter.alpha=alpha;
parameter.lam1=EigenVal1.M1;
parameter.lam2=EigenVal1.M2;
parameter.cl1=cl1;
parameter.cl2=cl2;
parameter.frq11=EigenVec1.rtheta1;
parameter.frq12=EigenVec2.rtheta1;
parameter.frq21=EigenVec1.rtheta2;
parameter.frq22=EigenVec2.rtheta2;
parameter.fqq11=EigenVec1.thetatheta1;
parameter.fqq12=EigenVec2.thetatheta1;
parameter.fqq21=EigenVec1.thetatheta2;
parameter.fqq22=EigenVec2.thetatheta2;
%% Calculating Phi
% This section calculates phi, the solution to the integral equations, and
% then uses krenk interpolation to find the calibration matrix entries
[phi1,phi2] = CalculatePhi(parameter);
% Krenk interpolation -  first we determine the phi_1 values, then apply
% the formulae set out in eq 21-24.
phi1x1=0;    % Initialising phi1x1
phi1y1=0;    % Initialising phi1y1
phi2x1=0;    % Initialising phi2x1
phi2y1=0;    % Initialising phi2y1
phi1x2=0;    % Initialising phi1x2
phi1y2=0;    % Initialising phi1y2
phi2x2=0;    % Initialising phi2x2
phi2y2=0;    % Initialising phi2y2
% Looping to find phi_1 values
for i=1:N
    phi1x1 = phi1x1 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi1(2*i-1,1);
    phi1y1 = phi1y1 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi1(2*i,1);
    phi2x1 = phi2x1 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi2(2*i-1,1);
    phi2y1 = phi2y1 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi2(2*i,1);
    phi1x2 = phi1x2 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi1(2*i-1,2);
    phi1y2 = phi1y2 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi1(2*i,2);
    phi2x2 = phi2x2 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi2(2*i-1,2);
    phi2y2 = phi2y2 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi2(2*i,2);
end
% Krenk interpolation formulae
a=pi*sqrt(2)/cl1^(parameter.lam1-1)*phi1y1;
b=pi*sqrt(2)/cl1^(parameter.lam2-1)*phi1y2;
c=pi*sqrt(2)/cl1^(parameter.lam1-1)*phi1x1;
d=pi*sqrt(2)/cl1^(parameter.lam2-1)*phi1x2;
e=pi*sqrt(2)/cl2^(parameter.lam1-1)*phi2y1;
f=pi*sqrt(2)/cl2^(parameter.lam2-1)*phi2y2;
g=pi*sqrt(2)/cl2^(parameter.lam1-1)*phi2x1;
h=pi*sqrt(2)/cl2^(parameter.lam2-1)*phi2x2;
end
%% Calculate Phi function
% This function calculates the unknown in a pair of SIE's, the structure
% parameter eighteen inputs, outlined above, detailing the state of stress in
% the uncracked body, the numerical parameters, and the properties of the
% crack. The output is the matrices phi1 and phi2, which provides the solution to the
% integral equations set out in equations 36-39.
function [phi1,phi2] = CalculatePhi(parameter)
% Read parameter - reading into local variables those brought into the
% function in the structure
n=parameter.n;
N=parameter.N;
lam1=parameter.lam1;
lam2=parameter.lam2;
l1=parameter.cl1;
l2=parameter.cl2;
lR=l2/l1;
frq11=parameter.frq11;
frq12=parameter.frq12;
frq21=parameter.frq21;
frq22=parameter.frq22;
fqq11=parameter.fqq11;
fqq12=parameter.fqq12;
fqq21=parameter.fqq21;
fqq22=parameter.fqq22;
phi1=parameter.theta1;
phi2=parameter.theta2;
% Numerical quadrature for notch faces - implementing the quadrature
% formulae set out by Hills et al.
Index=(1:n)';
t=cos(pi.*(2.*Index-1)./(2*n+1));
s=cos((pi.*2.*Index)./(2*n+1));
Wnotch=2.*(1-s)./(2*n+1);
parameter.t=t;
parameter.s=s;
parameter.Wnotch=Wnotch;
% Numerical quadrature for crack faces - implementing the quadrature
% formulae set out by Hills et al.
Idx=(1:N)';
v=cos(pi.*(2.*Idx)./(2*N+1));
u=cos((pi.*(2.*Idx-1))./(2*N+1));
W=2.*(1+u)./(2*N+1);
parameter.v=v;
parameter.u=u;
parameter.W=W;
%% Kernel functions
% This section calculatesd the kernels for all the pairs of positions we
% will use. This first section is for dislocations both placed and observed
% on crack 1, with angle theta1. We also calculate the kernel for indivdual
% dislocation, we will use these later.
parameter.c=l1;
parameter.faclR=1;
parameter.theta=phi1;
parameter.rho=phi1;
[Fxyy1,Fxxy1,Fyyy1,Fyxy1] = WedgeKernels(parameter);
[Krtt1,Krrt1,Kttt1,Ktrt1] = PolarKernelsLT(v,phi1,u',phi1); 
% The second section is for dislocations both placed and observed
% on crack 2, with angle theta2
parameter.c=l2;
parameter.faclR=1;
parameter.theta=phi2;
parameter.rho=phi2;
[Fxyy2,Fxxy2,Fyyy2,Fyxy2] = WedgeKernels(parameter);
[Krtt2,Krrt2,Kttt2,Ktrt2] = PolarKernelsLT(v,phi2,u',phi2);
% The third section is for dislocations both placed on crack 1 with angle theta1 and observed
% on crack 2, with angle theta2. Note the multiplier lR on the observation
% points
parameter.c=l1;
Rparameter=parameter;
Rparameter.faclR=lR;
Rparameter.rho=phi1;
[RFxyy,RFxxy,RFyyy,RFyxy] = WedgeKernels(Rparameter);
[RKrtt,RKrrt,RKttt,RKtrt] = PolarKernelsLT(lR*(v+1)-1,phi2,u',phi1); 
% The fourth section is for dislocations both placed on crack 2 with angle theta2 and observed
% on crack 1, with angle theta1. Note the multiplier 1/lR on the
% integration points
parameter.c=l2;
parameterR=parameter;
parameterR.faclR=1/lR;
parameterR.theta=phi1;
[FxyyR,FxxyR,FyyyR,FyxyR] = WedgeKernels(parameterR);
[KrttR,KrrtR,KtttR,KtrtR] = PolarKernelsLT(1/lR*(v+1)-1,phi1,u',phi2); 
%% Matrix assembly
% Initialisation of matrices
A=zeros(4*N,4*N);                                                          % Matrix A
F=zeros(4*N,2);                                                            % Right hand side
% 'A' Matrix assembly - this is the coefficeint matrix. We use the
% individual dislocation kernels as well here, for the 'singular' term, as
% it is not always singular when the observation point is on one crack and
% the dislocation on the other
A(1:4:end,1:4:end)=pi.*W'.*(Fxxy1+Krrt1);
A(1:4:end,2:4:end)=pi.*W'.*(Fyxy1+Ktrt1);
A(1:4:end,3:4:end)=pi.*W'.*(FxxyR+KrrtR);
A(1:4:end,4:4:end)=pi.*W'.*(FyxyR+KtrtR);
A(2:4:end,1:4:end)=pi.*W'.*(Fxyy1+Krtt1);
A(2:4:end,2:4:end)=pi.*W'.*(Fyyy1+Kttt1);
A(2:4:end,3:4:end)=pi.*W'.*(FxyyR+KrttR);
A(2:4:end,4:4:end)=pi.*W'.*(FyyyR+KtttR);
A(3:4:end,1:4:end)=pi.*W'.*(RFxxy+RKrrt);
A(3:4:end,2:4:end)=pi.*W'.*(RFyxy+RKtrt);
A(3:4:end,3:4:end)=pi.*W'.*(Fxxy2+Krrt2);
A(3:4:end,4:4:end)=pi.*W'.*(Fyxy2+Ktrt2);
A(4:4:end,1:4:end)=pi.*W'.*(RFxyy+RKrtt);
A(4:4:end,2:4:end)=pi.*W'.*(RFyyy+RKttt);
A(4:4:end,3:4:end)=pi.*W'.*(Fxyy2+Krtt2);
A(4:4:end,4:4:end)=pi.*W'.*(Fyyy2+Kttt2);
% This section creates the matrix for the bilateral solution for each case
F(1:4:end,1)=-frq11.*(l1.*(v+1)./2).^(lam1-1);                    % Right hand side mode 1 eq 1
F(2:4:end,1)=-fqq11.*(l1.*(v+1)./2).^(lam1-1);                    % Right hand side mode 1 eq 2
F(3:4:end,1)=-frq12.*(lR.*l1.*(v+1)./2).^(lam1-1);                % Right hand side mode 1 eq 3
F(4:4:end,1)=-fqq12.*(lR.*l1.*(v+1)./2).^(lam1-1);                % Right hand side mode 1 eq 4
F(1:4:end,2)=-frq21.*(l1.*(v+1)./2).^(lam2-1);                    % Right hand side mode 2 eq 1
F(2:4:end,2)=-fqq21.*(l1.*(v+1)./2).^(lam2-1);                    % Right hand side mode 2 eq 2
F(3:4:end,2)=-frq22.*(lR.*l1.*(v+1)./2).^(lam2-1);                % Right hand side mode 2 eq 3
F(4:4:end,2)=-fqq22.*(lR.*l1.*(v+1)./2).^(lam2-1);                % Right hand side mode 2 eq 4
%% Calculate phi values
% This section solves the equations that we have generated. It also
% separates the phi matrix into smaller pieces for more easy use later.
phi=A\F;
phi1=zeros(2*N,2);                  % Crack 1 phi
phi2=zeros(2*N,2);                  % Crack 2 phi
for i=1:N
    phi1(2*i-1,:)=phi(4*i-3,:);
    phi1(2*i,:)=phi(4*i-2,:);
    phi2(2*i-1,:)=phi(4*i-1,:);
    phi2(2*i,:)=phi(4*i,:);
end
end
%% Wedge kernels
% This function calculates the value of the kernels of a dislocation in an
% arbitrary angle wedge(2*alpha) observation angle theta, core angle rho, using the method set out in paper 1
function [Fxyy,Fxxy,Fyyy,Fyxy] = WedgeKernels(parameter)
% Read parameter - reading into local variables those brought into the
% function in the structure
N=parameter.N;
cl=parameter.c;
faclR=parameter.faclR;
theta=parameter.theta;
rho=parameter.rho;
n=parameter.n;
facl=parameter.facl;
alpha=parameter.alpha;
t=parameter.t;
s=parameter.s;
Wnotch=parameter.Wnotch;
v=parameter.v;
u=parameter.u;
% We now change the normalisation of the coordinates to normalise by the
% dislocation location for each dislocation
x=faclR*cl./2.*(v+1);      % Observation points
xi=cl./2.*(u+1);           % Dislocation core points
z=xi;                      % Radius of dislocation
l=facl*z;                  % Length of free surfaces
% Renormalising by dislocation location
v=1-2.^(1-2.*x'./l);       % Renormalised integration points
u=1-2.^(1-2.*z./l);        % Renormalised collocation points
% Assemble matrix A - consider the matrix in 4 vertical blocks, and n
% horizontal blocks. Each of the n blocks is a repeating pattern of 4
A=zeros(4*n,4*n);
% Top vertical block
[Krtt,Krrt,Kttt,Ktrt]=PolarKernelsMT(t,alpha,s',alpha);
A(1:n,1:4:end)=pi.*Wnotch'.*Krtt;
A(1:n,2:4:end)=pi.*Wnotch'.*Kttt;
A(n+1:2*n,1:4:end)=pi.*Wnotch'.*Krrt;
A(n+1:2*n,2:4:end)=pi.*Wnotch'.*Ktrt;
% Second vertical block
[Krtt,Krrt,Kttt,Ktrt]=PolarKernelsMT(t,alpha,s',-alpha);
A(1:n,3:4:end)=pi.*Wnotch'.*Krtt;
A(1:n,4:4:end)=pi.*Wnotch'.*Kttt;
A(n+1:2*n,3:4:end)=pi.*Wnotch'.*Krrt;
A(n+1:2*n,4:4:end)=pi.*Wnotch'.*Ktrt;
% Third vertical block
[Krtt,Krrt,Kttt,Ktrt]=PolarKernelsMT(t,-alpha,s',alpha);
A(2*n+1:3*n,1:4:end)=pi.*Wnotch'.*Krtt;
A(2*n+1:3*n,2:4:end)=pi.*Wnotch'.*Kttt;
A(3*n+1:4*n,1:4:end)=pi.*Wnotch'.*Krrt;
A(3*n+1:4*n,2:4:end)=pi.*Wnotch'.*Ktrt;
% Bottom vertical block
[Krtt,Krrt,Kttt,Ktrt]=PolarKernelsMT(t,-alpha,s',-alpha);
A(2*n+1:3*n,3:4:end)=pi.*Wnotch'.*Krtt;
A(2*n+1:3*n,4:4:end)=pi.*Wnotch'.*Kttt;
A(3*n+1:4*n,3:4:end)=pi.*Wnotch'.*Krrt;
A(3*n+1:4*n,4:4:end)=pi.*Wnotch'.*Ktrt;
% Right hand side for radial and tangential dislocation - in four vertical
% blocks
Fr=zeros(4*n,N);
Ft=zeros(4*n,N);
% Top two blocks
[Krtt,Krrt,Kttt,Ktrt] = PolarKernelsMT(t,alpha,u',rho);          % Kernel functions
Fr(1:n,:)=-Krtt./(l'./2./((1-u').*log(2)));                                           % Normal traction at positive half-line
Fr(n+1:2*n,:)=-Krrt./(l'./2./((1-u').*log(2)));                                       % Shear traction at positive half-line
Ft(1:n,:)=-Kttt./(l'./2./((1-u').*log(2)));                                           % Normal traction at positive half-line
Ft(n+1:2*n,:)=-Ktrt./(l'./2./((1-u').*log(2)));                                       % Shear traction at positive half-line
% Bottom two blocks
[Krtt,Krrt,Kttt,Ktrt] = PolarKernelsMT(t,-alpha,u',rho);         % Kernel functions
Fr(2*n+1:3*n,:)=-Krtt./(l'./2./((1-u').*log(2)));                                     % Normal traction at negative half-line
Fr(3*n+1:4*n,:)=-Krrt./(l'./2./((1-u').*log(2)));                                     % Shear traction at negative half-line
Ft(2*n+1:3*n,:)=-Kttt./(l'./2./((1-u').*log(2)));                                     % Normal traction at negative half-line
Ft(3*n+1:4*n,:)=-Ktrt./(l'./2./((1-u').*log(2)));                                     % Shear traction at negative half-line
%% Radial direction
phi = A\Fr;
% Extract different distributions of dislocations
phirp = phi(1:4:end,:);           % Density at positve free surface in radial direction
phitp = phi(2:4:end,:);           % Density at positive free surface in tangential direction
phirn = phi(3:4:end,:);           % Density at negative free surface in radial direction
phitn = phi(4:4:end,:);           % Density at negative free surface in tangential direction
% Calculating the kernels we need to construct the stress from phi
Fxyy = zeros(N,N);
Fxxy = zeros(N,N);
for i = 1:N
    % Kernels for dislocation on the free surfaces
    [Krttp,Krrtp,Ktttp,Ktrtp] = PolarKernelsMT(v(i,:)',theta,s',alpha);
    [Krttn,Krrtn,Ktttn,Ktrtn] = PolarKernelsMT(v(i,:)',theta,s',-alpha);
    % Stress build up, phi*weights*kernels*pi*c/2 on both half-lines
    Fxyy(:,i)=(pi.*Wnotch'*(phirp(:,i).*Krttp' + phitp(:,i).*Ktttp').*cl./2)'+(pi.*Wnotch'*(phirn(:,i).*Krttn' + phitn(:,i).*Ktttn').*cl./2)';
    Fxxy(:,i)=(pi.*Wnotch'*(phirp(:,i).*Krrtp' + phitp(:,i).*Ktrtp').*cl./2)'+(pi.*Wnotch'*(phirn(:,i).*Krrtn' + phitn(:,i).*Ktrtn').*cl./2)';
end
%% Tangential direction
phi = A\Ft;
% Extract different distributions of dislocations
phirp = phi(1:4:end,:);           % Density at positive free surface in radial direction
phitp = phi(2:4:end,:);           % Density at positive free surface in tangential direction
phirn = phi(3:4:end,:);           % Density at negative free surface in radial direction
phitn = phi(4:4:end,:);           % Density at negative free surface in tangential direction
% Calculating the kernels we need to construct the stress from phi
Fyyy = zeros(N,N);
Fyxy = zeros(N,N);
for i = 1:N
    % Kernels for dislocation onthe free surfaces
    [Krttp,Krrtp,Ktttp,Ktrtp] = PolarKernelsMT(v(i,:)',theta,s',alpha);
    [Krttn,Krrtn,Ktttn,Ktrtn] = PolarKernelsMT(v(i,:)',theta,s',-alpha);
    % Stress build up, phi*weights*kernels*pi*c/2 on both half-lines
    Fyyy(:,i)=(pi.*Wnotch'*(phirp(:,i).*Krttp' + phitp(:,i).*Ktttp').*cl./2)'+(pi.*Wnotch'*(phirn(:,i).*Krttn' + phitn(:,i).*Ktttn').*cl./2)';
    Fyxy(:,i)=(pi.*Wnotch'*(phirp(:,i).*Krrtp' + phitp(:,i).*Ktrtp').*cl./2)'+(pi.*Wnotch'*(phirn(:,i).*Krrtn' + phitn(:,i).*Ktrtn').*cl./2)';
end
end
%% Kernel function
% We have two slightly different functions here. The underlying formulae are the same, 
% but different transform, Mellin or linear, result in slightly different functions. 
% Firstly, the Mellin transform verion. This function defines the kernel of a dislocation at (s+1,phi) observed
% at (t+1,theta)
function [Krtt,Krrt,Kttt,Ktrt] = PolarKernelsMT(t,theta,s,phi)
    Krtt = (s-1).^(-1).*(log((-2).*(s-1).^(-1)).^2-2.*cos(theta-1.*phi).*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1))+log((-2).*(t-1).^(-1)).^2).^(-2).*((4+cos(2.*(theta-1.*phi)))...
      .*log((-2).*((-1)+s).^(-1)).^2.*log((-2).*(t-1).^(-1))+log((-2).*(t-1).^(-1)).^3-2.*cos(theta-1.*phi).*log((-2).*(s-1).^(-1)).*(log((-2).*(s-1).^(-1)).^2+2.*log((-2).*(t-1).^(-1)).^2)).*sin(theta-1.*phi);
    Krrt = (s-1).^(-1).*(log((-2).*(s-1).^(-1))-1.*cos(theta-phi).*log((-2).*(t-1).^(-1))).*(log((-2).*(s-1).^(-1)).^2-2.*cos(theta-phi).*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1))+...
        log((-2).*(t-1).^(-1)).^2).^(-2).*(cos(2.*(theta-phi)).*log((-2).*(s-1).^(-1)).^2+log((-2).*(t-1).^(-1)).*((-2).*cos(theta-phi).*log((-2).*(s-1).^(-1))+log((-2).*(t-1).^(-1))));
    Kttt = (1/2).*(s-1).^(-1).*(log((-2).*(s-1).^(-1)).^2-2.*cos(theta+(-1).*phi).*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1))+log((-2).*(t-1).^(-1)).^2).^(-2).*(2.*log((-2).*(s-1).^(-1)).^3+(-1).*(5.*cos(theta-phi)...
        +cos(3.*(theta-1.*phi))).*log((-2).*(s-1).^(-1)).^2.*log((-2).*(t-1).^(-1))+2.*(1+2.*cos(2.*(theta-1.*phi))).*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1)).^2-2.*cos(theta-1.*phi).*log((-2).*(t-1).^(-1)).^3);
    Ktrt = (-1).*(s-1).^(-1).*log((-2).*(t-1).^(-1)).*(log((-2).*((-1)+s).^(-1)).^2-2.*cos(theta-1.*phi).*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1))+log((-2).*(t-1).^(-1)).^2).^(-2).*(cos(2.*(theta-phi))...
        .*log((-2).*(s-1).^(-1)).^2+log((-2).*(t-1).^(-1)).*((-2).*cos(theta-1.*phi).*log((-2).*(s-1).^(-1))+log((-2).*(t-1).^(-1)))).*sin(theta-1.*phi);
end
% Secondly, the linear transform verion. This function defines the kernel of a dislocation at (s+1,phi) observed
% at (t+1,theta)
function [Krtt,Krrt,Kttt,Ktrt] = PolarKernelsLT(t,theta,s,phi)
% Temporary variables
COSV=cos(theta-phi);
SINV=sin(theta-phi);
COS2V=cos(2.*(theta-phi));
% Kernel formulae
Krtt=(-1).*(2+s.*(2+s)+t.*(2+t)+(-2).*(1+s).*(1+t).*COSV).^(-2).*((-2).*(1+s).*(3+s.*(2+s)+2.*t.*(2+t)).*COSV+(1+t).*...
    (5+4.*s.*(2+s)+t.*(2+t)+(1+s).^2.*COS2V)).*SINV;
Krrt=((-1)-s+(1+t).*COSV).*(2+s.*(2+s)+t.*(2+t)+(-2).*(1+s).*(1+t).*COSV).^(-2).*((1+t).*(1+t+(-2).*(1+s).*COSV)+(1+s).^2.*COS2V);
Kttt=(1/2).*((1+s).^2+(1+t).^2+(-2).*(1+s).*(1+t).*COSV).^(-2).*((1+t).*(7+5.*s.*(2+s)+2.*t.*(2+t)).*COSV-(1+s).*...
    (4+2.*s.*(2+s)+2.*t.*(2+t)+4.*(1+t).^2.*COS2V-(1+s).*(1+t).*cos(3.*(theta-phi))));
Ktrt=(1+t).*(2+s.*(2+s)+t.*(2+t)+(-2).*(1+s).*(1+t).*COSV).^(-2).*((1+t).*(1+t+(-2).*(1+s).*COSV)+(1+s).^2.*COS2V).*SINV;
end