% This script calculates the calibration matrix entries for two cracks of
% length cl1 and cl2, at angles theta1 and theta2 respectively, placed in a
% wedge of angle 2*alpha. This script allows the user to vary the ratio of
% the crack length
% 
% This script was written in MATLAB R2017b and requires TwoCrackCalib
%
% The results produced are displayed in figure 6
% 
% Daniel Riddoch, 2021-03-01
% University of Oxford, Department of Engineering Science
clc
clear
close all
%% Inputs
% This section contains all the inputs to the problem, both physical and
% numerical.
Pts=75;                                     % This is the number of points that will be calculated(default value is 75)
alpha=0.75*pi;                              % Wedge half-angle(radians)
theta1=pi/4;                                % Crack 1 angle(radians)
theta2=-pi/4;                               % Crack 2 angle(radians)
cl1=1;                                      % Crack 1 length
cl2=linspace(0.1,10,Pts);                     % Vector of crack 2 input length(crack 1 length may be varied instead if desired)
N=100;                                      % Number of points on the crack face(default value is 100)
n=500;                                      % Number of Gauss points for notch faces(default value is 500)
facl=200;                                   % Length of free surfaces for notch faces(default value is 200)
%% Calibration loop
% This section loops through the alpha variable to calculate all the
% calibration matrix entries
RawCals=zeros(Pts,8);       % Initialising matrix
for j=1:Pts
    [a,b,c,d,e,f,g,h]=TwoCrackCalib(alpha,theta1,theta2,cl1,cl2(j),n,facl,N);             % Calculating calibration
    RawCals(j,:)=[a,b,c,d,e,f,g,h]';                                                      % Storing calibration data
end
%% Figures
% This section plots the calibration matrix entries against the varying
% wedge angle
Font=36;                                        % Font size
AxisLim=[min(log(cl2./cl1))./pi max(log(cl2./cl1))./pi -2 2]; % Axis limits for the plot
figure
% Plotting a and e entries
subplot(2,2,1)
plot(log(cl2./cl1),RawCals(:,1),'-b','LineWidth',2)
hold on
plot(log(cl2./cl1),RawCals(:,5),'-r','LineWidth',2)
xlabel('log$\frac{l_2}{l_1}$','interpreter','latex','FontSize',Font)
ylabel('Matrix coefficient','interpreter','latex','FontSize',Font)
legend({'a','e'},'interpreter','latex','FontSize',Font)
set(gca,'FontSize',Font)
axis(AxisLim)
% Plotting b and f entries
subplot(2,2,2)
plot(log(cl2./cl1),RawCals(:,2),'-b','LineWidth',2)
hold on
plot(log(cl2./cl1),RawCals(:,6),'-r','LineWidth',2)
xlabel('log$\frac{l_2}{l_1}$','interpreter','latex','FontSize',Font)
ylabel('Matrix coefficient','interpreter','latex','FontSize',Font)
legend({'b','f'},'interpreter','latex','FontSize',Font)
set(gca,'FontSize',Font)
axis(AxisLim)
% Plotting c and g entries
subplot(2,2,3)
plot(log(cl2./cl1),RawCals(:,3),'-b','LineWidth',2)
hold on
plot(log(cl2./cl1),RawCals(:,7),'-r','LineWidth',2)
xlabel('log$\frac{l_2}{l_1}$','interpreter','latex','FontSize',Font)
ylabel('Matrix coefficient','interpreter','latex','FontSize',Font)
legend({'c','g'},'interpreter','latex','FontSize',Font)
set(gca,'FontSize',Font)
axis(AxisLim)
% Plotting d and h entries
subplot(2,2,4)
plot(log(cl2./cl1),RawCals(:,4),'-b','LineWidth',2)
hold on
plot(log(cl2./cl1),RawCals(:,8),'-r','LineWidth',2)
xlabel('log$\frac{l_2}{l_1}$','interpreter','latex','FontSize',Font)
ylabel('Matrix coefficient','interpreter','latex','FontSize',Font)
legend({'d','h'},'interpreter','latex','FontSize',Font)
set(gca,'FontSize',Font)
axis(AxisLim)