% 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 angle
% alpha
% 
% This script was written in MATLAB R2017b and requires TwoCrackCalib
%
% The results produced are displayed in figure 7
% 
% 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=25;                                     % This is the number of points that will be calculated(default value is 75)
alpha=linspace(0.75*pi,0.9*pi,Pts);         % Wedge half-angle(radians)
theta1=alpha-pi/2;                          % Crack 1 angle(radians)
theta2=pi/2-alpha;                          % Crack 2 angle(radians)
cl1=1;                                      % Crack 1 length
cl2=0.75;                                   % Crack 2 length
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(j),theta1(j),theta2(j),cl1,cl2,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(alpha)./pi max(alpha)./pi -2 2]; % Axis limits for the plot
figure
% Plotting a and e entries
subplot(2,2,1)
plot(alpha./pi,RawCals(:,1),'-b','LineWidth',2)
hold on
plot(alpha./pi,RawCals(:,5),'-r','LineWidth',2)
xlabel('$\frac{\alpha}{\pi}$','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(RawCals./pi,RawCals(:,2),'-b','LineWidth',2)
hold on
plot(RawCals./pi,RawCals(:,6),'-r','LineWidth',2)
xlabel('$\frac{\alpha}{\pi}$','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(RawCals./pi,RawCals(:,3),'-b','LineWidth',2)
hold on
plot(RawCals./pi,RawCals(:,7),'-r','LineWidth',2)
xlabel('$\frac{\alpha}{\pi}$','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(RawCals./pi,RawCals(:,4),'-b','LineWidth',2)
hold on
plot(RawCals./pi,RawCals(:,8),'-r','LineWidth',2)
xlabel('$\frac{\alpha}{\pi}$','interpreter','latex','FontSize',Font)
ylabel('Matrix coefficient','interpreter','latex','FontSize',Font)
legend({'d','h'},'interpreter','latex','FontSize',Font)
set(gca,'FontSize',Font)
axis(AxisLim)