function dpdt = odesealfcn(t, p, N, index)
%This function sets the ODE system for Sealed Volumes pressure evolution to
%be solved by ode15s integrating solver.
%INPUT:
% t : integrating variable of ode15s (in this case it is an angular position)
% p : variable to be integrated by ode15s (in this case it a length N
%     vector, each slot is associated with a sealed volume)
% N : Number of contemporary Sealed Volumes; this variable determine ODE
%     system dimension
%index : vector containing the indexes of angular positions where the Number 
%        of Sealed Volumes increases/decreases of 1
%OUTPUT:
%dpdt : function vector containing N governing equations for the N
%       contemporary Sealed Volumes (to be solved by ode15s)

%Global variable used in this function
global pcav 
global rhocav
global Z
global w
global omega
global TC
global LC
global DR
global pa
global pD
global pm
global BM
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]
VV = M(5); %Inter teeth Vane Volume [mm3]
MT = M(7); %Mean tooth thickness [mm]
M = readmatrix(filename,'Sheet', 2);
TR = M(:, 1);  %tooth profile radii and thickness vectors [mm]
filename = 'GPgeometry.xlsx'; 
M =  readmatrix(filename,'Sheet', 1, 'Range', 'A2');
TETA = M(:, 1); %sampling angular positions vector [deg]
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]

%Linear interpolation of thseal1 and thseal2: function is designed to have
%a neat discontinuity when the number of Sealed Volumes jumps from N to N+1
%and from N+1 to N (or from N to N-1 and from N-1 to N)
th1 = @(t) polyval(polyfit(TETA(1:index(1)-1), thseal1(1:index(1)-1), 1), t).*(t < TETA(index(1))) +...
           polyval(polyfit(TETA(index(1):index(2)-1), thseal1(index(1):index(2)-1), 1), t).*(t >= TETA(index(1)) & t < TETA(index(2))) +...
           polyval(polyfit(TETA(index(2):end), thseal1(index(2):end), 1), t).*(t >= TETA(index(2)));
th2 = @(t) polyval(polyfit(TETA(1:index(1)-1), thseal2(1:index(1)-1), 1), t).*(t < TETA(index(1))) +...
           polyval(polyfit(TETA(index(1):index(2)-1), thseal2(index(1):index(2)-1), 1), t).*(t >= TETA(index(1)) & t < TETA(index(2))) +...
           polyval(polyfit(TETA(index(2):end), thseal2(index(2):end), 1), t).*(t >= TETA(index(2)));

TH = TR(end)-TR(1); %tooth height [mm]
b = w;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINITION OF LEAKAGES BASED ON GEOMETRY (TIP, SIDE, DRAINAGE, BACKEND GROOVE)
Qtip = @(p1, p2, ht) b*ht.^3/(12*mu*TW_t)*(p1-p2) + .5*(b*ht*omega*TR(end));
                                                   %[mm3/s] !!!!!!!!!!!!!!
Qlat = @(p1, p2)  TH*LC.^3/(12*mu*MT)*(p1-p2) + .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] !!!!!!!!!!!!!!
                                                   
Qgroove = @(p, A) 0.65*A.*(1e3*sqrt(2*abs(pm - p)./density(pm,BM,pcav,rhocav))).*(pm>=p) + ...
                 (-0.65*A.*(1e3*sqrt(2*abs(pm - p)./density(pm,BM,pcav,rhocav)))).*(pm<p);
                                                   %[mm3/s] !!!!!!!!!!!!!!
DVDT = @(h1, h2)  TR(end)*b*(h1-h2);                                               
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

htip = tipGap(N, TC, BC,  th1(t), th2(t)); %radial clearances distribution
                                           %(if the numeb of Sealed Volumes
                                           %is N then htip is a length N+1
                                           %vector)
Agr = groove_leak_area(N, th1(t), th2(t), GArea, Gangle(1), Gangle(2));  
                                           %Length N vector containing
                                           %connection areas of each sealed
                                           %volume to the backend groove
                                           %(most actually have no
                                           %connection, so a lot of vector
                                           %components are zero)
dpdt = zeros(N, 1);

%ODE system setting; ordered from the first sealed volume (i.e. suction
%side) to the last sealed volume (i.e. delivery side)
for ii = 1 : N
  if N ~= 1
    if ii == 1   %first sealed volume (directly connected to suction volume 
                 %and second sealed volume)
        dpdt(ii) = pi/180 * BM/VV*( (Qtip(p(ii+1), p(ii), htip(ii+1)) -...
            Qtip(p(ii), pa, htip(ii)) + 2*Qlat(p(ii+1), p(ii)) -...
            2*Qlat(p(ii), pa) - 2*Qdrain(p(ii)) + 2*Qgroove(p(ii), Agr(ii)))/omega - DVDT(htip(ii+1), htip(ii)) );
    elseif ii == N   %last sealed volume (directly connected to second-last
                     %sealed volume and delivery volume
        dpdt(ii) = pi/180 * BM/VV*( (Qtip(pm, p(ii), htip(ii+1)) -...
            Qtip(p(ii), p(ii-1), htip(ii)) + 2*Qlat(pm, p(ii)) -...
            2*Qlat(p(ii), p(ii-1)) - 2*Qdrain(p(ii)) + 2*Qgroove(p(ii), Agr(ii)))/omega - DVDT(htip(ii+1), htip(ii)) );
    else            %generic in-between sealed volume
        dpdt(ii) = pi/180 * BM/VV*( (Qtip(p(ii+1), p(ii), htip(ii+1)) -...
            Qtip(p(ii), p(ii-1), htip(ii)) + 2*Qlat(p(ii+1), p(ii)) -...
            2*Qlat(p(ii), p(ii-1)) - 2*Qdrain(p(ii)) + 2*Qgroove(p(ii), Agr(ii)))/omega - DVDT(htip(ii+1), htip(ii)) );
    end

  elseif N == 1     %if there is only one contemporary sealed volume it is 
                    %directly connected to suction and delivery volumes
       dpdt(ii) = pi/180 * BM/VV*( (Qtip(pm, p(ii), htip(ii+1)) -...
            Qtip(p(ii), pa, htip(ii)) + 2*Qlat(pm, p(ii)) -...
            2*Qlat(p(ii), pa) - 2*Qdrain(p(ii)))/omega - DVDT(htip(ii+1), htip(ii)) );
  end
end

end
