def boundary(costheta, a=1, epsilon=0.1, nu=0):
    """Projected boundary of a prolate spheroid

    Compute the boundary according to equation (4) in
    :cite:`Boyde2009` with the addition of the
    Poisson's ratio of the object.

    .. math::

       B(\\theta) = a (1+\\epsilon)
       \\left[ (1+\\epsilon)^2 - \\epsilon (1+\\nu)
       (2+\\epsilon (1-\\nu)) \\cos^2 \\theta \\right]^{-1/2}

    This boundary function was derived for a prolate spheroid under
    the assumption that the semi-major axis :math:`a` and the
    semi-minor axes :math:`b=c` are defined as

    .. math::

       a = b \\cdot \\frac{1+ \\epsilon}{1- \\nu \\epsilon}

    The boundary function :math:`B(\\theta)` can be derived with
    the above relation using the equation for a prolate spheroid.

    Parameters
    ----------
    costheta: float or np.ndarray
        Cosine of polar coordinates :math:`\\theta`
        at which to compute the boundary.
    a: float
        Equatorial radii of prolate spheroid (semi-minor axis).
    epsilon: float
        Stretch ratio; defines size of semi-major axis:
        :math:`a = (1+\\epsilon) b`. Note that this is not
        the eccentricity of the prolate spheroid.
    nu: float
        Poisson's ratio :math:`\\nu` of the material.

    Returns
    -------
    B: 1d ndarray
        Radial object boundary in dependence of theta
        :math:`B(\\theta)`.

    Notes
    -----
    For :math:`\\nu=0`, the above equation becomes
    equation (4) in :cite:`Boyde2009`.
    """
    x = costheta
    B = a * (1 + epsilon) / ((1 + epsilon) ** 2 - epsilon * (1 + nu) * (2 +
        epsilon * (1 - nu)) * x ** 2) ** 0.5
    return B
