%the aim is to compute the ellipse sector instead of the circle sector

%one of the solution that I found is the following one: see the notebook,
%Must change the part of the code that has the hypotesis on the circle and
%substitute.

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

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

l= 0.016;           %[m] stroke length -> wear increment
L= 1000;            %[m] test duration - total sliding distance
N= round(L/l);      %overall number of duration spans -> number of cycles doubled

% e=0.95;       %Iset a guess value for the form of the elliptical area


%% 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('Hertz 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

s= zeros(N,1);          %initializing the incremental sliding distance vector
a= zeros(N,1);          %initializing the major contact axis vector
b= zeros(N,1);          %initializing the minor contact axis vector
delta= zeros(N,1);      %initializing the contact body convergence vector
d= zeros(N+1,1);        %initializing the wear track width vector
Aw= zeros(N,1);         %initializing the cross-sectional area of wear track vector
r_t= zeros(N,1);        %initializing the wear track profile radius vector
theta= zeros(N,1);      %initializing the wear track angular amplitude vector
arc= zeros(N,1);        %initializing the wear track arc profile length vector
hw= zeros(N,1);         %initializing the underfomed wear depth vector
hw_tot= zeros(N,1);     %initializing the wear depth vector + elastic deformation
corrections= zeros(N,2);    %initializing the correction flags matrix
%column 1 for r_t limited to r_s
%column 2 for d limited to linearized arc profile

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

a_t= zeros(N,1);
b_t= zeros(N,1);
delta_t= zeros(N,1);
hw_t= zeros(N,1);
hw_tot_t= zeros(N,1);
cos_tau=zeros(N,1);
Q=zeros(N,1);
M=zeros(N,1);
C=zeros(N,1);
r_eq= zeros(N,1);        %initializing the wear track profile radius vector

flag_1=zeros(N,1);
flag_2=zeros(N,1);
flag_3=zeros(N,1);

%A_0= pi*(a_0)^2; Hertzian contact area sphere-plane
%pm_0= F/A_0; Hertzian mean contact pressure


%%
% Wear increment loop
for j=1:N  %j -> stroke counter
    s(j)= j*l;  %cumulative sliding distance
    deltaAw= k(j)*F;
    
    % It is assumed that the wear volume is inscribed into a circular profile if elastic deformation is recovered
    % Circular segment area formula: A = 0.5 * r_t^2 * (2*theta - sen(2*theta))
    % Wear track radius: r_t= a_0 * 1/sen(theta)
    % Wear track angular amplitude -> theta
    % Wear track widht --> d
    % r_t, theta and d are unknown -> system of 2 equations for 3 unknowns: one must add an hypotesys on wear track width d
    
    if j==1 %First iteration (first stroke)
        a_0= ( 3*F/(4/d_s)*Er )^(1/3);      %Hertzian contact RADIUS sphere-plane at first contact (Shaeffler handbook formula)
        Aw(j)= deltaAw;                     %Cross sectional area of the wear track
        a(j)=a_0;
        d(j) =2*a_0;                        %I know the track amplitude from the first iteration -> the width of the sphere-plane hertzian contact
        f= @(x)( (6*Aw(j)/d(j)^2)*( 1 + (x^4)/36 - (x^2)/3 ) - x );     %equation to solve theta with: d(1)=2*a_0
        X0 = fzero(f, atan(a_0/r_s));                                   %solutions vector with 1 positive value only (I hope!!)
        theta(j) = X0(X0>0 && X0<pi);                                   %calculated wear track amplitude,  within its phisical limits
        r_t(j)= a_0 / sin(theta(j));
        flag_1(j)=1;
        
        Q(j)=a_0;
        M(j)=a_0;
        
        hw(j)=M(j);
        
        
        
    else %strokes else than the first
        
        Aw(j)= Aw(j-1) + deltaAw;                                     %Cumulative cross sectional area of the wear track
        
        
        %Numerical solution of the system r,theta
        syms q m
        X0=[Q(j-1)*1.001, M(j-1)/1.001];
        eqns = [(q^2-m^2)^(1/2) == e*q, Aw(j) == 3.14*q*m/2]; %numerical solution of the system of equations
        S = vpasolve(eqns,[q m], X0);        %S is the vector (?) of solutions
        
        Q(j)=S.q;                           %extracting the value of
        M(j)=S.m;                           %extracting the value of
        S.q;
        S.m;
        
        d(j)=2*Q(j);
        hw(j)=M(j);
        
        flag_1(j)=2;
    end
    
    %     %Wear track profile and wear depth calculation
    %     if  a(j)>1.001*r_s                              %Hertz thoery is assumed to hold
    %         hw(j)= m(j);          % wear depth
    %         hw_t(j)=r_t(j) * (1-cos(theta(j)));
    %         arc(j)= 2* r_t(j) * theta(j); % arc profile of the wear track
    %         r_eq(j)=sqrt(Q(j)*M(j));
    %         [a(j),b(j),delta(j)]= Hertz_contact(F,E_s,E_d,ni_s,ni_d,r_s,r_eq(j),r_s,inf); %Hertz contact calculation at the end of the stroke, pay attention on the right order in which
    %VERIFY IF CONSISTENT WITH VALUES FORM SCHAEFFLER HANDBOOK FORMULA --> AXIAL BALL BEARINGS
    %         hw_tot(j)= hw(j) + delta(j+1);
    
    %Hertz contact calculation at the end of the stroke from the Schaeffler handbook formulas
    %         [a(j),b(j),delta_t(j),cos_tau(j)]= Hertz_contact_table(F,E_s,E_d,ni_s,ni_d,r_s,r_t(j),r_s,inf,CosTauTable);
    %         hw_tot_t(j)= hw_t(j) + delta_t(j+1);
    %
    %         flag_2(j)=1;
    %
    %         if 2*a(j)<arc(j)
    %             d(j+1)= 2*a(j);     %HP: wear track amplitude of the following iteration is identified approximately as the major Hertz contact axis
    %             %but the theoretical major contact axis may expand too much if r_t(j)/r_s -> 1.001, according to Hertz theory
    %
    %             flag_3(j)=1;
    %
    %         else
    %             d(j+1)= arc(j);     %HP: wear track amplitude d is artificially limited to the linearized arc profile of the wear track
    %             corrections(j,2)=1; %flag to identify the correction for track amplitude
    %
    %
    %             flag_3(j)=2;
    %         end
    
    %     else %HP: wear track radius is never smaller that the sphere radius for phisical limits
    %         a(j)= 1.001*r_s;                              % Hertz theory does not hold anymore!
    %         corrections(j,1)=1;                             % flag to identified the correction for radius
    %         theta(j)=((3/2)*Aw(j)/Q(j)^2)^(1/3);          % the track radius is known a priori here --> function f is useless
    %         hw(j)= Q(j) * (1-cos(theta(j)));              % wear depth
    %         arc(j)= 2* M(j) * theta(j);                   % arc profile of the wear track
    %         d(j+1)= 2*M(j)*sin(theta(j));                 % wear track width is calculated here
    %         hw_tot(j)= hw(j) + delta(j);
    %
    %
    %         flag_2(j)=2;
    %     end
    
end


% i set a series of for cycle in order to extract the values of the
% coefficient of hertz contact and evaluate if they are correct or not

size_1=0;
size_2=0;

% for i=1:N
%
%     if corrections(i,1) == 0
%         size_1=size_1+1;
%     else
%         size_2=size_2+1;
%     end
% end

results=rand(size_1,2);  % initializing the resul matrix in the form [a, b, delta, a_t, b_t, delta_t]

for i=1:size_1
    results(i,1)=a(i);
    results(i,2)=b(i);
    %     results(i,3)=delta(i);
    %     results(i,4)=a_t(i);
    %     results(i,5)=b_t(i);
    %     results(i,6)=delta_t(i);
end


%% Results
close all
% figure()
% plot(s,(hw_tot*10^6),'DisplayName','Total wear penetration depth');
% title('Total wear penetration depth [um]');
% xlabel('Test duration [m]');

figure()
plot(s,(Aw*10^12),'DisplayName','Simulated cross-sectional area')
hold on
plot(WearTrackArea(:,1),WearTrackArea(:,2),'xr','DisplayName','Exprimental values')
title('Corss-sectional wear track area [um^2]')
xlabel('Test duration [m]')
legend('location','southeast')
hold off

% figure()
% loglog(s,(r_t*10^3),'DisplayName','Wear track profile radius')
% title('Wear track profile radius [mm]')
% xlabel('Test duration [m]')

figure()
plot(s(1:end-1),(d(2:end-1)*10^3),'DisplayName','Wear track amplitude')
title('Wear track width [mm]')
xlabel('Test duration [m]')

figure()
plot(s(2:end),(hw(2:end)*10^6),'DisplayName','Simulated wear track depth')
hold on
plot(WearTrackDepth(:,1),WearTrackDepth(:,2),'xr','DisplayName','Exprimental values')
title('Wear track maximum depth [um]')
xlabel('Test duration [m]')
legend('location','southeast')
hold off

% figure()
% plot(1:size_1, results(:,1)*10^6, 'DisplayName', 'Hamrock & Dawson')
% hold on
% plot(1:size_1, results(:,4)*10^6,'.','MarkerIndices',1:10:size_1, 'DisplayName', 'Hertz theory')
% title('Major contact axis [um]')
% xlabel('Stroke count')
% legend('location','northwest')
% hold off
%
% figure()
% plot(1:size_1, results(:,2)*10^6, 'DisplayName', 'Hamrock & Dawson')
% hold on
% plot(1:size_1, results(:,5)*10^6, '.','MarkerIndices',1:10:size_1, 'DisplayName', 'Hertz theory')
% title('Minor contact axis [um]')
% xlabel('Stroke count')
% legend('location','southwest')
% hold off

% figure()
% plot(1:size_1, results(:,3), 'DisplayName', 'Hamrock & Dawson')
% hold on
% plot(1:size_1, results(:,6), 'DisplayName', 'Hertz theory')
% title('contact body convergence')
% xlabel('Stroke count')
% legend('location','northeast')
% hold off

% figure()
% plot(1:size_1, hw_tot(1:size_1)*10^6, 'DisplayName', 'Hamrock & Dawson')
% hold on
% plot(1:size_1, hw_tot_t(1:size_1)*10^6,'r--', 'DisplayName', 'Hertz theory')
% title('Wear track maximum depth [um]')
% xlabel('Stroke count')
% legend('location','southwest')
% hold off
