function dpdt = odesealfcn(t, p, N, index)
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 from .xlsx FILE 
filename = 'Gprofile.xlsx'; 
M = readmatrix(filename,'Sheet', 1, 'Range','A4:G4');
TW_r = M(2); TW_t = M(3); VV = M(5); MR = M(6); MT = M(7);
M = readmatrix(filename,'Sheet', 2);
TR = M(:, 1); 
filename = 'GPgeometry.xlsx'; 
M =  readmatrix(filename,'Sheet', 1, 'Range', 'A2');
TETA = M(:, 1); thseal1 = M(:, 11); thseal2 = M(:, 12);

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);
b = w;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINITION OF LEAKAGES BASED ON GEOMETRY (TIP, SIDE, DRAINAGE)
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)); %tip coinvolti nei tip leakages
Agr = groove_leak_area(N, th1(t), th2(t), GArea, Gangle(1), Gangle(2)); %vettore di dimensione N

dpdt = zeros(N, 1);

for ii = 1 : N
  if N ~= 1
    if ii == 1 
        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
        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 
        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
       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
