from sage.interfaces.phc import phc
R7blk.<w,r12,r13,r23,m1,m2,m3> = MPolynomialRing(QQ,7,order = "lex(1),lex(3),degrevlex(3)")
ac3eqs = [-2*r12^3*r13^3*r23^3*m1 - 2*r12^3*r13^3*r23^3*m2 - 2*r12^3*r13^3*r23^3*m3 + r12^3*r13^3*m3 + r12^3*r23^3*m3 - r12*r13^5*m3 + r12*r13^3*r23^2*m3 + r12*r13^2*r23^3*m3 - r12*r23^5*m3 + 2*r13^3*r23^3*m1 + 2*r13^3*r23^3*m2, -r12^5*r13*m2 - 2*r12^3*r13^3*r23^3*m1 - 2*r12^3*r13^3*r23^3*m2 - 2*r12^3*r13^3*r23^3*m3 + r12^3*r13^3*m2 + r12^3*r13*r23^2*m2 + 2*r12^3*r23^3*m1 + 2*r12^3*r23^3*m3 + r12^2*r13*r23^3*m2 + r13^3*r23^3*m2 - r13*r23^5*m2, -r12^5*r23*m1 - 2*r12^3*r13^3*r23^3*m1 - 2*r12^3*r13^3*r23^3*m2 - 2*r12^3*r13^3*r23^3*m3 + 2*r12^3*r13^3*m2 + 2*r12^3*r13^3*m3 + r12^3*r13^2*r23*m1 + r12^3*r23^3*m1 + r12^2*r13^3*r23*m1 - r13^5*r23*m1 + r13^3*r23^3*m1]
ac3eqm = [x.subs({m1:1,m2:1,m3:1}).factor() for x in ac3eqs]
ThreeBodTest = phc.blackbox(ac3eqm,R7blk)
ThreeBodSols = ThreeBodTest.solutions(get_failures=True)
ThreeBodSols.sort()
ThreeBodTriangs = []
maxrsq = max([x[0]^2+x[1]^2 for x in sol for sol in ThreeBodSols])
minrsq = min([x[0]^2+x[1]^2 for x in sol for sol in ThreeBodSols])
for index in range(len(ThreeBodSols)):
    sol = ThreeBodSols[index]
    sol_maxrsq = max([x[0]^2+x[1]^2 for x in sol])
    hue_index = float(sol_maxrsq-minrsq)/(maxrsq-minrsq)
    ThreeBodTriangs.append(line([sol[0],sol[1],sol[2],sol[0]],hue=hue_index))
show(sum(ThreeBodTriangs)+circle((0,0),1),frame=True,figsize=[6,6])