r'''Returns the Number of Transfer Units of a heat exchanger at a specified 
    heat capacity rate, effectiveness, and configuration. The following
    configurations are supported:
        
        * Counterflow (ex. double-pipe)
        * Parallel (ex. double pipe inefficient configuration)
        * Shell and tube exchangers with even numbers of tube passes,
          one or more shells in series (TEMA E (one pass shell) only)
        * Crossflow, single pass, fluids unmixed
        * Crossflow, single pass, Cmax mixed, Cmin unmixed
        * Crossflow, single pass, Cmin mixed, Cmax unmixed
        * Boiler or condenser

    These situations are normally not those which occur in real heat exchangers,
    but are useful for academic purposes. More complicated expressions are 
    available for other methods. These equations are confirmed in [1]_, [2]_,
    and [3]_.
    
    For parallel flow heat exchangers:

    .. math::
        NTU = -\frac{\ln[1 - \epsilon(1+C_r)]}{1+C_r}
        
    For counterflow heat exchangers:

    .. math::
        NTU = \frac{1}{C_r-1}\ln\left(\frac{\epsilon-1}{\epsilon C_r-1}\right)
        
        NTU = \frac{\epsilon}{1-\epsilon} \text{ if } C_r = 1

    For TEMA E shell-and-tube heat exchangers with one shell pass, 2n tube 
    passes:

    .. math::
        (NTU)_1 = -(1 + C_r^2)^{-0.5}\ln\left(\frac{E-1}{E+1}\right)
        
        E = \frac{2/\epsilon_1 - (1 + C_r)}{(1 + C_r^2)^{0.5}}

    For TEMA E shell-and-tube heat exchangers with more than one shell pass, 2n  
    tube passes (this model assumes each exchanger has an equal share of the 
    overall NTU or said more plainly, the same UA):

    .. math::
        \epsilon_1 = \frac{F-1}{F-C_r}
        
        F = \left(\frac{\epsilon C_r-1}{\epsilon-1}\right)^{1/n}
        
        NTU = n(NTU)_1
        
    For cross-flow (single-pass) heat exchangers with both fluids unmixed, 
    there is no analytical solution. However, the function is monotonically
    increasing, and a closed-form solver is implemented as 'crossflow approximate',
    guaranteed to solve for :math:`10^{-7} < NTU < 10^5`. The exact solution
    for 'crossflow' uses the approximate solution's initial guess as a starting
    point for Newton's method. Some issues are noted at effectivenesses higher
    than 0.9 and very high NTUs, because the numerical integral term approaches
    1 too quickly.

    For cross-flow (single-pass) heat exchangers with Cmax mixed, Cmin unmixed:

    .. math::
        NTU = -\ln\left[1 + \frac{1}{C_r}\ln(1 - \epsilon C_r)\right]
        
    For cross-flow (single-pass) heat exchangers with Cmin mixed, Cmax unmixed:

    .. math::
        NTU = -\frac{1}{C_r}\ln[C_r\ln(1-\epsilon)+1]

    For cases where `Cr` = 0, as in an exchanger with latent heat exchange,
    flow arrangement does not matter: 

    .. math::
        NTU = -\ln(1-\epsilon)

    Parameters
    ----------
    effectiveness : float
        The thermal effectiveness of the heat exchanger, [-]
    Cr : float
        The heat capacity rate ratio, of the smaller fluid to the larger
        fluid, [-]
    subtype : str, optional
        The subtype of exchanger; one of 'counterflow', 'parallel', 'crossflow'
        'crossflow approximate', 'crossflow, mixed Cmin', 
        'crossflow, mixed Cmax', 'boiler', 'condenser', 'S&T', or 'nS&T' where 
        n is the number of shell and tube exchangers in a row.

    Returns
    -------
    NTU : float
        Thermal Number of Transfer Units [-]

    Notes
    -----
    Unlike :obj:`ht.hx.effectiveness_from_NTU`, not all inputs can 
    calculate the NTU - many exchanger types have effectiveness limits
    below 1 which depend on `Cr` and the number of shells in the case of
    heat exchangers. If an impossible input is given, an error will be raised
    and the maximum possible effectiveness will be printed.
    
    >>> NTU_from_effectiveness(.99, Cr=.7, subtype='5S&T')
    Traceback (most recent call last):
    Exception: The specified effectiveness is not physically possible for this configuration; the maximum effectiveness possible is 0.974122977755.

    Examples
    --------
    Worst case, parallel flow:
    
    >>> NTU_from_effectiveness(effectiveness=0.5881156068417585, Cr=0.7, subtype='parallel')
    5.000000000000012
    
    Crossflow, somewhat higher effectiveness:
        
    >>> NTU_from_effectiveness(effectiveness=0.8444821799748551, Cr=0.7, subtype='crossflow')
    5.000000000000859

    Counterflow, better than either crossflow or parallel flow:

    >>> NTU_from_effectiveness(effectiveness=0.9206703686051108, Cr=0.7, subtype='counterflow')
    5.0
    
    Shell and tube exchangers:
    
    >>> NTU_from_effectiveness(effectiveness=0.6834977044311439, Cr=0.7, subtype='S&T')
    5.000000000000071
    >>> NTU_from_effectiveness(effectiveness=0.9205058702789254, Cr=0.7, subtype='50S&T')
    4.999999999999996


    Overall case of rating an existing heat exchanger where a known flowrate
    of steam and oil are contacted in crossflow, with the steam side mixed,
    known inlet and outlet temperatures, and unknown UA
    (based on example 10-8 in [3]_):

    >>> Cp_oil = 1900 # J/kg/K
    >>> Cp_steam = 1860 # J/kg/K
    >>> m_steam = 5.2 # kg/s
    >>> m_oil = 1.45 # kg/s
    >>> Thi = 130 # °C
    >>> Tci = 15 # °C
    >>> Tco = 85 # °C # Design specification
    >>> Q = Cp_oil*m_oil*(Tci-Tco)
    >>> dTh = Q/(m_steam*Cp_steam)
    >>> Tho = Thi + dTh
    >>> Cmin = calc_Cmin(mh=m_steam, mc=m_oil, Cph=Cp_steam, Cpc=Cp_oil)
    >>> Cmax = calc_Cmax(mh=m_steam, mc=m_oil, Cph=Cp_steam, Cpc=Cp_oil)
    >>> Cr = calc_Cr(mh=m_steam, mc=m_oil, Cph=Cp_steam, Cpc=Cp_oil)
    >>> effectiveness = -Q/Cmin/(Thi-Tci)
    >>> NTU = NTU_from_effectiveness(effectiveness, Cr, subtype='crossflow, mixed Cmax')
    >>> UA = UA_from_NTU(NTU, Cmin)
    >>> U = 200 # Assume this was calculated; would actually need to be obtained iteratively as U depends on the exchanger geometry
    >>> A = UA/U
    >>> Tho, Cmin, Cmax, Cr
    (110.06100082712986, 2755.0, 9672.0, 0.2848428453267163)
    >>> effectiveness, NTU, UA, A
    (0.6086956521739131, 1.1040839095588, 3041.751170834494, 15.208755854172471)

    References
    ----------
    .. [1] Bergman, Theodore L., Adrienne S. Lavine, Frank P. Incropera, and
       David P. DeWitt. Introduction to Heat Transfer. 6E. Hoboken, NJ:
       Wiley, 2011.
    .. [2] Shah, Ramesh K., and Dusan P. Sekulic. Fundamentals of Heat 
       Exchanger Design. 1st edition. Hoboken, NJ: Wiley, 2002.
    .. [3] Holman, Jack. Heat Transfer. 10th edition. Boston: McGraw-Hill 
       Education, 2009.

def NTU_from_effectiveness(effectiveness, Cr, subtype='counterflow'):
    r'''Returns the Number of Transfer Units of a heat exchanger at a specified 
    heat capacity rate, effectiveness, and configuration. The following
    configurations are supported:
        
        * Counterflow (ex. double-pipe)
        * Parallel (ex. double pipe inefficient configuration)
        * Shell and tube exchangers with even numbers of tube passes,
          one or more shells in series (TEMA E (one pass shell) only)
        * Crossflow, single pass, fluids unmixed
        * Crossflow, single pass, Cmax mixed, Cmin unmixed
        * Crossflow, single pass, Cmin mixed, Cmax unmixed
        * Boiler or condenser

    These situations are normally not those which occur in real heat exchangers,
    but are useful for academic purposes. More complicated expressions are 
    available for other methods. These equations are confirmed in [1]_, [2]_,
    and [3]_.
    
    For parallel flow heat exchangers:

    .. math::
        NTU = -\frac{\ln[1 - \epsilon(1+C_r)]}{1+C_r}
        
    For counterflow heat exchangers:

    .. math::
        NTU = \frac{1}{C_r-1}\ln\left(\frac{\epsilon-1}{\epsilon C_r-1}\right)
        
        NTU = \frac{\epsilon}{1-\epsilon} \text{ if } C_r = 1

    For TEMA E shell-and-tube heat exchangers with one shell pass, 2n tube 
    passes:

    .. math::
        (NTU)_1 = -(1 + C_r^2)^{-0.5}\ln\left(\frac{E-1}{E+1}\right)
        
        E = \frac{2/\epsilon_1 - (1 + C_r)}{(1 + C_r^2)^{0.5}}

    For TEMA E shell-and-tube heat exchangers with more than one shell pass, 2n  
    tube passes (this model assumes each exchanger has an equal share of the 
    overall NTU or said more plainly, the same UA):

    .. math::
        \epsilon_1 = \frac{F-1}{F-C_r}
        
        F = \left(\frac{\epsilon C_r-1}{\epsilon-1}\right)^{1/n}
        
        NTU = n(NTU)_1
        
    For cross-flow (single-pass) heat exchangers with both fluids unmixed, 
    there is no analytical solution. However, the function is monotonically
    increasing, and a closed-form solver is implemented as 'crossflow approximate',
    guaranteed to solve for :math:`10^{-7} < NTU < 10^5`. The exact solution
    for 'crossflow' uses the approximate solution's initial guess as a starting
    point for Newton's method. Some issues are noted at effectivenesses higher
    than 0.9 and very high NTUs, because the numerical integral term approaches
    1 too quickly.

    For cross-flow (single-pass) heat exchangers with Cmax mixed, Cmin unmixed:

    .. math::
        NTU = -\ln\left[1 + \frac{1}{C_r}\ln(1 - \epsilon C_r)\right]
        
    For cross-flow (single-pass) heat exchangers with Cmin mixed, Cmax unmixed:

    .. math::
        NTU = -\frac{1}{C_r}\ln[C_r\ln(1-\epsilon)+1]

    For cases where `Cr` = 0, as in an exchanger with latent heat exchange,
    flow arrangement does not matter: 

    .. math::
        NTU = -\ln(1-\epsilon)

    Parameters
    ----------
    effectiveness : float
        The thermal effectiveness of the heat exchanger, [-]
    Cr : float
        The heat capacity rate ratio, of the smaller fluid to the larger
        fluid, [-]
    subtype : str, optional
        The subtype of exchanger; one of 'counterflow', 'parallel', 'crossflow'
        'crossflow approximate', 'crossflow, mixed Cmin', 
        'crossflow, mixed Cmax', 'boiler', 'condenser', 'S&T', or 'nS&T' where 
        n is the number of shell and tube exchangers in a row.

    Returns
    -------
    NTU : float
        Thermal Number of Transfer Units [-]

    Notes
    -----
    Unlike :obj:`ht.hx.effectiveness_from_NTU`, not all inputs can 
    calculate the NTU - many exchanger types have effectiveness limits
    below 1 which depend on `Cr` and the number of shells in the case of
    heat exchangers. If an impossible input is given, an error will be raised
    and the maximum possible effectiveness will be printed.
    
    >>> NTU_from_effectiveness(.99, Cr=.7, subtype='5S&T')
    Traceback (most recent call last):
    Exception: The specified effectiveness is not physically possible for this configuration; the maximum effectiveness possible is 0.974122977755.

    Examples
    --------
    Worst case, parallel flow:
    
    >>> NTU_from_effectiveness(effectiveness=0.5881156068417585, Cr=0.7, subtype='parallel')
    5.000000000000012
    
    Crossflow, somewhat higher effectiveness:
        
    >>> NTU_from_effectiveness(effectiveness=0.8444821799748551, Cr=0.7, subtype='crossflow')
    5.000000000000859

    Counterflow, better than either crossflow or parallel flow:

    >>> NTU_from_effectiveness(effectiveness=0.9206703686051108, Cr=0.7, subtype='counterflow')
    5.0
    
    Shell and tube exchangers:
    
    >>> NTU_from_effectiveness(effectiveness=0.6834977044311439, Cr=0.7, subtype='S&T')
    5.000000000000071
    >>> NTU_from_effectiveness(effectiveness=0.9205058702789254, Cr=0.7, subtype='50S&T')
    4.999999999999996


    Overall case of rating an existing heat exchanger where a known flowrate
    of steam and oil are contacted in crossflow, with the steam side mixed,
    known inlet and outlet temperatures, and unknown UA
    (based on example 10-8 in [3]_):

    >>> Cp_oil = 1900 # J/kg/K
    >>> Cp_steam = 1860 # J/kg/K
    >>> m_steam = 5.2 # kg/s
    >>> m_oil = 1.45 # kg/s
    >>> Thi = 130 # °C
    >>> Tci = 15 # °C
    >>> Tco = 85 # °C # Design specification
    >>> Q = Cp_oil*m_oil*(Tci-Tco)
    >>> dTh = Q/(m_steam*Cp_steam)
    >>> Tho = Thi + dTh
    >>> Cmin = calc_Cmin(mh=m_steam, mc=m_oil, Cph=Cp_steam, Cpc=Cp_oil)
    >>> Cmax = calc_Cmax(mh=m_steam, mc=m_oil, Cph=Cp_steam, Cpc=Cp_oil)
    >>> Cr = calc_Cr(mh=m_steam, mc=m_oil, Cph=Cp_steam, Cpc=Cp_oil)
    >>> effectiveness = -Q/Cmin/(Thi-Tci)
    >>> NTU = NTU_from_effectiveness(effectiveness, Cr, subtype='crossflow, mixed Cmax')
    >>> UA = UA_from_NTU(NTU, Cmin)
    >>> U = 200 # Assume this was calculated; would actually need to be obtained iteratively as U depends on the exchanger geometry
    >>> A = UA/U
    >>> Tho, Cmin, Cmax, Cr
    (110.06100082712986, 2755.0, 9672.0, 0.2848428453267163)
    >>> effectiveness, NTU, UA, A
    (0.6086956521739131, 1.1040839095588, 3041.751170834494, 15.208755854172471)

    References
    ----------
    .. [1] Bergman, Theodore L., Adrienne S. Lavine, Frank P. Incropera, and
       David P. DeWitt. Introduction to Heat Transfer. 6E. Hoboken, NJ:
       Wiley, 2011.
    .. [2] Shah, Ramesh K., and Dusan P. Sekulic. Fundamentals of Heat 
       Exchanger Design. 1st edition. Hoboken, NJ: Wiley, 2002.
    .. [3] Holman, Jack. Heat Transfer. 10th edition. Boston: McGraw-Hill 
       Education, 2009.
    '''
    if Cr > 1:
        raise Exception('Heat capacity rate must be less than 1 by definition.')

    if subtype == 'counterflow':
        # [2]_ gives the expression 1./(1-Cr)*log((1-Cr*eff)/(1-eff)), but
        # this is just the same equation rearranged differently.
        if Cr < 1:
            return 1./(Cr - 1.)*log((effectiveness - 1.)/(effectiveness*Cr - 1.))
        elif Cr == 1:
            return effectiveness/(1. - effectiveness)
    elif subtype == 'parallel':
        if effectiveness*(1. + Cr) > 1:
            raise Exception('The specified effectiveness is not physically \
possible for this configuration; the maximum effectiveness possible is %s.' % (1./(Cr + 1.)))
        return -log(1. - effectiveness*(1. + Cr))/(1. + Cr)
    elif 'S&T' in subtype:
        # [2]_ gives the expression
        # D = (1+Cr**2)**0.5
        # 1/D*log((2-eff*(1+Cr-D))/(2-eff*(1+Cr + D)))
        # This is confirmed numerically to be the same equation rearranged
        # differently
        str_shells = subtype.split('S&T')[0]
        shells = int(str_shells) if str_shells else 1
        
        F = ((effectiveness*Cr - 1.)/(effectiveness - 1.))**(1./shells)
        e1 = (F - 1.)/(F - Cr)
        E = (2./e1 - (1. + Cr))/(1. + Cr**2)**0.5
        
        if (E - 1.)/(E + 1.) <= 0:
            # Derived with SymPy
            max_effectiveness = (-((-Cr + sqrt(Cr**2 + 1) + 1)/(Cr + sqrt(Cr**2 + 1) - 1))**shells + 1)/(Cr - ((-Cr + sqrt(Cr**2 + 1) + 1)/(Cr + sqrt(Cr**2 + 1) - 1))**shells)
            raise Exception('The specified effectiveness is not physically \
possible for this configuration; the maximum effectiveness possible is %s.' % (max_effectiveness))
        
        NTU = -(1. + Cr*Cr)**-0.5*log((E - 1.)/(E + 1.))
        return shells*NTU
    elif subtype == 'crossflow':
        # Can't use a bisect solver here because at high NTU there's a derivative of 0
        # due to the integral term not changing when it's very near one
        guess = NTU_from_effectiveness(effectiveness, Cr, 'crossflow approximate')
        def to_solve(NTU, Cr, effectiveness):
            return effectiveness_from_NTU(NTU, Cr, subtype='crossflow') - effectiveness
        return newton(to_solve, guess, args=(Cr, effectiveness))
    elif subtype == 'crossflow approximate':
        # This will fail if NTU is more than 10,000 or less than 1E-7, but
        # this is extremely unlikely to occur in normal usage.
        # Maple and SymPy and Wolfram Alpha all failed to obtain an exact
        # analytical expression even with coefficients for 0.22 and 0.78 or 
        # with an explicit value for Cr. The function has been plotted,
        # and appears to be monotonic - there is only one solution.
        def to_solve(NTU, Cr, effectiveness):
            return (1. - exp(1./Cr*NTU**0.22*(exp(-Cr*NTU**0.78) - 1.))) - effectiveness
        return ridder(to_solve, a=1E-7, b=1E5, args=(Cr, effectiveness))
    
    elif subtype == 'crossflow, mixed Cmin':
        if Cr*log(1. - effectiveness) < -1:
            raise Exception('The specified effectiveness is not physically \
possible for this configuration; the maximum effectiveness possible is %s.' % (1. - exp(-1./Cr)))
        return -1./Cr*log(Cr*log(1. - effectiveness) + 1.)
    
    elif subtype ==  'crossflow, mixed Cmax':
        if 1./Cr*log(1. - effectiveness*Cr) < -1:
            raise Exception('The specified effectiveness is not physically \
possible for this configuration; the maximum effectiveness possible is %s.' % (((exp(Cr) - 1.0)*exp(-Cr)/Cr)))
        return -log(1. + 1./Cr*log(1. - effectiveness*Cr))
    
    elif subtype in ['boiler', 'condenser']:
        return -log(1. - effectiveness)
    else:
        raise Exception('Input heat exchanger type not recognized')