function [xu, yu, xl, yl, epsilon, th0] = contact0(D, De, PA, Z, m, tcx, tcy)
%This function evaluates the coordinates of the extreme first and second 
%gear contacts.
%It also evaluates the angular rotation needed to reach the positioning
%where a second contact just occurred (i.e. the trapped volume is just
%born), from the default meshed gears positioning.
% INPUT: 
% 1. D : pitch diameter [mm]
% 2. De : outer diameter (tip) [mm]
% 3. PA : working pressure angle [deg]
% 4. Z : teeth number
% 5. m : gear module
% 6. tcx, tcy : tooth profile coordinates 
%OUTPUT: 
% 1. xu, yu : extreme upper contact coordinates
% 2. xl, yl : extreme lower contact coordinates
% 3. epsilon : rapporto di condotta (1 < epsilon < 2)
% 4. th0 : angular rotation for correct positining of meshing gears at
%          starting angular sampling position [rad]

rb = .5*D*cosd(PA);  %base radius
L = sqrt( (.5*De)^2 - rb^2 ) - rb*tand(PA); %half length of the effective 
                                            %contact segment
xu = -L*sind(PA); yu = L*cosd(PA);  %corresponding upper contact coordinates
                                    %in the pump's reference frame
xl = -xu; yl = -yu;
epsilon = 2*L / (m*pi*cosd(PA)); %formula to calculate (rapporto di condotta) 

O1H = .5*D - abs(xu);   %upper contact coordinates in driver gear reference frame
phi = atan( yu / O1H );   %upper contact angular coord. in driver gear reference frame
r = sqrt( O1H^2 + yu^2 ); %upper contact radial coord. in driver gear reference frame

tcx = tcx(length(tcx)/2 + 1 : end); %first tooth profile coordinates (right side)
tcy = tcy(length(tcy)/2 + 1 : end);
rcx = sqrt( tcx.^2 + tcy.^2 );      %corresponding radii and angular position
thcx = atan( tcx./tcy );            %(driver gear reference frame)

%the following lines are for determining required angular rotation.
th0=0; flag=0;
for ii = length(rcx) : -1 : 1   %cycle for every tooth right side point
    if abs(rcx(ii) - r) < 1e-3 * r %if the radius of the point corresponds to the
                                   %radius of the upper contact then we
                                   %need to rotate the driver gear of a
                                   %certain angular value in order to make
                                   %said point match the extreme upper
                                   %contact point
      th0 = pi/2 - phi - thcx(ii); %required angular rotation is calculated here:
                                   %rotating the driver gear of th0
                                   %(clockwise) will make the tooth right
                                   %side touch the extreme upper contact
                                   %point.
      flag = 1;
      ii= length(rcx)+10;   %when wanted point is found cycle is stopped
    end
end

%ERROR MESSAGE if things go wrong for some reason
if flag == 0, msgbox('CANNOT ROTATE','ERROR'), end 
end
