def temperature_effectiveness_plate(R1, NTU1, Np1, Np2, counterflow=True, 
                                    passes_counterflow=True, reverse=False):
    r'''Returns the temperature effectiveness `P1` of side 1 of a plate heat 
    exchanger with a specified side 1 heat capacity ratio `R1`, side 1 number
    of transfer units `NTU1`, number of passes on sides 1 and 2 (respectively
    `Np1` and `Np2`). 
        
    For all cases, the function also takes as arguments whether the exchanger 
    is setup in an overall counter or parallel orientation `counterflow`, and 
    whether or not individual stream passes are themselves counterflow or
    parallel. 
    
    The 20 supported cases are as follows. (the first number of sides listed
    refers to side 1, and the second number refers to side 2):
        
    * 1 pass/1 pass parallelflow
    * 1 pass/1 pass counterflow
    * 1 pass/2 pass
    * 1 pass/3 pass or 3 pass/1 pass (with the two end passes in parallel)
    * 1 pass/3 pass or 3 pass/1 pass (with the two end passes in counterflow)
    * 1 pass/4 pass 
    * 2 pass/2 pass, overall parallelflow, individual passes in parallel 
    * 2 pass/2 pass, overall parallelflow, individual passes counterflow
    * 2 pass/2 pass, overall counterflow, individual passes parallelflow 
    * 2 pass/2 pass, overall counterflow, individual passes counterflow 
    * 2 pass/3 pass or 3 pass/2 pass, overall parallelflow 
    * 2 pass/3 pass or 3 pass/2 pass, overall counterflow
    * 2 pass/4 pass or 4 pass/2 pass, overall parallel flow
    * 2 pass/4 pass or 4 pass/2 pass, overall counterflow flow
    
    For all plate heat exchangers, there are two common formulas used by most
    of the expressions.
    
    .. math::
        P_p(x, y) = \frac{1 - \exp[-x(1 + y)]}{1 + y}
        
        P_c(x, y) = \frac{1 - \exp[-x(1 - y)]}{1 - y\exp[-x(1 - y)]}
        
    The main formulas used are as follows. Note that for some cases such as
    4 pass/2 pass, the formula is not shown because it is that of 2 pass/4 
    pass, but with R1, NTU1, and P1 conversions.
        
    For 1 pass/1 pass paralleflow (streams symmetric):
        
    .. math::
        P_1 = P_p(NTU_1, R_1)
        
    For 1 pass/1 pass counterflow (streams symmetric):
    
    .. math::
        P_1 = P_c(NTU_1, R_1)
            
    For 1 pass/2 pass (any of the four possible configurations):
        
    .. math::
        P_1 = 0.5(A + B - 0.5ABR_1)
        
        A = P_p(NTU_1, 0.5R_1)
        
        B = P_c(NTU_1, 0.5R_1)
        
    For 1 pass/3 pass (with the two end passes in parallel):
        
    .. math::
        P_1 = \frac{1}{3}\left[B + A\left(1 - \frac{R_1 B}{3}\right)\left(2 
        - \frac{R_1 A}{3}\right)\right]
        
        A = P_p\left(NTU_1, \frac{R_1}{3}\right)
        
        B = P_c\left(NTU_1, \frac{R_1}{3}\right)
        
    For 1 pass/3 pass (with the two end passes in counterflow):
        
    .. math::
        P_1 = \frac{1}{3}\left[A + B\left(1 - \frac{R_1 A}{3}\right)\left(2
        - \frac{R_1 B}{3}\right)\right]
            
        A = P_p\left(NTU_1, \frac{R_1}{3}\right)
        
        B = P_c\left(NTU_1, \frac{R_1}{3}\right)
        
    For 1 pass/4 pass (any of the four possible configurations):
    
    .. math::
        P_1 = \frac{1-Q}{R_1}
        
        Q = \left(1 - \frac{AR_1}{4}\right)^2\left(1 - \frac{BR_1}{4}\right)^2
        
        A = P_p\left(NTU_1, \frac{R_1}{4}\right)
        
        B = P_c\left(NTU_1, \frac{R_1}{4}\right)
        
    For 2 pass/2 pass, overall parallelflow, individual passes in parallel 
    (stream symmetric):
        
    .. math::
        P_1 = P_p(NTU_1, R_1)
        
    For 2 pass/2 pass, overall parallelflow, individual passes counterflow
    (stream symmetric):
        
    .. math::
        P_1 = B[2 - B(1 + R_1)]
        
        B = P_c\left(\frac{NTU_1}{2}, R_1\right)
        
    For 2 pass/2 pass, overall counterflow, individual passes parallelflow 
    (stream symmetric):
        
    .. math::
        P_1 = \frac{2A - A^2(1 + R_1)}{1 - R_1 A^2}
        
        A = P_p\left(\frac{NTU_1}{2}, R_1\right)
        
    For 2 pass/2 pass, overall counterflow and individual passes counterflow 
    (stream symmetric):
        
    .. math::
        P_1 = P_c(NTU_1, R_1)
        
    For 2 pass/3 pass, overall parallelflow:
        
    .. math::
        P_1 = A + B - \left(\frac{2}{9} + \frac{D}{3}\right)
        (A^2 + B^2) - \left(\frac{5}{9} + \frac{4D}{3}\right)AB
        + \frac{D(1+D)AB(A+B)}{3} - \frac{D^2A^2B^2}{9}
        
        A = P_p\left(\frac{NTU_1}{2}, D\right)
        
        B = P_c\left(\frac{NTU_1}{2}, D\right)
        
        D = \frac{2R_1}{3}
        
    For 2 pass/3 pass, overall counterflow:
        
    .. math::
        P_1 = \frac{A + 0.5B + 0.5C + D}{R_1}
        
        A = \frac{2R_1 EF^2 - 2EF + F - F^2}
        {2R_1 E^2 F^2 - E^2 - F^2 - 2EF + E + F}
        
        B = \frac{A(E-1)}{F}
        
        C = \frac{1 - A}{E}
        
        D = R_1 E^2 C - R_1 E + R_1 - \frac{C}{2}
        
        E = \frac{3}{2R_1 G}
        
        F = \frac{3}{2R_1 H}
        
        G = P_c\left(\frac{NTU_1}{2}, \frac{2R_1}{3}\right)
        
        H = P_p\left(\frac{NTU_1}{2}, \frac{2R_1}{3}\right)
        
    For 2 pass/4 pass, overall parallel flow:
        
    .. math::
        P_1 = 2D - (1 + R_1)D^2
        
        D = \frac{A + B - 0.5ABR_1}{2}
        
        A = P_p\left(\frac{NTU_1}{2}, \frac{R_1}{2}\right)
        
        B = P_c\left(\frac{NTU_1}{2}, \frac{R_1}{2}\right)
        
    For 2 pass/4 pass, overall counterflow flow:
        
    .. math::
        P_1 = \frac{2D - (1+R_1)D^2}{1 - D^2 R_1}
        
        D = \frac{A + B - 0.5ABR_1}{2}
        
        A = P_p\left(\frac{NTU_1}{2}, \frac{R_1}{2}\right)
        
        B = P_c\left(\frac{NTU_1}{2}, \frac{R_1}{2}\right)
                
    Parameters
    ----------
    R1 : float
        Heat capacity ratio of the heat exchanger in the P-NTU method,
        calculated with respect to stream 1 [-]
    NTU1 : float
        Thermal number of transfer units of the heat exchanger in the P-NTU 
        method, calculated with respect to stream 1 [-]
    Np1 : int
        Number of passes on side 1 [-]
    Np2 : int
        Number of passes on side 2 [-]
    counterflow : bool
        Whether or not the overall flow through the heat exchanger is in
        counterflow or parallel flow, [-]
    passes_counterflow : bool
        In addition to the overall flow direction, in some cases individual 
        passes may be in counter or parallel flow; this controlls that [-]
    reverse : bool
        Used **internally only** to allow cases like the 1-4 formula to work  
        for the 4-1 flow case, without having to duplicate the code [-]

    Returns
    -------
    P1 : float
        Thermal effectiveness of the heat exchanger in the P-NTU method,
        calculated with respect to stream 1 [-]

    Notes
    -----
    For diagrams of these heat exchangers, see [3]_.
    In all cases, each pass is assumed to be made up of an infinite number
    of plates. The fluid velocities must be uniform across the plate channels,
    and the flow must be uniformly distributed between the channels. The heat
    transfer coefficient is also assumed constant.
    
    The defaults of counterflow=True and passes_counterflow=True will always
    result in the most efficient heat exchanger option, normally what is
    desired.
    
    If a number of passes which is not supported is provided, an exception is
    raised.

    Examples
    --------
    Three passes on side 1; one pass on side 2; two end passes in counterflow
    orientation.
    
    >>> temperature_effectiveness_plate(R1=1/3., NTU1=1., Np1=3, Np2=1)
    0.5743514352720835
    
    If the same heat exchanger (in terms of NTU1 and R1) were operating with
    sides 1 and 2 switched, a slightly less efficient design results.
    
    >>> temperature_effectiveness_plate(R1=1/3., NTU1=1., Np1=1, Np2=3)
    0.5718726757657066
    
    References
    ----------
    .. [1] Shah, Ramesh K., and Dusan P. Sekulic. Fundamentals of Heat 
       Exchanger Design. 1st edition. Hoboken, NJ: Wiley, 2002.
    .. [2] Rohsenow, Warren and James Hartnett and Young Cho. Handbook of Heat
       Transfer, 3E. New York: McGraw-Hill, 1998.
    .. [3] Kandlikar, S. G., and R. K. Shah. "Asymptotic Effectiveness-NTU 
       Formulas for Multipass Plate Heat Exchangers." Journal of Heat Transfer 
       111, no. 2 (May 1, 1989): 314-21. doi:10.1115/1.3250679.
    .. [4] Kandlikar, S. G., and R. K. Shah. "Multipass Plate Heat Exchangers
       Effectiveness-NTU Results and Guidelines for Selecting Pass 
       Arrangements." Journal of Heat Transfer 111, no. 2 (May 1, 1989): 
       300-313. doi:10.1115/1.3250678.   
    '''
    if Np1 == 1 and Np2 == 1 and counterflow:
        return Pc(NTU1, R1)
    elif Np1 == 1 and Np2 == 1 and not counterflow:
        return Pp(NTU1, R1)
    elif Np1 == 1 and Np2 == 2:
        # There are four configurations but all have the same formula
        # They do behave different depending on the number of available plates
        # but this model assues infinity
        # There are four more arrangements that are equivalent as well
        A = Pp(NTU1, 0.5*R1)
        B = Pc(NTU1, 0.5*R1)
        return 0.5*(A + B - 0.5*A*B*R1)
    elif Np1 == 1 and Np2 == 3 and counterflow:
        # There are six configurations, two formulas
        # Each behaves differently though as a function of number of plates
        A = Pp(NTU1, R1/3.)
        B = Pc(NTU1, R1/3.)
        return 1/3.*(A + B*(1. - R1*A/3.)*(2. - R1*B/3.))
    elif Np1 == 1 and Np2 == 3 and not counterflow:
        A = Pp(NTU1, R1/3.)
        B = Pc(NTU1, R1/3.)
        return 1/3.*(B + A*(1. - R1*B/3.)*(2. - R1*A/3.))
    elif Np1 == 1 and Np2 == 4:
        # four configurations
        # Again a function of number of plates, but because expressions assume
        # infinity it gets ignored and they're the same
        A = Pp(NTU1, 0.25*R1)
        B = Pc(NTU1, 0.25*R1)
        t1 = (1. - 0.25*A*R1)
        t2 = (1. - 0.25*B*R1)
        t3 = t1*t2 # minor optimization
        return (1. - t3*t3)/R1
    elif Np1 == 2 and Np2 == 2:
        if counterflow and passes_counterflow:
            return Pc(NTU1, R1)
        elif counterflow and not passes_counterflow:
            A = Pp(0.5*NTU1, R1)
            return (2.*A - A*A*(1. + R1))/(1. - R1*A*A)
        elif not counterflow and passes_counterflow:
            B = Pc(0.5*NTU1, R1)
            return B*(2. - B*(1. + R1))
        elif not counterflow and not passes_counterflow:
            return temperature_effectiveness_plate(R1, NTU1, Np1=1, Np2=1, 
                                                   counterflow=False)
    elif Np1 == 2 and Np2 == 3:
        # One place says there are four configurations; no other discussion is
        # presented
        if counterflow:
            H = Pp(0.5*NTU1, 2./3.*R1)
            G = Pc(0.5*NTU1, 2./3.*R1)
            E = 1./(2./3.*R1*G)
            F = 1./(2./3.*R1*H)
            E2 = E*E
            F2 = F*F
            A = (2.*R1*E*F2 - 2.*E*F + F - F2)/(2.*R1*E2*F2 - E2 - F2 - 2.*E*F + E + F)
            C = (1. - A)/E
            D = R1*E*E*C - R1*E + R1 - 0.5*C
            B = A*(E - 1.)/F
            return (A + 0.5*B + 0.5*C + D)/R1
        elif not counterflow:
            D = 2*R1/3.
            A = Pp(NTU1/2, D)
            B = Pc(NTU1/2, D)
            return (A + B - (2/9. + D/3.)*(A*A + B*B)
                    -(5./9. + 4./3.*D)*A*B
                    + D*(1. + D)*A*B*(A + B)/3.
                    - D*D*A*A*B*B/9.)
    elif Np1 == 2 and Np2 == 4:
        # Both cases are correct for passes_counterflow=True or False
        if counterflow:
            A = Pp(0.5*NTU1, 0.5*R1)
            B = Pc(0.5*NTU1, 0.5*R1)
            D = 0.5*(A + B - 0.5*A*B*R1)
            return (2.*D - (1. + R1)*D*D)/(1. - D*D*R1)
        elif not counterflow:
            A = Pp(0.5*NTU1, 0.5*R1)
            B = Pc(0.5*NTU1, 0.5*R1)
            D = 0.5*(A + B - 0.5*A*B*R1)
            return 2.*D - ((1. + R1)*D*D)
    if not reverse:
        # only the assymetric cases will be able to solve by flipping things
        # Note that asymmetric performs differently depending on the arguments
        # The user still needs to input R1, NTU1 for side 1
        # so if they want to do a 3-1 instead of a 1-3 as is implemented here
        # They give R1 and NTU1 for "3 pass" side instead of the "1 pass" side 
        # and will get back P1 for the "3 pass" side.
        R2 = 1./R1
        NTU2 = NTU1*R1
        P2 = temperature_effectiveness_plate(R1=R2, NTU1=NTU2, Np1=Np2, Np2=Np1,
                                             counterflow=counterflow, 
                                             passes_counterflow=passes_counterflow, 
                                             reverse=True)
        P1 = P2*R2
        return P1
    
    raise Exception('Supported number of passes does not have a formula available')