function [nVol, Q_tip, Q_lat, Q_drain, Q_mesh, Q_rlf, Q_back, Q_leak] = etaVol(theta, p_seal, p_trap, Displ, rho)
%This function evaluates leakages trends along a complete angular pitch and
%pump volumetric efficiency.
%INPUT
%theta : vector containing angular sampling positions over a complete angular pitch
%p_seal : matrix containing sealed volumes pressures in corresponding
%         angular positions (#rows = maximum number of contemporary sealed volumes,
%         #columns = same length of theta vector)
%p_trap : vector containing trapped volume pressure in corresponding
%         angular positions (ony when trapped volume exists)
%Displ : theoretical pump displacement [mm3/rev]
% rho : fluid density (in delivery volume, @ delivery pressure)
%OUTPUT
%nVol : volumetric efficiency
% Q_tip, Q_lat, Q_drain, Q_mesh, Q_rlf, Q_back, Q_leak : vectors containing
%          leakages in angular sampling positions (same length of theta vector) [mm3/s]
%NOTE: if there is no backend groove, Q_back is an all null vector

%Global variables to be used in this function 
global Z
global w
global omega
global TC
global LC
global DR
global pa
global pD
global pm
global mu 
global BC
global GArea
global Gangle

%INPUT data from Excel FILE 
filename = 'Gprofile.xlsx'; 
M = readmatrix(filename,'Sheet', 1, 'Range','A4:G4');
TW_t = M(3); %tooth width at tip radius [mm]
MT = M(7); %tooth mean thickness [mm]
M = readmatrix(filename,'Sheet', 2);
TR = M(:, 1); %tooth radii vector [mm]

filename = 'GPgeometry.xlsx'; 
M =  readmatrix(filename,'Sheet', 1, 'Range', 'A2');
UDA = M(:, 4);  %Upper Discharge Area (connection area w/ discharge relief groove [mm2]
NoSV = M(:, 6); %Number of (contemporary) Sealed Volumes vector 
thcontact1 = M(:, 7);  %angular coord. of meshing contact #1 [deg] (driving gear ref. sys)
thseal1 = M(:, 11);  %angular coord. of last tooth in contact w/ casing [deg]
thseal2 = M(:, 12);  %angular coord. of last tooth in contact w/ casing [deg]

UDA = UDA(1 : length(p_trap)); %crop UDA vector to match p_trap vector length

TH = TR(end) - TR(1); %Tooth heigth (outer radius - root radius)

%Characteristic dimensions of meshing area leakage
bmesh = 0.5*(TR(end)+TR(1)) - TR(1); %meatus width [mm]
lmesh = MT;  %meatus length [mm]
hmesh = LC;  %meatus height [mm]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINITION OF LEAKAGES BASED ON GEOMETRY (TIP, SIDE, DRAINAGE,
% DELIVERY RELIEF GROOVE, BACKEND GROOVE)
Qtip = @(dp, ht) w*ht.^3/(12*mu*TW_t)*(dp) + .5*(w*ht*omega*TR(end));
                                                   %[mm3/s] !!!!!!!!!!!!!!
Qlat = @(dp)  TH*LC.^3/(12*mu*MT)*(dp) + .5*(TH*LC*omega*(.5*(TR(1)+TR(end))));
                                                   %[mm3/s] !!!!!!!!!!!!!!
Qdrain = @(p)  (.5*TR(1)*2*pi/Z)*LC.^3/(12*mu*(TR(1)-DR))*(p-pD);
                                                   %[mm3/s] !!!!!!!!!!!!!! 
Qrlf = @(Adel, p) (0.65*Adel.*(1e3*sqrt(2*abs(pm - p)./rho))).*(pm>=p) + ...
                   0.*(pm<p);
                                                   %[mm3/s] !!!!!!!!!!!!!!
Qbckgr = @(p, A) 0.65*A.*(1e3*sqrt(2*abs(pm - p)./rho)).*(pm>=p) + ...
                   0.*(pm<p);
                                                   %[mm3/s] !!!!!!!!!!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Initialize leakage vectors
Q_tip = zeros(1, length(theta)); Q_lat = Q_tip; Q_drain = Q_tip;
Q_mesh = Q_tip; Q_rlf = Q_tip; Q_back = Q_tip; Q_leak = Q_tip;

%Leakages evaluation at every sampling angular position
for ii = 1 : length(theta)

  %LAST SEALED TOOTH LEAKAGES (TIP AND SIDE) FOR BOTH WHEELS 
  if NoSV(ii) == max(NoSV) %in this case p_seal(end,ii)~=pm
      %compute radial clearances distribution
      ht = tipGap(NoSV(ii), TC, BC,  thseal1(ii), thseal2(ii));
      htip(ii) = ht(end); %last sealed tooth radial clearance

      %pressure delta across last sealed tooth
      Dp(ii) = pm - p_seal(end, ii);
      if Dp(ii) == 0 && ii<size(p_seal,2), Dp(ii) = pm - p_seal(end, ii+1); end
     
      %last tooth tip and side leakages [mm3/s]
      Q_tip(ii) = 2*Qtip(Dp(ii), htip(ii));
      Q_lat(ii) = 2*2*Qlat(Dp(ii)); %(*2 because both sides)
      if Q_tip(ii) < 0, Q_tip(ii) = 0; end
      if Q_lat(ii) < 0, Q_lat(ii) = 0; end
  else               %in this case p_seal(end,ii)==pm (i.e. the last sealed volume 
                     %opened at delivery volume); last actual sealed volume
                     %pressure is p_seal(end-1,ii)
      %compute radial clearances distribution
      ht = tipGap(NoSV(ii), TC, BC,  thseal1(ii), thseal2(ii));
      htip(ii) = ht(end);  %last sealed tooth radial clearance

      %pressure delta across last sealed tooth
      Dp(ii) = pm - p_seal(end-1, ii);
      if Dp(ii) == 0 && ii<size(p_seal,2), Dp(ii) = pm - p_seal(end-1, ii+1); end
      
      %last tooth tip and side leakages [mm3/s]
      Q_tip(ii) = 2*Qtip(Dp(ii), htip(ii));
      Q_lat(ii) = 2*2*Qlat(Dp(ii)); %(*2 because both sides)
      if Q_tip(ii) < 0, Q_tip(ii) = 0; end
      if Q_lat(ii) < 0, Q_lat(ii) = 0; end
  end

  %DRAINAGE LEAKAGES (SINGLE WHEEL)
  %how many delivery vanes are leaking to the drainage. Drainage leaking 
  %vanes are those between first gears contact point angular position and
  %last sealed tooth angular position
  if thcontact1(ii) <= 90, drain_open = thseal1(ii) - thcontact1(ii);
  elseif thcontact1(ii) >= 270, drain_open = thseal1(ii) + (360 - thcontact1(ii));
  end
  %TOTAL DRAINAGE LEAKAGE
  Q_drain(ii) = 2*2*ceil(drain_open/(360/Z))*Qdrain(pm); %(*2 because both sides)

  %MESHING ZONE LEAKAGE
  Q_mesh(ii) = 2*2*bmesh*hmesh.^3/(12*mu*lmesh)*(pm-pa); %[mm3/s] !!!!!

  %DELIVERY RELIEF GROOVE LEAKAGE
  if ii <= length(p_trap)  %when trapped volume exists
      Q_rlf(ii) = 2*Qrlf(UDA(ii), p_trap(ii)); %(*2 because both plates)
  else, Q_rlf(ii) = 0;
  end

  %BACKEND GROOVE LEAKAGE
  %connection area w/ sealed volumes vector
  Agr = groove_leak_area(NoSV(ii), thseal1(ii), thseal2(ii), GArea, Gangle(1), Gangle(2));
  %Backend leakage for every sealed volume (most of them have no connection
  %w/ backend groove so their leakage is actually zero)
  for kk = 1:length(Agr)
  Q_back_vec(kk) = Qbckgr(p_seal(kk, ii), Agr(kk));
  end
  Q_back(ii) = 2*2*sum(Q_back_vec); %(*2 because both plates, *2 because both wheels)

  %TOTAL LEAKAGE (both wheels)
  Q_leak(ii) = Q_tip(ii) + Q_lat(ii) + Q_drain(ii) + Q_mesh(ii) + Q_rlf(ii) + Q_back(ii); %[mm3/s] !!!!!
end

nVol = 1 - mean(Q_leak)/(Displ*omega/(2*pi)); % 1-(Mean leakage and mean theorical flow ratio)!

end