% pc_co2.m

% Kalles fra pc.m
% Returnerer TVNderivert, tar inn: tilstander, og wl,wv og Htot fra ventil
% samt Q.

function dy=pc_co2(y,u)
global b

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

Tl = y(1:b.n);          % Kolonnevektorer 
Vl = y(b.n+1:2*b.n);
Nl = y(2*b.n+1:3*b.n);
Tv = y(3*b.n+1:4*b.n);
Vv = b.V_pc(1)-Vl;
Nv = y(4*b.n+1:5*b.n);

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

o = NaN*zeros(b.n,1);
dUdTl=o; dUdVl=o; dUdNl=o; dUdTv=o; dUdVv=o; dUdNv=o;
dpdTl=o; dpdVl=o; dpdNl=o; dpdTv=o; dpdVv=o; dpdNv=o;
dmydTl=o; dmydVl=o; dmydNl=o; dmydTv=o; dmydVv=o; dmydNv=o;
Hl=o; Hv=o; pl=o; pv=o; myl=o; myv=o;
for i=1:b.n
    Al = co2_sw([Tl(i),Vl(i),Nl(i)]'); %eller x?
    Av = co2_sw([Tv(i),Vv(i),Nv(i)]');
    
    %Elementer til M-matrisen. (Blir radvektorer)
    dUdTl(i) = -Tl(i)*Al.H(1,1); % (= Cvliq)
    dUdVl(i) = Al.g(2) - Tl(i)*Al.H(1,2);
    dUdNl(i) = Al.g(3) - Tl(i)*Al.H(1,3);
    dUdTv(i) = -Tv(i)*Av.H(1,1); % (= Cvvap)
    dUdVv(i) = Av.g(2) - Tv(i)*Av.H(1,2);
    dUdNv(i) = Av.g(3) - Tv(i)*Av.H(1,3);
    
    dpdTl(i) = -Al.H(2,1);
    dpdVl(i) = -Al.H(2,2);
    dpdNl(i) = -Al.H(2,3);
    dpdTv(i) = -Av.H(2,1);
    dpdVv(i) = -Av.H(2,2);
    dpdNv(i) = -Av.H(2,3);

    dmydTl(i)= Al.H(3,1);
    dmydVl(i)= Al.H(3,2);
    dmydNl(i)= Al.H(3,3);
    dmydTv(i)= Av.H(3,1);
    dmydVv(i)= Av.H(3,2);
    dmydNv(i)= Av.H(3,3);
    
    % Elementer til x: (Høyre side av likn)
    Hl(i) = (-Al.g(1)*Tl(i) + Al.g(3)*Nl(i))/Nl(i); % H=TS+myN
    Hv(i) = (-Av.g(1)*Tv(i) + Av.g(3)*Nv(i))/Nv(i); % [J/mol]
    pl(i) = -Al.g(2);                               % [Pa]
    pv(i) = -Av.g(2); % trykk i hvert kontr.volum
    myl(i)= Al.g(3);                                % [J/mol]
    myv(i)= Av.g(3);
end
Dpl = pl - [pl(2:end); b.p_pc_ut(1)];
Dpv = pv - [pv(2:end); b.p_pc_ut(1)];

% 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 = [wl_ventil; wl]/b.MW_pc(1); % [molliq/s]
nv = [wv_ventil; wv]/b.MW_pc(1); % [molvap/s]

nl_inn = nl(1:end-1);
nl_ut  = nl(2:end);
nv_inn = nv(1:end-1);
nv_ut  = nv(2:end);

Hut = Hl.*nl_ut + Hv.*nv_ut;                % [J/s]
Hinn = [H_inn; Hut(1:end-1)];               % H_inn(1) gitt av u øverst


Udot = Hinn - Hut + Q;                      % [J/s]
Vdot = zeros(b.n,1);
Ndot = (nl_inn+nv_inn) - (nl_ut+nv_ut);     % [mol/s]
deltaT = zeros(b.n,1);
deltap = pl-pv;                             % [Pa]
deltamy= myl-myv;                           % [J/mol]

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

x = [Udot; Vdot; Ndot; deltaT; -deltap; -deltamy];
xx=[x(1),x(21),x(41),x(61),x(81),x(101)]';

M = [diag(dUdTl)  diag(dUdVl)  diag(dUdNl)  diag(dUdTv)   diag(dUdVv)   diag(dUdNv);  %Energibalanse
     zero           one         zero         zero          one           zero;    %Volumbal
     zero           zero        one          zero          zero          one;     %Massebal
     one            zero        zero         -one          zero          zero;    %Templikevekt
     diag(dpdTl)  diag(dpdVl)  diag(dpdNl)  -diag(dpdTv)  -diag(dpdVv)  -diag(dpdNv);  %Trykklikev
     diag(dmydTl) diag(dmydVl) diag(dmydNl) -diag(dmydTv) -diag(dmydVv) -diag(dmydNv)];%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)       % bytt denne til lengde nedover (mer korrekt)
    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*b.n); dy(5*b.n+1:6*b.n)];   % Ø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
