Source code for ddg.math.complex

"""Complex analysis for plane geometry.
"""

import numpy as np

# --- sampling ---
def complex_sample(size=None):
    """Wraps random_sample() to get corresponding number of complex valued samples
    where real and imaginary part are identically uniformly distributed on [0,1]"""
    return np.random.random_sample(size) + 1j * np.random.random_sample(size)


# --- conversion ---
def to_array(z):
    """Corresponding element in R^2"""
    return np.array([z.real, z.imag])

def to_complex(x):
    """Converts 2D-array to complex number"""
    return complex(x[0], x[1])

def homogeneous(z):
    """Homogeneous coordinates in CP^2 interpreting z as affine coordinates."""
    if np.isinf(z):
        return np.array([1., 0.])
    else:
        return np.array([z, 1.])


# --- plane geometry ---
def scalar_product(z1, z2):
    """Real scalar product of the corresponding elements in R^2"""
    return (z1.conjugate() * z2).real

def determinant(z1, z2):
    """Real determinant of the corresponding elements in R^2"""
    return (z1.conjugate() * z2).imag

def rel_angle(u, v):
    """Positive angle in [0,2*pi) between complex numbers u,v"""
    alpha = np.angle(v/u)
    return alpha if alpha >= 0 else 2*np.pi + alpha


# --- cross-ratios ---
def cr(z1,z2,z3,z4):
    """Complex cross-ratio (z1 - z2)/(z2 - z3) * (z3 - z4)/(z4 - z1)"""
    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. or d == 0.:
        return complex('infinity')
    else:
        return (a*c)/(b*d)

# def cross_ratio(z1, z2, z3, z4):
#   return ((z1-z2)*(z3-z4))/((z2-z3)*(z4-z1))

def fourth_point_from_cross_ratio(z, z1, z2, q):
  """Get z12 in cr(z, z1, z12, z2) = q"""
  a = (z-z1)/(z2-z)
  return (a*z2 + q*z1)/(a + q)


# --- circumcircles ---
def triangle_area_from_lengths(a, b, c):
  """Numerically stable version of Heron's formula
     for the area of a tringagle with side lengths a,b,c. """
  return np.sqrt( (a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)) ) / 4

def circumcircle_center(z1, z2, z3):
    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 np.array([u,v])

def circumcircle_radius(z1, z2, z3):
    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*triangle_area_from_lengths(l_xy, l_yz, l_zx) )


# --- quads ---
def diagonals_intersection(z1, z2, z3, z4):
    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)

def intersect_edges(z1, z2, z3, z4):
    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)])
    x = np.linalg.solve(A,b)
    return complex(*x)