function [angle] = meshing_gears(tcx, tcy, O1, O2, PA, Z)
%This function is used to rotate the driven gear (clockwise) and make it
%mesh with the driver gear: gears must have 1/2 contact point(s) and must
%not intersect.
% INPUT: 
% 1. tcx, tcy : profilo del dente (lato sx e dx)
% 2. O1, O2 : Xcoordinate of driver (1) and driven (2) gears in the pump's
%             reference frame
% 3. PA : working pressure angle [deg]
% 4. Z : teeth number 
%OUTPUT: 
% angle : rotation angle to be applied to driven gear (2)

%the following lines rotate the tooth profile (right side only) 90degree 
%clockwise for the driver gear (1) and 90degree counterclockwise for the
%driven gear (2), so that both teeth lies horizontally
[X1, Y1] = rotateZ(tcx(length(tcx)/2 + 1 : end), ...
                   tcy(length(tcy)/2 + 1 : end), -pi/2); %coordinates in driver gear 
                                                         %reference frame
X1 = X1 + O1;                         %coordinates in pump's
                                      %reference frame
[X2, Y2] = rotateZ(tcx(length(tcx)/2 + 1 : end), ...
                   tcy(length(tcy)/2 + 1 : end), pi/2);  %coordinates in driven gear
                                                         %reference frame
Xm = 0; Ym = 0; %contact point initialization

tol = abs(1e-3 * O1);

%Contact point coordinates evaluation: it lies on the Direct Line of Action
%(effective contacts segment)
for i = 1:length(X1)
    %the following lines are for detecting the driver gear tooth right
    %side point lying on the Direct Line of Action.
    %Exact matching is practically impossible as the tooth profile is given
    %by points (so it is not continuous) so a certain tolerance on relative
    %distance is set to achieve the identification.
    if ( abs(X1(i) + Y1(i)*tand(PA)) < tol )
        Ym = Y1(i);          %Xm, Ym store the coordinates of the 
        Xm = X1(i);          %corresponding contact point.
        i= length(X1)+10;    %When contact point is found cycle is stopped
    end
end
flag = 0;

%the following kines are for evaluating the required driven gear rotating
%angle for correct meshing. Angle is not known beforehand, but it is for
%sure between 0 and the angular pitch value (360/Z): 1000 intermediate
%angles are considered. 

teta = deg2rad(linspace(1, 360 / Z, 1000));

%The following 'for cycle k' tries every intermediate angle value until any
%of the points belonging to the tooth right side of the driven gear matches
%with the previously evaluated contact point (Xm,Ym)

for k = 1 : length(teta)
  [XX2, YY2] = rotateZ(X2, Y2, teta(k)); %rotation of driven gear tooth right
                                         %side (driven gear reference frame)
  XX2 = XX2 + O2;  %coordinates in pump's reference frame

  for i = 1 : length(X2)   %searching for the point matching contact point
                           %coordinates (Xm,Ym); when it is found the inner 
                           %cycle is stopped
    if abs(XX2(i) - Xm) < tol && abs(YY2(i) - Ym) < tol, flag = 1; break, end 
  end
    if flag == 1, break, end   %when point is found the outer for cycle is
                               %stopped too
end

%ERROR MESSAGE is displayed if something goes wrong
if flag == 0 ,  msgbox('CANNOT MESH','ERROR'), end

angle = teta(k); %the trial angle selected by the previous 'for cycle' is 
                 %returned 
end
