clc
clear all
close all

%% Geometric and material parameters
r_s= 0.005;         %Sphere radius [m]
d_s= 2*r_s;         %Sphere diameter [m]
E_s= 390*10^9;      %[Pa] Sphere elastic modulus
ni_s= 0.26;         %Sphere Poisson ration
r_lim=1.001*r_s;    %Cut off for "plane" [m]
r_cut=0.05;         %I set this value in order to reduce the amount of point in the model  [m]
E_d= 218*10^9;      %[Pa] Disc Elastic modulus
ni_d= 0.30;         %Disc Poisson ration

F= 10;              %[N] Applied vertical load
Er= 0.5*( (1-ni_s^2)/E_s + (1-ni_d^2)/E_d ); %Reduced elastic modulus of material pair

M=10000;

r_1=linspace(10000,10,M);
r_2=linspace(0.999*10,1,M);

r_plane= [inf, r_1, r_2, linspace(0.999,r_lim,M)];                    %Put together the two vectors
N=length(r_plane);                      %N is the number of elements of the total vector

error=0.99*1e-4;    % I set the error 99%of the precision of the table
%% Experimental wear law
% Archard coefficients from excel analysis of wear curves
k1= 8*10^(-15); % [m3/(m*N)] dimensional wear rate k = V/(s*F)- wear-in
L1= 23; %[m] end of wear-in
k2= 3*10^(-14); % [m3/(m*N)] dimensional wear rate k = V/(s*F)- wear phase I
L2= 175; %[m] end of wear phase I
k3= 5*10^(-14); % [m3/(m*N)] dimensional wear rate k = V/(s*F)- wear phase II

% k= zeros(N,1); %initializing the dimensional wear coefficient vector
% k( 1:round(L1/l) )= k1;
% k( (round(L1/l)+1) : round(L2/l) )= k2;
% k( (round(L2/l)+1) : end )= k3;

% Experimental data import
opts = delimitedTextImportOptions("NumVariables", 2);
opts.DataLines = [1, Inf]; % Specify range and delimiter
opts.Delimiter = "\t";
opts.VariableNames = ["VarName1", "VarName2"]; % Specify column names and types
opts.VariableTypes = ["double", "double"];
opts.ExtraColumnsRule = "ignore"; % Specify file level properties
opts.EmptyLineRule = "read";
opts = setvaropts(opts, ["VarName1", "VarName2"], "DecimalSeparator", ",");% Specify variable properties

% Import the data of wear track cross-sectional area
WearTrackArea = readtable('Wear track area.txt', opts);     %File must be in the same folder as the script
% Convert to numeric matrix
WearTrackArea = table2array(WearTrackArea);

% Import the data of wear track depth
WearTrackDepth = readtable('Wear track depth.txt', opts);   %File must be in the same folder as the script
% Convert to numeric matrix
WearTrackDepth = table2array(WearTrackDepth);

% Clear temporary variables
clear opts

% Table's data import, first have to change load options
opts = delimitedTextImportOptions("NumVariables", 5);
opts.DataLines = [1, Inf]; % Specify range and delimiter
opts.Delimiter = "\t";
opts.VariableNames = ["VarName1", "VarName2", "VarName3", "VarName4", "VarName5"]; % Specify column names and types
opts.VariableTypes = ["double", "double", "double", "double", "double"];
opts.ExtraColumnsRule = "ignore"; % Specify file level properties
opts.EmptyLineRule = "read";
opts = setvaropts(opts, ["VarName1", "VarName2", "VarName3", "VarName4", "VarName5"], "DecimalSeparator", ".");% Specify variable properties

% Import the table
CosTauTable = readtable('Schaeffler table.txt', opts);   %File must be in the same folder as the script
% Convert to numeric matrix
CosTauTable = table2array(CosTauTable);

% Clear temporary variables
clear opts

%% Initializing the vectors for interations containing the END-OF-STROKE pysical and geometrical quantities

a= zeros(N,1);              %initializing the major contact axis vector
b= zeros(N,1);              %initializing the minor contact axis vector
r_t= zeros(N,1);            %initializing the wear track profile radius vector

%initializing the major contact axis vector, as it's done for the Hamrock
%and dawson simplified formulae, it's done for the handbook's formulas
% _t is the subscript to indicate the values came from the Table of
% Schaeffler's handbook

a_t= zeros(N,1);
b_t= zeros(N,1);
cos_tau=zeros(N,1);
den=zeros(N,1);
m_zeta=zeros(N,1);
m_eta=zeros(N,1);
m_zeta_eta=zeros(N,1);
m_psi_zeta=zeros(N,1);
zeta=zeros(N,1);
eta=zeros(N,1);
zeta_eta=zeros(N,1);
psi_zeta=zeros(N,1);
k_sum=zeros(N,1);
curvatura=zeros(N,1);

%% Evaluation of the contact area

for j=1:N
           
    r_t(j)=r_plane(j);
    [a(j),b(j)]= Hertz_contact(F,E_s,E_d,ni_s,ni_d,r_s,r_t(j),r_s,inf);                                         %Hertz contact calculation, approximation Hamrock&Dawson
    
    %HERTZIAN POINT CONTACT solution according to tabulate value functions of
    %cos tau, from Schaeffler handbook
    %
    % I and II ar the subscripts for the two perpendicular planes of curvature (load direction belogns to them)
    % 1 and 2 are the subscripts for the two contacting bodies
    % Fn is the load normal to the contact
    
    % E1 and E2 are entered in [Pa]
    % r are entered in [m] WITH CONCAVE/CONVEX SIGNS !!!
    % type "inf" for r if the corresponding curvature is a plane
    
    %    E is the reduced elastic modulus of the contact pair
    % 1/E = 0.5*( ((1-v1^2)/E1)+((1-v2^2)/E2) ) = Einv
    
    % R values stand for curvature (inverse of radius) quatities [1/m]
    % r values stand for radius quantities [m]
    
    r_matrix= zeros(4,1); % matrix of input radii
    r_matrix(1)= r_s;
    r_matrix(2)= r_t(j);
    r_matrix(3)= r_s;
    r_matrix(4)= inf;
    R_matrix= zeros(4,1); % matrix of input curvature
        
    for w=1:4
        if  r_matrix(w)==inf
            R_matrix(w)= 0;
        else
            R_matrix(w)= 1/r_matrix(w);
        end
    end
    
    R= R_matrix(1) + R_matrix(3) + R_matrix(2) + R_matrix(4);    % inverse of the overall reduced radius of curvature
    R_I= R_matrix(1) + R_matrix(2);                              % inverse of the reduced radius of curvature in plane I
    R_II= R_matrix(3) + R_matrix(4);                             % inverse of the reduced radius of curvature in plane II
    %  cos_tau = abs(R_matrix(1) - R_matrix(3) + R_matrix(2) - R_matrix(4))/R ;
    %GENERAL FORM
    
    %     cos_tau=1/(2*r2_I/r1_I-1);
    %     k=2/r1_I*(2-(0.5*r1_I/r2_I));
    
    %Identification of planes I and II as conventional x,y planes
    
    if R_I >= R_II
        Rx= R_I;
        Ry= R_II;
    else
        Rx= R_II;
        Ry= R_I;
    end
    
    r= 1/R;     %[m] overall reduced radius of curvature
    rx= 1/Rx;   %[m] reduced radius of curvature in plane x
    ry= 1/Ry;   %[m] reduced radius of curvature in plane y
    
    Einv= 0.5*( ((1-ni_s^2)/E_s)+((1-ni_d^2)/E_d) );     %Done it before
    E= 1/Einv;                                     %reduced elastic modulus
    
    cos_tau_exact=CosTauTable(:,1);
    % Find the right or the nearest value of cos tau from the table, then
    % extract the correct values zeta, eta , ecc.
    
    cos_tau(j)=1/((2*r_t(j)/r_s)-1);            %cos tau compute for an axial ball bearing
    k_sum(j)=2/r_s*(2-(r_s/r_t(j)));                      %sum of k compute for an axial ball bearing
    
    difference=cos_tau(j)-cos_tau_exact;
    
%     [MIN, i]= min(abs(difference));    % provare con un'indicizzazione tipo (>0 && )
    
    true=difference>=0;
    S=sum(true);
    
   % the table has this form, 1 column cos tau, 2 column zeta, 3 eta, 4 zeta*eta 
   % 5 psi/zeta
   
   
   % Now the table is from 0 --> 1. Must evaluate the right choices for the
   % interpolation between two cells.
   
   if abs(difference(S))<error
       
       zeta(j)=CosTauTable(S,2);        %right value of zeta
       eta(j)=CosTauTable(S,3);          %eta correct
       zeta_eta(j)=CosTauTable(S,4);    % zeta*eta correct
       psi_zeta(j)=CosTauTable(S,5);    % psi/zeta correct
%        disp('if')
       
   else
       
       den(j)=(CosTauTable(S+1,1)-CosTauTable(S,1));              %extracting the RIGHT coefficient from the table
       m_zeta(j)=(CosTauTable(S+1,2)-CosTauTable(S,2))/den(j);           %I compute the slope of the linear function
       m_eta(j)=(CosTauTable(S+1,3)-CosTauTable(S,3))/den(j);
       m_zeta_eta(j)=(CosTauTable(S+1,4)-CosTauTable(S,4))/den(j);
       m_psi_zeta(j)=(CosTauTable(S+1,5)-CosTauTable(S,5))/den(j);
       
       %Now I have to compute the right value
       
       zeta(j)=m_zeta(j)*(cos_tau(j)-CosTauTable(S,1))+CosTauTable(S,2);        %right value of zeta    MUST EVALUATED THIS POINT
       eta(j)=m_eta(j)*(cos_tau(j)-CosTauTable(S,1))+CosTauTable(S,3);          %eta correct
       zeta_eta(j)=m_zeta_eta(j)*(cos_tau(j)-CosTauTable(S,1))+CosTauTable(S,4);    % zeta*eta correct
       psi_zeta(j)=m_psi_zeta(j)*(cos_tau(j)-CosTauTable(S,1))+CosTauTable(S,5);    % psi/zeta correct
       
%        disp('else')
   end
    
    
    a_t(j)= zeta(j)*(3*F/(E*k_sum(j)))^(1/3);                      %major axis of ellipse, i've removed the k_sum
    b_t(j)= eta(j)*(3*F/(E*k_sum(j)))^(1/3);                       %minor axis of ellipse
%     delta_t(j)= psi_zeta(j)*(9*F^2*k(j)/(8*E^2))^(1/3);      %contact body convergence
    % p_med= Fn/(pi*a*b);
    % p_max= 3/2*p_med;
    p_max(j)=(1/zeta_eta(j))*(3*F*(E*k_sum(j))^2)^(1/3);
    cos_tau;
        
    curvatura(j)=1/r_t(j);
end



%% Results

figure()
plot(curvatura, a*10^6, 'DisplayName', 'Hamrock & Dawson')
hold on
plot(curvatura, a_t*10^6, 'DisplayName', 'Hertz theory')
title('Major contact axis [um]')
xlabel('Curvature')
legend('location','best')
hold off

figure()
plot(curvatura, b*10^6, 'DisplayName', 'Hamrock & Dawson')
hold on
plot(curvatura, b_t*10^6, 'DisplayName', 'Hertz theory')
title('Minor contact axis [um]')
xlabel('Curvature')
legend('location','best')
hold off



figure()
plot(curvatura(4:end-3), a(4:end-3)*10^6, 'DisplayName', 'Hamrock & Dawson')
hold on
plot(curvatura(4:end-3), a_t(4:end-3)*10^6, 'DisplayName', 'Hertz theory')
title('Major contact axis [um]')
xlabel('Curvature')
legend('location','best')
hold off

figure()
plot(curvatura(4:end-3), b(4:end-3)*10^6, 'DisplayName', 'Hamrock & Dawson')
hold on
plot(curvatura(4:end-3), b_t(4:end-3)*10^6, 'DisplayName', 'Hertz theory')
title('Minor contact axis [um]')
xlabel('Curvature')
legend('location','best')
hold off



















%Another solution for explain the result, I use loglog in order to see if
%the results can be closer


% figure()
% loglog(r_t(end:-1:1), 10^6*a(end:-1:1), 'DisplayName', 'Hamrock & Dawson')
% hold on
% loglog(r_t(end:-1:1), 10^6*a_t(end:-1:1), 'DisplayName', 'Hertz theory')
% title('Major contact axis [um]')
% xlabel('Plane radius [m]')
% legend('location','best')
% hold off
% 
% figure()
% loglog(r_t(end:-1:1), 10^6*b(end:-1:1), 'DisplayName', 'Hamrock & Dawson')
% hold on
% loglog(r_t(end:-1:1), 10^6*b_t(end:-1:1), 'DisplayName', 'Hertz theory')
% title('Minor contact axis [um]')
% xlabel('Plane radius [m]')
% ylabel('a, a_t [m]')
% legend('location','best')
% hold off
