clear all; close all; clc;

% :::: Hertz Contact: Universal point contact ::::
% Created by:   Maarten Dekker
% Date:         2013-04-09
% Version:      #01

% Parameters body 1
E1 = 210e9; % [N/m^2] Youngs modulus
v1 = 0.3;   % [-] Poisson ratio
r1 = 5e-3;  % [m] small raidus
R1 = 5e-3;  % [m] large radius

% Parameters body 2
E2 = 210e9; % [N/m^2] Youngs modulus
v2 = 0.3;   % [-] Poisson ratio
r2 = 5e-3;  % [m] small raidus
R2 = 5e-3; % [m] large radius

% Load
F = 10;     % [N]   

% Computations
E = ((1-v1^2)/2/E1 + (1-v2^2)/2/E2)^(-1); % equivalent Youngs Modulus [N/m^2]
R = R1*R2/(R1 + R2);                      % Equivalent large radius [m]
r = r1*r2/(r1 + r2);                      % Equivalent small radius [m]
w = R/r;                                  % check
if w < 1
    disp('Error: Large radii are not larger than small radii')
    clear all; pause
elseif w >= 1 && w <= 25
    a = 2^(-1/3)*w^(11/24)*(3*F*r/E)^(1/3); % ellipse half long axis [m]
elseif w > 25 && w <1e5
    a = 2^(5/3)/pi*w^(8/24)*(3*F*r/E)^(1/3); % ellipse half long axis [m]
else
    disp('Error: Large radii too large; use other Hertz contact')
    clear all; pause
end
b = 2^(-1/3)*w^(-4/24)*(3*F*r/E)^(1/3);   % ellipse half short axis [m]
d = a^2/2/R + b^2/2/r;                    % approach of bodies [m]
sigma_max = 1.5*F/pi/a/b;                 % maximum stress [N/mm^2]
sigma_av  = F/pi/a/b;
C=F/d;

x = -a:2*a/1e3:a;
y = b.*sqrt(1-(x/a).^2);
plot(x*1e6,y*1e6,'r',x*1e6,-y*1e6,'r'); axis equal
xlabel(['x = 2a = 2 * ',num2str(round(a*1e9)/1e3),' [\mum]']); 
ylabel(['y = 2b = 2 * ',num2str(round(b*1e9)/1e3),' [\mum]']); 
title(['\delta = ',num2str(round(d*1e9)/1e3),' \mum, \sigma_m_a_x = ',num2str(round(sigma_max/1e6)),' MPa, C = ',num2str(round(C)),' ^N/_m'])
