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

% Simple Heat Exchanger
% Written by Julian Straus 20150730

% This function is a simple Heat Exchanger, in which the output temperature
% is calculated using the effectiveness - NTU method.  

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

% 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)+273.15;
hot.xin = x(4:itg.k);

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

par.cond = par.cond*1e6;%scl.W;

% Calculation of the heat capacity ratio
minQc = min(cold.Qin*par.cpCold,hot.Qin*par.cpHot);
maxQc = max(cold.Qin*par.cpCold,hot.Qin*par.cpHot);
Cstar = minQc/maxQc;

% Calculation of the number of transfer units
NTU = par.UA/minQc;

% Calculation of the effectiveness
expv = exp(-NTU*(1-Cstar));
num = 1-expv;
den = 1/minQc-1/maxQc*expv;
E = num/den;

% Calculation of the energy flowrate
Q = E*(hot.Tin-cold.Tin);   % = k*u/v

% Mass balances
Ceq = zeros(k2,1);
Ceq([1 k1+1]) = [hot.Qin - hot.Qout...
                 cold.Qin - cold.Qout];

% Pressure balances
Ceq([2 k1+2]) = [par.kHot^2*(hot.pin^2-hot.pout^2)-hot.Qin^2;...
                 par.kCold^2*(cold.pin^2-cold.pout^2)-cold.Qin^2];
% Ceq([2 k1+2]) = [par.kHot*sqrt(hot.pin^2-hot.pout^2)-hot.Qin;...
%                  par.kCold*sqrt(cold.pin^2-cold.pout^2)-cold.Qin];

% Temperature balances
Ceq([3 k1+3]) = [Q-hot.Qin*par.cpHot*(hot.Tin-hot.Tout)-par.cond;...
                 Q-cold.Qin*par.cpCold*(cold.Tout-cold.Tin)+par.cond];
   
% Component mass balances
Ceq([4:k1 k1+4:k2]) = [hot.Qin*hot.xin   - hot.Qout*hot.xout...
                       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
        On14 = ones(1,4);
        On1j = ones(1,itg.j);
        
        % Definition of values used for the Jacobian of the mass balance
        R.Q = [1  1    k1+1  k1+1];
        C.Q = [1  k1+1 k2+1  k3+1];
        V.Q = [1 -1    1    -1];
        
        % Definition of values used for the Jacobian of the energy balance
        R.T = [3 3    3    k1+3 k1+3 k1+3 3 3    k1+3 k1+3];
        C.T = [3 k1+3 k2+3 3    k2+3 k3+3 1 k2+1 1    k2+1]; 
        V.T = zeros(1,10);
        V.T(1:6) = [E-hot.Qin*par.cpHot,...             Derivative of Ceq(3) w.r.t. hot.Tin
                    hot.Qin*par.cpHot,...               Derivative of Ceq(3) w.r.t. hot.Tout
                   -E,...                               Derivative of Ceq(3) w.r.t. cold.Tin
                    E,...                           	Derivative of Ceq(k1+3) w.r.t. hot.Tin
                   -E+cold.Qin*par.cpCold,...      		Derivative of Ceq(k1+3) w.r.t. cold.Tin
                   -cold.Qin*par.cpCold];...            Derivative of Ceq(k1+3) w.r.t. cold.Tout

        if cold.Qin*par.cpCold > hot.Qin*par.cpHot
            dexpmax = -par.UA/(maxQc*cold.Qin)*expv;
            dexpmin =  par.UA/(minQc*hot.Qin)*expv;
            Qpmax = (hot.Tin-cold.Tin)*...                                      k
                    (-dexpmax*den - ...                                         u'*v
                    (-(1/maxQc*dexpmax+(-1/(maxQc*cold.Qin))*expv)*num))/...    v'*u
                    den^2;                                                    % v^2
            Qpmin = (hot.Tin-cold.Tin)*...                                      k
                    (-dexpmin*den - ...                                         u'*v
                    (-1/(minQc*hot.Qin)-1/maxQc*dexpmin)*num)/...               v'*u
                    den^2;
            V.T(7:10) = [Qpmin - par.cpHot*(hot.Tin-hot.Tout),...           Derivative of Ceq(3) w.r.t. hot.Qin
                         Qpmax,...                                          Derivative of Ceq(3) w.r.t. cold.Qin
                         Qpmin,...                                          Derivative of Ceq(k1+3) w.r.t. hot.Qin
                         Qpmax - par.cpCold*(cold.Tout-cold.Tin)];...       Derivative of Ceq(k1+3) w.r.t. cold.Qin
        else
            dexpmax = -par.UA/(maxQc*hot.Qin)*expv;
            dexpmin =  par.UA/(minQc*cold.Qin)*expv;
            Qpmax = (hot.Tin-cold.Tin)*...                                      k
                    (-dexpmax*den - ...                                         u'*v
                    (-(1/maxQc*dexpmax+(-1/(maxQc*hot.Qin))*expv)*num))/...     v'*u
                    den^2;                                                    % v^2
            Qpmin = (hot.Tin-cold.Tin)*...                                      k
                    (-dexpmin*den - ...                                         u'*v
                    (-1/(minQc*cold.Qin)-1/maxQc*dexpmin)*num)/...              v'*u
                    den^2;
            V.T(7:10) = [Qpmax - par.cpHot*(hot.Tin-hot.Tout),...           Derivative of Ceq(3) w.r.t. hot.Qin
                         Qpmin,...                                          Derivative of Ceq(3) w.r.t. cold.Qin
                         Qpmax,...                                          Derivative of Ceq(k1+3) w.r.t. hot.Qin
                         Qpmin - par.cpCold*(cold.Tout-cold.Tin)];...       Derivative of Ceq(k1+3) w.r.t. cold.Qin
        end
        
        % Definition of values used for the Jacobian of the pressure balance
        R.p = [2*On14  (k1+2)*On14];
        C.p = [1 2 k1+2 k4+1 k2+1 k2+2 k3+2 k4+2];
        V.p = [-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
               -2*cold.Qin,...                                      Derivative of Ceq(k1+2) w.r.t. cold.Qin
                2*par.kCold^2*cold.pin,...                          Derivative of Ceq(k1+2) w.r.t. cold.pin
               -2*par.kCold^2*cold.pout,...                         Derivative of Ceq(k1+2) w.r.t. cold.pount
                2*par.kCold*(cold.pin^2-cold.pout^2)];...           Derivative of Ceq(k1+2) w.r.t. kCold
%         V.p = [-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(k1+2) w.r.t. cold.Qin
%                 par.kCold*cold.pin/sqrt(cold.pin^2-cold.pout^2),... Derivative of Ceq(k1+2) w.r.t. cold.pin
%                -par.kCold*cold.pout/sqrt(cold.pin^2-cold.pout^2),...Derivative of Ceq(k1+2) w.r.t. cold.pount
%                 sqrt(cold.pin^2-cold.pout^2)];...                   Derivative of Ceq(k1+2) w.r.t. kCold
        
        % Definition of values used for the Jacobian of the component 
        % mass balance
        R.x = [(4:k1)'*On14 (k1+4:k2)'*On14];
        R.x = R.x(:)';
        C.x = [On1j 4:k1 On1j*(k1+1) k1+4:k2 On1j*(k2+1) k2+4:k3 On1j*(k3+1) k3+4:k4];
        V.x = [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
               cold.xin',...                Derivative of Ceq(k1+4:k2) w.r.t. cold.Qin
               cold.Qin*On1j,...            Derivative of Ceq(k1+4:k2) w.r.t. cold.xin
              -cold.xout',...               Derivative of Ceq(k1+4:k2) w.r.t. cold.Qout
              -cold.Qout*On1j];...          Derivative of Ceq(k1+4:k2) w.r.t. cold.xout
       
        % Definition of the combined vectors
        R.tot = [R.Q R.p R.T R.x];
        C.tot = [C.Q C.p C.T C.x];
        V.tot = [V.Q V.p V.T V.x]...
                .*scl.Eq(R.tot)'.*scl.Var(C.tot);
            
       	% Definition of the sparse Jacobian matrix
        dCeqdx = sparse(R.tot,C.tot,V.tot,k2,k4+2);
end