% Script for calculating the transferred heat and outlet temperatures in a
% heat exchanger using the NTU method.

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.00001 && beta >= 0.99999
        
        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.00001 || beta < 0.99999
    
        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