clc
clear all
close all force

set(gca, 'color', [0 0.1 0.1])
set(gca, 'XColor', [0 0.5 0.5]);
set(gca, 'YColor', [0 0.5 0.5]);
axis equal;
grid on; hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This main FILE is dedicated to geometry. It:
%1.Creates gears profiles from user input
%2.Writes Excel FILE containing main data on gears (radii, tooth thickness,
%  ecc.)
%3.Meshes the gears together
%4.Evaluates meshing contact points, trapped volume, relief groove connection
%            areas, number of sealed volumes over fixed step angular
%            sampling positions
%5.Plots graphic results
%6.Writes Excel FILE containing main data on meshing contact points,  
%  trapped volume, relief groove connection areas, number of sealed volumes

%Please, run this before going to 'pressure_field.m'. 'pressure_field.m'
%reads 'geometry.m' output Excel FILES as input!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% USER INPUT (GEAR PROFILE) %%%%%%%%%%%%%%%%%%%%%%%%%%
Z = 10; %define number of teeth
PA = 20; %define Pressure angle [deg]
w = 33; %define Tooth Width [mm]
m = 4.69; %define gear module
ACDFac = 2.15;  %Addendum Circle Diameter Factor (STANDARD = 2) 
                %(see function 'MakeGear1.m')
DCDFac = 2.2;  %Dedendum Circle Diameter Factor (STANDARD = 2) 
               %(see function 'MakeGear1.m')

lres = 0.01; % Define resolution of gear profile drawing (best to keep it  
             % between 0.1 - 0.01 for decent accuracy + low running time
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% USER INPUT (CASING PROFILE) %%%%%%%%%%%%%%%%%%%%%%%%%%
Tup_case = 115; %define casing angular extension (delivery side) [deg]
Tdown = 107; %define casing angular extension (suction side) [deg]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% USER INPUT (PERIPHERAL GROOVE ANGULAR EXTENSION) %%%%%%%%%%%%%
%This is used only to visualize the peripheral groove extension on the
%meshing animation
Tgroove = 20; %define backend groove extension from horizontal [deg] (+)if going down
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% USER INPUT (RELIEF GROOVES PROFILES) %%%%%%%%%%%%%%%%%%%%%%%%%
Dup = 6.0; %define delivery relief groove vertical distance from centerline [mm]
           %(towards delivery side) [mm]    (please insert a positive value)
Ddown = 5.0; %define suction relief groove vertical distance from centerline
             %(towards suction side) [mm]   (please insert a positive value)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% USER INPUT (ANGULAR STEP FOR GEOMETRY EVALUATION) %%%%%%%%%%%%%% 
teta_step = 0.5; %angular step for meshing sampling [deg], keep it small to
                 %improve accuracy (teta_step =0.5 [deg] is advised)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% USER INPUT (TOLERANCE ON TRAPPED VOLUME SHAPE RECOGNITION) %%%%%%%% 
Tol_trap = m * 1e-02; %if trapped volume computation is aborted, 
                      %try a less strict tolerance (Tol_trap = m * 1e-02 is
                      %advised as default, seems to work in most cases)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Some minor adjustment to input values to make them coherent with the
%reference frame
Tgroove = Tgroove + Tup_case;
Tup = Tup_case - Tgroove;


%The following lines make two identical gears (driver and driven) with evolving
%teeth profiles using function 'MakeGear1.m' (see 'MakeGear1.m' for more
%details).
%Complete gear coordinates are stored in vectors g1x, g1y
%Tooth profile coordinates are stored in vectors tcx, tcy
%PCD1 : prime circle diam. [mm]
%ACD1 : addendum circle diam. [mm]
%DCD1 : dedendum circle diam. [mm]
%TT : tooth (mean) thickness [mm]
[g1x, g1y, tcx, tcy, PCD1, ACD1, DCD1, TT] = MakeGear1(Z, PA, m, lres, ACDFac, DCDFac);
g2x=g1x; g2y=g1y; PCD2=PCD1; ACD2=ACD1; DCD2=DCD1; %driven gear (2) is identical to driver gear (1)

figure(1)
plot( tcx, tcy, 'linewidth', 1.5) %plot single tooth profile for visual check

% EVALUATE GEOMETRIC DIMENSIONS
TH = .5*(ACD1-DCD1); %tooth height [mm]
TW = -2*tcx(1:length(tcx)/2); %vector containing tooth thickness along height [mm]
                              %(from root to tip)
TR = tcy(1:length(tcy)/2); %vector containing corresponding radial coordinates [mm]
                           %(from root to tip)
TW_r = TW(1); %tooth thickness @ root [mm]
TW_t = TW(end); %tooth thickness @ tip [mm]
GAP = (pi*DCD1 - Z*TW_r) / Z; % @ root [mm]
       %Calculated as [(Perimeter of the root diameter circle)-(root tooth 
       %thickness * # of teeth)]/(# of interteeth vanes). 
       %(# of interteeth vanes == # of teeth)
VA = (pi*(ACD1/2)^2 - polyarea(g1x, g1y)) / Z; %Interteeth vane area [mm^2]
       %Calculated as [(Area of the external gear diameter circle)-(area of
       %the solid gear)]/(# of teeth) . Area of the solid gear is
       %calculated w/ MATLAB function 'polyarea'
VV = VA*w; %Interteeth Vane Volume [mm^3]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The following lines are for identifying all Z tooth tip edges (X,Y) 
%coordinates (clockwise) (driver gear)
% xLTE, yLTE vectors contain coordinates of the left side edges
% xRTE, yRTE vectors contain coordinates of the right side edges
%This operation will be useful for tracking sealed volumes!!!

for ii=1:Z
   %identify left and right edges of tooth profile (contained in tcx, tcy)
   %and rotate them clockwise of consecutive angular pitches (Z times) 
   %to get all teeth tip edges
  [xLTE(ii), yLTE(ii)] = rotateZ(tcx(length(tcx)/2), ...
                                 tcy(length(tcy)/2), -2*pi*(ii-1)/Z); 
  [xRTE(ii), yRTE(ii)] = rotateZ(tcx(length(tcx)/2 + 1), ...
                                 tcy(length(tcy)/2 + 1), -2*pi*(ii-1)/Z);

%De-comment if you want to plot all gears tooth tip edges for visual check
%   figure(2)
%   plot(xLTE(ii), yLTE(ii), '.r', xRTE(ii), yRTE(ii), '.b' )
%   hold on 
%   pause(.5)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %de-activate some annoying MATLAB default warning messages when polyshapes
 %will be created
warning('off', 'MATLAB:polyshape:repairedBySimplify')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%The following lines are for identifying all right tooth side coordinates (clockwise)
%This will be useful for the following meshing of the driver and driven
%gears (right tooth sides of driver gear will be meshing w/ right tooth
%sides of driven gear !!!)

is_RS = zeros(length(g1x) - 1, 1); %initialize gera right sides identifiers
T = fix((length(g1x)-1) / Z) ; %T is the length of the portion of the gear coordinates vector
                               %corresponding to a whole angular pitch.
                               %fix(X) rounds the elements of X to the nearest integers
                               %towards zero. It is used to avoid tr being
                               %not an integer
%identifier vector is_RS is set at 1 if the corresponding gear coordinates vectors 
%are associated with one of the Z right tooth sides. (is_RS = 0 otherwise)
for ii = 1 : Z
    tr = (ii-1)*T;
    is_RS( length(tcx)/2 + 1 + tr: length(tcx) + tr ) = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
WB = 1.00*(PCD1/2 + PCD2/2); %wheelbase [mm]  %USER INPUT
PA = acosd(0.5*PCD1*cosd(PA)/(.5*WB)); %working pressure angle given
                                           %wheelbase [deg]
O1 = - WB / 2; O2 = -O1; %driver (1) and driven (2) gear axis x-axis coordinates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%WRITE OUTPUT .xlsx FILE 'Gprofile.xlsx'
%'Sheet 1' contains:
%teeth number, gear thickness, working pressure angle, pitch radius, base radius;
%tooth height, root and tip thickness, root clearance, Vane volume, tooth mean thickness
%'Sheet 2' contains:
%tooth profile radii and corresponding thicknesses
filename = 'Gprofile.xlsx'; 
                            % please change file name if you want to create
                            % different geometries (or delete the existing
                            % .xslx file)
string = ["Teeth number" "Gear thick." "Working pressure angle" "Pitch Radius" "Base radius"];
writematrix(string, filename, 'Sheet', 1, 'Range', 'A1')
              %                                          Rbase = Rp*cosd(alpha_w); %[mm]
writematrix([Z, w, PA, -O1/2, -O1/2*cosd(PA)], filename, 'Sheet', 1, 'Range', 'A2')
string = ["Tooth height" "Rooth thick." "Tip thick." "Root clearance" "Vane volume", "Mean Thick"];
writematrix(string, filename, 'Sheet', 1, 'Range', 'A3')
writematrix([TH, TW_r, TW_t, GAP, VV, TT], filename, 'Sheet', 1, 'Range', 'A4')
string = ["Radius" "Tooth thick"];
writematrix(string, filename, 'Sheet', 2, 'Range', 'A1')
writematrix([TR', TW'], filename, 'Sheet', 2, 'Range', 'A2')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The following lines are for adjusting the driven gear to be in meshing
%contact w/ driver gear. 
%Driver gear remains fixed here. Phase angle is evaluated w/ function
%'meshing_gears.m' and the driven gear (2) is then rotated according to the
%phase angle (counterclockwise) in order to be in contact w/ driver gear (1)
phase = meshing_gears(tcx, tcy, O1, O2, PA, Z);
[g2x, g2y] = rotateZ(g2x, g2y, phase); 

%The following lines are for rotating both meshed gears in order to be in
%the angular position where a second contact point occurs (i.e.
%trapped volume is just born): this will be our starting angular position
%for the following evaluations.
%teta0 is the angular rotation to be applied to the driver gear (clockwise)
%and to the driven gear (counterclockwise) in order to achieve such
%positioning.
%(See function 'contact0.m' for more details) 
[UCx, UCy, LCx, LCy, condotta, teta0] = contact0(2*abs(O1), ACD1, PA, Z, m, tcx, tcy);
%Warning message if (Rapporto di condotta) is <1 (i.e. at least 1 contact
%point is not guaranteed: in this case please reconsider wheelbase value
%before proceeding)
if condotta < 1, msgbox('CHECK YOUR INPUT PARAMETERS! GEAR MESHING IS NOT GARANTEED, \eps < 1','WARNING'), end
%rotate both driver (1) and driven (2) gears to the starting sampling angular position 
[g1x, g1y] = rotateZ(g1x, g1y, -teta0);
[g2x, g2y] = rotateZ(g2x, g2y, teta0);
[xLTE, yLTE] = rotateZ(xLTE, yLTE, -teta0);
[xRTE, yRTE] = rotateZ(xRTE, yRTE, -teta0);

%Gears are now correctly positioned at the starting sampling positioning!!!

%VISUAL PLOTTING OF CASING: casing (X,Y) coordinates for plotting (both
%left and right sides)
[Xcase_s, Ycase_s, Xcase_d, Ycase_d] = casing(Tup_case, Tdown, O1 , .5*ACD1);
%VISUAL PLOTTING OF PERIPHERAL GROOVE: peripheral groove (X,Y) coordinates for plotting
%(both left and right sides)
%NOTE: as the peripheral groove is a circumference arc, function 'casing.m' is used again
[Xgro_s, Ygro_s, Xgro_d, Ygro_d] = casing(Tup_case, -Tup, O1 , .5*ACD1*.9);
%VISUAL PLOTTING OF RELIEF GROOVES: delivery (up) and suction (down) relief
%grooves (X,Y) coordinates
[Xdis_up, Ydis_up] = discharge_shape(Dup, abs(O1), abs(O1)); %discharge relief groove
[Xdis_down, Ydis_down] = discharge_shape(-Ddown, -abs(O1), abs(O1)); %suction relief grooves

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %de-activate some annoying MATLAB default warning messages when polyshapes
 %will be created
warning('off', 'MATLAB:polyshape:repairedBySimplify')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%The following 'for cycle ii' is for evaluating relevant geometry
%properties during gear meshing. Sampling is done from starting angular
%position and for a whole angular pitch (properties are the same every
%angular pitch)

th1 = -deg2rad(teta_step); %angular sampling step form user input 
                           %is converted from [deg] to [rad]
for ii = 1 : ((360 / Z) / teta_step )
    TETA(ii) = (ii-1)*teta_step;    %current angular position w/ respect to start
    [g1x, g1y] = rotateZ(g1x, g1y, th1); %gears are rotated of an angular step
    [g2x, g2y] = rotateZ(g2x, g2y, -th1); %gears are rotated of an angular step
    G1x = g1x + O1; G1y = g1y; %gears Xcoordinates in pump's reference frame
    G2x = g2x + O2; G2y = g2y; %gears Xcoordinates in pump's reference frame

    [xLTE, yLTE] = rotateZ(xLTE, yLTE, th1);
    [xRTE, yRTE] = rotateZ(xRTE, yRTE, th1);
    
    %identify trapped volume polygon and build its coordinate vectors
    %Xtrap, Ytrap. Track contact points coordinates in vectors Xc, Yc.
    %(see 'create_polygon.m' for details)
    %Driver gear (1) and driven gear (2) coordinates are inputs.
    [Xtrap, Ytrap, Xc, Yc] = create_polygon( G1x, g1y, G2x, g2y, is_RS, PA, UCx, UCy, LCx, LCy, .1*length(tcx)/2, Tol_trap );
    [TA(ii) , TV(ii)]= evaluateArea( Xtrap, Ytrap, w );  %evaluate trapped area
                                                         %and corresponding
                                                         %volume given polygon
                                                         %coordinate vectors
   
    %the following lines are for evaluating the radial and angular
    %coordinates of the gear contact points
    if length(Xc) == 1  %if there is only one contact point
         %angular coord. of contact 1 in driver gear reference frame
         angle_contact(ii, 1) = rad2deg(mod(atan2(Yc, Xc-O1), 2*pi));
         %angular coord. of contact 2 does not exist (set to 10000)
         angle_contact(ii, 2) = 1e4;
         %radial coord. of contact 1 in driver gear reference frame
         r_contact(ii, 1) = sqrt(Yc^2 + (Xc-O1)^2);
         %radial coord. of contact 2 does not exist (set to 10000)
         r_contact(ii, 2) = 1e4;
         
         %radial coord. of contact 1 in driven gear reference frame
         %(evaluated as the module of the vectorial difference between the
         %vector from the driver wheel center to contact 1 and
         %the vector from the driver wheel center to the driven wheel
         %center
         vec1 = [-2*O1, 0]; vec2 = [Xc-O1, Yc];
         %radial coord. of contact 1 in driven gear reference frame
         r_contact2(ii, 1) = norm(vec2 - vec1 , 2);
         %radial coord. of contact 2 does not exist (set to 10000)
         r_contact2(ii, 2) = 1e4;
    else                               %if there are two contact points
         %angular coord. of contact 1 in driver gear reference frame
         angle_contact(ii, 1) = rad2deg(mod(atan2(Yc(1), Xc(1)-O1), 2*pi));
         %angular coord. of contact 2 in driver gear reference frame
         angle_contact(ii, 2) = rad2deg(mod(atan2(Yc(2), Xc(2)-O1), 2*pi));
         %radial coord. of contact 1 in driver gear reference frame
         r_contact(ii, 1) = sqrt(Yc(1)^2 + (Xc(1)-O1)^2);
         %radial coord. of contact 2 in driver gear reference frame
         r_contact(ii, 2) = sqrt(Yc(2)^2 + (Xc(2)-O1)^2);

         %radial coord. of contact 1 in driven gear reference frame
         %(evaluated as the module of the vectorial difference between the
         %vector from the driver wheel center to contact 1 and
         %the vector from the driver wheel center to the driven wheel
         %center
         vec1 = [-2*O1, 0]; vec2 = [Xc(1)-O1, Yc(1)]; vec3 = [Xc(2)-O1, Yc(2)];
         %radial coord. of contact 1 in driven gear reference frame
         r_contact2(ii, 1) = norm(vec2 - vec1 , 2);
         %radial coord. of contact 2 in driven gear reference frame
         r_contact2(ii, 2) = norm(vec3 - vec1 , 2);
    end
    
    %the following lines are to evaluate upper discharge area and lower 
    %discharge area
    if ii>1 && TV(ii-1) == 0  %if trapped volume does not exist
        TA(ii) = 0; TV(ii) = 0; UDA(ii) = 0; LDA(ii) = 0;
    else                      %if trapped volume exists
        [polyup, polydown, UDA(ii), LDA(ii)] = discharge_Area( Xtrap, Ytrap, ...
                                   Xdis_up, Ydis_up,Xdis_down, Ydis_down );
           %Function 'discharge_Area.m' returns the connection areas shapes
           %as polygons (polyup, polydown)
    end
    
    %the following lines are for tracking the coordinates of the
    %last left tooth edge in contact w/ the casing (xxl, yyl) and last
    %right tooth edge in contact w/ the casing (xxr, yyr) (both w/ respect 
    %to the driver gear reference frame) and the number of contemporary 
    %Sealed Volumes
    [NoSV(ii), xxl, yyl, xxr, yyr] = sealedVolumes(Tup_case, Tdown, xLTE, yLTE, xRTE, yRTE, Z); 
                                                                                    
    angle_seal(ii, 1) = rad2deg(mod(atan2(yyl, xxl), 2*pi));
             %from XY coordinates to angular position ( w/ respect to the 
             %driver gear reference frame), last sealed tooth
    angle_seal(ii, 2) = rad2deg(mod(atan2(yyr, xxr), 2*pi));
             %from XY coordinates to angular position ( w/ respect to the 
             %driver gear reference frame), first sealed tooth 

    %Plotting of the gears meshing; shape of trapped volume and delivery
    %and suction relief grooves are coloured (green for trapped volume,
    %yellow for relief grooves)
    figure(3)
    clf;
    set(gca, 'color', [0 0.1 0.1])
    set(gca, 'XColor', [0 0.5 0.5]);
    set(gca, 'YColor', [0 0.5 0.5]);
    axis equal;
    grid on; hold on;
    plot(G1x, g1y, 'w');
    plot(G2x, g2y, 'w');
    plot([PCD1, -PCD1]*tand(PA),  [-PCD1, PCD1], '--y', ...
          [UCx, LCx], [UCy, LCy], '-r')
    plot( Xcase_s, Ycase_s, '-w', Xcase_d, Ycase_d, '-w')
    plot(Xgro_s, Ygro_s, ':w', Xgro_d, Ygro_d, ':w')
    plot(Xc, Yc,  '.c', 'LineWidth', 1.25)
    plot(xxr +O1, yyr, 'om', xxl +O1, yyl, 'om', 'linewidth', 1.5)
    plot(Xdis_up, Ydis_up, '-b', Xdis_down, Ydis_down, '-b' )
    if length(Xtrap) > 1,  fill(Xtrap, Ytrap, 'g'), end
    if UDA(ii) > 0, fill(polyup.Vertices(:,1),polyup.Vertices(:,2), 'y') , end
    if LDA(ii) > 0, fill(polydown.Vertices(:,1),polydown.Vertices(:,2), 'y') , end
    title(['Number of Sealed Vane: ',num2str(NoSV(ii))]) 
    xlim([-1.5*WB, 1.5*WB])
    ylim([-PCD1, PCD1])
    xlabel('[mm]')
    ylabel('[mm]')
    pause(0.1)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%EVALUATE TRAPPED VOLUME AND TRAPPED VOLUME'S DERIVATIVES WITH RESPECT TO 
%ANGULAR POSITION (starting from first contact point)
VOL = zeros(1,  length(TV)); 
if  sum(TV~=0) ~= 0
    for ii = 1 : length(TV)
        if TV(ii) == 0, VOL = TV(1 : ii - 1); break, end
    end
end
%the following lines are simply for evaluationg the trapped volume's derivative
%(simple application of centered finite differences method)
dVol = zeros(1,  length(VOL));
dVol(1) = (VOL(2) - VOL(1)) / (-th1);
dVol(end) = (VOL(end) - VOL(end - 1)) / (-th1);
for ii = 2 : length(VOL) - 1
    dVol(ii) = (VOL(ii + 1) - VOL(ii - 1)) / (-2*th1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plotting trapped volume and its derivative (derivative is plotted for 
%angular positions where trapped volume exists only)
figure(4)
plot(TETA(1 : length(dVol)), dVol, TETA, TV)
legend('dV_{trap}/d\theta [mm^3/rad]', 'V_{trap} [mm^3]')
xlabel('\theta [deg]')
%Plotting delivery connection area and suction connection area along
%meshing (for angular positions where trapped volume exists only)
figure(5)
plot(TETA(1 : length(dVol)), UDA(1 : length(dVol)), ...
     TETA(1 : length(dVol)), LDA(1 : length(dVol)), ...
     TETA(1 : length(dVol)), TA(1 : length(dVol)) )
legend('Delivery connection Area [mm^2]', 'Suction connection Area [mm^2]', 'Trapped Area [mm^2]')
xlabel('\theta [deg]')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%WRITE OUTPUT .xlsx FILE 'GPgeometry.xlsx'
%It contains:
%1.sampling angular positions
%2.trapped volume trend
%3.its derivative
%4.delivery relief groove connection area
%5.suction relief groove connection area
%6.,7. angular coordinates of both gear contact points (w/ respect to
%      driver gear reference frame)
%8.,9. radial coordinates of both gear contact points (w/ respect to 
%      driver gear reference frame)
%10.,11. angular coordinates of last (Thseal1) and first (Thseal2) sealed
%        teeth tips (w/ respect to driver gear reference frame)
%12.,13. radial coordinates of both gear contact points (w/ respect to  
%        driven gear reference frame)
filename = 'GPgeometry.xlsx';
                              % please change file name if you want to
                              % create different geometries (or delete the
                              % existing  .xslx file)
string = ["THETA" "TrapVol" "dVOL/dTH" "DelivArea" "SuctArea"...
          "SealedVanes" "ThContact1" "ThContact2" "rContact1" "rContact2" "ThSeal1" "ThSeal2" "r2Contact1" "r2Contact2"];
writematrix(string, filename, 'Sheet', 1, 'Range', 'A1:N1')
writematrix(TETA',filename,'Sheet',1,'Range','A2')
writematrix(TV',filename,'Sheet',1,'Range','B2')
writematrix(dVol',filename,'Sheet',1,'Range','C2')
writematrix(UDA(1:length(dVol))',filename,'Sheet',1,'Range','D2')
writematrix(LDA(1:length(dVol))',filename,'Sheet',1,'Range','E2')
writematrix(NoSV',filename,'Sheet',1,'Range','F2')
writematrix([angle_contact, r_contact, angle_seal, r_contact2],filename,'Sheet',1,'Range','G2')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
