% pc_co2.m

% Kalles fra pc.m
% Returnerer TVNderivert, tar inn: tilstander, og pl,pv,wl,wv og Htot fra
% forrige delvolum, samt Q.

function dy=pc_co2(y,u)
global b

pl_neste = u(1);
pv_neste = u(2);
nl_inn = u(3)/b.MW_pc(1);
nv_inn = u(4)/b.MW_pc(1);       %[mol/s]
H_inn = u(5);                   %[J/s] Totalt for l og v 
Q = u(6);                       %[J/s] Beregnet i pc.m. For hvert kontrollvolum

Tl = y(1);                      
Vl = y(2);
Nl = y(3);
Tv = y(4);
Vv = b.V_pc(1)-Vl;
Nv = y(5);

if Tl<0
    disp('pc_co2.m, Tl=')
    disp(Tl)
end

Al = co2_sw([Tl,Vl,Nl]');
Av = co2_sw([Tv,Vv,Nv]');

%Elementer til M-matrisen.
dUdTl = -Tl*Al.H(1,1); % (= Cvliq)
dUdVl = Al.g(2) - Tl*Al.H(1,2);
dUdNl = Al.g(3) - Tl*Al.H(1,3);
dUdTv = -Tv*Av.H(1,1); % (= Cvvap)
dUdVv = Av.g(2) - Tv*Av.H(1,2);
dUdNv = Av.g(3) - Tv*Av.H(1,3);
    
% Elementer til x: (Høyre side av likn)
Hl = (-Al.g(1)*Tl + Al.g(3)*Nl)/Nl;         % H=TS+myN
Hv = (-Av.g(1)*Tv + Av.g(3)*Nv)/Nv;         % [J/mol]
pl = -Al.g(2);                              % [Pa]
pv = -Av.g(2); % trykk i hvert kontr.volum
myl= Al.g(3);                               % [J/mol]
myv= Av.g(3);

Dpl = pl - pl_neste;
Dpv = pv - pv_neste;

I = find(Dpl<0);
Dpl(I) = 0;
I = find(Dpv<0);
Dpv(I) = 0;

rhol = Nl/Vl*b.MW_pc(1);% [kg/m3]
rhov = Nv/Vv*b.MW_pc(1);
rhotot = (Nl+Nv)/b.V_pc(1)*b.MW_pc(1);

wl = b.Cv_pc(1)*Vl/b.V_pc(1)*rhol*sqrt(abs(Dpl)/rhotot)*sign(Dpl);%[kgliq/s]
wv = b.Cv_pc(1)*Vv/b.V_pc(1)*rhov*sqrt(abs(Dpv)/rhotot)*sign(Dpv);%[kgvap/s]

nl_ut = wl/b.MW_pc(1); % [molliq/s]
nv_ut = wv/b.MW_pc(1); % [molvap/s]

H_ut = Hl*nl_ut + Hv*nv_ut;                 % [J/s]

Udot = H_inn - H_ut + Q;                    % [J/s]
Vdot = 0;
Ndot = (nl_inn+nv_inn) - (nl_ut+nv_ut);     % [mol/s]
deltaT = 0;
deltap = pl-pv;                             % [Pa]
deltamy= myl-myv;                           % [J/mol]


x = [Udot; Vdot; Ndot; -deltaT; -deltap; -deltamy];

M = [dUdTl  dUdVl  dUdNl   dUdTv  dUdVv   dUdNv;    %Energibalanse
      0       1      0       0      1      0;       %Volumbal
      0       0      1       0      0      1;       %Massebal
      1       0      0      -1      0      0;       %Templikevekt
     -Al.H(2,1:3)  Av.H(2,1:3);                     %Trykklikev
     Al.H(3,1:3)  -Av.H(3,1:3)];                    %Potensiallikev

 
M = sparse(M);

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

B=NaN*zeros(size(M));
for i=1:length(M)
    B(i,:) = b.D1(i)*M(i,:).*b.D2';
end
B = sparse(B);
q = b.D1.*x;

% z er egentlig (skalert) dydt, dvs ydot
z = inv(B)*q; % T,V,N
dy = z.*b.D2;
dy = [dy(1:4); dy(6)];   % Ønsker ikke å returnere Vv

g=isreal(Udot); if g==0, disp('kompl elem i U'), end
g=isreal(Vdot); if g==0, disp('kompl elem i V'), end
g=isreal(Ndot); if g==0, disp('kompl elem i N'), end
g=isreal(deltaT); if g==0, disp('kompl elem i T'), end
g=isreal(deltap); if g==0, disp('kompl elem i p'), end
g=isreal(deltamy); if g==0, disp('kompl elem i my'), end