function [xu, yu, xl, yl, epsilon, th0] = contact0(D, De, PA, Z, m, tcx, tcy)
%% Questa funzione calcola la posizione angolare della ruota 1 (motrice) 
%  quando si genera il volume di ingranamento
% INPUT: 
% 1. D : Diametro primitivo [mm]
% 2. De : Diametro di testa (tip) [mm]
% 3. PA : angolo di pressione [deg]
% 4. Z : numero di denti 
% 5. m : modulo ruota
% 6. tcx, tcy : profilo del dente (lato sx e dx)
%OUTPUT: 
% 1. xu, yu : coord contatto superiore
% 2. xl, yl : coord contatto inferiore
% 3. epsilon : rapporto di condotta (1 < epsilon < 2)
% 4. th0 : angolo @t0 [rad], si chiude il volume di ingranamento

rb = .5*D*cosd(PA);
L = sqrt( (.5*De)^2 - rb^2 ) - rb*tand(PA);
xu = -L*sind(PA); yu = L*cosd(PA);
xl = -xu; yl = -yu;
epsilon = 2*L / (m*pi*cosd(PA));

O1H = .5*D - abs(xu);
phi = atan( yu / O1H );
r = sqrt( O1H^2 + yu^2 );

tcx = tcx(length(tcx)/2 + 1 : end);
tcy = tcy(length(tcy)/2 + 1 : end);
rcx = sqrt( tcx.^2 + tcy.^2 ); 
thcx = atan( tcx./tcy );
% figure(100)
% plot( thcx, rcx, '*', [0, max(thcx)], [r,r])

th0=0; flag=0;
for ii = length(rcx) : -1 : 1
    if abs(rcx(ii) - r) < 1e-3 * r
      th0 = pi/2 - phi - thcx(ii);
      flag = 1;
      ii= length(rcx)+10;
    end
end

if flag == 0, msgbox('CANNOT ROTATE','ERROR'), end
end