# Draw curves in the Hesse pencil x^3 + y^3 + z^3 = 3kxyz
#
# To show off the curve symmetry we use symmetric projective coordinates:
# lines x=0, y=0, and z=0 form an equilateral triangle
# with 0:0:1, 0:1:0, and 1:0:0 as vertices and 1:1:1 at its center.
# For this view the line at infinity is x+y+z = 0.

from PADS.SVG import *
from math import cos,sin,pi
import sys

o = 200+200j
bbox = 400+400j
r = 45
a = o + 0.5j*r + 3**0.5*r/2
b = o + 0.5j*r - 3**0.5*r/2
c = o + -1j*r

def p2c(p): # projective to complex, unused
    px,py,pz = p
    return (px*a+py*b+pz*c)/(px+py+pz)

def c2p(q): # complex to projective
    q = (q-o)/r     # normalize
    z = (-q.imag * 2 + 1)/3
    x = q.real/3**0.5 + (1-z)/2
    y = -q.real/3**0.5 + (1-z)/2
    return (x,y,z)

def between(p,q,s): return (1-s)*p+s*q

def root(p,q,f): # find root of function by binary search
    fp = f(p)
    fq = f(q)
    if fp > fq:
        p,q,fp,fq = q,p,fq,fp
    if fp > 0: return p
    if fq < 0: return q
    for i in range(40):
        mid = (p+q)/2
        if f(mid) < 0:
            p = mid
        else: q = mid
    return between(p,q,mid)

svg = SVG(bbox,sys.stdout)

def hesse(k):
    def thefun(p):
        x,y,z = c2p(p)
        return x**3+y**3+z**3-3*k*x*y*z
    return thefun

def onside(p,q,n):
    return [between(p,q,j/n) for j in range(n)]

def beyond(p,q,n):
    return [between(p,q,1-5**(j/n)) for j in range(1,n)]

# degenerate case, k=infinity, three black lines
svg.group(fill=colors.none,stroke=colors.black)
for p,q in [(a,b),(b,c),(c,a)]:
    svg.segment(6*p-5*q,6*q-5*p)
svg.ungroup()

# other degenerate case, k=1, isolated point (+ line @ infty not shown)
svg.circle(o,3,fill=colors.black,stroke=colors.none)

# case k > 1 has two lobes, the outer one with three real components
def innerlobe(k,pc):
    h = hesse(k)
    curve = []
    for i in 0,1,2:
        p,q = [a,b,c][i-1], [a,b,c][i]
        for tripoint in onside(p,q,9):
            curve.append(root(o,tripoint,h))
    curve = curve + curve[:3]
    svg.polycurve(curve,opacity=pc)

def outerlobe(k,pc):
    h = hesse(k)
    for i in 0,1,2:
        p,q,r = [a,b,c][i],[a,b,c][i-1],[a,b,c][i-2]
        onlobe = root(p,8*p-7*o,h)
        focus = 2*onlobe-p
        wedge = list(reversed(beyond(p,q,20)))+[p]+beyond(p,r,20)
        svg.polycurve([root(focus,w,h) for w in wedge],opacity=pc)

# k >> 1, draw both lobes
svg.group(fill=colors.none,stroke=colors.blue)
# outerlobe(1,"100%")
inners = [3,4,6,12,36]
opacities = ["100%", "85%","70%","55%", "40%"]
for i in range(5):
    innerlobe(inners[i],opacities[i])
    outerlobe(inners[i],opacities[i])
svg.ungroup()

# k > 1 but close to 1, only inner lobe visibl
svg.group(fill=colors.none,stroke=colors.green)
for i in range(4):
    innerlobe([2,1.5,1.25,1.1][i],["40%","60%","80%","100%"][i])
svg.ungroup()

# k < 1, one lobe
def singlelobe(k,pc):
    h = hesse(k)
    for i in 0,1,2:
        p,q,r = [a,b,c][i],[a,b,c][i-1],[a,b,c][i-2]
        midpoint = (p+q)/2
        focus = 15*midpoint-14*o
        border = list(reversed(beyond(p,r,12)))+onside(p,q,12)+[q]+beyond(q,r,12)
        svg.polycurve([root(focus,w,h) for w in border],opacity=pc)

svg.group(fill=colors.none,stroke=colors.red)
outers = [0,-0.5,-1,-2,-8]
opacities = ["100%", "85%","70%","55%", "40%"]
for i in range(5):
    singlelobe(outers[i],opacities[i])
svg.ungroup()
svg.close()