% This function takes as its inputs the angle of the wedge alpha, the angle
% and length of the crack, theta and cl, and 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 13, and the eigenvector
% structure.
% 
% This script was written in MATLAB R2017b.
%
% The results produced are not displayed
% 
% Daniel Riddoch, 2021-03-01
% University of Oxford, Department of Engineering Science
function [a,b,c,d,EigenVec]=SingleCrackCalib(alpha,theta,cl,n,facl,N)
%% Eigensolution
% This section calculates the eigensolution for the state of stress in the
% uncracked body
[EigenVal,EigenVec]=WilliamsEigenSol(alpha,theta);  % Williams eigensolution
grq1=EigenVec.rtheta1/EigenVec.thetatheta1;         % Realigned eigenvector
grq2=EigenVec.rtheta2/EigenVec.thetatheta2;         % Realigned eigenvector
%% Creating parameter structure
% This section folds the variables needed by other functions into a single
% structure for ease.
parameter.N=N;
parameter.c=cl;
parameter.theta=theta;
parameter.n=n;
parameter.facl=facl;
parameter.alpha=alpha;
parameter.lam1=EigenVal.M1;
parameter.lam2=EigenVal.M2;
parameter.g1=grq1;
parameter.g2=grq2;
%% Calculating Phi
% This section calculates phi, the solution to the integral equations, and
% then uses krenk interpolation to find the calibration matrix entries
[phi] = CalculatePhi(parameter);
% Krenk interpolation -  first we determine the phi_1 values, then apply
% the formulae set out in eq 11 and 12.
phix1=0;    % Initialising phix1
phiy1=0;    % Initialising phiy1
phix2=0;    % Initialising phix2
phiy2=0;    % Initialising phiy2
% Looping to find phi_1 values
for i=1:N
    phix1 = phix1 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi(2*i-1,1);
    phiy1 = phiy1 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi(2*i,1);
    phix2 = phix2 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi(2*i-1,2);
    phiy2 = phiy2 + 2/(2*N+1)*cot(pi/2*(2*i-1)/(2*N+1))*sin(pi*N*(2*i-1)/(2*N+1))*phi(2*i,2);
end
% Krenk interpolation formulae
a=pi*sqrt(2)*phiy1;
b=pi*sqrt(2)*phiy2;
c=pi*sqrt(2)*phix1;
d=pi*sqrt(2)*phix2;
end
%% Calculate Phi function
% This function calculates the unknown in a pair of SIE's, the structure
% parameter ten 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 matrix phi, which provides the solution to the
% integral equations set out in equations 34 and 35
function [phi] = 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;
g1=parameter.g1;
g2=parameter.g2;
% 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;
% Initialisation of matrices
A=zeros(2*N,2*N);                                                     % Matrix A - the coefficient matrix containing the kernels
F=zeros(2*N,1);                                                       % Right hand side - the bilateral tractions
[Fxyy,Fxxy,Fyyy,Fyxy] = WedgeKernels(parameter);                      % Kernel functions - calculating the kernels for the pairs of quadrature points
% Matrix A has four distinct parts, each calculated here
A(1:2:end,1:2:end)=pi.*W'.*Fxyy;
A(1:2:end,2:2:end)=pi.*W'.*(Fyyy+1./(v'-u)');
A(2:2:end,1:2:end)=pi.*W'.*(Fxxy+1./(v'-u)');
A(2:2:end,2:2:end)=pi.*W'.*Fyxy;
% Matrix A has four distinct parts, each calculated here
F(1:2:end,1)=-((v+1)./2).^(lam1-1);
F(2:2:end,1)=-g1.*((v+1)./2).^(lam1-1);
F(1:2:end,2)=-((v+1)./2).^(lam2-1);
F(2:2:end,2)=-g2.*((v+1)./2).^(lam2-1);
% Calculate phi values
phi=A\F;
end
%% Kernel function
% This function applies 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;
theta=parameter.theta;
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=cl./2.*(v+1);      % Observation points
xi=cl./2.*(u+1);     % Dislocation core points
l=facl*cl;           % Length of free surfaces
% Renormalising by dislocation location
v = 1-2.^(1-2.*x./l);
u = 1-2.^(1-2.*xi./l);
% 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]=PolarKernels(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]=PolarKernels(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]=PolarKernels(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]=PolarKernels(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);
Multip=(l./2./((1-u').*log(2)));
% Top two blocks
[Krtt,Krrt,Kttt,Ktrt] = PolarKernels(t,alpha,u',theta);          % Kernel functions
Fr(1:n,:)=-Krtt./Multip;                                           % Normal traction at positive half-line
Fr(n+1:2*n,:)=-Krrt./Multip;                                       % Shear traction at positive half-line
Ft(1:n,:)=-Kttt./Multip;                                           % Normal traction at positive half-line
Ft(n+1:2*n,:)=-Ktrt./Multip;                                       % Shear traction at positive half-line
% Bottom two blocks
[Krtt,Krrt,Kttt,Ktrt] = PolarKernels(t,-alpha,u',theta);         % Kernel functions
Fr(2*n+1:3*n,:)=-Krtt./Multip;                                     % Normal traction at negative half-line
Fr(3*n+1:4*n,:)=-Krrt./Multip;                                     % Shear traction at negative half-line
Ft(2*n+1:3*n,:)=-Kttt./Multip;                                     % Normal traction at negative half-line
Ft(3*n+1:4*n,:)=-Ktrt./Multip;                                     % Shear traction at negative half-line
% Having set up the left hand side and two right hand side matrices, we
% first consider the 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
% Initialising matrices
Fxyy = zeros(N,N);
Fxxy = zeros(N,N);
% Calculating kernels at the observation points
[Krttp,Krrtp,Ktttp,Ktrtp] = PolarKernels(v,theta,s',alpha);
[Krttn,Krrtn,Ktttn,Ktrtn] = PolarKernels(v,theta,s',-alpha);
for i = 1:N
    % Stress build up, phi*weights*kernels*pi*c/2 on both half-lines
    Fxyy(i,:)=(pi.*Wnotch'*(phirp.*Krttp(i,:)' + phitp.*Ktttp(i,:)').*cl./2)+(pi.*Wnotch'*(phirn.*Krttn(i,:)' + phitn.*Ktttn(i,:)').*cl./2);
    Fxxy(i,:)=(pi.*Wnotch'*(phirp.*Krrtp(i,:)' + phitp.*Ktrtp(i,:)').*cl./2)+(pi.*Wnotch'*(phirn.*Krrtn(i,:)' + phitn.*Ktrtn(i,:)').*cl./2);
end
% Next we consider the tangential direction
phi = A\Ft;
% Extract different distributions of dislocations
phirp = phi(1:4:end,:);           % Density at pos free surf in rad direction
phitp = phi(2:4:end,:);           % Density at pos free surf in tan direction
phirn = phi(3:4:end,:);           % Density at neg free surf in rad direction
phitn = phi(4:4:end,:);           % Density at neg free surf in tan direction
% Initialising matrices
Fyyy = zeros(N,N);
Fyxy = zeros(N,N);
% Calculating kernels at the observation points
[Krttp,Krrtp,Ktttp,Ktrtp] = PolarKernels(v,theta,s',alpha);
[Krttn,Krrtn,Ktttn,Ktrtn] = PolarKernels(v,theta,s',-alpha);
for i = 1:N
    % Stress build up, phi*weights*kernels*pi*c/2 on both half-lines
    Fyyy(i,:)=(pi.*Wnotch'*(phirp.*Krttp(i,:)' + phitp.*Ktttp(i,:)').*cl./2)+(pi.*Wnotch'*(phirn.*Krttn(i,:)' + phitn.*Ktttn(i,:)').*cl./2);
    Fyxy(i,:)=(pi.*Wnotch'*(phirp.*Krrtp(i,:)' + phitp.*Ktrtp(i,:)').*cl./2)+(pi.*Wnotch'*(phirn.*Krrtn(i,:)' + phitn.*Ktrtn(i,:)').*cl./2);
end
end
%% Kernel function
% This function defines the kernel of a dislocation at (s+1,phi) observed
% at (t+1,theta)
function [Krtt,Krrt,Kttt,Ktrt] = PolarKernels(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