function [xx_up, yy_up, Xc1, Yc1] = create_polygon( X1, Y1, X2, Y2, is_right1, PA, xu, yu, xl, yl, S, tol)
%This function builds the trapped volume shape as a CLOSED polygon and tracks
%the contact point(s).
%INPUT:
% X1, Y1, X2, Y2 : driver and driven gears coordinates in pump's reference frame
% is_right1 : vector containing identifiers of points belonging to teeth
%             right sides (is_right1==1) and not (is_right1==0)
% PA : pressure angle [deg]
% xu, yu : extreme upper contact point coord.
% xl, yl : extreme lower contact point coord.
% S : 'spacer' for contact points searching (see below)
% tol : tolerance to be applied when searching for contact point lying on
%       the Direct Line of Action of the gears
%OUTPUT: 
% xx_up, yy_up : polygon shape coordinates
% Xc1, Yc1 : upper and lower contact point coordinates 
%           (length Xc1 = 2, length Yc1 = 2 if there are 2 contact points, 
%            length Xc1 = 1, length Yc1 = 1 if there is only 1 contact point)

Xc1 = zeros(1,2); Yc1 = zeros(1,2); Xc2 = zeros(1,2); Yc2 = zeros(1,2);
coord1 = zeros(1,2); coord2 = zeros(1,2);

jj = 1; ii = 1;
while ii <= length(X1)
    %the following lines search for the point(s) belonging to any of the teeth
    %right sides and lying on the Direct Line of Action (DLA) (i.e the
    %effective contacts segment).
    %Exact matching is practically impossible as the gear profile is given
    %by points (so it is not continuous) so a certain tolerance on relative
    %distance is set to achieve the identification.
    if (Y1(ii) <= yu) && (Y1(ii) >= yl) && ( abs(X1(ii) + Y1(ii)*tand(PA)) < tol ) && (is_right1(ii) == 1)
               %when any tooth right side point satisfies the relative
               %tolerance its coordinates are loaded in vectros Xc1, Yc1
        coord1(jj) = ii;
        Xc1(jj) = X1(ii);
        Yc1(jj) = Y1(ii);
        jj = jj + 1;
        ii = ii + round(S); %S serves as a 'spacer' in order to go search for
                            %another contact point far enough from the one
                            %just found: this is because otherwise the function
                            %would identify as another contact the
                            %point immediately after the first one (being
                            %close enough to satisfy the same tolerance,
                            %tol)
    end
    if jj == 3, break ,end  %when two distinct contact points are found cycle is stopped
    ii = ii + 1;
end

Xc1 = Xc1(Xc1~=0);
Yc1 = Yc1(Yc1~=0);

 if length(Xc1) == 1, xx_up = 0; yy_up = 0; %if there is only one contact point
                                            %then polygon shape does not exist
 elseif length(Xc1) == 2 %there are 2 contact points and trapped volume exists
     for jj = 1:2
         %the following lines are for detecting the point belonging to the
         %teeth right sides of the driven gear (2) which are in contact
         %with the two contact points (again, exact matching is practically
         %impossible as the gear profile is given by points: points of driven
         %gear achiving minimum distance from expected upper and lower
         %contact points are chosen).
       for ii = 1 : length(X2)
           if (Y2(ii) <= 1.1*yu) && (Y2(ii) >= 1.1*yl) && abs(X2(ii) - Xc1(jj)) <= 1e-2 * Xc1(jj) && abs(Y2(ii) - Yc1(jj)) <= 1e-2 * Yc1(jj)
               coord2(jj) = ii;
               Xc2(jj) = X2(ii);
               Yc2(jj) = Y2(ii);
               break
           else dist(ii) = sqrt( (X2(ii) - Xc1(jj))^2 + (Y2(ii) - Yc1(jj))^2 );
           end
          [M,Ind] = min(dist); %points of driven gear achieving minimum distance
          coord2(jj) =  Ind;   %from expected contact points are chosen 
       end
     end
    
    %Build polygon contour: select all driver gear (1) points belonging in
    %between the upper and lower contact points
    if Y1(coord1(1)) < Y1(coord1(2)), coord1 = flip(coord1); end
    if coord1(1) < coord1(2)
        c1x = X1(coord1(1) : coord1(2));
        c1y = Y1(coord1(1) : coord1(2));
    else 
        c1x = [X1(coord1(1) : end), X1(1 : coord1(2))];
        c1y = [Y1(coord1(1) : end), Y1(1 : coord1(2))];
    end
    %Build polygon contour: select all driven gear (2) points belonging in
    %between the upper and lower contact points
    if Y2(coord2(2)) < Y2(coord2(1)), coord2 = flip(coord2); end
    if coord2(1) < coord2(2)
        c2x = X2(coord2(1) : coord2(2));
        c2y = Y1(coord2(1) : coord2(2));
    else 
        c2x = [X2(coord2(1) : end), X2(1 : coord2(2))];
        c2y = [Y2(coord2(1) : end), Y2(1 : coord2(2))];
    end
    %Finally build complete polygon shape by seaming driver and driven gear
    %portions together
    xx_up = [c1x, c2x];
    yy_up = [c1y, c2y];
   
end

end