Coverage for local/lib/python2.7/site-packages/sage/geometry/hyperbolic_space/hyperbolic_geodesic.py : 93%

Point in UHP 1/2*(sqrt(2)*e^(1/2*arccosh(3)) - sqrt(2) + (I - 1)*e^(1/2*arccosh(3)) + I - 1)/((1/4*I - 1/4)*sqrt(2)*e^(1/2*arccosh(3)) - (1/4*I - 1/4)*sqrt(2) + 1/2*e^(1/2*arccosh(3)) + 1/2)

Check that floating points remain floating points

in :meth:`midpoint` ::

sage: UHP = HyperbolicPlane().UHP()

sage: g = UHP.get_geodesic(CC(0,1), CC(2,2))

sage: g.midpoint()

Point in UHP 0.666666666666667 + 1.69967317119760*I

sage: parent(g.midpoint().coordinates())

Complex Field with 53 bits of precision

"""

from sage.matrix.matrix_symbolic_dense import Matrix_symbolic_dense

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)

# If the matrix is symbolic then needs to be simplified in order to

# make the calculations easier for the symbolic calculus module.

if isinstance(S, Matrix_symbolic_dense):

S = S.simplify_full().simplify_full()

S_1 = S.inverse()

T = matrix([[exp(d), 0], [0, 1]])

M = S_1 * 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

P_3 = moebius_transform(M, end_p)

return self._model.get_point(P_3)

def angle(self, other): # UHP

r"""

Return the angle between any two given completed geodesics if

they intersect.

INPUT:

- ``other`` -- a hyperbolic geodesic in the UHP model

OUTPUT:

- the angle in radians between the two given geodesics

EXAMPLES::

sage: UHP = HyperbolicPlane().UHP()

sage: g = UHP.get_geodesic(2, 4)

sage: h = UHP.get_geodesic(3, 3 + I)

sage: g.angle(h)

1/2*pi

sage: numerical_approx(g.angle(h))

1.57079632679490

.. PLOT::

UHP = HyperbolicPlane().UHP()

g = UHP.get_geodesic(2, 4)

h = UHP.get_geodesic(3, 3 + I)

sphinx_plot(g.plot()+h.plot())

If the geodesics are identical, return angle 0::

sage: g.angle(g)

0

It is an error to ask for the angle of two geodesics that do not

intersect::

sage: g = UHP.get_geodesic(2, 4)

sage: h = UHP.get_geodesic(5, 7)

sage: g.angle(h)

Traceback (most recent call last):

...

ValueError: geodesics do not intersect

"""

if self.is_parallel(other):

raise ValueError("geodesics do not intersect")

# Make sure the segments are complete or intersect

if (not(self.is_complete() and other.is_complete()) and

not self.intersection(other)):

print("Warning: Geodesic segments do not intersect. "

"The angle between them is not defined.\n"

"Returning the angle between their completions.")

# Make sure both are in the same model

if other._model is not self._model:

other = other.to_model(self._model)

(p1, p2) = sorted([k.coordinates()

for k in self.ideal_endpoints()], key=str)

(q1, q2) = sorted([k.coordinates()

for k in other.ideal_endpoints()], key=str)

# if the geodesics are equal, the angle between them is 0

if (abs(p1 - q1) < EPSILON and

abs(p2 - q2) < EPSILON):

return 0

elif p2 != infinity: # geodesic not a straight line

# So we send it to the geodesic with endpoints [0, oo]

T = HyperbolicGeodesicUHP._crossratio_matrix(p1, (p1 + p2) / 2, p2)

else:

# geodesic is a straight line, so we send it to the geodesic

# with endpoints [0,oo]

T = HyperbolicGeodesicUHP._crossratio_matrix(p1, p1 + 1, p2)

# b1 and b2 are the endpoints of the image of other

b1, b2 = [moebius_transform(T, k) for k in [q1, q2]]

# If other is now a straight line...

if (b1 == infinity or b2 == infinity):

# then since they intersect, they are equal

return 0

return real(arccos((b1 + b2) / abs(b2 - b1)))

##################

# Helper methods #

##################

@staticmethod

def _get_B(a):

r"""

Helper function to get an appropriate matrix transforming

(0,1,inf) -> (0,I,inf) based on the type of a

INPUT:

- ``a`` -- an element to identify the class of the resulting matrix.

EXAMPLES::

sage: UHP = HyperbolicPlane().UHP()

sage: g = UHP.random_geodesic()

sage: B = g._get_B(CDF.an_element()); B

[ 1.0 0.0]

[ 0.0 -1.0*I]

sage: type(B)

::

sage: B = g._get_B(SR(1)); B

[ 1 0]

[ 0 -I]

sage: type(B)

::

sage: B = g._get_B(complex(1)); B

[ 1.0 0.0]

[ 0.0 -1.0*I]

sage: type(B)

::

sage: B = g._get_B(QQbar(1+I)); B

[ 1 0]

[ 0 -I]

sage: type(B[1,1])

sage: type(B)

"""

from sage.structure.element import Element

from sage.symbolic.expression import Expression

from sage.rings.complex_double import CDF

if isinstance(a, (int, float, complex)): # Python number

a = CDF(a)

if isinstance(a, Expression): # symbolic

P = SR

zero = SR.zero()

one = SR.one()

I = SR("I")

elif isinstance(a, Element): # Sage number

P = a.parent()

zero = P.zero()

one = P.one()

I = P.gen()

if I.is_one() or (I*I).is_one() or not (-I*I).is_one():

raise ValueError("invalid number")

else:

raise ValueError("not a complex number")

return matrix(P, 2, [one, zero, zero, -I])

def _to_std_geod(self, p):

r"""

Given the coordinates of a geodesic in hyperbolic space, return the

hyperbolic isometry that sends that geodesic to the geodesic

through 0 and infinity that also sends the point ``p`` to `i`.

INPUT:

- ``p`` -- the coordinates of the

EXAMPLES::

sage: UHP = HyperbolicPlane().UHP()

sage: (p1, p2) = [UHP.random_point() for k in range(2)]

sage: g = UHP.get_geodesic(p1, p2)

sage: m = g.midpoint()

sage: A = g._to_std_geod(m.coordinates()) # Send midpoint to I.

sage: A = UHP.get_isometry(A)

sage: [s, e]= g.complete().endpoints()

sage: bool(abs(A(s).coordinates()) < 10**-9)

True

sage: bool(abs(A(m).coordinates() - I) < 10**-9)

True

sage: bool(abs(A(e).coordinates()) > 10**9)

True

TESTS:

Check that floating points remain floating points through this method::

sage: H = HyperbolicPlane()

sage: g = H.UHP().get_geodesic(CC(0,1), CC(2,2))

sage: gc = g.complete()

sage: parent(gc._to_std_geod(g.start().coordinates()))

Full MatrixSpace of 2 by 2 dense matrices over Complex Field

with 53 bits of precision

"""

[s, e] = [k.coordinates() for k in self.complete().endpoints()]

B = HyperbolicGeodesicUHP._get_B(p)

# outmat below will be returned after we normalize the determinant.

outmat = B * HyperbolicGeodesicUHP._crossratio_matrix(s, p, e)

outmat = outmat / outmat.det().sqrt()

if (outmat - outmat.conjugate()).norm(1) < 10**-9:

# Small imaginary part.

outmat = (outmat + outmat.conjugate()) / 2

# Set it equal to its real part.

return outmat

@staticmethod

def _crossratio_matrix(p0, p1, p2): # UHP

r"""

Given three points (the list `p`) in `\mathbb{CP}^{1}` in affine

coordinates, return the linear fractional transformation taking

the elements of `p` to `0`, `1`, and `\infty`.

INPUT:

- a list of three distinct elements

of `\mathbb{CP}^1` in affine coordinates; that is, each element

must be a complex number, `\infty`, or symbolic.

OUTPUT:

- an element of `\GL(2,\CC)`

EXAMPLES::

sage: from sage.geometry.hyperbolic_space.hyperbolic_geodesic \

....: import HyperbolicGeodesicUHP

sage: from sage.geometry.hyperbolic_space.hyperbolic_isometry \

....: import moebius_transform

sage: UHP = HyperbolicPlane().UHP()

sage: (p1, p2, p3) = [UHP.random_point().coordinates()

....: for k in range(3)]

sage: A = HyperbolicGeodesicUHP._crossratio_matrix(p1, p2, p3)

sage: bool(abs(moebius_transform(A, p1)) < 10**-9)

True

sage: bool(abs(moebius_transform(A, p2) - 1) < 10**-9)

True

sage: bool(moebius_transform(A, p3) == infinity)

True

sage: (x,y,z) = var('x,y,z');

sage: HyperbolicGeodesicUHP._crossratio_matrix(x,y,z)

[ y - z -x*(y - z)]

[ -x + y (x - y)*z]

"""

if p0 == infinity:

return matrix([[0, -(p1 - p2)], [-1, p2]])

elif p1 == infinity:

return matrix([[1, -p0], [1, -p2]])

elif p2 == infinity:

return matrix([[1, -p0], [0, p1 - p0]])

return matrix([[p1 - p2, (p1 - p2)*(-p0)],

[p1 - p0, (p1 - p0)*(-p2)]])

#***********************************************************************

# Other geodesics

#***********************************************************************

class HyperbolicGeodesicPD(HyperbolicGeodesic):

r"""

A geodesic in the Poincaré disk model.

Geodesics in this model are represented by segments of circles contained

within the unit disk that are orthogonal to the boundary of the disk,

plus all diameters of the disk.

INPUT:

- ``start`` -- a :class:`HyperbolicPoint` in hyperbolic space

representing the start of the geodesic

- ``end`` -- a :class:`HyperbolicPoint` in hyperbolic space

representing the end of the geodesic

EXAMPLES::

sage: PD = HyperbolicPlane().PD()

sage: g = PD.get_geodesic(PD.get_point(I), PD.get_point(-I/2))

sage: g = PD.get_geodesic(I,-I/2)

sage: h = PD.get_geodesic(-1/2+I/2,1/2+I/2)

.. PLOT::

PD = HyperbolicPlane().PD()

g = PD.get_geodesic(I,-I/2)

h = PD.get_geodesic(-0.5+I*0.5,0.5+I*0.5)

sphinx_plot(g.plot()+h.plot(color='green'))

"""

def plot(self, boundary=True, **options):

r"""

Plot ``self``.

EXAMPLES:

First some lines::

sage: PD = HyperbolicPlane().PD()

sage: PD.get_geodesic(0, 1).plot()

Graphics object consisting of 2 graphics primitives

.. PLOT::

sphinx_plot(HyperbolicPlane().PD().get_geodesic(0, 1).plot())

::

sage: PD.get_geodesic(0, 0.3+0.8*I).plot()

Graphics object consisting of 2 graphics primitives

.. PLOT::

PD = HyperbolicPlane().PD()

sphinx_plot(PD.get_geodesic(0, 0.3+0.8*I).plot())

Then some generic geodesics::

sage: PD.get_geodesic(-0.5, 0.3+0.4*I).plot()

Graphics object consisting of 2 graphics primitives

sage: g = PD.get_geodesic(-1, exp(3*I*pi/7))

sage: G = g.plot(linestyle="dashed",color="red"); G

Graphics object consisting of 2 graphics primitives

sage: h = PD.get_geodesic(exp(2*I*pi/11), exp(1*I*pi/11))

sage: H = h.plot(thickness=6, color="orange"); H

Graphics object consisting of 2 graphics primitives

sage: show(G+H)

.. PLOT::

PD = HyperbolicPlane().PD()

PD.get_geodesic(-0.5, 0.3+0.4*I).plot()

g = PD.get_geodesic(-1, exp(3*I*pi/7))

G = g.plot(linestyle="dashed",color="red")

h = PD.get_geodesic(exp(2*I*pi/11), exp(1*I*pi/11))

H = h.plot(thickness=6, color="orange")

sphinx_plot(G+H)

"""

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

show = deprecated_function_alias(20530, plot)

class HyperbolicGeodesicKM(HyperbolicGeodesic):

r"""

A geodesic in the Klein disk model.

Geodesics are represented by the chords, straight line segments with ideal

endpoints on the boundary circle.

INPUT:

- ``start`` -- a :class:`HyperbolicPoint` in hyperbolic space

representing the start of the geodesic

- ``end`` -- a :class:`HyperbolicPoint` in hyperbolic space

representing the end of the geodesic

EXAMPLES::

sage: KM = HyperbolicPlane().KM()

sage: g = KM.get_geodesic(KM.get_point((0.1,0.9)), KM.get_point((-0.1,-0.9)))

sage: g = KM.get_geodesic((0.1,0.9),(-0.1,-0.9))

sage: h = KM.get_geodesic((-0.707106781,-0.707106781),(0.707106781,-0.707106781))

sage: P = g.plot(color='orange')+h.plot(); P

Graphics object consisting of 4 graphics primitives

.. PLOT::

KM = HyperbolicPlane().KM()

g = KM.get_geodesic((0.1,0.9),

(-0.1,-0.9))

h = KM.get_geodesic((-0.707106781,-0.707106781),

(0.707106781,-0.707106781))

sphinx_plot(g.plot(color='orange')+h.plot())

"""

def plot(self, boundary=True, **options):

r"""

Plot ``self``.

EXAMPLES::

sage: HyperbolicPlane().KM().get_geodesic((0,0), (1,0)).plot()

Graphics object consisting of 2 graphics primitives

.. PLOT::

KM = HyperbolicPlane().KM()

sphinx_plot(KM.get_geodesic((0,0), (1,0)).plot())

"""

opts = {'axes': False, 'aspect_ratio': 1}

opts.update(self.graphics_options())

opts.update(options)

pic = line([k.coordinates() for k in self.endpoints()], **opts)

if boundary:

pic += self._model.get_background_graphic()

return pic

show = deprecated_function_alias(20530, plot)

class HyperbolicGeodesicHM(HyperbolicGeodesic):

r"""

A geodesic in the hyperboloid model.

Valid points in the hyperboloid model satisfy :math:`x^2 + y^2 - z^2 = -1`

INPUT:

- ``start`` -- a :class:`HyperbolicPoint` in hyperbolic space

representing the start of the geodesic

- ``end`` -- a :class:`HyperbolicPoint` in hyperbolic space

representing the end of the geodesic

EXAMPLES::

sage: from sage.geometry.hyperbolic_space.hyperbolic_geodesic import *

sage: HM = HyperbolicPlane().HM()

sage: p1 = HM.get_point((4, -4, sqrt(33)))

sage: p2 = HM.get_point((-3,-3,sqrt(19)))

sage: g = HM.get_geodesic(p1, p2)

sage: g = HM.get_geodesic((4, -4, sqrt(33)), (-3, -3, sqrt(19)))

.. PLOT::

HM = HyperbolicPlane().HM()

p1 = HM.get_point((4, -4, sqrt(33)))

p2 = HM.get_point((-3,-3,sqrt(19)))

g = HM.get_geodesic(p1, p2)

sphinx_plot(g.plot(color='blue'))

"""

def plot(self, show_hyperboloid=True, **graphics_options):

r"""

Plot ``self``.

EXAMPLES::

sage: from sage.geometry.hyperbolic_space.hyperbolic_geodesic \

....: import *

sage: g = HyperbolicPlane().HM().random_geodesic()

sage: g.plot()

Graphics3d Object

.. PLOT::

sphinx_plot(HyperbolicPlane().HM().random_geodesic().plot())

"""

x = SR.var('x')

opts = self.graphics_options()

opts.update(graphics_options)

v1, u2 = [vector(k.coordinates()) for k in self.endpoints()]

# Lorentzian Gram Shmidt. The original vectors will be

# u1, u2 and the orthogonal ones will be v1, v2. Except

# v1 = u1, and I don't want to declare another variable,

# hence the odd naming convention above.

# We need the Lorentz dot product of v1 and u2.

v1_ldot_u2 = u2[0]*v1[0] + u2[1]*v1[1] - u2[2]*v1[2]

v2 = u2 + v1_ldot_u2 * v1

v2_norm = sqrt(v2[0]**2 + v2[1]**2 - v2[2]**2)

v2 = v2 / v2_norm

v2_ldot_u2 = u2[0]*v2[0] + u2[1]*v2[1] - u2[2]*v2[2]

# Now v1 and v2 are Lorentz orthogonal, and |v1| = -1, |v2|=1

# That is, v1 is unit timelike and v2 is unit spacelike.

# This means that cosh(x)*v1 + sinh(x)*v2 is unit timelike.

hyperbola = cosh(x)*v1 + sinh(x)*v2

endtime = arcsinh(v2_ldot_u2)

from sage.plot.plot3d.all import parametric_plot3d

pic = parametric_plot3d(hyperbola, (x, 0, endtime), **graphics_options)

if show_hyperboloid:

pic += self._model.get_background_graphic()

return pic

show = deprecated_function_alias(20530, plot)

« index coverage.py v4.5.1, created at 2018-05-01 10:21