def ut_haider(dp, mu, phi, rhog, rhos):
    """
    Calculate terminal velocity of a particle as discussed in the Haider and
    Levenspiel [5]_. Valid for particle sphericities of 0.5 to 1. Particle
    diameter should be an equivalent spherical diameter, such as the diameter
    of a sphere having same volume as the particle.

    To determine the terminal velocity for a range of particle sphericities, Haider
    and Levenspiel first define two dimensionless quantities

    .. math::

       d_{*} = d_p \\left[ \\frac{g\\, \\rho_g (\\rho_s - \\rho_g)}{\\mu^2} \\right]^{1/3} \\
       u_* = \\left[ \\frac{18}{d{_*}^2} + \\frac{2.3348 - 1.7439\\, \\phi}{d{_*}^{0.5}} \\right]^{-1}

    where :math:`0.5 \\leq \\phi \\leq 1` and particle diameter :math:`d_p` is an
    equivalent spherical diameter, the diameter of a sphere having the same volume
    as the particle. The relationship between :math:`u_*` and :math:`u_t` is given
    by

    .. math:: u_* = u_t \\left[ \\frac{\\rho{_g}^2}{g\\, \\mu\\, (\\rho_s - \\rho_g)} \\right]^{1/3}

    The terminal velocity of the particle can finally be determined by rearranging
    the above equation such that

    .. math:: u_t = u_* \\left[ \\frac{g\\, \\mu\\, (\\rho_s - \\rho_g)}{\\rho{_g}^2} \\right]^{1/3}

    Parameters
    ----------
    dp : float
        Diameter of particle [m]
    mu : float
        Viscosity of gas [kg/(m s)]
    phi : float
        Sphericity of particle [-]
    rhog : float
        Density of gas [kg/m^3]
    rhos : float
        Density of particle [kg/m^3]

    Returns
    -------
    ut : float
        Terminal velocity of a particle [m/s]

    Example
    -------
    >>> ut_haider(0.00016, 1.8e-5, 0.67, 1.2, 2600)
    0.8857

    References
    ----------
    .. [5] A. Haider and O. Levenspiel. Drag coefficient and terminal velocity
           of spherical and nonspherical particles. Powder Technology,
           58:63–70, 1989.
    """
    if phi > 1.0 or phi < 0.5:
        raise ValueError('Sphericity must be 0.5 <= phi <= 1.0')

    d_star = dp * ((9.81 * rhog * (rhos - rhog)) / (mu**2))**(1 / 3)
    u_star = (18 / (d_star**2) + ((2.3348 - 1.7439 * phi) / (d_star**0.5)))**-1
    ut = u_star * ((9.81 * (rhos - rhog) * mu) / rhog**2)**(1 / 3)
    return ut