#!/usr/bin/python

import sys
import optparse
from math import *

usage='usage: %prog [--smooth] <level>\nDraw levels 0 through <level> of the dragon curve as an .obj and write the result to stdout.'
parser=optparse.OptionParser(usage)
parser.add_option('-s','--smooth',type='int',help='subdivide to given level')
options,args=parser.parse_args()
if len(args)!=1:
    parser.error('expected <level> argument')
else:
    toplevel=int(args[0])

s=[]
patterns={}
for i in range(toplevel):
    patterns[i]=s
    s=s+[True]+[not f for f in reversed(s)]

if options.smooth:
    corner_shift=0.
else:
    corner_shift=1./20

def Vertices(level):
    x,y=0.,0.
    ds=.5**(level/2.)
    dx,dy=ds,0.
    theta=-2*pi*level/8
    dx,dy=cos(theta)*dx-sin(theta)*dy,sin(theta)*dx+cos(theta)*dy
    vertices=[(x,y)]
    for f in patterns[level]:
        x+=dx
        y+=dy
        odx,ody=dx,dy
        if f:
            dx,dy=-dy,dx
        else:
            dx,dy=dy,-dx
        vertices.append((x+corner_shift*(dx-odx),y+corner_shift*(dy-ody)))
    vertices.append((x+dx,y+dy))
    return vertices

# heights
z=0.
heights=[]
for level in range(toplevel):
    heights.append(z)
    z+=.25*.8**level

vertices=[]
for level in range(toplevel):
    vertices.append(Vertices(level))

# write obj header
print '''\
# simple obj file format:
#   # vertex at coordinates (x,y,z)
#   v x y z
#   # triangle with vertices a,b,c
#   f a b c
#   # vertices are indexed starting from 1'''

if not options.smooth: # simple version drawing sharp corners

    # generate vertices
    bases={0:0}
    for level in range(toplevel):
        v=vertices[level]
        z=heights[level]
        for k in range(len(v)):
            x,y=v[k]
            print 'v %g %g %g'%(x,y,z)
            if k<len(v)-1:
                xn,yn=v[k+1]
                print 'v %g %g %g'%((x+xn)/2,(y+yn)/2,z)
        bases[level+1]=bases[level]+2*len(v)-1

    # generate triangles
    for level in range(toplevel-1):
        n=len(vertices[level+1])-1
        lo=bases[level]
        hi=bases[level+1]
        for k in range(n):
            a=lo+k+1
            b=hi+k+k+1
            print 'f %d %d %d %d'%(a,a+1,b+2,b)

else: # complicated version drawing smooth subdivided surface
    X,Y=[[[p[i] for p in v] for v in vertices] for i in 0,1]

    def Subdivide(x):
        result=[x[0]]
        result.append(.5*(x[0]+x[1]))
        for i in range(1,len(x)-1):
            result.append(.75*x[i]+.125*(x[i-1]+x[i+1]))
            result.append(.5*(x[i]+x[i+1]))
        result.append(x[-1])
        return result

    # subdivide until each level has the same number of vertices
    t_subdivisions=(len(X[-1])-1)*2**options.smooth
    for level in range(toplevel):
        for _ in range(options.smooth+toplevel-level-1):
            X[level]=Subdivide(X[level])
            Y[level]=Subdivide(Y[level])
        assert len(X[level])==t_subdivisions+1

    # generate vertices by subdividing in the height direction for each strip
    s_subdivisions=(toplevel-1)*2**options.smooth
    for _ in range(options.smooth):
        heights=Subdivide(heights) 
    assert len(heights)==s_subdivisions+1
    for i in range(t_subdivisions+1):
        x=[X[k][i] for k in range(toplevel)]
        y=[Y[k][i] for k in range(toplevel)]
        for _ in range(options.smooth):
            x=Subdivide(x)
            y=Subdivide(y)
        assert len(x)==s_subdivisions+1
        for j in range(s_subdivisions+1):
            print 'v %g %g %g'%(x[j],y[j],heights[j])

    # generate quads
    for i in range(t_subdivisions):
        for j in range(s_subdivisions):
            a=(s_subdivisions+1)*i+j+1
            b,c,d=a+1,a+s_subdivisions+2,a+s_subdivisions+1
            print 'f %d %d %d %d'%(a,b,c,d)