function sys=sepmod(x,u,mode,ID,L,rho_L,RTmg,orientation,FB)

persistent hvec VLvec

mL=x(1);
mG=x(2);
if ~FB
    qLin=u(1);
    qGin=u(2);
    qLout=u(3);
    qGout=u(4);
elseif FB
    mLin=u(1);
    mGin=u(2);
    mLout=u(3);
    mGout=u(4);
end
A=pi*ID^2/4;
V_T=A*L;
V_L=mL/rho_L;
if orientation
    if isempty(hvec)
        hvec=0:ID/100:ID;
        phi=acos((ID-2*hvec)/ID);
        AL=A*phi/pi-ID^2*cos(phi).*sin(phi)/4;
        VLvec=AL*L;
    end
    h=interp1(VLvec,hvec,V_L);
else
    h=V_L/A;
end
V_G=V_T-V_L;
rho_G=mG/V_G;
P=mG*RTmg/V_G;
if isequal(mode,'derivatives')
    if ~FB
        sys=[(qLin-qLout)*rho_L;(qGin-qGout)*rho_G];
    elseif FB
        sys=[mLin-mLout;mGin-mGout];
    end
        
elseif isequal(mode,'measurements')
    sys=[h;P;rho_G];
end

