% s_ngventil.m
function [sys,x0,str,ts] = s_pc(t,x,u,flag,ni)
% u = pl pv wng Hng fra neste, png wl wv Htot fra forrige, 

global b
switch flag

  % INITALIZATION
  % -------------
  case 0,      
        [sys,x0,str,ts]=Initialize(ni);
    
  % DERIVATIVES
  % -----------
  case 1,
    dy = pc(t,x,u);   % x er tilstandene TVN og p (y el dydt i pc.m)
    sys= dy;
    
  % OUTPUTS 
  % -------
  case 3,
    Tl = x(1);              % Kolonnevektorer (nedover)
    Vl = x(2);
    Nl = x(3);
    Tv = x(4);
    Vv = b.V_pc(1)-Vl;
    Nv = x(5);
    Tng= x(6);
    Vng= x(7);
    Nng= x(8);
    Tvegg= x(9);
    
    xv = Nv/(Nv+Nl);
    
    Al = co2_sw([Tl,Vl,Nl]');
    Av = co2_sw([Tv,Vv,Nv]');
    Ang = ng_pr([Tng,Vng,Nng,Nng*1e-6,Nng*1e-6]');
    pl  = -Al.g(2); % trykk i hvert kontr.volum
    pv  = -Av.g(2);
    png = -Ang.g(2);
    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]
    Hng = (-Ang.g(1)*Tng + Ang.g(3)*Nng)/Nng; % [J/mol]

    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);
    rhong = Nng/Vng*b.MW_pc(2);
    Dpl = pl-u(1);
    Dpv = pv-u(2);
    Dpng= png-u(5);
    
    I = find(Dpl<0);
    Dpl(I) = 0;
    I = find(Dpv<0);
    Dpv(I) = 0;
    I = find(Dpng<0);
    Dpng(I) = 0;
    
    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]
    wng = b.Cv_pc(2)*rhong*sqrt(abs(Dpng)/rhong)*sign(Dpng);%[kgNG/s]
    nl = wl/b.MW_pc(1); % [molliq/s]
    nv = wv/b.MW_pc(1); % [molvap/s]
    nng= wng/b.MW_pc(2);
    
    Htot = Hl*nl + Hv*nv;             % [J/s]
    Hng  = Hng*nng;
    
    sys = [pl pv wng Hng png wl wv Htot xv Tl Vl Nl Tv Nv Tng Vng Nng Tvegg]';
    
  case {2,4,9}
    sys=[];
    
  % UNEXPECTED FLAGS
  % ----------------
  otherwise
    error(['Unhandled flag = ',num2str(flag)]);
    
end


function [sys,x0,str,ts]=Initialize(ni)
global b 

sizes = simsizes;

sizes.NumContStates  = 9; % TVN(liq) TN(vap) TVN(ng) T(vegg)
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 18; % pl pv wng Hng png wl wv Htot T(liq),vl,T(ng)
sizes.NumInputs      = 8; % pl pv wng Hng png wl wv Htot 
sizes.DirFeedthrough = 1; % 0 hvis utgang ikke direkte funksjon av inngang, ellers 1
sizes.NumSampleTimes = 1; % at least one sample time is needed

sys = simsizes(sizes);

load initialB
% initialize the initial conditions
x0 = initB(ni).y0; % TVN(liq) TN(vap) TVN(ng) T(vegg)

% str is always an empty matrix
str = [];

% initialize the array of sample times
ts  = [0 0];
%data=initialize;
% end mdlInitializeSizes

