% pc_ng.m

% Kalles fra pc.m
% Returnerer TVNderivert, tar inn: tilstander, og w og H fra forrgie
% delvolum, samt Q.

function dy=pc_ng(y,u)
global b

n_inn = u(1)/b.MW_pc(2);        %[mol/s]
H_inn = u(2);                   %[J/s]
p_neste = u(3);
Q = u(4);                       %[J/s]

T = y(1);
V = y(2);
N = y(3);

A = ng_pr([T,V,N,N*1e-6,N*1e-6]');
    
%Elementer til M-matrisen.
dUdT = -T*A.H(1,1); % (= Cv)
dUdV = A.g(2) - T*A.H(1,2);
dUdN = A.g(3) - T*A.H(1,3);
    
% Elementer til x: (Høyre side av likn)
H = (-A.g(1)*T + A.g(3)*N)/N; % H=TS+myN [J/mol]
p = -A.g(2); % trykk i hvert kontr.volum
my= A.g(3);                                     % [J/mol]


Dp = p - p_neste;
I = find(Dp<0);
Dp(I) = 0;

rho = N/V*b.MW_pc(2);% [kg/m3]

w = b.Cv_pc(2)*rho*sqrt(abs(Dp)/rho)*sign(Dp);%[m3/s]
n_ut = w/b.MW_pc(2);       % [mol/s]

H_ut = H*n_ut;             % [J/s]

Udot = H_inn - H_ut + Q;
Vdot = 0;
Ndot = n_inn - n_ut;

x = [Udot; Vdot; Ndot];

M = [dUdT  dUdV  dUdN;  %Energibalanse
      0      1     0 ;    %Volumbal
      0      0     1];     %Massebal
M = sparse(M);

l=isreal(M);
if l==0
   disp('Komplekse elementer i M_ng')
end

% z er egentlig dydt, dvs ydot
dy = inv(M)*x; % T,V,N


