%NTU_CALC calculates the heat duty [W] and outlet temperatures [°C] in a
%heat exchanger as a function of inlet temperatures and flow rates. The 
%syntax is
%
%   [X] = NTU_calc(u)
%
%where u is a vector containing:
%
%   u(1) - cold side inlet temperature [°C]
%   u(2) - cold side flow rate [kg/s]
%   u(3) - hot side inlet temperature [°C]
%   u(4) - hot side flow rate [kg/s]
%   u(5) - UA number [W/K]
%   u(6) - hot side medium specific heat capacity [J/kg*K]
%   u(7) - cold side medium specific heat capacity [J/kg*K]
%
%The output X contains:
%
%   X(1) - hot side outlet temperature [°C]
%   X(2) - cold side outlet temperature [°C]
%   X(3) - heat duty [W]
%
%Helge Smedsrud
%IKP, NTNU
%2008

function [X] = NTU_calc(Tic,wc,Tih,wh,UA,cph,cpc)

if wc <= 0 || wh <= 0 || Tic < 0 || Tih < 0

    X = [0 0 0];
    
else

    beta = (wc*cpc)/(wh*cph);
    
    if beta <= 1.0001 && beta >= 0.9999
        
        delta = beta - 1;
        
        ed = UA*delta/(wc*cpc);
    
        S = 1 - ed + ((-ed)^2)/2 + ((-ed)^3)/factorial(3) + ...
            ((-ed)^4)/factorial(4) + ((-ed)^5)/factorial(5) - 1 - delta;
    
        gamma = -ed*(1+S)+1;
    
        es = UA*(1+S)/(wc*cpc);
    
        D(1,1) = gamma/(1+es);
        D(1,2) = es*beta/(1+es);
        D(2,1) = es/(1+es);
        D(2,2) = 1/(1+es);
        
    elseif beta > 1.0001 || beta < 0.9999
    
        epsilon = -1; % (-1 for counter-current flow; 1 for co-current flow)
    
        alpha = UA/(wh*cph) + (epsilon*UA/(wc*cpc));
    
        gamma = exp(-alpha);
    
        D(1,1) = gamma*(1-beta)/(gamma-beta);
        D(1,2) = beta*(gamma-1)/(gamma-beta);
        D(2,1) = (gamma-1)/(gamma-beta);
        D(2,2) = (1-beta)/(gamma-beta);
    
    elseif isnan(beta) == 1
        
        D = [0 0;0 0];
        
    end

    To = D*[Tih;Tic];
    q = wh*cph*(Tih-To(1));
        
        if isnan(To(1)) == 1 || isnan(To(2)) == 1 || isnan(q) == 1
    
            To(1) = 0;
            To(2) = 0;
            q = 0;
        
        end

    X = [To(1) To(2) q];

end