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

% Simple Compressor
% Written by Julian Straus 20160317

% This function defines a simple compressor based on polytopic efficiency
% calculations. It calculates the nonlinear equality constraints for
% optimization.

% 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):    Compressor Speed in [RPM]
%       x(k1+2):    Compressor Duty in [W]

%       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:
%                   .W          Power of compressor in [MW]
%                   .polyeff    Polytropic efficiency of the compressor
%                   .n          Polytropic factor

%       scl:        Scaling factors of the system, must include:
%                   .Q          Scaling factor for the flow
%                   .p          Scaling factor for the pressure
%                   .T          Scaling factor for the temperature
%                   .W          Scaling factor for the power

%       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 as
%       Ceq         Equality Constraint for the otpimizer/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';

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

S   = x(k1+1);
W   = x(k1+2);

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

% Calculation of the values for the polytropic efficiency and exponent 
% based on the inlet conditions
ValVec = [Qin pin Tin S];
eta = par.polyeff.eval(ValVec);
n = par.n.eval(ValVec);

% Predefinition of the vector for the equality constraints and n/(n-1)
Ceq = zeros(9,1);

% Definition of auxiliary values
ratio = (n-1)/n;
Fract = (eta*W)/(con.R*Tin*Qin);
FractRatio = Fract*ratio;

% Mass balance
Ceq(1) = Qin - Qout;

% Pressure calculation
Ceq(2) = pout - pin*(FractRatio+1)^(1/ratio);

% Temperature calculation
Ceq3aux = Tin*FractRatio;
Ceq(3) = Tout-Ceq3aux-Tin;

% Calculation of the values for the power based on the inlet conditions
Ceq(4) =  W - par.W.eval(ValVec);

% Component mass balance
Ceq(5:end) = Qin*xin - Qout*xout;

% Scaling of the Equality Constraints
Ceq = Ceq.*scl.Eq;

% Jacobian of the functions
switch nargout
    case 2
        % Definition of auxilliary values
        On1j = ones(1,itg.j);
        
        % Definition of the indices for the rows (R) and columns (C)
        R.Q = [1 1];
        R.p = [2 2 2 2 2 2];
        R.T = [3 3 3 3 3 3];
        R.x = (5:k1+1)'*ones(1,4);
        R.W = [4 4 4 4 4];
        R.x = R.x(:)';
        
        C.Q = [1 k1+2+1];
        C.p = [1:3 k1+1 k1+2+2 k1+2];
        C.T = [1:3 k1+1 k1+2+3 k1+2];
        C.W = [1 2 3 k1+1 k1+2];
        C.x = [On1j On1j*(k1+2+1) 4:k1 k1+2+4:k2+2];
        
        % Evaluation of the Jacobian of the compressor parameters w.r.t. 
        % the inlet conditions.
        deta = par.polyeff.eval_jacobian(ValVec);
        dn = par.n.eval_jacobian(ValVec);
        dW = par.W.eval_jacobian(ValVec);
        
        % Definition of the partial derivatives of equation 2 w.r.t. the
        % compressor parameters
        Ceq2aux1 = pin*(FractRatio+1)^(1/ratio-1)*Fract;
        dCeq2deta = -Ceq2aux1/eta;
        dCeq2dn = -pin*(FractRatio+1)^(1/ratio)*...
                 (n*(FractRatio/(n-1)-FractRatio/n)/((1+FractRatio)*(n-1))+...
                 (-n/(n-1)^2+1/(n-1))*log(1+FractRatio));
        
        % Definition of the partial derivatives of equation 2
        dCeq2dQin = Ceq2aux1/Qin              + dCeq2dn*dn(1) + dCeq2deta*deta(1);
        dCeq2dpin = -(FractRatio+1)^(1/ratio) + dCeq2dn*dn(2) + dCeq2deta*deta(2);
        dCeq2dTin = Ceq2aux1/Tin              + dCeq2dn*dn(3) + dCeq2deta*deta(3);
        dCeq2dS   =                           + dCeq2dn*dn(4) + dCeq2deta*deta(4);
        
        % Definition of the partial derivatives of equation 3 w.r.t the
        % compressor parameters
        dCeq3deta = -Ceq3aux/eta;
        dCeq3dn = -Tin*Fract/n^2;
        
        % Definition of the partial derivatives of equation 3
        dCeq3dQin = Ceq3aux/Qin + dCeq3dn*dn(1) + dCeq3deta*deta(1);
        dCeq3dpin =               dCeq3dn*dn(2) + dCeq3deta*deta(2);
        dCeq3dTin = -1 +          dCeq3dn*dn(3) + dCeq3deta*deta(3);
        dCeq3dS   =               dCeq3dn*dn(4) + dCeq3deta*deta(4);
        
        % Definition of the values
        V.Q = [1 -1];...                    Derivative of Ceq(1) w.r.t. Qin and Qout
              
        V.p = [dCeq2dQin,...                Derivative of Ceq(2) w.r.t. Qin
               dCeq2dpin,...                Derivative of Ceq(2) w.r.t. pin
               dCeq2dTin,...                Derivative of Ceq(2) w.r.t. Tin
               dCeq2dS,...                  Derivative of Ceq(2) w.r.t. S
               1,...                        Derivative of Ceq(2) w.r.t. pout
              -Ceq2aux1/W];...              Derivative of Ceq(2) w.r.t. W
           
        V.T = [dCeq3dQin,...                Derivative of Ceq(3) w.r.t. Qin
               dCeq3dpin,...                Derivative of Ceq(3) w.r.t. pin
               dCeq3dTin,...                Derivative of Ceq(3) w.r.t. Tin
               dCeq3dS,...                  Derivative of Ceq(3) w.r.t. S
               1,...                        Derivative of Ceq(3) w.r.t. Tout
              -Ceq3aux/W];...               Derivative of Ceq(3) w.r.t. W
              
        V.W = [-dW,...                      Derivative of Ceq(4) w.r.t. Qin, pin, Tin, and Sin
               1];...                       Derivative of Ceq(4) w.r.t. W
                
       	V.x = [xin',...                     Derivative of Ceq(5:9) w.r.t. Qin
              -xout'...                     Derivative of Ceq(5:9) w.r.t. Qout
               Qin*On1j,...                 Derivative of Ceq(5:9) w.r.t. xin
              -Qout*On1j];...               Derivative of Ceq(5:9) w.r.t. xout
       
        % Definition of the combined vectors
        R.tot = [R.Q R.p R.T R.W R.x];
        C.tot = [C.Q C.p C.T C.W C.x];
        V.tot = [V.Q V.p V.T V.W V.x]...
                .*scl.Eq(R.tot)'.*scl.Var(C.tot);
         
        % Definition of the sparse matrix
        dCeqdx = sparse(R.tot,C.tot,V.tot,k1+1,k2+2);
end