#mobius.py
import cmath, numpy
class circle_involution:
def __init__(self, center, radius):
self.center = center
self.radius = radius
def transform(self):
c = self.center
r = self.radius
return translation(conj(c))*self.scaling(r)*involution*self.inv_scaling(1.0/r)*self.inv_translation(-c)
def of(self, z):
"Apply to homogenous coordinate z."
w = self.transform()*z
return conj(w)
def of_circle(self, center, radius):
"Apply to circle involution."
if radius == 0:
return 0, 0
# We find the involution of center, radius
c = self.center
r = self.radius
# We conjugate it with our own involution
C = translation(conj(c))*scaling(r)*involution*scaling(1.0/r)*translation(-c)
D = translation(center)*scaling(radius)*involution*scaling(1.0/radius)*translation(conj(-center))
E = translation(conj(c))*scaling(r)*involution*scaling(1.0/r)*translation(-c)
image_transform = C*D*E
# The result is itself an involution whose center/radius we return
a = image_transform.item(0)
b = image_transform.item(1)
c = image_transform.item(2)
d = image_transform.item(3)
a /= c
b /= c
d /= c
c = 1.0
new_center = -d
new_radius = abs(b + (new_center*conj(new_center)))**.5
return new_center, new_radius
def translation(t):
return numpy.matrix([[1.0, t], [0.0, 1.0]])
def scaling(r):
return numpy.matrix([[r, 0.0], [0.0, 1.0]])
def conj(z):
return z - 2j*z.imag
def normalize(z):
"Normalize homogenous coordinates so last is 1.0"
a = z.item(0)
b = z.item(1)
if b != 0.0 + 0j:
return numpy.matrix([[a/b], [1.0 + 0j]])
else:
return numpy.matrix([[1.0 + 0j], [0.0 + 0j]])
involution = numpy.matrix([[0.0, 1.0], [1.0, 0.0]])
