function [Ceq, dCeqdx] = HExSimpH(x,par,scl,con,~,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 hot stream.

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

%       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 [-]

%       x(k4+1):    Transferred energy 

% 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+2)';

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

Tref    = x(k1+1)+273.15;

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

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

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

Qtrans    = x(k4+1+1);

% Mass balance
Ceq = zeros(k3+1,1);
Ceq(1) = (hot.Qin - hot.Qout);

% Pressure balance
Ceq(2) = (par.kHot^2*(hot.pin^2-hot.pout^2)-hot.Qin^2);
% Ceq(2) = par.kHot*sqrt(host.pin^2-hot.pout^2)-hot.Qin;

% Temperature definition
Ceq(3) = (Tref - hot.Tout);

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

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

% Transferred energy mass balance
mHInd = [2 3 5:k1];
mH = (par.Hin.eval(hot.in(mHInd))-par.Hout.eval(hot.out(mHInd)));
Ceq(k3+1) = Qtrans-hot.Qin*(con.M*hot.xin)*mH;


% 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.QH = [1 1];
        C.QH = [1 k1+1+1];
        
        R.pTH = [2 2 2      2      3      3];
        C.pTH = [1 2 k1+2+1 k4+2+1 k1+3+1 k1+1];
        
        R.xH = [4:k1 4:k1 4:k1 4:k1];
        C.xH = [On1j 4:k1 On1j*(k1+1+1) k1+4+1:k2+1];
        
        R.C = k1+1:k3;
        C.C = R.C+k1+1;
        
        C.Qt = [1:k1 k1+2+1 k1+3+1 k1+5+1:k2+1 k4+1+1];
        R.Qt = (k3+1)*ones(size(C.Qt));
        
        % Definition of the values at the positions
        V.QH = [1 -1];...                       Derivative of Ceq(1) w.r.t. hot.Qin and hot.Qout
            
        V.pTH = [-2*hot.Qin,...                         Derivative of Ceq(2) w.r.t. hot.Qin 
                  2*par.kHot^2*hot.pin,...              Derivative of Ceq(2) w.r.t. hot.pin 
                 -2*par.kHot^2*hot.pout,...             Derivative of Ceq(2) w.r.t. hot.pout 
                  2*par.kHot*(hot.pin^2-hot.pout^2),...	Derivative of Ceq(2) w.r.t. kHot
                 -1,...                                 Derivative of Ceq(3) w.r.t. hot.Tout
                  1];...                                Derivative of Ceq(3) w.r.t. Tref
                 
%         V.pTH = [-1,...                                             Derivative of Ceq(2) w.r.t. hot.Qin
%                   par.kHot*hot.pin/sqrt(hot.pin^2-hot.pout^2),...   Derivative of Ceq(2) w.r.t. hot.pin
%                  -par.kHot*hot.pout/sqrt(hot.pin^2-hot.pout^2),...  Derivative of Ceq(2) w.r.t. hot.pout
%                   sqrt(hot.pin^2-hot.pout^2),...                    Derivative of Ceq(2) w.r.t. kHot
%                  -1];...                                            Derivative of Ceq(3) w.r.t. hot.Tout
                
        V.xH = [hot.xin',...                 Derivative of Ceq(4:k1) w.r.t. hot.Qin
                hot.Qin*On1j,...             Derivative of Ceq(4:k1) w.r.t. hot.xin
               -hot.xout',...                Derivative of Ceq(4:k1) w.r.t. hot.Qout
               -hot.Qout*On1j];...           Derivative of Ceq(4:k1) w.r.t. hot.xout
              
        V.C = ones(1,k2);...                 Derivative of Ceq(k1+1:k3) w.r.t. cold outlet variables
            
        % Calculation of the values for the transferred heat equation
        dHin  = par.Hin.eval_jacobian(hot.in(mHInd));
        dHout = par.Hout.eval_jacobian(hot.out(mHInd));
        dCeqk31x = hot.Qin*con.M*mH + hot.Qin*(con.M*hot.xin)*[0 dHin(3:end)];
        
        V.Qt = [-(hot.xin'*con.M')*mH,...                   Derivative of Ceq(k3+1) w.r.t. hot.Qin
                -hot.Qin*(hot.xin'*con.M')*dHin(1:2),...    Derivative of Ceq(k3+1) w.r.t. hot.pin and hot.Tin
                -dCeqk31x,...                               Derivative of Ceq(k3+1) w.r.t. hot.out
                 hot.Qin*(hot.xin'*con.M')*dHout,...        Derivative of Ceq(k3+1) w.r.t. hot.out
                 1];...                                     Derivative of Ceq(k3+1) w.r.t. Qtrans
                           
        % Definition of the combined vectors
        R.tot = [R.QH R.pTH R.xH R.C R.Qt];
        C.tot = [C.QH C.pTH C.xH C.C C.Qt];
        V.tot = [V.QH V.pTH V.xH V.C V.Qt]...
                .*scl.Eq(R.tot)'.*scl.Var(C.tot);
            
        % Definition of the sparse matrix
        dCeqdx = sparse(R.tot,C.tot,V.tot,k3+1,k4+3);
end