#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]])