#!/usr/bin/env python3

# https://math.stackexchange.com/a/4500458/746312

from math import sqrt, sin, cos, acos, pi

class Point:
    def __init__ (self, x, y):
        self.x, self.y = x, y

    def __add__ (self, Q):
        return Point (self.x + Q.x, self.y + Q.y)

    def __sub__ (self, Q):
        return Point (self.x - Q.x, self.y - Q.y)

    def __truediv__ (self, Q):
        return Point (self.x / Q, self.y / Q)

    def __mul__ (self, Q):
        if isinstance (Q, Point):
            return self.x * Q.x + self.y * Q.y
        return Point (self.x * Q, self.y * Q)

    def __rmul__ (self, Q):
        return self * Q

    def abs2 (self):
        return self * self

    def abs (self):
        return sqrt (self.abs2())

    def bot (self):
        return Point (self.y, -self.x)

    def cos (self, Q):
        return (self * Q) / (self.abs() * Q.abs())

    def __repr__ (self):
        return "[%s, %s]" % (self.x, self.y)

    def rot (self, phi):
        s, c = sin (phi), cos(phi)
        return Point (c * self.x + s * self.y,
                      c * self.y - s * self.x)

def center (A, B, phi):
    """Get the center M of a circle through points A and B such that a
       triangle ABC has an angle of PHI at C.  This means C lies on the
       circumcircle around M through A and B."""

    rotAB = (B-A).rot (pi/2 - phi)
    lam = (0.5 * (B*B + A*A) - A*B) / (rotAB * (B-A))
    return A + lam * rotAB

def intersect (m1, m2, B):
    """Let Ck be a circle around mk through B.  Return the intersection
       of C1 and C2 that is not B."""
    # Reflect B at line m1---m2.
    m = m2 - m1
    lam =  ((B - m1) * m) / (m * m)
    beta = ((B - m1) * m.bot()) / (m * m)
    return m1 + lam * m - beta * m.bot()

def show_solution (A, B, C, alpha, theta):
    M1 = center (A, B, alpha)
    M2 = center (B, C, theta)
    print ("M1 =", M1)
    print ("M2 =", M2)

    U = intersect (M1, M2, B)
    print ("U =", U)

    w1 = (U-A).cos(U-B)
    w2 = (U-B).cos(U-C)
    print (acos(w1), "= pi *", acos(w1) / pi)
    print (acos(w2), "= pi *", acos(w2) / pi)

A = Point (2, 3)
B = Point (6, 5)
C = Point (8, 5)

alpha = 45 * pi/180
theta = 55 * pi/180

show_solution (A, B, C, alpha, theta)
show_solution (C, B, A, theta, alpha)