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

% Simple Valve
% Written by Julian Straus 20160317

% This function is a pressure drop equation

% The state variables are defined as
%       x(1):       Inlet flow in [kmol/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 [kmol/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

%       con:        Natural constants, must include:
%                   .R      Ideal gas constant in [J/k/mol]

%       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);

% Definition of auxiliary variables
f = 0.003;
M = par.M'*xin;
K = par.D*(pi/4*par.D^2)^2/(4*M*f*par.L*con.R*Tin);

% Mass balances
Ceq = zeros(k1,1);
Ceq(1) = (Qin - Qout);
    
% Momentum balance
Ceq(2) = (1/K*Qin^2 - (pin^2-pout^2));

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

% 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 auxiliary values for the derivative of K
        dKdTin = -K/Tin;
        dKdxin = -K/M*par.M';
        On1j = ones(1,itg.j);
        
        % Definition of the indices for the rows (R) and columns (C) of the
        % momentum and energy balance
        R.pT = [2 2 2    2 On1j*2 3 3];
        C.pT = [1 2 k1+2 3 4:k1   3 k1+3];
        
        % Definition of the values of the momentum and energy balance at 
        % the positions
        V.pT = [2*Qin,...                   Derivative of Ceq(2) w.r.t. Qin
               -2*K*pin,...                 Derivative of Ceq(2) w.r.t. pin
                2*K*pout,...                Derivative of Ceq(2) w.r.t. pout
                dKdTin*(pin^2-pout^2),...   Derivative of Ceq(2) w.r.t. Tin
                dKdxin*(pin^2-pout^2),...	Derivative of Ceq(2) w.r.t. xin
                1,...                       Derivative of Ceq(3) w.r.t. Tin
               -1];...                      Derivative of Ceq(3) w.r.t. Tout

        V.pT = [2*Qin/K,...                 Derivative of Ceq(2) w.r.t. Qin
               -2*pin,...                   Derivative of Ceq(2) w.r.t. pin
                2*pout,...                  Derivative of Ceq(2) w.r.t. pout
                Qin^2/K/Tin,...             Derivative of Ceq(2) w.r.t. Tin
               -Qin^2/K*M./par.M',...       Derivative of Ceq(2) w.r.t. xin
                1,...                       Derivative of Ceq(3) w.r.t. Tin
               -1];...                      Derivative of Ceq(3) w.r.t. Tout
               
        % Definition of the indices for the rows (R), columns (C), and 
        % values of the component mass balance
      	R.x = (4:k1)'*ones(1,4);
        R.x = R.x(:)';
        C.x = [On1j 4:k1 On1j*(k1+1) k1+4:k2];
        V.x = [xin',...         Derivative of Ceq(7:11) w.r.t. hot.Qin
               Qin*On1j,...     Derivative of Ceq(7:11) w.r.t. hot.xin
              -xout',...        Derivative of Ceq(7:11) w.r.t. hot.Qout
              -Qout*On1j];...   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);
         
end