% This function determines the Williams solution eigenvectors and
% eigenvalues for a 2 alpha angle wedge. The eigenvectors are evaluated at
% the angle phi, all angles inputs are assumed to be in radians. 
% The method followed is outlined by Barber in his book
% Elasticity. The function outputs are the Williams eigenvalues, in the
% structure EigenVal, and the eigenvectors in the EigenVec structure.
% 
% This script was written in MATLAB R2017b.
%
% The results produced are not displayed
% 
% Daniel Riddoch, 2021-02-15
% University of Oxford, Department of Engineering Science
function [EigenVal,EigenVec]=WilliamsEigenSol(alpha,phi)
%% Finding eigenvalues
% This section finds the solutions of the characteristic equations
if alpha>2.25   % Finding a pair of singular solutions
    EigenVal.M1=fzero(@(x) x.*sin(2.*alpha)+sin(2.*alpha.*x),0.5);
    EigenVal.M2=fzero(@(x) x.*sin(2.*alpha)-sin(2*alpha.*x),0.6);
else            % Finding a singular and a bounded solution
    EigenVal.M1=fzero(@(x) x.*sin(2.*alpha)+sin(2.*alpha.*x),0.5);
    EigenVal.M2=fzero(@(x) x.*sin(2.*alpha)-sin(2.*alpha.*x),1.5);
end
%% Temporary variables
% In this section, we introduce a few variables which will help make the
% rest of this function more readable
COSPA1=cos((EigenVal.M1+1).*alpha);
COSPP1=cos((EigenVal.M1+1).*phi);
COSNA1=cos((EigenVal.M1-1).*alpha);
COSNP1=cos((EigenVal.M1-1).*phi);
SINPA1=sin((EigenVal.M1+1).*alpha);
SINPP1=sin((EigenVal.M1+1).*phi);
SINNA1=sin((EigenVal.M1-1).*alpha);
SINNP1=sin((EigenVal.M1-1).*phi);
COSPA2=cos((EigenVal.M2+1).*alpha);
COSPP2=cos((EigenVal.M2+1).*phi);
COSNA2=cos((EigenVal.M2-1).*alpha);
COSNP2=cos((EigenVal.M2-1).*phi);
SINPA2=sin((EigenVal.M2+1).*alpha);
SINPP2=sin((EigenVal.M2+1).*phi);
SINNA2=sin((EigenVal.M2-1).*alpha);
SINNP2=sin((EigenVal.M2-1).*phi);
%% Finding eigenvalues
% In this section we give the formulae presented by Barber and use them to
% calculate the eigenvectors of the solution
EigenVec.rr1=(COSNA1.*COSPP1-((EigenVal.M1-3)./(EigenVal.M1+1)).*COSPA1.*COSNP1)./(COSPA1-COSNA1);
EigenVec.rr2=(SINNA2.*SINPP2-((EigenVal.M2-3)./(EigenVal.M2+1)).*SINPA2.*SINNP2)./(SINNA2-((EigenVal.M2-1)./(EigenVal.M2+1)).*SINPA2);
EigenVec.rtheta1=(SINNA1.*SINPP1-SINPA1.*SINNP1)./(SINNA1-((EigenVal.M1+1)./(EigenVal.M1-1)).*SINPA1);
EigenVec.rtheta2=(COSNA2.*COSPP2-COSPA2.*COSNP2)./(COSNA2-COSPA2);
EigenVec.thetatheta1=(COSNA1.*COSPP1-COSPA1.*COSNP1)./(COSNA1-COSPA1);
EigenVec.thetatheta2=(SINNA2.*SINPP2-SINPA2.*SINNP2)./(((EigenVal.M2-1)./(EigenVal.M2+1)).*SINPA2-SINNA2);
end