isogonal.py

from sympy import expand, poly, simplify
from cartesian import *

def foot(P1, P2, P3):
    x, y = symbols('x, y')
    P1P2 = line(P1, P2)
    P3H = Eq((P1[0] - P2[0])*(x - P3[0]) + (P1[1] - P2[1])*(y - P3[1]), 0)
    return intersect(P1P2, P3H)

def circle(P1, P2, P3):
    # return F(x, y) such that F(x, y) = 0 is the circle's equation
    D, E, F, x, y = symbols('D, E, F, x, y')
    circle_eq = Eq(x**2 + y**2 + D*x + E*y + F, 0)
    circle_eqs = []
    circle_eqs.append(circle_eq.subs(x, P1[0]).subs(y, P1[1]))
    circle_eqs.append(circle_eq.subs(x, P2[0]).subs(y, P2[1]))
    circle_eqs.append(circle_eq.subs(x, P3[0]).subs(y, P3[1]))
    s = solve(circle_eqs, (D, E, F))
    return x**2 + y**2 + s[D]*x + s[E]*y + s[F]

def center(c):
    x, y = symbols('x, y')
    return cancel(-poly(c, x).nth(1)/2), cancel(-poly(c, y).nth(1)/2)

def cos2eq(a2, b2, c2, d2, e2, f2):
    left_num = fraction(cancel(a2 + b2 - c2))
    left_den1 = fraction(cancel(a2))
    left_den2 = fraction(cancel(b2))
    right_num = fraction(cancel(d2 + e2 - f2))
    right_den1 = fraction(cancel(d2))
    right_den2 = fraction(cancel(e2))
    left = left_num[0]**2*left_den1[1]*left_den2[1]*right_num[1]**2*right_den1[0]*right_den2[0]
    right = right_num[0]**2*right_den1[1]*right_den2[1]*left_num[1]**2*left_den1[0]*left_den2[0]
    return expand(left - right) == 0

def main():
    a, b, c, d, e, x, y = symbols('a, b, c, d, e, x, y')
    A, B, C, P = (0, a), (b, 0), (c, 0), (d, e)
    D, E, F = foot(B, C, P), foot(C, A, P), foot(A, B, P)
    print('D:', D)
    print('E:', E)
    print('F:', F)
    pedal_circle = circle(D, E, F)
    print('Pedal Circle:', pedal_circle, '= 0')
    M = center(pedal_circle)
    print('M:', M)
    Q = (2*M[0] - P[0], 2*M[1] - P[1])
    print('Q:', Q)
    G, H, J = foot(B, C, Q), foot(C, A, Q), foot(A, B, Q)
    print('G:', G)
    print('H:', H)
    print('J:', J)
    print('Is G on Pedal Circle?', simplify(pedal_circle.subs(x, G[0]).subs(y, G[1])) == 0)
    print('Is H on Pedal Circle?', simplify(pedal_circle.subs(x, H[0]).subs(y, H[1])) == 0)
    print('Is J on Pedal Circle?', simplify(pedal_circle.subs(x, J[0]).subs(y, J[1])) == 0)
    AB2, AC2, BC2 = dist2(A, B), dist2(A, C), dist2(B, C)
    PA2, PB2, PC2 = dist2(P, A), dist2(P, B), dist2(P, C)
    QA2, QB2, QC2 = dist2(Q, A), dist2(Q, B), dist2(Q, C)
    print('cos²∠PAB = cos²∠QAC?', cos2eq(PA2, AB2, PB2, QA2, AC2, QC2))
    print('cos²∠PAC = cos²∠QAB?', cos2eq(PA2, AC2, PC2, QA2, AB2, QB2))
    print('cos²∠PBA = cos²∠QBC?', cos2eq(PB2, AB2, PA2, QB2, BC2, QC2))
    print('cos²∠PBC = cos²∠QBA?', cos2eq(PB2, BC2, PC2, QB2, AB2, QA2))
    print('cos²∠PCA = cos²∠QCB?', cos2eq(PC2, AC2, PA2, QC2, BC2, QB2))
    print('cos²∠PCB = cos²∠QCA?', cos2eq(PC2, BC2, PB2, QC2, AC2, QA2))

if __name__ == '__main__':
    main()