% This fuction calculates the slip stick boundary for a complete contact.
% This is done by optimising to minimise the error found by reintroducing
% the discarded equation from the bounded-bounded quadrature. The
% simultaneous and sequential dislocation methods are used in both cases
% with constants N and n(=5*N) being varied ecach time. The results are 
% plotted to show which method converges first.
%
% This code requires WilliamsEigenSol, FindTolSim, CalculatePhiSeq, 
% FindStresses and FindTolSeq to run and was written in Matlab R2017b.
%
% Daniel Riddoch, 03/02/23
% Department of Engineering Science, University of Oxford
% close all
clear
clc
%% Parameter specification
% Stress intensity factors
K1    = -1;
K2    = 1;
% Contact properties
parameter.alpha = 0.75*pi;                                  % Half angle of notch
parameter.fc    = 0.25;                                      % Coefficient of friction
% Crack properties
parameter.theta = pi-parameter.alpha;                       % Angle of crack
% Input for notch faces
parameter.facl  = 180;                                      % Length of free surfaces
parameter.norm  = 'mellin2';                                % Normalisation for free notch faces
%% Williams solution
[EigenVal,EigenVec]=WilliamsEigenSol(parameter.alpha,parameter.theta);  % Williams eigensolution
parameter.lam1 = EigenVal.M1;                                 % Williams eigenvalue - mode I   
parameter.lam2 = EigenVal.M2;                                 % Williams eigenvalue - mode II
parameter.g1   = EigenVec.rtheta1/EigenVec.thetatheta1;       % Realigned eigenvector
parameter.g2   = EigenVec.rtheta2/EigenVec.thetatheta2;       % Realigned eigenvector
parameter.d0   = ((EigenVec.thetatheta1*K1)/(EigenVec.thetatheta2*K2))^(1/(parameter.lam2-parameter.lam1));    % Characteristic length dimension
%% Varying N
Steps=25;                                            % Number of different values of coefficient of friction(Default value is 25)
Nmin=100;                                            % Minimum permitted coefficient of friction(Default is 100)
Nmax=2500;                                           % Maximum permitted coefficient of friction(Default is 2500)
NV=linspace(Nmin,Nmax,Steps);                         % Vector of coefficient of friction values
for i=1:Steps
    c0d0 = (-(parameter.fc-parameter.g1)/(parameter.fc-parameter.g2))^(1/(parameter.lam2-parameter.lam1));% Normalised slip zone length by Williams solution
    N     = round(NV(i));
    ID    = round(N/2);                                 % The number of the equation to exclude
    ID    = logical([zeros(ID-1,1);1;zeros(N-ID+1,1)]); % Logical vector
    %% Parameter
    % Creating parameter structure to transfer variables
    parameter.ID    = ID;
    parameter.N     = N;
    parameter.n     = N*5;
    %% Find slip zone size - Simultaneous method
    Options.Display = 'iter';
    [c,fval,exitflag(i),output] = fsolve(@(d) FindTolSim(parameter,d),2.5*c0d0,Options);
    cRatSim(i)  = c./c0d0
    %% Find slip zone size - Simultaneous method
    parameter.f=parameter.fc;
    Options.Display = 'iter';
    [c,fval,exitflag(i),output] = fsolve(@(d) FindTolSeq(parameter,d),2.5*c0d0,Options);
    cRatSeq(i)  = c./c0d0
end
%% Plotting
figure
hold on
plot(NV,cRatSim,'LineWidth',2)
plot(NV,cRatSeq,'LineWidth',2)
xlabel(strcat('$N$'),'interpreter','latex','FontSize',16)
ylabel(strcat('$\frac{c}{c_0}$'),'interpreter','latex','FontSize',16)
legend({'Simultaneous','Sequential'},'interpreter','latex','FontSize',16)
ax = gca;
set(ax, ...
  'Box'         , 'on'     , ...
  'TickDir'     , 'in'     , ...
  'TickLength'  , [.02 .02] , ...
  'TickLabelInterpreter'  , 'latex' , ...
  'FontSize'    , 16 , ...
  'XMinorTick'  , 'off'      , ...
  'YMinorTick'  , 'off'      , ...
  'XGrid'       , 'on'      , ...
  'YGrid'       , 'on'      , ...
  'XColor'      , [0 0 0], ...
  'YColor'      , [0 0 0], ...
  'LineWidth'   , 1.0         );