
function sys=momentbal(t,u,x,D1,D2,h,l_o,l_g,blocks_o,blocks_g,p1_out,p3_out,ventiltype_g,ventiltype_o)
%funksjonen momentbal tar for seg beregningene av trykkforandringen i rørene etter ventilene

%Beskrivlese av parametre: p1=trykket i gassrøret	v1=hastigheten i gassrøret		
%									p3=trykket i oljerøret	v2=hastigheten i oljerøret

global M z T rho_o Cv_o Cv_g Y my1 my2 e c; 

u_o=u(1); %pådraget fra regulatoren som reg. totalnivået
u_g=u(2); %pådraget fra regulatoren som reg. trykket
p=u(3);   %trykket i separatoren
h_t=u(4); %totalt væskenivå

R=8314;   % [J/kmol*K]
g=9.81;	 % gravitasjon  [m/s2]

rho_g_t=p*M/(z*R*T); %tetthet i separatortanken

p1=x(blocks_g+1:2*blocks_g); %trykket i gassrøret
rho_g=p1*M/(z*R*T);  %Tilstandsligning. Finner tettheten i røret
v1=x(1:blocks_g)./rho_g(1:blocks_g); %hastighetene i gassrøret
v2=x(2*blocks_g+1:2*blocks_g+blocks_o); %hastighetene i oljerøret
p3=x(2*blocks_g+blocks_o+1:2*blocks_g+2*blocks_o); %trykket i oljerøret

alfa_g=asin(h/l_g); %vinkel på helling av gassrøret 
alfa_o=asin(h/l_o); %vinkel på helling av oljerøret

delta_z_g=l_g/blocks_g;  %[m] lengde på hver blokk i gassrøret
delta_z_o=l_o/blocks_o;  %[m] lengde på hver blokk i oljerøret

%lamda=0.005*(1+((2*10^4*e/Dh)+(10^6/Re))^(1/3))
%Re=D*v*rho/(my*10^(-3))


%Beregninger i gassrøret

if ventiltype_g== 1  %lineær ventilkarakteristikk
   wut_g=Y*Cv_g*u_g*sqrt((p-p1(1))*p*M/(z*R*T));  %gassmengde ut av separatoren
   
elseif ventiltype_g== 2  %equal percentage ventilkarakteristikk
   wut_g=Y*Cv_g*(40^(u_g-1))*sqrt((p-p1(1))*p*M/(z*R*T));
   
else
   wut_g=0;
   
end

rho_g_v1_inn=4*wut_g/(pi*D1^2);   %Grensebetingelse  -impuls ut av gassventilen   

Dh_g=D1;   %Hydraulisk diameter -rørets indre diameter
lamda_g=(0.005*(1+((2*10^4*e/Dh_g)+(10^6*my1*10^(-3)/(D1*rho_g.*v1))).^(1/3)))';  %Friksjonsfaktoren

%d(rho_g*v1)/dt=-d(rho_g*v1^2+p1)/delta_z-lamda*0.5*rho_g*v1^2*(1/D1)+rho_g*g_z   Momentbalansen
sys(1)=(((rho_g_v1_inn.^2/rho_g_t)+p1(1))-(rho_g(1).*v1(1).^2+p1(2)))/delta_z_g-0.5*lamda_g(1).*rho_g(1).*v1(1).^2*(1/D1)-rho_g(1)*g*sin(alfa_g);

sys(2:blocks_g-1)=((rho_g(1:blocks_g-2).*v1(1:blocks_g-2).^2+p1(2:blocks_g-1))-(rho_g(2:blocks_g-1).*v1(2:blocks_g-1).^2+p1(3:blocks_g)))/delta_z_g-0.5*lamda_g(2:blocks_g-1).*rho_g(2:blocks_g-1).*v1(2:blocks_g-1).^2*(1/D1)-rho_g(2:blocks_g-1)*g*sin(alfa_g);

sys(blocks_g)=((rho_g(blocks_g-1).*v1(blocks_g-1).^2+p1(blocks_g))-(rho_g(blocks_g).*v1(blocks_g).^2+p1_out))/delta_z_g-0.5*lamda_g(blocks_g).*rho_g(blocks_g).*v1(blocks_g).^2*(1/D1)-rho_g(blocks_g)*g*sin(alfa_g); 

%dp1/dt=-(d(rho_g*v1)/delta_z)*(z*R*T/M)  Massebalansen
sys(blocks_g+1)=(((rho_g_v1_inn)-(rho_g(1).*v1(1)))/delta_z_g)*(z*R*T/M);

sys(blocks_g+2:2*blocks_g-1)=(((rho_g(1:blocks_g-2).*v1(1:blocks_g-2))-(rho_g(2:blocks_g-1).*v1(2:blocks_g-1)))/delta_z_g)*(z*R*T/M);

sys(2*blocks_g)=(((rho_g(blocks_g-1).*v1(blocks_g-1))-(rho_g(blocks_g).*v1(blocks_g)))/delta_z_g)*(z*R*T/M);  

%Beregninger i oljerøret

if ventiltype_g== 1  %lineær ventilkarakteristikk
   qut_o=Cv_o*u_o*sqrt((p+rho_o*g*h_t-p3(1))/rho_o); %mengde olje ut av separatoren
   
elseif ventiltype_g== 2   %equal percentage ventilkarakteristikk
   
   qut_o=Cv_o*(40^(u_o-1))*sqrt((p+rho_o*g*h_t-p3(1))/rho_o);
   
else
   qut_o=0;
   
end

   
v2_inn=4*qut_o/(pi*D2^2);   %Grensebetingelse -hastighet ut av oljeventilen 

Dh_o=D2;  %hydraulisk diameter  -rørets indre diameter
lamda_o=(0.005*(1+((2*10^4*e/Dh_o)+(10^6*my2*10^(-3)/(D2*v2*rho_o))).^(1/3)))';  %Frisksjonsfaktoren

%dv2/dt=-d(v2^2)/delta_z-(1/rho_o)*(dp3/delta_z)-0.5*lamda*v2^2*(1/D2)+g_z  Momentbalansen
sys(2*blocks_g+1)=((v2_inn.^2)-(v2(1).^2))/delta_z_o+(1/rho_o)*((p3(1)-p3(2))/delta_z_o)-0.5*lamda_o(1).*v2(1).^2*(1/D2)-g*sin(alfa_o);

sys(2*blocks_g+2:2*blocks_g+blocks_o-1)=((v2(1:blocks_o-2).^2)-(v2(2:blocks_o-1).^2))/delta_z_o+(1/rho_o)*((p3(2:blocks_o-1)-p3(3:blocks_o))/delta_z_o)-0.5*lamda_o(2:blocks_o-1).*v2(2:blocks_o-1).^2*(1/D2)-g*sin(alfa_o);

sys(2*blocks_g+blocks_o)=((v2(blocks_o-1).^2)-(v2(blocks_o).^2))/delta_z_o+(1/rho_o)*((p3(blocks_o)-p3_out)/delta_z_o)-0.5*lamda_o(blocks_o).*v2(blocks_o).^2*(1/D2)-g*sin(alfa_o);

%dp3/dt=-c^2*(d(rho_o*v2))/delta_z   Massebalansen
sys(2*blocks_g+blocks_o+1)=c^2*rho_o*(v2_inn-v2(1))/delta_z_o;

sys(2*blocks_g+blocks_o+2:2*blocks_g+2*blocks_o-1)=c^2*rho_o*((v2(1:blocks_o-2))-(v2(2:blocks_o-1)))/delta_z_o;

sys(2*blocks_g+2*blocks_o)=c^2*rho_o*(v2(blocks_o-1)-v2(blocks_o))/delta_z_o;

t;



