def _parametric_plot3d_curve(f, urange, plot_points, **kwds):
    r"""
    This function is used internally by the
    ``parametric_plot3d`` command.
    """
    from sage.plot.misc import setup_for_eval_on_grid
    g, ranges = setup_for_eval_on_grid(f, [urange], plot_points)
    f_x,f_y,f_z = g
    w = [(f_x(u), f_y(u), f_z(u)) for u in xsrange(*ranges[0], include_endpoint=True)]
    return line3d(w, **kwds)

def _parametric_plot3d_curve(f, urange, plot_points, **kwds):
    r"""
    Return a parametric three-dimensional space curve.
    This function is used internally by the
    :func:`parametric_plot3d` command.

    There are two ways this function is invoked by
    :func:`parametric_plot3d`.

    - ``parametric_plot3d([f_x, f_y, f_z], (u_min,
      u_max))``:
      `f_x, f_y, f_z` are three functions and
      `u_{\min}` and `u_{\max}` are real numbers

    - ``parametric_plot3d([f_x, f_y, f_z], (u, u_min,
      u_max))``:
      `f_x, f_y, f_z` can be viewed as functions of
      `u`

    INPUT:

    - ``f`` - a 3-tuple of functions or expressions, or vector of size 3

    - ``urange`` - a 2-tuple (u_min, u_max) or a 3-tuple
      (u, u_min, u_max)

    - ``plot_points`` - (default: "automatic", which is 75) initial
      number of sample points in each parameter; an integer.

    EXAMPLES:

    We demonstrate each of the two ways of calling this.  See
    :func:`parametric_plot3d` for many more examples.

    We do the first one with a lambda function, which creates a
    callable Python function that sends `u` to `u/10`::

        sage: parametric_plot3d((sin, cos, lambda u: u/10), (0,20)) # indirect doctest
        Graphics3d Object

    Now we do the same thing with symbolic expressions::

        sage: u = var('u')
        sage: parametric_plot3d((sin(u), cos(u), u/10), (u,0,20))
        Graphics3d Object

    """
    from sage.plot.misc import setup_for_eval_on_grid
    g, ranges = setup_for_eval_on_grid(f, [urange], plot_points)
    f_x, f_y, f_z = g
    w = [(f_x(u), f_y(u), f_z(u)) for u in xsrange(*ranges[0], include_endpoint=True)]
    return line3d(w, **kwds)

def _parametric_plot3d_surface(f, urange, vrange, plot_points, boundary_style, **kwds):
    r"""
    This function is used internally by the
    ``parametric_plot3d`` command.
    """
    from sage.plot.misc import setup_for_eval_on_grid
    g, ranges = setup_for_eval_on_grid(f, [urange,vrange], plot_points)
    urange = srange(*ranges[0], include_endpoint=True)
    vrange = srange(*ranges[1], include_endpoint=True)
    G = ParametricSurface(g, (urange, vrange), **kwds)
    
    if boundary_style is not None:
        for u in (urange[0], urange[-1]):
            G += line3d([(g[0](u,v), g[1](u,v), g[2](u,v)) for v in vrange], **boundary_style)
        for v in (vrange[0], vrange[-1]):
            G += line3d([(g[0](u,v), g[1](u,v), g[2](u,v)) for u in urange], **boundary_style)
    return G

#.........这里部分代码省略.........
      ``(x,xmin,xmax)``

    - ``(ymin, ymax)`` -- 2-tuple, the range of ``y`` values OR 3-tuple
      ``(y,ymin,ymax)``

    The following inputs must all be passed in as named parameters:

    - ``plot_points`` -- integer (default: 25); number of points to plot
      in each direction of the grid

    - ``cmap`` -- a colormap (type ``cmap_help()`` for more information).

    - ``interpolation`` -- string (default: ``'catrom'``), the interpolation
      method to use: ``'bilinear'``, ``'bicubic'``, ``'spline16'``,
      ``'spline36'``, ``'quadric'``, ``'gaussian'``, ``'sinc'``,
      ``'bessel'``, ``'mitchell'``, ``'lanczos'``, ``'catrom'``,
      ``'hermite'``, ``'hanning'``, ``'hamming'``, ``'kaiser'``


    EXAMPLES:

    Here we plot a simple function of two variables.  Note that
    since the input function is an expression, we need to explicitly
    declare the variables in 3-tuples for the range::

        sage: x,y = var('x,y')
        sage: density_plot(sin(x)*sin(y), (x, -2, 2), (y, -2, 2))
        Graphics object consisting of 1 graphics primitive


    Here we change the ranges and add some options; note that here
    ``f`` is callable (has variables declared), so we can use 2-tuple ranges::

        sage: x,y = var('x,y')
        sage: f(x,y) = x^2*cos(x*y)
        sage: density_plot(f, (x,-10,5), (y, -5,5), interpolation='sinc', plot_points=100)
        Graphics object consisting of 1 graphics primitive

    An even more complicated plot::

        sage: x,y = var('x,y')
        sage: density_plot(sin(x^2 + y^2)*cos(x)*sin(y), (x, -4, 4), (y, -4, 4), cmap='jet', plot_points=100)
        Graphics object consisting of 1 graphics primitive

    This should show a "spotlight" right on the origin::

        sage: x,y = var('x,y')
        sage: density_plot(1/(x^10+y^10), (x, -10, 10), (y, -10, 10))
        Graphics object consisting of 1 graphics primitive

    Some elliptic curves, but with symbolic endpoints.  In the first
    example, the plot is rotated 90 degrees because we switch the
    variables `x`, `y`::

        sage: density_plot(y^2 + 1 - x^3 - x, (y,-pi,pi), (x,-pi,pi))
        Graphics object consisting of 1 graphics primitive

    ::

        sage: density_plot(y^2 + 1 - x^3 - x, (x,-pi,pi), (y,-pi,pi))
        Graphics object consisting of 1 graphics primitive

    Extra options will get passed on to show(), as long as they are valid::

        sage: density_plot(log(x) + log(y), (x, 1, 10), (y, 1, 10), dpi=20)
        Graphics object consisting of 1 graphics primitive

    ::

        sage: density_plot(log(x) + log(y), (x, 1, 10), (y, 1, 10)).show(dpi=20) # These are equivalent

    TESTS:

    Check that :trac:`15315` is fixed, i.e., density_plot respects the
    ``aspect_ratio`` parameter. Without the fix, it looks like a thin line
    of width a few mm. With the fix it should look like a nice fat layered
    image::

        sage: density_plot((x*y)^(1/2), (x,0,3), (y,0,500), aspect_ratio=.01)
        Graphics object consisting of 1 graphics primitive

    Default ``aspect_ratio`` is ``"automatic"``, and that should work too::

        sage: density_plot((x*y)^(1/2), (x,0,3), (y,0,500))
        Graphics object consisting of 1 graphics primitive

    """
    from sage.plot.all import Graphics
    from sage.plot.misc import setup_for_eval_on_grid
    g, ranges = setup_for_eval_on_grid([f], [xrange, yrange], options['plot_points'])
    g = g[0]
    xrange,yrange=[r[:2] for r in ranges]

    xy_data_array = [[g(x, y) for x in xsrange(*ranges[0], include_endpoint=True)]
                              for y in xsrange(*ranges[1], include_endpoint=True)]

    g = Graphics()
    g._set_extra_kwds(Graphics._extract_kwds_for_show(options, ignore=['xmin', 'xmax']))
    g.add_primitive(DensityPlot(xy_data_array, xrange, yrange, options))
    return g

#.........这里部分代码省略.........
    Here we plot a simple function of two variables::

        sage: x,y = var('x,y')
        sage: region_plot(cos(x^2+y^2) <= 0, (x, -3, 3), (y, -3, 3))
         
    Here we play with the colors::

        sage: region_plot(x^2+y^3 < 2, (x, -2, 2), (y, -2, 2), incol='lightblue', bordercol='gray')
        
    An even more complicated plot, with dashed borders::

        sage: region_plot(sin(x)*sin(y) >= 1/4, (x,-10,10), (y,-10,10), incol='yellow', bordercol='black', borderstyle='dashed', plot_points=250)

    A disk centered at the origin::

        sage: region_plot(x^2+y^2<1, (x,-1,1), (y,-1,1))

    A plot with more than one condition (all conditions must be true for the statement to be true)::

        sage: region_plot([x^2+y^2<1, x<y], (x,-2,2), (y,-2,2))

    Since it doesn't look very good, let's increase ``plot_points``::

        sage: region_plot([x^2+y^2<1, x<y], (x,-2,2), (y,-2,2), plot_points=400)

    To get plots where only one condition needs to be true, use a function.
    Using lambda functions, we definitely need the extra ``plot_points``::

        sage: region_plot(lambda x,y: x^2+y^2<1 or x<y, (x,-2,2), (y,-2,2), plot_points=400)
    
    The first quadrant of the unit circle::

        sage: region_plot([y>0, x>0, x^2+y^2<1], (x,-1.1, 1.1), (y,-1.1, 1.1), plot_points = 400)

    Here is another plot, with a huge border::

        sage: region_plot(x*(x-1)*(x+1)+y^2<0, (x, -3, 2), (y, -3, 3), incol='lightblue', bordercol='gray', borderwidth=10, plot_points=50)

    If we want to keep only the region where x is positive::

        sage: region_plot([x*(x-1)*(x+1)+y^2<0, x>-1], (x, -3, 2), (y, -3, 3), incol='lightblue', plot_points=50)

    Here we have a cut circle::

        sage: region_plot([x^2+y^2<4, x>-1], (x, -2, 2), (y, -2, 2), incol='lightblue', bordercol='gray', plot_points=200)

    The first variable range corresponds to the horizontal axis and
    the second variable range corresponds to the vertical axis::

        sage: s,t=var('s,t')
        sage: region_plot(s>0,(t,-2,2),(s,-2,2))

    ::

        sage: region_plot(s>0,(s,-2,2),(t,-2,2))

    """

    from sage.plot.plot import Graphics
    from sage.plot.misc import setup_for_eval_on_grid
    import numpy

    if not isinstance(f, (list, tuple)):
        f = [f]

    f = [equify(g) for g in f]

    g, ranges = setup_for_eval_on_grid(f, [xrange, yrange], plot_points)
    xrange,yrange=[r[:2] for r in ranges]

    xy_data_arrays = numpy.asarray([[[func(x, y) for x in xsrange(*ranges[0], include_endpoint=True)]
                                     for y in xsrange(*ranges[1], include_endpoint=True)]
                                    for func in g],dtype=float)
    xy_data_array=numpy.abs(xy_data_arrays.prod(axis=0))
    # Now we need to set entries to negative iff all
    # functions were negative at that point.
    neg_indices = (xy_data_arrays<0).all(axis=0)
    xy_data_array[neg_indices]=-xy_data_array[neg_indices]

    from matplotlib.colors import ListedColormap
    incol = rgbcolor(incol)
    outcol = rgbcolor(outcol)
    cmap = ListedColormap([incol, outcol])
    cmap.set_over(outcol)
    cmap.set_under(incol)
    
    g = Graphics()
    g._set_extra_kwds(Graphics._extract_kwds_for_show(options, ignore=['xmin', 'xmax']))
    g.add_primitive(ContourPlot(xy_data_array, xrange,yrange, 
                                dict(contours=[-1e307, 0, 1e307], cmap=cmap, fill=True, **options)))

    if bordercol or borderstyle or borderwidth:
        cmap = [rgbcolor(bordercol)] if bordercol else ['black']
        linestyles = [borderstyle] if borderstyle else None
        linewidths = [borderwidth] if borderwidth else None
        g.add_primitive(ContourPlot(xy_data_array, xrange, yrange, 
                                    dict(linestyles=linestyles, linewidths=linewidths,
                                         contours=[0], cmap=[bordercol], fill=False, **options)))
    
    return g

#.........这里部分代码省略.........
        sage: P

    ::

        sage: P=contour_plot(y^2 + 1 - x^3 - x, (x,-pi,pi), (y,-pi,pi), fill=False, cmap='hsv',labels=True,\
        ...    contours=[-4,0,4],  label_fmt=lambda x: "$z=%s$"%x, label_colors='black', label_inline=True, \
        ...    label_fontsize=12)
        sage: P

    ::

        sage: P=contour_plot(y^2 + 1 - x^3 - x, (x,-pi,pi), (y,-pi,pi), \
        ...    fill=False, cmap='hsv', labels=True, label_fontsize=18)
        sage: P

    ::

        sage: P=contour_plot(y^2 + 1 - x^3 - x, (x,-pi,pi), (y,-pi,pi), \
        ...    fill=False, cmap='hsv', labels=True, label_inline_spacing=1)
        sage: P

    ::

        sage: P= contour_plot(y^2 + 1 - x^3 - x, (x,-pi,pi), (y,-pi,pi), \
        ...    fill=False, cmap='hsv', labels=True, label_inline=False)
        sage: P
    
    We can change the color of the labels if so desired::

        sage: contour_plot(f, (-2,2), (-2,2), labels=True, label_colors='red')

    We can add a colorbar as well::

        sage: f(x,y)=x^2-y^2
        sage: contour_plot(f, (x,-3,3), (y,-3,3), colorbar=True)

    ::

        sage: contour_plot(f, (x,-3,3), (y,-3,3), colorbar=True,colorbar_orientation='horizontal')

    ::

        sage: contour_plot(f, (x,-3,3), (y,-3,3), contours=[-2,-1,4],colorbar=True)

    ::

        sage: contour_plot(f, (x,-3,3), (y,-3,3), contours=[-2,-1,4],colorbar=True,colorbar_spacing='uniform')

    ::

        sage: contour_plot(f, (x,-3,3), (y,-3,3), contours=[0,2,3,6],colorbar=True,colorbar_format='%.3f')

    ::

        sage: contour_plot(f, (x,-3,3), (y,-3,3), labels=True,label_colors='red',contours=[0,2,3,6],colorbar=True)

    ::

        sage: contour_plot(f, (x,-3,3), (y,-3,3), cmap='winter', contours=20, fill=False, colorbar=True)

    This should plot concentric circles centered at the origin::

        sage: x,y = var('x,y')
        sage: contour_plot(x^2+y^2-2,(x,-1,1), (y,-1,1))

    Extra options will get passed on to show(), as long as they are valid::

        sage: f(x, y) = cos(x) + sin(y)
        sage: contour_plot(f, (0, pi), (0, pi), axes=True)

    ::

        sage: contour_plot(f, (0, pi), (0, pi)).show(axes=True) # These are equivalent

    Note that with ``fill=False`` and grayscale contours, there is the 
    possibility of confusion between the contours and the axes, so use
    ``fill=False`` together with ``axes=True`` with caution::

        sage: contour_plot(f, (-pi, pi), (-pi, pi), fill=False, axes=True) 

    TESTS:

    To check that ticket 5221 is fixed, note that this has three curves, not two::

        sage: x,y = var('x,y')
        sage: contour_plot(x-y^2,(x,-5,5),(y,-3,3),contours=[-4,-2,0], fill=False)
    """
    from sage.plot.plot import Graphics
    from sage.plot.misc import setup_for_eval_on_grid
    g, ranges = setup_for_eval_on_grid([f], [xrange, yrange], options['plot_points'])
    g = g[0]
    xrange,yrange=[r[:2] for r in ranges]
    
    xy_data_array = [[g(x, y) for x in xsrange(*ranges[0], include_endpoint=True)]
                              for y in xsrange(*ranges[1], include_endpoint=True)]

    g = Graphics()
    g._set_extra_kwds(Graphics._extract_kwds_for_show(options, ignore=['xmin', 'xmax']))
    g.add_primitive(ContourPlot(xy_data_array, xrange, yrange, options))
    return g        

        sage: x,y = var('x,y')
        sage: plot_vector_field( (-x/sqrt(x^2+y^2), -y/sqrt(x^2+y^2)), (x, -10, 10), (y, -10, 10))

    ::

        sage: x,y = var('x,y')
        sage: plot_vector_field( (-x/sqrt(x+y), -y/sqrt(x+y)), (x, -10, 10), (y, -10, 10))

    Extra options will get passed on to show(), as long as they are valid::

        sage: plot_vector_field((x, y), (x, -2, 2), (y, -2, 2), xmax=10)
        sage: plot_vector_field((x, y), (x, -2, 2), (y, -2, 2)).show(xmax=10) # These are equivalent
    """
    from sage.plot.plot import Graphics
    from sage.plot.misc import setup_for_eval_on_grid
    z, ranges = setup_for_eval_on_grid([f,g], [xrange, yrange], options['plot_points'])
    f,g = z

    xpos_array, ypos_array, xvec_array, yvec_array = [],[],[],[]
    for x in xsrange(*ranges[0], include_endpoint=True):
        for y in xsrange(*ranges[1], include_endpoint=True):
            xpos_array.append(x)
            ypos_array.append(y)
            xvec_array.append(f(x,y))
            yvec_array.append(g(x,y))

    import numpy
    xvec_array = numpy.ma.masked_invalid(numpy.array(xvec_array, dtype=float))
    yvec_array = numpy.ma.masked_invalid(numpy.array(yvec_array, dtype=float))
    g = Graphics()
    g._set_extra_kwds(Graphics._extract_kwds_for_show(options))

def plot_vector_field3d(functions, xrange, yrange, zrange, plot_points=5, colors="jet", center_arrows=False, **kwds):
    r"""
    Plot a 3d vector field

    INPUT:

    - ``functions`` - a list of three functions, representing the x-,
      y-, and z-coordinates of a vector

    - ``xrange``, ``yrange``, and ``zrange`` - three tuples of the
      form (var, start, stop), giving the variables and ranges for each axis

    - ``plot_points`` (default 5) - either a number or list of three
      numbers, specifying how many points to plot for each axis

    - ``colors`` (default 'jet') - a color, list of colors (which are
      interpolated between), or matplotlib colormap name, giving the coloring
      of the arrows.  If a list of colors or a colormap is given,
      coloring is done as a function of length of the vector

    - ``center_arrows`` (default False) - If True, draw the arrows
      centered on the points; otherwise, draw the arrows with the tail
      at the point

    - any other keywords are passed on to the plot command for each arrow

    EXAMPLES::

        sage: x,y,z=var('x y z')
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi))
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),colors=['red','green','blue'])
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),colors='red')
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),plot_points=4)
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),plot_points=[3,5,7])
        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),center_arrows=True)

    TESTS:

    This tests that :trac:`2100` is fixed in a way compatible with this command::

        sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),center_arrows=True,aspect_ratio=(1,2,1))
    """
    (ff, gg, hh), ranges = setup_for_eval_on_grid(functions, [xrange, yrange, zrange], plot_points)
    xpoints, ypoints, zpoints = [srange(*r, include_endpoint=True) for r in ranges]
    points = [vector((i, j, k)) for i in xpoints for j in ypoints for k in zpoints]
    vectors = [vector((ff(*point), gg(*point), hh(*point))) for point in points]

    try:
        from matplotlib.cm import get_cmap

        cm = get_cmap(colors)
    except (TypeError, ValueError):
        cm = None
    if cm is None:
        if isinstance(colors, (list, tuple)):
            from matplotlib.colors import LinearSegmentedColormap

            cm = LinearSegmentedColormap.from_list("mymap", colors)
        else:
            cm = lambda x: colors

    max_len = max(v.norm() for v in vectors)
    scaled_vectors = [v / max_len for v in vectors]

    if center_arrows:
        return sum([plot(v, color=cm(v.norm()), **kwds).translate(p - v / 2) for v, p in zip(scaled_vectors, points)])
    else:
        return sum([plot(v, color=cm(v.norm()), **kwds).translate(p) for v, p in zip(scaled_vectors, points)])

def _parametric_plot3d_surface(f, urange, vrange, plot_points, boundary_style, **kwds):
    r"""
    Return a parametric three-dimensional space surface.
    This function is used internally by the
    :func:`parametric_plot3d` command.

    There are two ways this function is invoked by
    :func:`parametric_plot3d`.

    - ``parametric_plot3d([f_x, f_y, f_z], (u_min, u_max),
      (v_min, v_max))``:
      `f_x, f_y, f_z` are each functions of two variables

    - ``parametric_plot3d([f_x, f_y, f_z], (u, u_min,
      u_max), (v, v_min, v_max))``:
      `f_x, f_y, f_z` can be viewed as functions of
      `u` and `v`

    INPUT:

    - ``f`` - a 3-tuple of functions or expressions, or vector of size 3

    - ``urange`` - a 2-tuple (u_min, u_max) or a 3-tuple
      (u, u_min, u_max)

    - ``vrange`` - a 2-tuple (v_min, v_max) or a 3-tuple
      (v, v_min, v_max)

    - ``plot_points`` - (default: "automatic", which is [40,40]
      for surfaces) initial number of sample points in each parameter;
      a pair of integers.

    - ``boundary_style`` - (default: None, no boundary) a dict that describes
      how to draw the boundaries of regions by giving options that are passed
      to the line3d command.

    EXAMPLES:

    We demonstrate each of the two ways of calling this.  See
    :func:`parametric_plot3d` for many more examples.

    We do the first one with lambda functions::

        sage: f = (lambda u,v: cos(u), lambda u,v: sin(u)+cos(v), lambda u,v: sin(v))
        sage: parametric_plot3d(f, (0, 2*pi), (-pi, pi)) # indirect doctest
        Graphics3d Object

    Now we do the same thing with symbolic expressions::

        sage: u, v = var('u,v')
        sage: parametric_plot3d((cos(u), sin(u) + cos(v), sin(v)), (u, 0, 2*pi), (v, -pi, pi), mesh=True)
        Graphics3d Object
    """
    from sage.plot.misc import setup_for_eval_on_grid
    g, ranges = setup_for_eval_on_grid(f, [urange, vrange], plot_points)
    urange = srange(*ranges[0], include_endpoint=True)
    vrange = srange(*ranges[1], include_endpoint=True)
    G = ParametricSurface(g, (urange, vrange), **kwds)

    if boundary_style is not None:
        for u in (urange[0], urange[-1]):
            G += line3d([(g[0](u,v), g[1](u,v), g[2](u,v)) for v in vrange], **boundary_style)
        for v in (vrange[0], vrange[-1]):
            G += line3d([(g[0](u,v), g[1](u,v), g[2](u,v)) for u in urange], **boundary_style)
    return G

def plot3d_adaptive(f, x_range, y_range, color="automatic",
                    grad_f=None,
                    max_bend=.5, max_depth=5, initial_depth=4, num_colors=128, **kwds):
    r"""
    Adaptive 3d plotting of a function of two variables.

    This is used internally by the plot3d command when the option
    ``adaptive=True`` is given.

    INPUT:


    -  ``f`` - a symbolic function or a Python function of
       3 variables.

    -  ``x_range`` - x range of values: 2-tuple (xmin,
       xmax) or 3-tuple (x,xmin,xmax)

    -  ``y_range`` - y range of values: 2-tuple (ymin,
       ymax) or 3-tuple (y,ymin,ymax)

    -  ``grad_f`` - gradient of f as a Python function

    -  ``color`` - "automatic" - a rainbow of num_colors
       colors

    -  ``num_colors`` - (default: 128) number of colors to
       use with default color

    -  ``max_bend`` - (default: 0.5)

    -  ``max_depth`` - (default: 5)

    -  ``initial_depth`` - (default: 4)

    -  ``**kwds`` - standard graphics parameters


    EXAMPLES:

    We plot `\sin(xy)`::

        sage: from sage.plot.plot3d.plot3d import plot3d_adaptive
        sage: x,y=var('x,y'); plot3d_adaptive(sin(x*y), (x,-pi,pi), (y,-pi,pi), initial_depth=5)
        Graphics3d Object
        
    .. PLOT::
        
        from sage.plot.plot3d.plot3d import plot3d_adaptive
        x,y=var('x,y')
        sphinx_plot(plot3d_adaptive(sin(x*y), (x,-pi,pi), (y,-pi,pi), initial_depth=5))

    """
    if initial_depth >= max_depth:
        max_depth = initial_depth

    from sage.plot.misc import setup_for_eval_on_grid
    g, ranges = setup_for_eval_on_grid(f, [x_range,y_range], plot_points=2)
    xmin,xmax = ranges[0][:2]
    ymin,ymax = ranges[1][:2]

    opacity = kwds.get('opacity',1)

    if color == "automatic":
        texture = rainbow(num_colors, 'rgbtuple')
    else:
        if isinstance(color, list):
            texture = color
        else:
            kwds['color'] = color
            texture = Texture(kwds)

    factory = TrivialTriangleFactory()
    plot = TrianglePlot(factory, g, (xmin, xmax), (ymin, ymax), g = grad_f,
                        min_depth=initial_depth, max_depth=max_depth,
                        max_bend=max_bend, num_colors = None)

    P = IndexFaceSet(plot._objects)
    if isinstance(texture, (list, tuple)):
        if len(texture) == 2:
            # do a grid coloring
            xticks = (xmax - xmin)/2**initial_depth
            yticks = (ymax - ymin)/2**initial_depth
            parts = P.partition(lambda x,y,z: (int((x-xmin)/xticks) + int((y-ymin)/yticks)) % 2)
        else:
            # do a topo coloring
            bounds = P.bounding_box()
            min_z = bounds[0][2]
            max_z = bounds[1][2]
            if max_z == min_z:
                span = 0
            else:
                span = (len(texture)-1) / (max_z - min_z)    # max to avoid dividing by 0
            parts = P.partition(lambda x,y,z: int((z-min_z)*span))
        all = []
        for k, G in parts.iteritems():
            G.set_texture(texture[k], opacity=opacity)
            all.append(G)
        P = Graphics3dGroup(all)
    else:
#.........这里部分代码省略.........

def plot_vector_field(f_g, xrange, yrange, **options):
    r"""
    ``plot_vector_field`` takes two functions of two variables xvar and yvar
    (for instance, if the variables are `x` and `y`, take `(f(x,y), g(x,y))`)
    and plots vector arrows of the function over the specified ranges, with
    xrange being of xvar between xmin and xmax, and yrange similarly (see below).

    ``plot_vector_field((f, g), (xvar, xmin, xmax), (yvar, ymin, ymax))``

    EXAMPLES:

    Plot some vector fields involving sin and cos::

        sage: x,y = var('x y')
        sage: plot_vector_field((sin(x), cos(y)), (x,-3,3), (y,-3,3))
        Graphics object consisting of 1 graphics primitive

    ::

        sage: plot_vector_field(( y, (cos(x)-2)*sin(x)), (x,-pi,pi), (y,-pi,pi))
        Graphics object consisting of 1 graphics primitive

    Plot a gradient field::

        sage: u,v = var('u v')
        sage: f = exp(-(u^2+v^2))
        sage: plot_vector_field(f.gradient(), (u,-2,2), (v,-2,2), color='blue')
        Graphics object consisting of 1 graphics primitive

    Plot two orthogonal vector fields::

        sage: x,y = var('x,y')
        sage: a=plot_vector_field((x,y), (x,-3,3),(y,-3,3),color='blue')
        sage: b=plot_vector_field((y,-x),(x,-3,3),(y,-3,3),color='red')
        sage: show(a+b)

    We ignore function values that are infinite or NaN::

        sage: x,y = var('x,y')
        sage: plot_vector_field( (-x/sqrt(x^2+y^2), -y/sqrt(x^2+y^2)), (x, -10, 10), (y, -10, 10))
        Graphics object consisting of 1 graphics primitive

    ::

        sage: x,y = var('x,y')
        sage: plot_vector_field( (-x/sqrt(x+y), -y/sqrt(x+y)), (x, -10, 10), (y, -10, 10))
        Graphics object consisting of 1 graphics primitive

    Extra options will get passed on to show(), as long as they are valid::

        sage: plot_vector_field((x, y), (x, -2, 2), (y, -2, 2), xmax=10)
        Graphics object consisting of 1 graphics primitive
        sage: plot_vector_field((x, y), (x, -2, 2), (y, -2, 2)).show(xmax=10) # These are equivalent
    """
    (f, g) = f_g
    from sage.plot.all import Graphics
    from sage.plot.misc import setup_for_eval_on_grid
    z, ranges = setup_for_eval_on_grid([f,g], [xrange, yrange], options['plot_points'])
    f,g = z

    xpos_array, ypos_array, xvec_array, yvec_array = [],[],[],[]
    for x in xsrange(*ranges[0], include_endpoint=True):
        for y in xsrange(*ranges[1], include_endpoint=True):
            xpos_array.append(x)
            ypos_array.append(y)
            xvec_array.append(f(x,y))
            yvec_array.append(g(x,y))

    import numpy
    xvec_array = numpy.ma.masked_invalid(numpy.array(xvec_array, dtype=float))
    yvec_array = numpy.ma.masked_invalid(numpy.array(yvec_array, dtype=float))
    g = Graphics()
    g._set_extra_kwds(Graphics._extract_kwds_for_show(options))
    g.add_primitive(PlotField(xpos_array, ypos_array, xvec_array, yvec_array, options))
    return g

from sage.misc.decorators import options
@options(plot_points=20)
def plot_vector_field_on_curve( (xf, yf), (x, y), range, **options ):
    r"""Plot values of a vector-values function along points of a curve
    in the plane.

    Note this function doesn't plot the curve itself."""
    from sage.plot.all import Graphics
    from sage.misc.misc import xsrange
    from sage.plot.plot_field import PlotField
    from sage.plot.misc import setup_for_eval_on_grid
    zz, rangez = setup_for_eval_on_grid( (x, y, xf, yf), [ range ], options['plot_points'] )
    #print 'setup: ', zz, rangez
    x, y, xf, yf = zz
    xpos_array, ypos_array, xvec_array, yvec_array = [],[],[],[]
    for t in xsrange( *rangez[0], include_endpoint=True ):
       xpos_array.append( x(t) )
       ypos_array.append( y(t) )
       xvec_array.append( xf(t) )
       yvec_array.append( yf(t) )
    import numpy
    xvec_array = numpy.ma.masked_invalid(numpy.array(xvec_array, dtype=float))
    yvec_array = numpy.ma.masked_invalid(numpy.array(yvec_array, dtype=float))
    g = Graphics()
    g._set_extra_kwds(Graphics._extract_kwds_for_show(options))
    g.add_primitive(PlotField(xpos_array, ypos_array, xvec_array, yvec_array, options))
    return g

#.........这里部分代码省略.........
        sage: region_plot(s>0,(t,-2,2),(s,-2,2))
        Graphics object consisting of 1 graphics primitive

    ::

        sage: region_plot(s>0,(s,-2,2),(t,-2,2))
        Graphics object consisting of 1 graphics primitive

    An example of a region plot in 'loglog' scale::

        sage: region_plot(x^2+y^2<100, (x,1,10), (y,1,10), scale='loglog')
        Graphics object consisting of 1 graphics primitive

    TESTS:

    To check that :trac:`16907` is fixed::

        sage: x, y = var('x, y')
        sage: disc1 = region_plot(x^2+y^2 < 1, (x, -1, 1), (y, -1, 1), alpha=0.5)
        sage: disc2 = region_plot((x-0.7)^2+(y-0.7)^2 < 0.5, (x, -2, 2), (y, -2, 2), incol='red', alpha=0.5)
        sage: disc1 + disc2
        Graphics object consisting of 2 graphics primitives

    To check that :trac:`18286` is fixed::
        sage: x, y = var('x, y')
        sage: region_plot([x == 0], (x, -1, 1), (y, -1, 1))
        Graphics object consisting of 1 graphics primitive
        sage: region_plot([x^2+y^2==1, x<y], (x, -1, 1), (y, -1, 1))
        Graphics object consisting of 1 graphics primitive

    """

    from sage.plot.all import Graphics
    from sage.plot.misc import setup_for_eval_on_grid
    from sage.symbolic.expression import is_Expression
    from warnings import warn
    import numpy

    if not isinstance(f, (list, tuple)):
        f = [f]

    feqs = [equify(g) for g in f if is_Expression(g) and g.operator() is operator.eq and not equify(g).is_zero()]
    f = [equify(g) for g in f if not (is_Expression(g) and g.operator() is operator.eq)]
    neqs = len(feqs)
    if neqs > 1:
        warn("There are at least 2 equations; If the region is degenerated to points, plotting might show nothing.")
        feqs = [sum([fn**2 for fn in feqs])]
        neqs = 1
    if neqs and not bordercol:
        bordercol = incol
    if not f:
        return implicit_plot(feqs[0], xrange, yrange, plot_points=plot_points, fill=False, \
                             linewidth=borderwidth, linestyle=borderstyle, color=bordercol, **options)
    f_all, ranges = setup_for_eval_on_grid(feqs + f, [xrange, yrange], plot_points)
    xrange,yrange=[r[:2] for r in ranges]

    xy_data_arrays = numpy.asarray([[[func(x, y) for x in xsrange(*ranges[0], include_endpoint=True)]
                                     for y in xsrange(*ranges[1], include_endpoint=True)]
                                    for func in f_all[neqs::]],dtype=float)
    xy_data_array=numpy.abs(xy_data_arrays.prod(axis=0))
    # Now we need to set entries to negative iff all
    # functions were negative at that point.
    neg_indices = (xy_data_arrays<0).all(axis=0)
    xy_data_array[neg_indices]=-xy_data_array[neg_indices]

    from matplotlib.colors import ListedColormap

#.........这里部分代码省略.........

        sage: streamline_plot((x, y), (x,-2,2), (y,-2,2), xmax=10)
        Graphics object consisting of 1 graphics primitive
        sage: streamline_plot((x, y), (x,-2,2), (y,-2,2)).show(xmax=10) # These are equivalent

    .. PLOT::

        x, y = var('x y')
        g = streamline_plot((x, y), (x,-2,2), (y,-2,2), xmax=10)
        sphinx_plot(g)

    We can also construct streamlines in a slope field::

        sage: x, y = var('x y')
        sage: streamline_plot((x + y) / sqrt(x^2 + y^2), (x,-3,3), (y,-3,3))
        Graphics object consisting of 1 graphics primitive

    .. PLOT::

        x, y = var('x y')
        g = streamline_plot((x + y) / sqrt(x**2 + y**2), (x,-3,3), (y,-3,3))
        sphinx_plot(g)

    We choose some particular points the streamlines pass through::

        sage: pts = [[1, 1], [-2, 2], [1, -3/2]]
        sage: g = streamline_plot((x + y) / sqrt(x^2 + y^2), (x,-3,3), (y,-3,3), start_points=pts)
        sage: g += point(pts, color='red')
        sage: g
        Graphics object consisting of 2 graphics primitives

    .. PLOT::

        x, y = var('x y')
        pts = [[1, 1], [-2, 2], [1, -3/2]]
        g = streamline_plot((x + y) / sqrt(x**2 + y**2), (x,-3,3), (y,-3,3), start_points=pts)
        g += point(pts, color='red')
        sphinx_plot(g)

    .. NOTE::

        Streamlines currently pass close to ``start_points`` but do
        not necessarily pass directly through them. That is part of
        the behavior of matplotlib, not an error on your part.

    """
    # Parse the function input
    if isinstance(f_g, (list, tuple)):
        (f,g) = f_g
    else:
        from sage.functions.all import sqrt
        from inspect import isfunction
        if isfunction(f_g):
            f = lambda x,y: 1 / sqrt(f_g(x, y)**2 + 1)
            g = lambda x,y: f_g(x, y) * f(x, y)
        else:
            f = 1 / sqrt(f_g**2 + 1)
            g = f_g * f

    from sage.plot.all import Graphics
    from sage.plot.misc import setup_for_eval_on_grid
    z, ranges = setup_for_eval_on_grid([f,g], [xrange,yrange], options['plot_points'])
    f, g = z

    # The density values must be floats
    if isinstance(options['density'], (list, tuple)):
        options['density'] = [float(x) for x in options['density']]
    else:
        options['density'] = float(options['density'])

    xpos_array, ypos_array, xvec_array, yvec_array = [], [], [], []
    for x in xsrange(*ranges[0], include_endpoint=True):
        xpos_array.append(x)
    for y in xsrange(*ranges[1], include_endpoint=True):
        ypos_array.append(y)
        xvec_row, yvec_row = [], []
        for x in xsrange(*ranges[0], include_endpoint=True):
            xvec_row.append(f(x, y))
            yvec_row.append(g(x, y))
        xvec_array.append(xvec_row)
        yvec_array.append(yvec_row)

    import numpy
    xpos_array = numpy.array(xpos_array, dtype=float)
    ypos_array = numpy.array(ypos_array, dtype=float)
    xvec_array = numpy.ma.masked_invalid(numpy.array(xvec_array, dtype=float))
    yvec_array = numpy.ma.masked_invalid(numpy.array(yvec_array, dtype=float))

    if 'start_points' in options:
        xstart_array, ystart_array = [], []
        for point in options['start_points']:
            xstart_array.append(point[0])
            ystart_array.append(point[1])
        options['start_points'] = numpy.array([xstart_array, ystart_array]).T

    g = Graphics()
    g._set_extra_kwds(Graphics._extract_kwds_for_show(options))
    g.add_primitive(StreamlinePlot(xpos_array, ypos_array,
                                   xvec_array, yvec_array, options))
    return g