function [a,b,delta] = Hertz_contact(Fn,E1,E2,v1,v2,r1_I,r2_I,r1_II,r2_II)

%HERTZIAN POINT CONTACT solution according to Hamrock & Dowson approximated approach
%  
% I and II are 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)= r1_I;
r_matrix(2)= r2_I;
r_matrix(3)= r1_II;
r_matrix(4)= r2_II;
R_matrix= zeros(4,1); % matrix of input curvature 

for j=1:4
    if  r_matrix(j)==inf
        R_matrix(j)= 0;
    else
        R_matrix(j)= 1/r_matrix(j);
    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

%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-v1^2)/E1)+((1-v2^2)/E2) );     %inverse of the reduced elastic modulus
E= 1/Einv;                                     %reduced elastic modulus
k= 1.0339*(ry/rx)^0.636;                    %H&D ellipticity parameter
eta= 1.0003 + 0.5968*rx/ry;                 %H&D simplified elliptical integral
zeta= 1.5277 + 0.6023*log(ry/rx);           %H&D simplified elliptical integral

a= ((6*k^2*eta*Fn*r)/(pi*E))^(1/3);     %major axis of ellipse
b= ((6*eta*Fn*r)/(k*pi*E))^(1/3);       %minor axis of ellipse
delta= zeta*((4.5/(eta*r))*(Fn/(k*pi*E))^2)^(1/3);

%p_med= Fn/(pi*a*b);
%p_max= 3/2*p_med;

end

