% pc_ng.m

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

function dy=pc_ng(y,u)
global b

w_inn = u(1);      %[kg/m3]
H_inn = u(2);       %[J/s] Totalt for l og v (beregnet i tank.m)
Q = u(3:end);      %[J/s] Beregnet i pc.m. For hvert kontrollvolum

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

Dp = p - [b.p_pc_ut(2); p(1:end-1)];
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);%[kg/s]
w0ng = w;
save massestrom w0ng

n = [w; w_inn]/b.MW_pc(2);       % [mol/s]

n_inn = n(2:end);
n_ut  = n(1:end-1);

Hut = H.*n_ut;             % [J/s]
Hinn = [Hut(2:end); H_inn];% Hinn(1) gitt av u øverst

H0ng = Hut;
save entalpi H0ng;

% Q = Qhatt.*n_ut;
Udot = Hinn - Hut + Q;
Vdot = zeros(b.n,1);
Ndot = n_inn - n_ut;

one = eye(b.n); %diagonalmatrise med enere
zero = zeros(b.n);

x = [Udot; Vdot; Ndot];

M = [diag(dUdT)  diag(dUdV)  diag(dUdN);  %Energibalanse
     zero           one         zero   ;    %Volumbal
     zero           zero        one    ];     %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


