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


% Simple Heat Exchanger
% Written by Julian Straus 20150730

% This function is a simple Heat Exchanger, in which the output is defined
% with a constant temperature. No Heat Transfer calculations have to be
% performed. It is used, if we would like to specify the outlet temperature
% of the cold stream.

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

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

%       x(k2+1):    Inlet flow cold stream in [mol/s]
%       x(k2+2):    Inlet pressure cold stream in [bar]
%       x(k2+3):    Inlet temperature cold stream in [C]
%       x(k2+4:k3): Inlet composition cold stream in [-]

%       x(k3+1):    Outlet flow cold stream in [mol/s]
%       x(k3+2):    Outlet pressure cold stream in [bar]
%       x(k3+3):    Outlet temperature cold stream in [C]
%       x(k3+4:k4): Outlet composition cold stream in [-]

% The necessary structures are defined as
%       par:        Parameters ot the unit, must include:
%                   .UA:        Heat transfer coefficient [W/K/s]
%                   .cpHot:     Molar heat capacity hot stream[J/mol/K] 
%                   .cpCold:    Molar heat capacity cold stream[J/mol/K] 
%               	.kHot:      Pressure drop factor hot stream
%                   .kCold:     Pressure drop factor cold stream
%                   .cond:      Energy adjusting for condensation

%       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 condensation duty

%       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;
k3 = 3*itg.k;
k4 = 4*itg.k;

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

% Rearrangement of the variables
hot.Qin = x(1);
hot.pin = x(2);
hot.Tin = x(3);
hot.xin = x(4:itg.k)+273.15;

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

cold.Qin = x(k2+1);
cold.pin = x(k2+2);
cold.Tin = x(k2+3)+273.15;
cold.xin = x(k2+4:k3);

cold.Qout = x(k3+1);
cold.pout = x(k3+2);
cold.Tout = x(k3+3)+273.15;
cold.xout = x(k3+4:k4);

% Zero balance for unused Cold Variables
Ceq = zeros(k3,1);
Ceq(1:k2) = [hot.Qin;...
             hot.pin;...
             hot.Tin-273.15;...
             hot.xin;...
             hot.Qout;...
             hot.pout;...
             hot.Tout-273.15;...
             hot.xout];
            
% Mass balances
Ceq(k2+1) = (cold.Qin - cold.Qout);

% Pressure balance
Q2 = scl.Q^2;
Ceq(k2+2) = (par.kCold^2*(cold.pin^2-cold.pout^2)-cold.Qin^2);
% Ceq(k2+2) = (par.kCold*sqrt(cold.pin^2-cold.pout^2)-cold.Qin);

% Temperature definition
Ceq(k2+3) = (par.Tref+273.15 - cold.Tout);

% Component mass balance
Ceq(k2+4:k3) = (cold.Qin*cold.xin-cold.Qout*cold.xout);

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

% Gradient of the function
switch nargout
    case 2
        % Definition of auxiliary values
        On1j = ones(1,itg.j);
        
        % Definition of the indices for the rows (R) and columns (C)
        R.QC = [k2+1 k2+1];
        C.QC = [k2+1 k3+1];
        
        R.pTC = [k2+2 k2+2 k2+2 k2+2 k2+3];
        C.pTC = [k2+1 k2+2 k3+2 k4+2 k3+3];
        
        R.xC = [k2+4:k3 k2+4:k3 k2+4:k3 k2+4:k3];
        C.xC = [On1j*(k2+1) k2+4:k3 On1j*(k3+1) k3+4:k4];
        
        R.H = 1:k2;
        C.H = 1:k2;
        
        % Definition of the values at the positions
        V.QC = [1 -1];...                       Derivative of Ceq(k2+1) w.r.t. cold.Qin and cold.Qout
            
        V.pTC = [-2*cold.Qim,...                            Derivative of Ceq(k2+2) w.r.t. cold.Qin 
                  2*par.kCold^2*cold.pin,...                Derivative of Ceq(k2+2) w.r.t. cold.pin 
                 -2*par.kCold^2*cold.pout,...               Derivative of Ceq(k2+2) w.r.t. cold.pout 
                  2*par.kCold*(cold.pin^2-cold.pout^2),...	Derivative of Ceq(k2+2) w.r.t. kCold
                 -1];...                                    Derivative of Ceq(k2+3) w.r.t. cold.Tout
                 
%         V.pTC = [-1,...                                                 Derivative of Ceq(k2+2) w.r.t. cold.Qin
%                   par.kCold*cold.pin/sqrt(cold.pin^2-cold.pout^2),...   Derivative of Ceq(k2+2) w.r.t. cold.pin
%                  -par.kCold*cold.pout/sqrt(cold.pin^2-cold.pout^2),...  Derivative of Ceq(k2+2) w.r.t. cold.pount
%                   sqrt(cold.pin^2-cold.pout^2),...                      Derivative of Ceq(k2+2) w.r.t. kCold
%                  -1];...                                                Derivative of Ceq(k2+3) w.r.t. cold.Tout
                 
        V.xC = [cold.xin',...                   Derivative of Ceq(k2+4:k3) w.r.t. cold.Qin
                cold.Qin*On1j,...               Derivative of Ceq(k2+4:k3) w.r.t. cold.xin
               -cold.xout',...                  Derivative of Ceq(k2+4:k3) w.r.t. cold.Qout
               -cold.Qout*On1j];...             Derivative of Ceq(k2+4:k3) w.r.t. cold.xout
        
        V.H = ones(1,k2);...                    Derivative of Ceq(1:k2) w.r.t. hot outlet variables
           
        % Definition of the combined vectors
        R.tot = [R.QC R.pTC R.xC R.H];
        C.tot = [C.QC C.pTC C.xC C.H];
        V.tot = [V.QC V.pTC V.xC V.H]...
                .*scl.Eq(R.tot)'.*scl.Var(C.tot);
            
        % Definition of the sparse matrix
        dCeqdx = sparse(R.tot,C.tot,V.tot,k3,k4+2);
end