"""
Utility functions for complex analysis of plane geometry.
Uses and generates (numpy) complex numbers.
"""
import numpy as np
import ddg.math.euclidean2d as euclidean2d
###########
# sampling
###########
[docs]def complex_sample(size=None, seed=None):
"""
Get a sample of random complex numbers (or vectors).
Parameters
----------
size : int or tuple of ints (default=None)
size of sample of complex numbers
or shape of numpy.ndarray
seed : int (default=None)
seed for obtaining a deterministic function
Returns
-------
complex or numpy.ndarray of numpy.complex128
Notes
-----
Real and imaginary part are standard normally distributed
"""
gen = np.random.default_rng(seed=seed)
re = gen.standard_normal(size=size)
im = gen.standard_normal(size=size)
return re + 1j * im
#############
# conversion
#############
[docs]def to_array(z):
"""
Convert complex numbers to 2D arrays of real and imaginary part.
Parameters
----------
z : complex
Returns
-------
numpy.ndarray of shape (2,)
"""
return np.array([z.real, z.imag])
[docs]def to_complex(x):
"""
Convert 2D array to complex number.
Parameters
----------
x : sequence of length 2
Returns
-------
complex
"""
return complex(x[0], x[1])
[docs]def homogeneous(z):
"""
Get complex number in homogeneous coordinates.
z is interpreted as the affine coordinate of an element in CP^1.
Parameters
----------
z : complex
Returns
-------
numpy.ndarray of shape (2,)
"""
if np.isinf(z):
return np.array([1.0, 0.0])
else:
return np.array([z, 1.0])
#################
# plane geometry
#################
[docs]def scalar_product(z1, z2):
"""
For two complex numbers get scalar product of corresponding elements in R^2.
Parameters
----------
z1,z2 : complex
Returns
-------
float
"""
return (z1.conjugate() * z2).real
[docs]def determinant(z1, z2):
"""
For two complex numbers get determinant of matrix of corresponding elements in R^2.
Parameters
----------
z1,z2 : complex
Returns
-------
float
"""
return (z1.conjugate() * z2).imag
[docs]def rel_angle(z1, z2):
"""Get angle in [0,2*pi) between complex numbers z1,z2.
Parameters
----------
z1,z2 : complex
Returns
-------
float
"""
alpha = np.angle(z2 / z1)
return alpha if alpha >= 0 else 2 * np.pi + alpha
###############
# cross-ratios
###############
[docs]def cr(z1, z2, z3, z4):
"""Compute complex cross ratio.
cr(z1,z2,z3,z4) = (z1 - z2)/(z2 - z3) * (z3 - z4)/(z4 - z1)
Parameters
----------
z1,z2,z3,z4 : complex
Returns
-------
complex
"""
z1, z2, z3, z4 = map(homogeneous, [z1, z2, z3, z4])
a = np.linalg.det(np.array([z1, z2]))
b = np.linalg.det(np.array([z2, z3]))
c = np.linalg.det(np.array([z3, z4]))
d = np.linalg.det(np.array([z4, z1]))
if b == 0.0 or d == 0.0:
return complex("infinity")
else:
return (a * c) / (b * d)
[docs]def fourth_point_from_cross_ratio(z, z1, z2, q):
"""
Get z12 such that cr(z, z1, z12, z2) = q
For a quadrilateral of the form ::
z2_____z12
| |
| |
z_____z1
with given z,z1,z2 and q, compute z12 .
Parameters
----------
z,z1,z2,q : complex
Returns
-------
complex
"""
a = (z - z1) / (z2 - z)
return (a * z2 + q * z1) / (a + q)
################
# circumcircles
################
[docs]def circumcenter(z1, z2, z3):
"""
Get center of circle through three complex numbers.
Parameters
----------
z1,z2,z3 : complex
Returns
-------
numpy.complex128
"""
x, y, z = map(lambda x: np.array([x.real, x.imag]), [z1, z2, z3])
D = 2 * (x[0] * (y[1] - z[1]) + y[0] * (z[1] - x[1]) + z[0] * (x[1] - y[1]))
s_x = x[0] ** 2 + x[1] ** 2
s_y = y[0] ** 2 + y[1] ** 2
s_z = z[0] ** 2 + z[1] ** 2
u = (s_x * (y[1] - z[1]) + s_y * (z[1] - x[1]) + s_z * (x[1] - y[1])) / D
v = (s_x * (z[0] - y[0]) + s_y * (x[0] - z[0]) + s_z * (y[0] - x[0])) / D
return u + 1j * v
[docs]def circumradius(z1, z2, z3):
"""
Get radius of circle through three complex numbers.
Parameters
----------
z1,z2,z3 : complex
Returns
-------
float
"""
x, y, z = map(lambda x: np.array([x.real, x.imag]), [z1, z2, z3])
l_xy, l_yz, l_zx = map(np.linalg.norm, (x - y, y - z, z - x))
return (l_xy * l_yz * l_zx) / (4 * euclidean2d.triangle_area(l_xy, l_yz, l_zx))
########
# quads
########
[docs]def intersect_diags(z1, z2, z3, z4):
r"""
For a complex quadrilateral (z1,z2,z3,z4)
get the intersection point of the diagonals.
The four complex points are given as an argument in positive cyclic order. ::
z4-----z3
|\ /|
| x |
|/ \|
z1-----z2
Parameters
----------
z1,z2,z3,z4 : complex
Returns
-------
complex
"""
n1 = 1j * (z3 - z1)
n2 = 1j * (z4 - z2)
A = np.array([[n1.real, n1.imag], [n2.real, n2.imag]])
b = np.array([scalar_product(n1, z1), scalar_product(n2, z2)])
x = np.linalg.solve(A, b)
return complex(*x)
[docs]def intersect_edges(z1, z2, z3, z4):
r"""
For a complex quadrilateral (z1,z2,z3,z4)
get the intersection of the two opposite edges (z1,z2) and (z3,z4).
Will return complex infinity if intersection does not exist. ::
z4-----z3--
| | \
| | x
| | /
z1-----z2--
Parameters
----------
z1,z2,z3,z4 : complex
Returns
-------
complex
"""
n1 = 1j * (z2 - z1)
n2 = 1j * (z4 - z3)
A = np.array([[n1.real, n1.imag], [n2.real, n2.imag]])
b = np.array([scalar_product(n1, z1), scalar_product(n2, z3)])
if np.linalg.cond(A) < 1 / np.finfo(float).eps:
x = np.linalg.solve(A, b)
else:
return complex("infinity")
return complex(*x)