Python other.imag函数代码示例

    def cmp_ir(self,z):
        """ 
        returns -1 for left, 0 for in, and 1 for right from initial region
        cut line is on the north ray from pfp.

        Works only for real x0.
        """
        pfp = self.pfp
        x0 = self.x0
        if x0 > 0.5:
            print z,abs(z)
            if real(z) >= real(pfp) and abs(z) < abs(pfp):
                return 0
            if real(z) < real(pfp):
                return -1
            if real(z) > real(pfp):
                return 1
        else:
            if imag(z) > imag(pfp):
                if real(z) > real(pfp):
                    return 1
                if real(z) < real(pfp):
                    return -1
            if real(z) < real(pfp) and real(z) > log(real(pfp)) + log(sqrt(1+tan(imag(z))**2)):
                return 0
            if real(z) > real(pfp):
                return 1
            if real(z) < real(pfp):
                return -1

 def cmp_ir(self,z):
     """ 
     returns -1 for left, 0 for in, and 1 for right from initial region
     cut line is on the north ray from L.
     """
     L = self.L
     x0 = self.x0
     if x0 > 0.5:
         if real(z) > real(L) and abs(z) < abs(L):
             return 0
         if real(z) < real(L):
             return -1
         if real(z) > real(L):
             return 1
     else:
         if imag(z) > imag(L):
             if real(z) > real(L):
                 return 1
             if real(z) < real(L):
                 return -1
         if real(z) < real(L) and real(z) > log(real(L)) + log(sqrt(1+tan(imag(z))**2)):
             return 0
         if real(z) > real(L):
             return 1
         if real(z) < real(L):
             return -1

    def show(self, boundary=True, **options):
        r"""
        Plot ``self``.

        EXAMPLES::

            sage: HyperbolicPlane().UHP().get_geodesic(0, 1).show()
            Graphics object consisting of 2 graphics primitives
        """
        opts = {'axes': False, 'aspect_ratio': 1}
        opts.update(self.graphics_options())
        opts.update(options)
        end_1, end_2 = [CC(k.coordinates()) for k in self.endpoints()]
        bd_1, bd_2 = [CC(k.coordinates()) for k in self.ideal_endpoints()]
        if (abs(real(end_1) - real(end_2)) < EPSILON) \
                or CC(infinity) in [end_1, end_2]: #on same vertical line
            # If one of the endpoints is infinity, we replace it with a
            # large finite  point
            if end_1 == CC(infinity):
                end_1 = (real(end_2), (imag(end_2) + 10))
                end_2 = (real(end_2), imag(end_2))
            elif end_2 == CC(infinity):
                end_2 = (real(end_1), (imag(end_1) + 10))
                end_1 = (real(end_1), imag(end_1))
            from sage.plot.line import line
            pic = line((end_1, end_2), **opts)
            if boundary:
                cent = min(bd_1, bd_2)
                bd_dict = {'bd_min': cent - 3, 'bd_max': cent + 3}
                bd_pic = self._model.get_background_graphic(**bd_dict)
                pic = bd_pic + pic
                return pic
        else:
            center = (bd_1 + bd_2)/2 # Circle center
            radius = abs(bd_1 - bd_2)/2
            theta1 = CC(end_1 - center).arg()
            theta2 = CC(end_2 - center).arg()
            if abs(theta1 - theta2) < EPSILON:
                theta2 += pi
            [theta1, theta2] = sorted([theta1, theta2])
            from sage.calculus.var import var
            from sage.plot.plot import parametric_plot
            x = var('x')
            pic = parametric_plot((radius*cos(x) + real(center),
                                   radius*sin(x) + imag(center)),
                                  (x, theta1, theta2), **opts)
            if boundary:
                # We want to draw a segment of the real line.  The
                # computations below compute the projection of the
                # geodesic to the real line, and then draw a little
                # to the left and right of the projection.
                shadow_1, shadow_2 = [real(k) for k in [end_1, end_2]]
                midpoint = (shadow_1 + shadow_2)/2
                length = abs(shadow_1 - shadow_2)
                bd_dict = {'bd_min': midpoint - length, 'bd_max': midpoint +
                           length}
                bd_pic = self._model.get_background_graphic(**bd_dict)
                pic = bd_pic + pic
            return pic

    def show(self, boundary=True, **options):
        r"""
        Plot ``self``.

        EXAMPLES:

        First some lines::

            sage: PD = HyperbolicPlane().PD()
            sage: PD.get_geodesic(0, 1).show()
            Graphics object consisting of 2 graphics primitives
            sage: PD.get_geodesic(0, 0.3+0.8*I).show()
            Graphics object consisting of 2 graphics primitives

        Then some generic geodesics::

            sage: PD.get_geodesic(-0.5, 0.3+0.4*I).show()
            Graphics object consisting of 2 graphics primitives
            sage: PD.get_geodesic(-1, exp(3*I*pi/7)).show(linestyle="dashed", color="red")
            Graphics object consisting of 2 graphics primitives
            sage: PD.get_geodesic(exp(2*I*pi/11), exp(1*I*pi/11)).show(thickness=6, color="orange")
            Graphics object consisting of 2 graphics primitives
        """
        opts = {'axes': False, 'aspect_ratio': 1}
        opts.update(self.graphics_options())
        opts.update(options)
        end_1, end_2 = [CC(k.coordinates()) for k in self.endpoints()]
        bd_1, bd_2 = [CC(k.coordinates()) for k in self.ideal_endpoints()]
        # Check to see if it's a line
        if abs(bd_1 + bd_2) < EPSILON:
            pic = line([end_1, end_2], **opts)
        else:
            # If we are here, we know it's not a line
            # So we compute the center and radius of the circle
            invdet = RR.one() / (real(bd_1)*imag(bd_2) - real(bd_2)*imag(bd_1))
            centerx = (imag(bd_2) - imag(bd_1)) * invdet
            centery = (real(bd_1) - real(bd_2)) * invdet
            center = centerx + I * centery
            radius = RR(abs(bd_1 - center))
            # Now we calculate the angles for the arc
            theta1 = CC(end_1 - center).arg()
            theta2 = CC(end_2 - center).arg()
            theta1, theta2 = sorted([theta1, theta2])
            # Make sure the sector is inside the disk
            if theta2 - theta1 > pi:
                theta1 += 2 * pi
            pic = arc((centerx, centery), radius,
                      sector=(theta1, theta2), **opts)
        if boundary:
            pic += self._model.get_background_graphic()
        return pic

    def isometry_from_fixed_points(self, repel, attract):
        r"""
        Given two fixed points ``repel`` and ``attract`` as complex
        numbers return a hyperbolic isometry with ``repel`` as repelling
        fixed point and ``attract`` as attracting fixed point.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: UHP.isometry_from_fixed_points(2 + I, 3 + I)
            Traceback (most recent call last):
            ...
            ValueError: fixed points of hyperbolic elements must be ideal

            sage: UHP.isometry_from_fixed_points(2, 0)
            Isometry in UHP
            [  -1    0]
            [-1/3 -1/3]

        TESTS::

            sage: UHP = HyperbolicPlane().UHP()
            sage: UHP.isometry_from_fixed_points(0, 4)
            Isometry in UHP
            [  -1    0]
            [-1/5 -1/5]
            sage: UHP.isometry_from_fixed_points(UHP.get_point(0), UHP.get_point(4))
            Isometry in UHP
            [  -1    0]
            [-1/5 -1/5]
        """
        if isinstance(repel, HyperbolicPoint):
            repel = repel._coordinates
        if isinstance(attract, HyperbolicPoint):
            attract = attract._coordinates

        if imag(repel) + imag(attract) > EPSILON:
            raise ValueError("fixed points of hyperbolic elements must be ideal")
        repel = real(repel)
        attract = real(attract)
        if repel == infinity:
            A = self._moebius_sending([infinity, attract, attract + 1],
                                     [infinity, attract, attract + 2])
        elif attract == infinity:
            A = self._moebius_sending([repel, infinity, repel + 1],
                                     [repel, infinity, repel + 2])
        else:
            A = self._moebius_sending([repel, attract, infinity],
                                     [repel, attract, max(repel, attract) + 1])
        return self.get_isometry(A)

    def image_coordinates(self, x):
        """
        Return the image of the coordinates of the hyperbolic point ``x``
        under ``self``.

        EXAMPLES::

            sage: PD = HyperbolicPlane().PD()
            sage: KM = HyperbolicPlane().KM()
            sage: phi = KM.coerce_map_from(PD)
            sage: phi.image_coordinates(0.5+0.5*I)
            (0.666666666666667, 0.666666666666667)
        """
        return (2*real(x)/(Integer(1) + real(x)**2 + imag(x)**2),
                2*imag(x)/(Integer(1) + real(x)**2 + imag(x)**2))

    def boundary_point_in_model(self, p):
        r"""
        Check whether a complex number is a real number or ``\infty``.
        In the ``UHP.model_name_name``, this is the ideal boundary of
        hyperbolic space.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: UHP.boundary_point_in_model(1 + I)
            False
            sage: UHP.boundary_point_in_model(infinity)
            True
            sage: UHP.boundary_point_in_model(CC(infinity))
            True
            sage: UHP.boundary_point_in_model(RR(infinity))
            True
            sage: UHP.boundary_point_in_model(1)
            True
            sage: UHP.boundary_point_in_model(12)
            True
            sage: UHP.boundary_point_in_model(1 - I)
            False
            sage: UHP.boundary_point_in_model(-2*I)
            False
            sage: UHP.boundary_point_in_model(0)
            True
            sage: UHP.boundary_point_in_model(I)
            False
        """
        if isinstance(p, HyperbolicPoint):
            return p.is_boundary()
        im = abs(imag(CC(p)).n())
        return (im < EPSILON) or bool(p == infinity)

    def _evalf_(self, z, parent=None, algorithm=None):
        r"""
        EXAMPLES:

        If the imaginary part of `z` obeys `-\pi < z \leq \pi`, then
        `\operatorname{exp\_polar}(z)` is evaluated as `\exp(z)`::

            sage: exp_polar(1.0 + 2.0*I)
            -1.13120438375681 + 2.47172667200482*I

        If the imaginary part of `z` is outside of that interval the
        expression is left unevaluated::

            sage: exp_polar(-5.0 + 8.0*I)
            exp_polar(-5.00000000000000 + 8.00000000000000*I)

        An attempt to numerically evaluate such an expression raises an error::

            sage: exp_polar(-5.0 + 8.0*I).n()
            Traceback (most recent call last):
            ...
            ValueError: invalid attempt to numerically evaluate exp_polar()

        """
        from sage.functions.other import imag

        if (not isinstance(z, Expression)
            and bool(-const_pi < imag(z) <= const_pi)):
            return exp(z)
        else:
            raise ValueError("invalid attempt to numerically evaluate exp_polar()")

    def point_in_model(self, p):
        r"""
        Check whether a complex number lies in the open upper half plane.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: UHP.point_in_model(1 + I)
            True
            sage: UHP.point_in_model(infinity)
            False
            sage: UHP.point_in_model(CC(infinity))
            False
            sage: UHP.point_in_model(RR(infinity))
            False
            sage: UHP.point_in_model(1)
            False
            sage: UHP.point_in_model(12)
            False
            sage: UHP.point_in_model(1 - I)
            False
            sage: UHP.point_in_model(-2*I)
            False
            sage: UHP.point_in_model(I)
            True
            sage: UHP.point_in_model(0) # Not interior point
            False
        """
        if isinstance(p, HyperbolicPoint):
            return p.is_boundary()
        return bool(imag(CC(p)) > 0)

    def _dist_points(self, p1, p2):
        r"""
        Compute the distance between two points in the Upper Half Plane
        using the hyperbolic metric.

        INPUT:

        - ``p1``, ``p2`` -- the coordinates of the points

        EXAMPLES::

           sage: HyperbolicPlane().UHP()._dist_points(4.0*I, I)
           1.38629436111989
        """
        num = (real(p2) - real(p1))**2 + (imag(p2) - imag(p1))**2
        denom = 2 * imag(p1) * imag(p2)
        if denom == 0:
            return infinity
        return arccosh(1 + num/denom)

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: 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

    def _evalf_(self, x, **kwargs):
        """
        EXAMPLES::

            sage: from sage.functions.airy import airy_bi_simple
            sage: airy_bi_simple(0.0)
            0.614926627446001
            sage: airy_bi_simple(1.0 * I)
            0.648858208330395 + 0.344958634768048*I

        We can use several methods for numerical evaluation::

            sage: airy_bi_simple(3).n(algorithm='mpmath')
            14.0373289637302
            sage: airy_bi_simple(3).n(algorithm='mpmath', prec=100)
            14.037328963730232031740267314
            sage: airy_bi_simple(3).n(algorithm='scipy')  # rel tol 1e-10
            14.037328963730136
            sage: airy_bi_simple(I).n(algorithm='scipy')  # rel tol 1e-10
            0.648858208330395 + 0.34495863476804844*I
            
        TESTS::

            sage: parent(airy_bi_simple(3).n(algorithm='scipy'))                                          
            Real Field with 53 bits of precision
            sage: airy_bi_simple(3).n(algorithm='scipy', prec=200)
            Traceback (most recent call last):
            ...
            NotImplementedError: airy_bi not implemented for precision > 53
        """
        algorithm = kwargs.get('algorithm', 'mpmath') or 'mpmath'
        parent = kwargs.get('parent', None)
        if algorithm == 'scipy':
            if hasattr(parent, 'prec') and parent.prec() > 53:
                raise NotImplementedError("%s not implemented for precision > 53"%self.name())
            from sage.rings.all import RR, CC
            from sage.functions.other import real,imag
            from scipy.special import airy as airy
            if x in RR:
                y = airy(real(x))[2]
                if parent is None:
                    return RR(y)
            else:
                y = airy(complex(real(x),imag(x)))[2]
                if parent is None:
                    return CC(y)
            return parent(y)
        elif algorithm == 'mpmath':
            import mpmath
            from sage.libs.mpmath import utils as mpmath_utils
            return mpmath_utils.call(mpmath.airybi, x, parent=parent)
        else:
            raise ValueError("unknown algorithm '%s'" % algorithm)

    def symmetry_involution(self):
        r"""
        Return the involutory isometry fixing the given point.

        EXAMPLES::

            sage: HyperbolicPlane().UHP().get_point(3 + 2*I).symmetry_involution()
            Isometry in UHP
            [  3/2 -13/2]
            [  1/2  -3/2]
        """
        p = self._coordinates
        x, y = real(p), imag(p)
        if y > 0:
            M = matrix([[x/y, -(x**2/y) - y], [1/y, -(x/y)]])
            return self.parent().get_isometry(M)
        raise ValueError("cannot determine the isometry of a boundary point")

    def ideal_endpoints(self):
        r"""
        Determine the ideal (boundary) endpoints of the complete
        hyperbolic geodesic corresponding to ``self``.

        OUTPUT:

        - a list of 2 boundary points

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: UHP.get_geodesic(I, 2*I).ideal_endpoints()
            [Boundary point in UHP 0,
             Boundary point in UHP +Infinity]
            sage: UHP.get_geodesic(1 + I, 2 + 4*I).ideal_endpoints()
            [Boundary point in UHP -sqrt(65) + 9,
             Boundary point in UHP sqrt(65) + 9]
        """
        start = self._start.coordinates()
        end = self._end.coordinates()
        [x1, x2] = [real(k) for k in [start, end]]
        [y1, y2] = [imag(k) for k in [start, end]]
        M = self._model
        # infinity is the first endpoint, so the other ideal endpoint
        # is just the real part of the second coordinate
        if start == infinity:
            return [M.get_point(start), M.get_point(x2)]
        # Same idea as above
        if end == infinity:
            return [M.get_point(x1), M.get_point(end)]
        # We could also have a vertical line with two interior points
        if x1 == x2:
            return [M.get_point(x1), M.get_point(infinity)]
        # Otherwise, we have a semicircular arc in the UHP
        c = ((x1+x2)*(x2-x1) + (y1+y2)*(y2-y1)) / (2*(x2-x1))
        r = sqrt((c - x1)**2 + y1**2)
        return [M.get_point(c - r), M.get_point(c + r)]

    def midpoint(self):  # UHP
        r"""
        Return the (hyperbolic) midpoint of ``self`` if it exists.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: g = UHP.random_geodesic()
            sage: m = g.midpoint()
            sage: d1 = UHP.dist(m, g.start())
            sage: d2 = UHP.dist(m, g.end())
            sage: bool(abs(d1 - d2) < 10**-9)
            True

        Infinite geodesics have no midpoint::

            sage: UHP.get_geodesic(0, 2).midpoint()
            Traceback (most recent call last):
            ...
            ValueError: the length must be finite
        """
        if self.length() == infinity:
            raise ValueError("the length must be finite")

        start = self._start.coordinates()
        end = self._end.coordinates()
        d = self._model._dist_points(start, end) / 2
        S = self.complete()._to_std_geod(start)
        T = matrix([[exp(d), 0], [0, 1]])
        M = S.inverse() * T * S
        if ((real(start - end) < EPSILON)
                or (abs(real(start - end)) < EPSILON
                    and imag(start - end) < EPSILON)):
            end_p = start
        else:
            end_p = end
        return self._model.get_point(mobius_transform(M, end_p))