% 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
% theta2
% 
% 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=75;                                     % This is the number of points that will be calculated(default value is 75&
alpha=3*pi/4;                               % Wedge half-angle(radians)
theta1=pi/4;                                % Crack 1 angle(radians)
theta2=linspace(-pi/2,0,Pts);               % Crack 2 angle(radians)
cl1=1;                                      % Crack 1 length
cl2=1;                                      % 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 is value 500)
facl=200;                                   % Length of free surfaces for notch faces(default value is 200)
%% Calibration loop
% This section loops through the theta 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(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
% crack angle
Font=36;                                        % Font size
AxisLim=[min(theta2)./pi max(theta2)./pi -2 2]; % Axis limits for the plot
figure
% Plotting a and e entries
subplot(2,2,1)
plot(theta2./pi,RawCals(:,1),'-b','LineWidth',2)
hold on
plot(theta2./pi,RawCals(:,5),'-r','LineWidth',2)
xlabel('$\frac{\theta_2}{\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(theta2./pi,RawCals(:,2),'-b','LineWidth',2)
hold on
plot(theta2./pi,RawCals(:,6),'-r','LineWidth',2)
xlabel('$\frac{\theta_2}{\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(theta2./pi,RawCals(:,3),'-b','LineWidth',2)
hold on
plot(theta2./pi,RawCals(:,7),'-r','LineWidth',2)
xlabel('$\frac{\theta_2}{\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(theta2./pi,RawCals(:,4),'-b','LineWidth',2)
hold on
plot(theta2./pi,RawCals(:,8),'-r','LineWidth',2)
xlabel('$\frac{\theta_2}{\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)