def surface_density_toy_model(r, phi, param):
    """

    Surface density of the toy model matter blob.

    INPUT:

    - ``r`` -- Boyer-Lindquist radial coordinate `\\bar{r}` in the matter blob
    - ``phi`` -- Boyer-Lindquist azimuthal coordinate `\\bar{\\phi}` in the
      matter blob
    - ``param`` -- list of parameters defining the position and extent of the
      matter blob:

      - ``param[0]``: mean radius `r_0` (Boyer-Lindquist coordinate)
      - ``param[1]``: mean azimuthal angle `\\phi_0` (Boyer-Lindquist
        coordinate)
      - ``param[2]``: radial extent `\\lambda`
      - ``param[3]``: opening angle `\\Delta\\phi`
      - ``param[4]`` (optional): amplitude `\\Sigma_0`; if not provided,
        then `\\Sigma_0=1` is used

    OUTPUT:

    - surface density `\\Sigma(\\bar{r}, \\bar{\\phi})`

    EXAMPLES:

    Use with the default amplitude (`\\Sigma_0=1`)::

        sage: from kerrgeodesic_gw import surface_density_toy_model
        sage: param = [6.5, 0, 0.6, 0.1]
        sage: surface_density_toy_model(6.5, 0, param)
        1.0
        sage: surface_density_toy_model(6.7, -0.04, param)
        1.0
        sage: surface_density_toy_model(6.3, 0.04, param)
        1.0
        sage: surface_density_toy_model(7, -0.06, param)
        0.0
        sage: surface_density_toy_model(5., 0.06, param)
        0.0

    3D representation: `z=\\Sigma(\\bar{r}, \\bar{\\phi})` in terms of
    `x:=\\bar{r}\\cos\\bar\\phi` and `y:=\\bar{r}\\sin\\bar\\phi`::

        sage: s_plot = lambda r, phi: surface_density_toy_model(r, phi, param)
        sage: r, phi, z = var('r phi z')
        sage: plot3d(s_plot, (r, 6, 8), (phi, -0.4, 0.4),
        ....:        transformation=(r*cos(phi), r*sin(phi), z))
        Graphics3d Object

    .. PLOT::

        from kerrgeodesic_gw import surface_density_toy_model
        param = [6.5, 0, 0.6, 0.1]
        s_plot = lambda r, phi: surface_density_toy_model(r, phi, param)
        r, phi, z = var('r phi z')
        g = plot3d(s_plot, (r, 6, 8), (phi, -0.4, 0.4), \\
                   transformation=(r*cos(phi), r*sin(phi), z))
        sphinx_plot(g)

    Use with a non-default amplitude (`\\Sigma_0=2`)::

        sage: sigma0 = 2.
        sage: param = [6.5, 0, 0.6, 0.1, sigma0]
        sage: surface_density_toy_model(6.5, 0, param)
        2.0

    """
    r0, phi0 = param[0], param[1]
    Dr, Dphi = param[2], param[3]
    Sigma0 = float(param[4]) if len(param) == 5 else float(1)
    if abs(r - r0) < Dr / float(2) and abs(phi - phi0) < Dphi / float(2):
        return Sigma0
    return float(0)
