Ticket #4102: trac_symbolic_bessel_deprecation.patch

sage/functions/bessel.py

@@ -172 +172 @@
 #                  http://www.gnu.org/licenses/
 #*****************************************************************************
 
-from sage.rings.all import RR, Integer
-from sage.misc.latex import latex
-from sage.symbolic.function import BuiltinFunction, is_inexact
-from sage.symbolic.expression import Expression
-from sage.structure.coerce import parent
-import sage.structure.element
-from sage.libs.mpmath import utils as mpmath_utils
 from sage.functions.other import sqrt
 from sage.functions.log import exp
 from sage.functions.hyperbolic import sinh, cosh
+from sage.libs.mpmath import utils as mpmath_utils
+from sage.misc.latex import latex
+from sage.rings.all import RR, Integer
+from sage.structure.coerce import parent
+from sage.structure.element import get_coercion_model
 from sage.symbolic.constants import pi
+from sage.symbolic.function import BuiltinFunction, is_inexact
+from sage.symbolic.expression import Expression
+
+# remove after deprecation period
+from sage.calculus.calculus import maxima
+from sage.functions.other import real, imag
+from sage.functions.special import maxima_function
+from sage.misc.sage_eval import sage_eval
+from sage.rings.real_mpfr import RealField
+from sage.plot.plot import plot
+from sage.rings.all import ZZ
 
 
 class Function_Bessel_J(BuiltinFunction):
@@ -283 +292 @@
                                  conversions=dict(maxima='bessel_j',
                                                   mathematica='BesselJ'))
 
+    # remove after deprecation period
+    def __call__(self, *args, **kwds):
+        """
+        Custom `__call__` method which uses the old Bessel function code if the
+        `algorithm` or `prec` arguments are used. This should be removed after
+        the deprecation period.
+
+        EXAMPLES::
+
+            sage: bessel_J(0, 1.0, "maxima", 53)
+            doctest:1: DeprecationWarning: precision argument is deprecated; algorithm argument is currently deprecated, but will be available as a named keyword in the future
+            See http://trac.sagemath.org/4102 for details.
+            0.765197686558
+        """
+        if len(args) > 2 or len(kwds) > 0:
+            from sage.misc.superseded import deprecation
+            deprecation(4102, 'precision argument is deprecated; algorithm '
+                              'argument is currently deprecated, but will be '
+                              'available as a named keyword in the future')
+            return _bessel_J(*args, **kwds)
+        else:
+            return super(BuiltinFunction, self).__call__(*args, **kwds)
+
     def _eval_(self, n, x):
         """
         EXAMPLES::
@@ -296 +328 @@
         if (not isinstance(n, Expression) and
                 not isinstance(x, Expression) and
                 (is_inexact(n) or is_inexact(x))):
-            coercion_model = sage.structure.element.get_coercion_model()
+            coercion_model = get_coercion_model()
             n, x = coercion_model.canonical_coercion(n, x)
             return self._evalf_(n, x, parent(n))
 
@@ -420 +452 @@
                                  conversions=dict(maxima='bessel_y',
                                                   mathematica='BesselY'))
 
+    # remove after deprecation period
+    def __call__(self, *args, **kwds):
+        """
+        Custom `__call__` method which uses the old Bessel function code if the
+        `algorithm` or `prec` arguments are used. This should be removed after
+        the deprecation period.
+
+        EXAMPLES::
+
+            sage: bessel_Y(0, 1, "maxima", 53)
+            doctest:1: DeprecationWarning: precision argument is deprecated; algorithm argument is currently deprecated, but will be available as a named keyword in the future
+            See http://trac.sagemath.org/4102 for details.
+            0.0882569642156769
+        """
+        if len(args) > 2 or len(kwds) > 0:
+            from sage.misc.superseded import deprecation
+            deprecation(4102, 'precision argument is deprecated; algorithm '
+                              'argument is currently deprecated, but will be '
+                              'available as a named keyword in the future')
+            return _bessel_Y(*args, **kwds)
+        else:
+            return super(BuiltinFunction, self).__call__(*args, **kwds)
+
     def _eval_(self, n, x):
         """
         EXAMPLES::
@@ -432 +487 @@
         """
         if (not isinstance(n, Expression) and not isinstance(x, Expression) and
                 (is_inexact(n) or is_inexact(x))):
-            coercion_model = sage.structure.element.get_coercion_model()
+            coercion_model = get_coercion_model()
             n, x = coercion_model.canonical_coercion(n, x)
             return self._evalf_(n, x, parent(n))
 
@@ -560 +615 @@
                                  conversions=dict(maxima='bessel_i',
                                                   mathematica='BesselI'))
 
+    # remove after deprecation period
+    def __call__(self, *args, **kwds):
+        """
+        Custom `__call__` method which uses the old Bessel function code if the
+        `algorithm` or `prec` arguments are used. This should be removed after
+        the deprecation period.
+
+        EXAMPLES::
+
+            sage: bessel_I(0, 1, "maxima", 53)
+            doctest:1: DeprecationWarning: precision argument is deprecated; algorithm argument is currently deprecated, but will be available as a named keyword in the future
+            See http://trac.sagemath.org/4102 for details.
+            1.266065877752009
+        """
+        if len(args) > 2 or len(kwds) > 0:
+            from sage.misc.superseded import deprecation
+            deprecation(4102, 'precision argument is deprecated; algorithm '
+                              'argument is currently deprecated, but will be '
+                              'available as a named keyword in the future')
+            return _bessel_I(*args, **kwds)
+        else:
+            return super(BuiltinFunction, self).__call__(*args, **kwds)
+
     def _eval_(self, n, x):
         """
         EXAMPLES::
@@ -576 +654 @@
         """
         if (not isinstance(n, Expression) and not isinstance(x, Expression) and
                 (is_inexact(n) or is_inexact(x))):
-            coercion_model = sage.structure.element.get_coercion_model()
+            coercion_model = get_coercion_model()
             n, x = coercion_model.canonical_coercion(n, x)
             return self._evalf_(n, x, parent(n))
 
@@ -725 +803 @@
                                  conversions=dict(maxima='bessel_k',
                                                   mathematica='BesselK'))
 
+    # remove after deprecation period
+    def __call__(self, *args, **kwds):
+        """
+        Custom `__call__` method which uses the old Bessel function code if the
+        `algorithm` or `prec` arguments are used. This should be removed after
+        the deprecation period.
+
+        EXAMPLES::
+
+            sage: bessel_K(0, 1, "maxima", 53)
+            doctest:1: DeprecationWarning: precision argument is deprecated; algorithm argument is currently deprecated, but will be available as a named keyword in the future
+            See http://trac.sagemath.org/4102 for details.
+            0.0882569642156769
+        """
+        if len(args) > 2 or len(kwds) > 0:
+            from sage.misc.superseded import deprecation
+            deprecation(4102, 'precision argument is deprecated; algorithm '
+                              'argument is currently deprecated, but will be '
+                              'available as a named keyword in the future')
+            return _bessel_Y(*args, **kwds)
+        else:
+            return super(BuiltinFunction, self).__call__(*args, **kwds)
+
     def _eval_(self, n, x):
         """
         EXAMPLES::
@@ -738 +839 @@
         """
         if (not isinstance(n, Expression) and not isinstance(x, Expression) and
                 (is_inexact(n) or is_inexact(x))):
-            coercion_model = sage.structure.element.get_coercion_model()
+            coercion_model = get_coercion_model()
             n, x = coercion_model.canonical_coercion(n, x)
             return self._evalf_(n, x, parent(n))
 
@@ -925 +1026 @@
         _type = args[1]
         _nargs = 1
     else:
-        raise TypeError("at most two positional arguments may be specified, "
-                        + "see the docstring for Bessel")
+        from sage.misc.superseded import deprecation
+        deprecation(4102, 'precision argument is deprecated; algorithm '
+                          'argument is currently deprecated, but will be '
+                          'available as a named keyword in the future')
+        return _Bessel(*args, **kwds)
+
     # check for type inconsistency
     if _type is not None and 'typ' in kwds and _type != kwds['typ']:
         raise ValueError("inconsistent types given")
@@ -957 +1062 @@
         return lambda x: _f(_order, x)
     else:
         return _f
+
+####################################################
+###  to be removed after the deprecation period  ###
+####################################################
+
+
+def _bessel_I(nu,z,algorithm = "pari",prec=53):
+    r"""
+    Implements the "I-Bessel function", or "modified Bessel function,
+    1st kind", with index (or "order") nu and argument z.
+
+    INPUT:
+
+
+    -  ``nu`` - a real (or complex, for pari) number
+
+    -  ``z`` - a real (positive) algorithm - "pari" or
+       "maxima" or "scipy" prec - real precision (for PARI only)
+
+
+    DEFINITION::
+
+            Maxima:
+                             inf
+                            ====   - nu - 2 k  nu + 2 k
+                            \     2          z
+                             >    -------------------
+                            /     k! Gamma(nu + k + 1)
+                            ====
+                            k = 0
+
+            PARI:
+
+                             inf
+                            ====   - 2 k  2 k
+                            \     2      z    Gamma(nu + 1)
+                             >    -----------------------
+                            /       k! Gamma(nu + k + 1)
+                            ====
+                            k = 0
+
+
+
+    Sometimes ``bessel_I(nu,z)`` is denoted
+    ``I_nu(z)`` in the literature.
+
+    .. warning::
+
+       In Maxima (the manual says) i0 is deprecated but
+       ``bessel_i(0,*)`` is broken. (Was fixed in recent CVS patch
+       though.)
+
+    EXAMPLES::
+
+        sage: from sage.functions.bessel import _bessel_I
+        sage: _bessel_I(1,1,"pari",500)
+        0.565159103992485027207696027609863307328899621621092009480294489479255640964371134092664997766814410064677886055526302676857637684917179812041131208121
+        sage: _bessel_I(1,1)
+        0.565159103992485
+        sage: _bessel_I(2,1.1,"maxima")
+        0.16708949925104...
+        sage: _bessel_I(0,1.1,"maxima")
+        1.32616018371265...
+        sage: _bessel_I(0,1,"maxima")
+        1.2660658777520...
+        sage: _bessel_I(1,1,"scipy")
+        0.565159103992...
+
+    Check whether the return value is real whenever the argument is real (#10251)::
+
+        sage: _bessel_I(5, 1.5, algorithm='scipy') in RR
+        True
+
+    """
+    if algorithm=="pari":
+        from sage.libs.pari.all import pari
+        try:
+            R = RealField(prec)
+            nu = R(nu)
+            z = R(z)
+        except TypeError:
+            C = ComplexField(prec)
+            nu = C(nu)
+            z = C(z)
+            K = C
+        K = z.parent()
+        return K(pari(nu).besseli(z, precision=prec))
+    elif algorithm=="scipy":
+        if prec != 53:
+            raise ValueError, "for the scipy algorithm the precision must be 53"
+        import scipy.special
+        ans = str(scipy.special.iv(float(nu),complex(real(z),imag(z))))
+        ans = ans.replace("(","")
+        ans = ans.replace(")","")
+        ans = ans.replace("j","*I")
+        ans = sage_eval(ans)
+        return real(ans) if z in RR else ans # Return real value when arg is real
+    elif algorithm == "maxima":
+        if prec != 53:
+            raise ValueError, "for the maxima algorithm the precision must be 53"
+        return sage_eval(maxima.eval("bessel_i(%s,%s)"%(float(nu),float(z))))
+    else:
+        raise ValueError, "unknown algorithm '%s'"%algorithm
+
+def _bessel_J(nu,z,algorithm="pari",prec=53):
+    r"""
+    Return value of the "J-Bessel function", or "Bessel function, 1st
+    kind", with index (or "order") nu and argument z.
+
+    ::
+
+            Defn:
+            Maxima:
+                             inf
+                            ====          - nu - 2 k  nu + 2 k
+                            \     (-1)^k 2           z
+                             >    -------------------------
+                            /        k! Gamma(nu + k + 1)
+                            ====
+                            k = 0
+
+            PARI:
+
+                             inf
+                            ====          - 2k    2k
+                            \     (-1)^k 2      z    Gamma(nu + 1)
+                             >    ----------------------------
+                            /         k! Gamma(nu + k + 1)
+                            ====
+                            k = 0
+
+
+    Sometimes bessel_J(nu,z) is denoted J_nu(z) in the literature.
+
+    .. warning::
+
+       Inaccurate for small values of z.
+
+    EXAMPLES::
+
+        sage: from sage.functions.bessel import _bessel_J
+        sage: _bessel_J(2,1.1)
+        0.136564153956658
+        sage: _bessel_J(0,1.1)
+        0.719622018527511
+        sage: _bessel_J(0,1)
+        0.765197686557967
+        sage: _bessel_J(0,0)
+        1.00000000000000
+        sage: _bessel_J(0.1,0.1)
+        0.777264368097005
+
+    We check consistency of PARI and Maxima::
+
+        sage: n(_bessel_J(3,10,"maxima"))
+        0.0583793793051...
+        sage: n(_bessel_J(3,10,"pari"))
+        0.0583793793051868
+        sage: _bessel_J(3,10,"scipy")
+        0.0583793793052...
+
+    Check whether the return value is real whenever the argument is real (#10251)::
+        sage: _bessel_J(5, 1.5, algorithm='scipy') in RR
+        True
+    """
+
+    if algorithm=="pari":
+        from sage.libs.pari.all import pari
+        try:
+            R = RealField(prec)
+            nu = R(nu)
+            z = R(z)
+        except TypeError:
+            C = ComplexField(prec)
+            nu = C(nu)
+            z = C(z)
+            K = C
+        if nu == 0:
+            nu = ZZ(0)
+        K = z.parent()
+        return K(pari(nu).besselj(z, precision=prec))
+    elif algorithm=="scipy":
+        if prec != 53:
+            raise ValueError, "for the scipy algorithm the precision must be 53"
+        import scipy.special
+        ans = str(scipy.special.jv(float(nu),complex(real(z),imag(z))))
+        ans = ans.replace("(","")
+        ans = ans.replace(")","")
+        ans = ans.replace("j","*I")
+        ans = sage_eval(ans)
+        return real(ans) if z in RR else ans
+    elif algorithm == "maxima":
+        if prec != 53:
+            raise ValueError, "for the maxima algorithm the precision must be 53"
+        return maxima_function("bessel_j")(nu, z)
+    else:
+        raise ValueError, "unknown algorithm '%s'"%algorithm
+
+def _bessel_K(nu,z,algorithm="pari",prec=53):
+    r"""
+    Implements the "K-Bessel function", or "modified Bessel function,
+    2nd kind", with index (or "order") nu and argument z. Defn::
+
+                    pi*(bessel_I(-nu, z) - bessel_I(nu, z))
+                   ----------------------------------------
+                                2*sin(pi*nu)
+
+
+    if nu is not an integer and by taking a limit otherwise.
+
+    Sometimes bessel_K(nu,z) is denoted K_nu(z) in the literature. In
+    PARI, nu can be complex and z must be real and positive.
+
+    EXAMPLES::
+
+        sage: from sage.functions.bessel import _bessel_K
+        sage: _bessel_K(3,2,"scipy")
+        0.64738539094...
+        sage: _bessel_K(3,2)
+        0.64738539094...
+        sage: _bessel_K(1,1)
+        0.60190723019...
+        sage: _bessel_K(1,1,"pari",10)
+        0.60
+        sage: _bessel_K(1,1,"pari",100)
+        0.60190723019723457473754000154
+
+    TESTS::
+
+        sage: _bessel_K(2,1.1, algorithm="maxima")
+        Traceback (most recent call last):
+        ...
+        NotImplementedError: The K-Bessel function is only implemented for the pari and scipy algorithms
+
+        Check whether the return value is real whenever the argument is real (#10251)::
+
+        sage: _bessel_K(5, 1.5, algorithm='scipy') in RR
+        True
+
+    """
+    if algorithm=="scipy":
+        if prec != 53:
+            raise ValueError, "for the scipy algorithm the precision must be 53"
+        import scipy.special
+        ans = str(scipy.special.kv(float(nu),float(z)))
+        ans = ans.replace("(","")
+        ans = ans.replace(")","")
+        ans = ans.replace("j","*I")
+        ans = sage_eval(ans)
+        return real(ans) if z in RR else ans
+    elif algorithm == 'pari':
+        from sage.libs.pari.all import pari
+        try:
+            R = RealField(prec)
+            nu = R(nu)
+            z = R(z)
+        except TypeError:
+            C = ComplexField(prec)
+            nu = C(nu)
+            z = C(z)
+            K = C
+        K = z.parent()
+        return K(pari(nu).besselk(z, precision=prec))
+    elif algorithm == 'maxima':
+        raise NotImplementedError, "The K-Bessel function is only implemented for the pari and scipy algorithms"
+    else:
+        raise ValueError, "unknown algorithm '%s'"%algorithm
+
+
+def _bessel_Y(nu,z,algorithm="maxima", prec=53):
+    r"""
+    Implements the "Y-Bessel function", or "Bessel function of the 2nd
+    kind", with index (or "order") nu and argument z.
+
+    .. note::
+
+       Currently only prec=53 is supported.
+
+    Defn::
+
+                    cos(pi n)*bessel_J(nu, z) - bessel_J(-nu, z)
+                   -------------------------------------------------
+                                     sin(nu*pi)
+
+    if nu is not an integer and by taking a limit otherwise.
+
+    Sometimes bessel_Y(n,z) is denoted Y_n(z) in the literature.
+
+    This is computed using Maxima by default.
+
+    EXAMPLES::
+
+        sage: from sage.functions.bessel import _bessel_Y
+        sage: _bessel_Y(2,1.1,"scipy")
+        -1.4314714939...
+        sage: _bessel_Y(2,1.1)
+        -1.4314714939590...
+        sage: _bessel_Y(3.001,2.1)
+        -1.0299574976424...
+
+    TESTS::
+
+        sage: _bessel_Y(2,1.1, algorithm="pari")
+        Traceback (most recent call last):
+        ...
+        NotImplementedError: The Y-Bessel function is only implemented for the maxima and scipy algorithms
+    """
+    if algorithm=="scipy":
+        if prec != 53:
+            raise ValueError, "for the scipy algorithm the precision must be 53"
+        import scipy.special
+        ans = str(scipy.special.yv(float(nu),complex(real(z),imag(z))))
+        ans = ans.replace("(","")
+        ans = ans.replace(")","")
+        ans = ans.replace("j","*I")
+        ans = sage_eval(ans)
+        return real(ans) if z in RR else ans
+    elif algorithm == "maxima":
+        if prec != 53:
+            raise ValueError, "for the maxima algorithm the precision must be 53"
+        return RR(maxima.eval("bessel_y(%s,%s)"%(float(nu),float(z))))
+    elif algorithm == "pari":
+        raise NotImplementedError, "The Y-Bessel function is only implemented for the maxima and scipy algorithms"
+    else:
+        raise ValueError, "unknown algorithm '%s'"%algorithm
+
+class _Bessel():
+    """
+    A class implementing the I, J, K, and Y Bessel functions.
+
+    EXAMPLES::
+
+        sage: from sage.functions.bessel import _Bessel
+        sage: g = _Bessel(2); g
+        J_{2}
+        sage: print g
+        J-Bessel function of order 2
+        sage: g.plot(0,10)
+
+    ::
+
+        sage: _Bessel(2, typ='I')(pi)
+        2.61849485263445
+        sage: _Bessel(2, typ='J')(pi)
+        0.485433932631509
+        sage: _Bessel(2, typ='K')(pi)
+        0.0510986902537926
+        sage: _Bessel(2, typ='Y')(pi)
+        -0.0999007139289404
+    """
+    def __init__(self, nu, typ = "J", algorithm = None, prec = 53):
+        """
+        Initializes new instance of the Bessel class.
+
+        INPUT:
+
+         - ``typ`` -- (default: J) the type of Bessel function: 'I', 'J', 'K'
+           or 'Y'.
+
+         - ``algorithm`` -- (default: maxima for type Y, pari for other types)
+           algorithm to use to compute the Bessel function: 'pari', 'maxima' or
+           'scipy'.  Note that type K is not implemented in Maxima and type Y
+           is not implemented in PARI.
+
+         - ``prec`` -- (default: 53) precision in bits of the Bessel function.
+           Only supported for the PARI algorithm.
+
+        EXAMPLES::
+
+            sage: from sage.functions.bessel import _Bessel
+            sage: g = _Bessel(2); g
+            J_{2}
+            sage: _Bessel(1,'I')
+            I_{1}
+            sage: _Bessel(6, prec=120)(pi)
+            0.014545966982505560573660369604001804
+            sage: _Bessel(6, algorithm="pari")(pi)
+            0.0145459669825056
+
+        For the Bessel J-function, Maxima returns a symbolic result.  For
+        types I and Y, we always get a numeric result::
+
+            sage: b = _Bessel(6, algorithm="maxima")(pi); b
+            bessel_j(6, pi)
+            sage: b.n(53)
+            0.0145459669825056
+            sage: _Bessel(6, typ='I', algorithm="maxima")(pi)
+            0.0294619840059568
+            sage: _Bessel(6, typ='Y', algorithm="maxima")(pi)
+            -4.33932818939038
+
+        SciPy usually gives less precise results::
+
+            sage: _Bessel(6, algorithm="scipy")(pi)
+            0.0145459669825000...
+
+        TESTS::
+
+            sage: _Bessel(1,'Z')
+            Traceback (most recent call last):
+            ...
+            ValueError: typ must be one of I, J, K, Y
+        """
+        if not (typ in ['I', 'J', 'K', 'Y']):
+            raise ValueError, "typ must be one of I, J, K, Y"
+
+        # Did the user ask for the default algorithm?
+        if algorithm is None:
+            if typ == 'Y':
+                algorithm = 'maxima'
+            else:
+                algorithm = 'pari'
+
+        self._system = algorithm
+        self._order = nu
+        self._type = typ
+        prec = int(prec)
+        if prec < 0:
+            raise ValueError, "prec must be a positive integer"
+        self._prec = int(prec)
+
+    def __str__(self):
+        """
+        Returns a string representation of this Bessel object.
+
+        TEST::
+
+            sage: from sage.functions.bessel import _Bessel
+            sage: a = _Bessel(1,'I')
+            sage: str(a)
+            'I-Bessel function of order 1'
+        """
+        return self.type()+"-Bessel function of order "+str(self.order())
+
+    def __repr__(self):
+        """
+        Returns a string representation of this Bessel object.
+
+        TESTS::
+
+            sage: from sage.functions.bessel import _Bessel
+            sage: _Bessel(1,'I')
+            I_{1}
+        """
+        return self.type()+"_{"+str(self.order())+"}"
+
+    def type(self):
+        """
+        Returns the type of this Bessel object.
+
+        TEST::
+
+            sage: from sage.functions.bessel import _Bessel
+            sage: a = _Bessel(3,'K')
+            sage: a.type()
+            'K'
+        """
+        return self._type
+
+    def prec(self):
+        """
+        Returns the precision (in number of bits) used to represent this
+        Bessel function.
+
+        TESTS::
+
+            sage: from sage.functions.bessel import _Bessel
+            sage: a = _Bessel(3,'K')
+            sage: a.prec()
+            53
+            sage: B = _Bessel(20,prec=100); B
+            J_{20}
+            sage: B.prec()
+            100
+        """
+        return self._prec
+
+    def order(self):
+        """
+        Returns the order of this Bessel function.
+
+        TEST::
+
+            sage: from sage.functions.bessel import _Bessel
+            sage: a = _Bessel(3,'K')
+            sage: a.order()
+            3
+        """
+        return self._order
+
+    def system(self):
+        """
+        Returns the package used, e.g. Maxima, PARI, or SciPy, to compute with
+        this Bessel function.
+
+        TESTS::
+
+            sage: from sage.functions.bessel import _Bessel
+            sage: _Bessel(20,algorithm='maxima').system()
+            'maxima'
+            sage: _Bessel(20,prec=100).system()
+            'pari'
+        """
+        return self._system
+
+    def __call__(self,z):
+        """
+        Implements evaluation of all the Bessel functions directly
+        from the Bessel class. This essentially allows a Bessel object to
+        behave like a function that can be invoked.
+
+        TESTS::
+
+            sage: from sage.functions.bessel import _Bessel
+            sage: _Bessel(3,'K')(5.0)
+            0.00829176841523093
+            sage: _Bessel(20,algorithm='maxima')(5.0)
+            2.77033005213e-11
+            sage: _Bessel(20,prec=100)(5.0101010101010101)
+            2.8809188227195382093062257967e-11
+            sage: B = _Bessel(2,'Y',algorithm='scipy',prec=50)
+            sage: B(2.0)
+            Traceback (most recent call last):
+            ...
+            ValueError: for the scipy algorithm the precision must be 53
+        """
+        nu = self.order()
+        t = self.type()
+        s = self.system()
+        p = self.prec()
+        if t == "I":
+            return _bessel_I(nu,z,algorithm=s,prec=p)
+        if t == "J":
+            return _bessel_J(nu,z,algorithm=s,prec=p)
+        if t == "K":
+            return _bessel_K(nu,z,algorithm=s,prec=p)
+        if t == "Y":
+            return _bessel_Y(nu,z,algorithm=s,prec=p)
+
+    def plot(self,a,b):
+        """
+        Enables easy plotting of all the Bessel functions directly
+        from the Bessel class.
+
+        TESTS::
+
+            sage: from sage.functions.bessel import _Bessel
+            sage: plot(_Bessel(2),3,4)
+            sage: _Bessel(2).plot(3,4)
+            sage: P = _Bessel(2,'I').plot(1,5)
+            sage: P += _Bessel(2,'J').plot(1,5)
+            sage: P += _Bessel(2,'K').plot(1,5)
+            sage: P += _Bessel(2,'Y').plot(1,5)
+            sage: show(P)
+        """
+        nu = self.order()
+        s = self.system()
+        t = self.type()
+        if t == "I":
+            f = lambda z: _bessel_I(nu,z,s)
+            P = plot(f,a,b)
+        if t == "J":
+            f = lambda z: _bessel_J(nu,z,s)
+            P = plot(f,a,b)
+        if t == "K":
+            f = lambda z: _bessel_K(nu,z,s)
+            P = plot(f,a,b)
+        if t == "Y":
+            f = lambda z: _bessel_Y(nu,z,s)
+            P = plot(f,a,b)
+        return P