function [Ceq, dCeqdx] = ValveSimp(x,par,scl,~,~,itg)

% Simple Valve
% Written by Julian Straus 20160317

% This function is a simple valve, which calculates the pressure difference
% based on a fixed k value, or varies the k-value based on a fixed pressure
% drop.

% The state variables are defined as
%       x(1):       Inlet flow in [mol/s]
%       x(2):       Inlet pressure in [bar]
%       x(3):       Inlet temperature in [C]
%       x(4:k1):    Inlet composition in [-]

%       x(k1+1):    Outlet flow in [mol/s]
%       x(k1+2):    Outlet pressure in [bar]
%       x(k1+3):    Outlet temperatur in [C]
%       x(k1+4:k2): Outlet composition in [-]

% The necessary structures are defined as
%       par:        Parameters ot the unit, must include:
%                   .K:        Pressure drop coefficient [mol/s/P]

%       scl:        Scaling factors of the system

%       itg:        Intiger values, which should include
%                   .k    	Total number of variables per stream
%                   .j     	Total number of components in each stream

% The output is defined
%   	Ceq:        Equality Constraint for the optimizer/solver
%       dCeqdx:     Jacobian of the equality constraints

% Definition of multiples of itg.k
k1 = itg.k;
k2 = 2*itg.k;

% Rescaling of the variables
x = x.*scl.Var(1:k2)';

% Redefinition of the state variables
Qin = x(1);
pin = x(2);
Tin = x(3)+273.15;
xin = x(4:k1);

Qout = x(k1+1);
pout = x(k1+2);
Tout = x(k1+3)+273.15;
xout = x(k1+4:k2);

% Mass balances
Ceq = zeros(8,1);
Ceq(1) = (Qin - Qout);
    
% Momentum balance
Ceq(2) = (Qin^2 - par.K^2*(pin^2-pout^2));
% Ceq(2) = (Qin - par.K*sqrt(pin^2-pout^2));

% Energy balance;
Ceq(3) = (Tin/pin - Tout/pout);

% Component mass balances
Ceq(4:k1) = (Qin*xin - Qout*xout);

Ceq = Ceq.*scl.Eq;
    
% Gradient of the function
switch nargout
    case 2
        % Definition of the indices for the rows (R), columns (C), and 
        % values of the mass balance
        R.Q = [1  1];
        C.Q = [1  k1+1];
        V.Q = [1 -1];
        
        % Definition of the indices for the rows (R) and columns (C) of the
        % momentum and energy balance
        R.pT = [2 2 2    2    3 3 3    3];
        C.pT = [1 2 k1+2 k2+1 2 3 k1+2 k1+3];
        
        % Definition of the values of the momentum and energy balance at 
        % the positions
        V.pT = [2*Qin,...                   Derivative of Ceq(1) w.r.t. Qin
               -2*par.K^2*pin,...           Derivative of Ceq(1) w.r.t. pin
                2*par.K^2*pout,...          Derivative of Ceq(1) w.r.t. pout
               -2*par.K*(pin^2-pout^2),...	Derivative of Ceq(1) w.r.t. K
               -Tin/pin^2,...               Derivative of Ceq(2) w.r.t. pin
                1/pin,...                   Derivative of Ceq(2) w.r.t. Tin
                Tout/pout^2,...             Derivative of Ceq(2) w.r.t. pout
               -1/pout];...                 Derivative of Ceq(2) w.r.t. Tout

%         V.pT = [1,...                       Derivative of Ceq(1) w.r.t. Qin
%                -par.K*pin/sqrt(pin^2-pout^2),...           Derivative of Ceq(1) w.r.t. pin
%                 par.K*pout/sqrt(pin^2-pout^2),...          Derivative of Ceq(1) w.r.t. pout
%                -sqrt(pin^2-pout^2),...      Derivative of Ceq(1) w.r.t. K
%                -Tin/pin^2,...               Derivative of Ceq(2) w.r.t. pin
%                 1/pin,...                   Derivative of Ceq(2) w.r.t. Tin
%                 Tout/pout^2,...             Derivative of Ceq(2) w.r.t. pout
%                -1/pout];...                 Derivative of Ceq(2) w.r.t. Tout
               
        % Definition of the indices for the rows (R), columns (C), and 
        % values of the component mass balance
        On15 = ones(1,5);
      	R.x = (4:k1)'*ones(1,4);
        R.x = R.x(:)';
        C.x = [On15 4:k1 On15*(k1+1) k1+4:k2];
        V.x = [xin',...         Derivative of Ceq(7:11) w.r.t. hot.Qin
               Qin*On15,...     Derivative of Ceq(7:11) w.r.t. hot.xin
              -xout',...        Derivative of Ceq(7:11) w.r.t. hot.Qout
              -Qout*On15];...   Derivative of Ceq(7:11) w.r.t. hot.xout
           
        % Definition of the combined vectors
        R.tot = [R.Q R.pT R.x];
        C.tot = [C.Q C.pT C.x];
        V.tot = [V.Q V.pT V.x]...
                .*scl.Eq(R.tot)'.*scl.Var(C.tot);
        
        % Definition of the sparse matrix
        dCeqdx = sparse(R.tot,C.tot,V.tot,k1,k2+1);
         
end