#!/usr/bin/python
#########################################################################
# Author: Kai Ren
# Created Time: 2012-10-07 23:24:40
#########################################################################
import numpy as np

def gen_koch_surface(s0, x, y, z, l, depth):
    if depth == 0:
        # Finish iterations, and then output the current surface
        xf = []
        yf = []
        zf = []
        l *= 3
        for k in [0, 1, 3, 2, 0]:
            i = k & 1
            j = (k & 10) >> 1
            pt = s0 + (i * l) * x + (j * l) * y
            xf.append(pt[0])
            yf.append(pt[1])
            zf.append(pt[2])
        surfaces.append((xf[:-1], yf[:-1], zf[:-1]))
        return

    #Generate new facets from the surface (s0, x, y, z)
    for i in range(3):
        xl = float(i) * l
        for j in range(3):
            yl = float(j) * l
            if i == 1 and j == 1:
                new_s0 = s0 + xl * x  + yl * y + l * z
            else:
                new_s0 = s0 + xl * x + yl * y
            gen_koch_surface(new_s0, x, y, z, l/3.0, depth-1)
   
    gen_koch_surface(s0+l*x+l*y, x, z, -y, l/3.0, depth-1)
    gen_koch_surface(s0+2.0*l*x+l*y, y, z, x, l/3.0, depth-1)
    gen_koch_surface(s0+2.0*l*x+2.0*l*y, -x, z, y, l/3.0, depth-1)
    gen_koch_surface(s0+l*x+2.0*l*y, -y, z, x, l/3.0, depth-1)

koch_no = 4  # Number of iterations to generate Koch surface

#The intial plane for Koch surface
s0 = np.array([0., 0., 0.])
xa = np.array([1., 0., 0.])
ya = np.array([0., 1., 0.])
za = np.array([0., 0., 1.])
surfaces = []
gen_koch_surface(s0, xa, ya, za, 1/3.0, koch_no)

#Output the final surface
for surf in surfaces:
    for i in range(len(surf[0])):
        print surf[0][i], surf[1][i], surf[2][i],
    print