## {{{ http://code.activestate.com/recipes/577715/ (r3)
# Random 2D Slice Of 4D Mandelbrot Fractal
# FB - 201105231
"""Random 2D Slice Of 4D Mandelbrot Fractal.
                          Modified by Symion 2011
    Now works with Visual Python v5.40.
    Produce 2D slice of 4D Mandelbrot Fractal and Map it in 3D!

Visual Python Controls:
    Click Left Mouse Key = Navigate
    Press Right Mouse Key = Spin
    Press Both Keys = Zoom
    q = Point size - 1
    w = Point size + 1
    e = Point shape + 1
    r = Scene Alignment
"""
vp_flag = 0

if vp_flag == 0:
    from visual import *
else:
    import random
    from math import *
    from PIL import Image
imgx = 512
imgy = 512
forge = 0
pcc = 0
pii = 0
psize = 3
pshape = 0

if vp_flag == 0:
    scene.width=imgx
    scene.height=imgy
    image = [points(size=psize, shape="square")]
    forward = vector(scene.forward)
    print "Number of Points objects: {0}".format(len(image))
else:
    image = Image.new("RGB", (imgx, imgy))
#
print __doc__
# drawing area (xa < xb & ya < yb)
xa = -2.0
xb = 2.0
ya = -2.0
yb = 2.0
maxIt = 32 # max number of iterations allowed
maxit = maxIt / 2.0
if True:
    # random rotation angles to convert 2d plane to 4d plane
    xy = random.random() * 2.0 * pi
    xz = random.random() * 2.0 * pi
    xw = random.random() * 2.0 * pi
    yz = random.random() * 2.0 * pi
    yw = random.random() * 2.0 * pi
    zw = random.random() * 2.0 * pi

else:
    # default rotation angles
    xy=1.3536589728
    xz=2.30808965705
    xw=3.50029464114
    yz=3.37449518258
    yw=4.23401560176
    zw=2.44695022478

sxy = sin(xy)
cxy = cos(xy)
sxz = sin(xz)
cxz = cos(xz)
sxw = sin(xw)
cxw = cos(xw)
syz = sin(yz)
cyz = cos(yz)
syw = sin(yw)
cyw = cos(yw)
szw = sin(zw)
czw = cos(zw)

origx = (xa + xb) / 2.0
origy = (ya + yb) / 2.0

for ky in range(imgy):
    b = ky * (yb - ya) / (imgy - 1)  + ya
    for kx in range(imgx):
        a = kx * (xb - xa) / (imgx - 1)  + xa
        x = a
        y = b
        z = 0 # c = 0
        w = 0 # d = 0
        # 4d rotation around center of the plane
        x = x - origx
        y = y - origy
        x0 = x * cxy - y * sxy
        y = x * sxy + y * cxy
        x = x0 # xy-plane rotation
        x0 = x * cxz - z * sxz
        z = x * sxz + z * cxz
        x = x0 # xz-plane rotation 
        x0 = x * cxw - z * sxw
        w = x * sxw + z * cxw
        x = x0 # xw-plane rotation
        y0 = y * cyz - z * syz
        z = y * syz + z * cyz
        y = y0 # yz-plane rotation
        y0 = y * cyw - w * syw
        w = y * syw + w * cyw
        y = y0 # yw-plane rotation
        z0 = z * czw - w * szw
        w = z * szw + w * czw
        z = z0 # zw-plane rotation
        x = x + origx
        y = y + origy
        if forge:
            for i in range(maxIt):
                # iteration using quaternion numbers
                x0 = x * x - y * y - z * z - w * w + a
                y = 2.0 * x * y + b
                z = 2.0 * x * z
                w = 2.0 * x * w
                x = x0
                s = x * x + y * y + z * z + w * w # 4d absolute value squared
                if s > 4.0:
                    break
        else:
            for i in range(maxIt):
                # iteration using hyper-complex numbers
                x0 = x * x - y * y - z * z - w * w + a
                y0 = 2.0 * x * y - 2.0 * z * w + b
                z0 = 2.0 * x * z - 2.0 * y * w
                w = 2.0 * x * w + 2.0 * z * y
                x = x0
                y = y0
                z = z0
                s = x * x + y * y + z * z + w * w # 4d absolute value squared
                if s > 4.0:
                    break
        pcc += i
        pii += 1
        if (i%maxIt) != 0:
            if vp_flag == 0:
                if len(image[-1].pos) > 9999:
                    image.append(points(size=psize, shape="square"))
                    print "Number of Points objects: {0}".format(len(image))
                c = (i/maxit, i/maxit, i/maxit)
                image[-1].append((kx-256, ky-256, i/maxit))
                image[-1].color[-1] = c
            else:
                image.putpixel((kx, ky), (i % 4 * 64, i % 8 * 32, i % 16 * 16))
if vp_flag == 0:
    scene.autoscale=False
    scene.range=mag(scene.mouse.camera)/sqrt(3)
    scene.visible = True
else:
    image.save("4D_Mandelbrot_Fractal.png", "PNG")
print "Finished!"
print "{0} / {1} = {2}".format(pcc, pii, pcc / pii)
mess = "Base Set:\nxy={0}, xz={1}, xw={2}, yz={3}, yw={4}, zw={5}"
print mess.format(xy, xz, xw, yz, yw, zw)
if vp_flag == 0:
    while 1:
        if scene.mouse.events>0:
            mk = scene.mouse.getevent()
            if mk.release == "left":
                scene.center = mk.pos
        elif scene.kb.keys:
            km = scene.kb.getkey()
            if km in ["x", "X"]:
                break
            elif km in ["w"]:
                psize = psize%50
                psize += 1
                for a in image:
                    a.size = psize
            elif km in ["q"]:
                psize -= 1
                if psize<1:
                    psize = 50
                for a in image:
                    a.size = psize
            elif km in ["e"]:
                pshape = (pshape+1)%2
                for a in image:
                    a.shape = ["square","round"][pshape]
            elif km in ["r"]:
                scene.forward = forward