######################################################
#
#   Fractal Creator
#
#  (c)2005, Carmen DiMichele
#  carmen@dimichele.net
#
######################################################

import Blender
from Blender import NMesh
from Blender.Draw import *
from Blender.BGL import *

import math
from math import *


# Create globals

# Fractal Types:
#	[ sizeX, sizeY, sizeZ, boxX0, boxY0, boxX1, boxY1, maxN, cRe, cIm ]
fsTypes = "Mandelbrot|Julia|Julia of sin(z)"
fTypes = [
	 [6.0, 4.0, 2.0, -2.1, -1.1, 0.6, 1.1, 100, 0.0, 0.0],			# Mandelbrot
	 [6.0, 4.0, 2.0, -1.6, -1.0, 1.6, 1.0, 100, -0.74543, 0.11301],	# Julia
	 [6.0, 4.0, 2.0, pi, -2, 2*pi, 2, 50, 1.0, 0.2]					# Julia, sin z
	]

# Fractal variables
fType = Create(1)
sizeX = Create(fTypes[fType.val-1][0]); sizeY = Create(fTypes[fType.val-1][1]); sizeZ = Create(fTypes[fType.val-1][2])
boxX0 = Create(fTypes[fType.val-1][3]); boxY0 = Create(fTypes[fType.val-1][4])
boxX1 = Create(fTypes[fType.val-1][5]); boxY1 = Create(fTypes[fType.val-1][6])
maxN = Create(fTypes[fType.val-1][7])
zRe = Create(fTypes[fType.val-1][8]); zIm = Create(fTypes[fType.val-1][9])
linkScale = [sizeX.val/(boxX1.val - boxX0.val), sizeY.val/(boxY1.val - boxY0.val), sizeZ.val/maxN.val]
sizeLink = Create(1)
resX = Create(80); resY = Create(60)
smooth = Create(0);
filtBase = Create(0);
filtFlats = Create(1);
invert = Create(0);

# Preview variables
prevMaxWidth = 80; prevMaxHeight = 50
prevX0 = -1; prevY0 = -1; prevWidth = -1; prevHeight = -1
prevOffsX = -1; prevOffsY = -1


######################################################
#
#    Object Math
#
######################################################

# Fractal Type Functions #############################

def iter1(maxN, x, y):
	z = c = complex(x,y)
	for n in range(maxN):
		if abs(z) <= 2.0:
			z = z**2 + c
		else:
			break
	return n

 
def iter2(maxN, c, x, y):
	z = complex(x,y)
	for n in range(maxN):
		if abs(z) <= 2.0:
			z = z**2 + c
		else:
			break
	return n

 
def iter3(maxN, c, x, y):
	z = complex(x,y)
	for n in range(maxN):
		if abs(z.imag) < 50.0:
			z = c*complex(sin(z.real)*cosh(z.imag), cos(z.real)*sinh(z.imag))
		else:
			break
	return n

######################################################

def smoothing(me, res):

	me2 = NMesh.New()
  
	w = res[0]
	h = res[1]
	max = w*h
  
	for y in range(h):
		for x in range(w):
			sum = 0.0
			cells = 0
			a = x + y * w

			for dy in range(-1,2):
				for dx in range(-1,2):
					a1 = a + dx + dy * w
					if a1 >= 0 and a1 < max:
						sum = sum + me.verts[a1].co.z
						cells = cells + 1
      
			me2.verts.append(NMesh.Vert(me.verts[a].co.x, me.verts[a].co.y, sum/cells))

	return me2  

	
def surfFunc(me, res):
	global	filtBase, filtFlats

	for y in range(res[1] - 1):
		for x in range(res[0] - 1):
			a = x + y * res[0]

			#find mid-height of face
			midz = (me.verts[a].co.z + me.verts[a + res[0] + 1].co.z) / 2.0

			#find z-slope of face
			sloz = max(
				abs(me.verts[a].co.z - me.verts[a + res[0] + 1].co.z),
				abs(me.verts[a + res[0]].co.z - me.verts[a + 1].co.z))

			# filter out base
			if (filtBase.val == 1) | (((midz > 5) & (sloz > 0.001)) | (midz > 20)):

				# filter out flats
				if (filtFlats.val == 1) | (sloz > 0.001):

					# draw the face
					f = NMesh.Face()
					f.v.append(me.verts[a])
					f.v.append(me.verts[a + res[0]])
					f.v.append(me.verts[a + res[0] + 1])
					f.v.append(me.verts[a + 1])
					f.smooth=1
					me.faces.append(f)


def fractal():
	global	sizeX, sizeY, sizeZ
	global	resX, resY, smooth
	global	maxN
	global	boxX0, boxY0, boxX1, boxY1
	global	zRe, zIm

	me = NMesh.New()

	dx = (boxX1.val - boxX0.val) / resX.val
	dy = (boxY0.val - boxY1.val) / resY.val
	y = boxY1.val
	
	xmid = (boxX1.val - boxX0.val)/2.0 + boxX0.val
	ymid = (boxY1.val - boxY0.val)/2.0 + boxY0.val

	for yn in range(resY.val):
		x = boxX0.val
		for xn in range(resX.val):

			if fType.val == 1:
				h = iter1(maxN.val, x, y)
			elif fType.val == 2:
				h = iter2(maxN.val, complex(zRe.val, zIm.val), x, y)
			elif fType.val == 3:
				h = iter3(maxN.val, complex(zRe.val, zIm.val), x, y)

			if invert.val: h = maxN.val - h
			me.verts.append(NMesh.Vert(x-xmid, y-ymid, h))
			x = x + dx
		y = y + dy

	for i in range(smooth.val):
		me = smoothing(me, [resX.val, resY.val])
		
	surfFunc(me, [resX.val, resY.val])
	
	obj = NMesh.PutRaw(me)

	#scaling
	obj.SizeX = sizeX.val/(boxX1.val - boxX0.val)
	obj.SizeY = sizeY.val/(boxY1.val - boxY0.val)
	obj.SizeZ = sizeZ.val/maxN.val
		
	obj.setLocation(0.0,0.0,0.0)
	
 

#########################################################
#
#    Preview
#
#########################################################

def updatePrevBox():
	global prevMaxWidth, prevMaxHeight, prevWidth, prevHeight
	global prevOffsX, prevOffsY

	prevHeight = prevMaxHeight
	prevWidth = int((boxX1.val - boxX0.val)*float(prevMaxHeight)/(boxY1.val - boxY0.val))
	if prevWidth<1: prevWidth = 1
	if prevWidth>prevMaxWidth:
		prevWidth = prevMaxWidth
		prevHeight = int((boxY1.val - boxY0.val)*float(prevMaxWidth)/(boxX1.val - boxX0.val))
		if prevHeight<1: prevHeight = 1
	prevOffsX = prevX0 + prevMaxWidth - prevWidth
	prevOffsY = prevY0 + prevMaxHeight

def preview():
	global	boxX0, boxY0, boxX1, boxY1
	global	zRe, zIm
	global	maxN
	global	prevWidth, prevHeight
	global	fType

	glBegin(GL_POINTS)

	dx = (boxX1.val - boxX0.val) / prevWidth
	dy = (boxY0.val - boxY1.val) / prevHeight
	y = boxY1.val
	
	for yn in range(prevHeight):
		x = boxX0.val
		for xn in range(prevWidth):

			if fType.val == 1:
				c = float(iter1(maxN.val, x, y))/float(maxN.val)
			elif fType.val == 2:
				c = float(iter2(maxN.val, complex(zRe.val, zIm.val), x, y))/float(maxN.val)
			elif fType.val == 3:
				c = float(iter3(maxN.val, complex(zRe.val, zIm.val), x, y))/float(maxN.val)
			
			if invert.val: c = 1.0 - c
			glColor3f(c, c, c)
			glVertex2f(prevOffsX + xn, prevOffsY - yn)
			x = x + dx
		y = y + dy

	glEnd()


#########################################################
#
#    User Interface
#
#########################################################

def rect(x, y, w, h, r, g, b):
	glColor3f(r,g,b)
	glBegin(GL_POLYGON)
	glVertex2i(x, y)
	glVertex2i(x + w, y)
	glVertex2i(x + w, y - h)
	glVertex2i(x, y - h)
	glEnd()


def gui():
	global	sizeX, sizeY, sizeZ
	global	resX, resY, smooth
	global	maxN
	global	boxX0, boxY0, boxX1, boxY1
	global	zRe, zIm
	global	fType
	global	sizeLink
	global	prevX0, prevY0, prevMaxHeight
	global	filtBase, filtFlats, invert

	glClear(GL_COLOR_BUFFER_BIT)

	y = 462
	rect(0,y+14,248,y+14,0.3,0.3,0.5)

	glColor3f(0.9,0.9,0.9)
	glRasterPos2d(5, y); Text("Fractal Creator v1.0")
	y = y - 14; glRasterPos2d(5, y); Text("(c)2005,Carmen DiMichele")

	y = y - prevMaxHeight - 4; prevX0 = 165; prevY0 = y-4;
	updatePrevBox()
	
	y = y + 6; fType = Menu(fsTypes, 11, 4, y+22, 154, 18, fType.val, "Fractal type")
	filtBase = Toggle("base", 10, 4, y+5, 50, 15, filtBase.val, "Leave base surface in mesh")
	filtFlats = Toggle("flats", 10, 56, y+5, 50, 15, filtFlats.val, "Leave flat surfaces in mesh")
	invert = Toggle("invert", 10, 108, y+5, 50, 15, invert.val, "Invert height of mesh")

	rect(0,y-19,248,70,0.7,0.5,0.7); glColor3f(0.9,0.9,0.9)
	y = y - 23; glRasterPos2d(4, y); Text("Overall Size")
	y = y - 20; sizeX = Slider("X-size: ", 14, 4, y, 240, 18, sizeX.val, 1, 20, 0, "Overall x-axis size")
	y = y - 20; sizeY = Slider("Y-size: ", 14, 4, y, 240, 18, sizeY.val, 1, 20, 0, "Overall y-axis size")
	y = y - 20; sizeZ = Slider("Z-size: ", 14, 4, y, 240, 18, sizeZ.val, 1, 20, 0, "Overall z-axis size")

	rect(0,y-19,248,113,0.5,0.7,0.7); glColor3f(0.9,0.9,0.9)
	y = y - 23; glRasterPos2d(4, y); Text("Calculation Box")
	sizeLink = Toggle("link", 10, 194, y+3, 50, 14, sizeLink.val, "Link overall size settings to calculation box changes")
	y = y - 20; boxX0 = Slider("Lower X: ", 20, 4, y, 240, 18, boxX0.val, -10, 10, 0, "Lower x calculation box coordinate")
	y = y - 20; boxY0 = Slider("Lower Y: ", 21, 4, y, 240, 18, boxY0.val, -10, 10, 0, "Lower y calculation box coordinate")
	y = y - 20; boxX1 = Slider("Upper X: ", 22, 4, y, 240, 18, boxX1.val, -10, 10, 0, "Upper x calculation box coordinate")
	y = y - 20; boxY1 = Slider("Upper Y: ", 23, 4, y, 240, 18, boxY1.val, -10, 10, 0, "Upper y calculation box coordinate")
	y = y - 23; maxN = Slider("max N: ", 24, 4, y, 240, 18, maxN.val, 1, 1000, 0, "Maximum iterations allowed (height along z-axis, calculation points)")

	rect(0,y-19,248,70,0.5,0.5,0.7); glColor3f(0.9,0.9,0.9)
	y = y - 23; glRasterPos2d(4, y); Text("Resolution")
	y = y - 20; resX = Slider("X-res: ", 12, 4, y, 240, 18, resX.val, 1, 1000, 0, "X-axis resolution (calculation points)")
	y = y - 20; resY = Slider("Y-res: ", 12, 4, y, 240, 18, resY.val, 1, 1000, 0, "Y-axis resolution (calculation points)")
	y = y - 20; smooth = Slider("Smoothing: ", 12, 4, y, 240, 18, smooth.val, 0, 10, 0, "3x3 averaging smoothing iterations")

	rect(0,y-19,248,50,0.7,0.7,0.5); glColor3f(0.9,0.9,0.9)
	y = y - 23; glRasterPos2d(4, y); Text("c Factor")
	y = y - 20; zRe = Slider("Real:      ", 13, 4, y, 240, 18, zRe.val, -10, 10, 0, "Real part of complex c factor")
	y = y - 20; zIm = Slider("Imaginary: ", 13, 4, y, 240, 18, zIm.val, -10, 10, 0, "Imaginary part of complex c factor")

	y = y - 40; Button("Create", 3, 4, y, 80, 19)
	y = y - 0; Button("Exit", 1, 164, y, 80, 19)

	preview()


def zoom(win, mous, zm):
	global prevOffsX, prevOffsY, prevWidth, prevHeight
	global prevWidth, prevHeight
	global boxX0, boxX1, boxY0, boxY1

	x = mous[0]-win[0]-prevOffsX
	if 0 <= x <= prevWidth:
		y = mous[1]-win[1]-prevOffsY+prevHeight
		if 0 <= y <= prevHeight:
			w = (boxX1.val - boxX0.val) / zm
			h = (boxY1.val - boxY0.val) / zm
			boxX0.val = boxX0.val + (float(x)/float(prevWidth) * (boxX1.val - boxX0.val)) - w/2
			boxX1.val = boxX0.val + w
			boxY0.val = boxY0.val + (float(y)/float(prevHeight) * (boxY1.val - boxY0.val)) - h/2
			boxY1.val = boxY0.val + h
	

def event(evt, val):
	if evt==ESCKEY and not val:
		Exit()
	
	elif evt==LEFTMOUSE and val:
		zoom( Blender.Window.GetScreenInfo(14)[0].get('vertices'),
			  Blender.Window.GetMouseCoords(),
			  1.25)
		Redraw(1)
		return
				
	elif evt==RIGHTMOUSE and val:
		zoom( Blender.Window.GetScreenInfo(14)[0].get('vertices'),
			  Blender.Window.GetMouseCoords(),
			  1.0/1.25)
		Redraw(1)
		return


def bevent(evt):
	if   evt==1: Exit()
	elif evt==3: fractal()

	# Fractal Type Menu
	elif evt==11:
		sizeX.val = fTypes[fType.val-1][0]; sizeY.val = fTypes[fType.val-1][1]; sizeZ.val = fTypes[fType.val-1][2]
		boxX0.val = fTypes[fType.val-1][3]; boxY0.val = fTypes[fType.val-1][4]
		boxX1.val = fTypes[fType.val-1][5]; boxY1.val = fTypes[fType.val-1][6]
		maxN.val = fTypes[fType.val-1][7]
		zRe.val = fTypes[fType.val-1][8]; zIm.val = fTypes[fType.val-1][9]
		
	# Calculation Box Sliders                
	elif evt==20:
		if boxX0.val > boxX1.val: boxX1.val = boxX0.val + 0.01
		if sizeLink.val == 1:
			sizeX.val = linkScale[0] * (boxX1.val - boxX0.val)
		updatePrevBox();

	elif evt==21:
		if boxY0.val > boxY1.val: boxY1.val = boxY0.val + 0.01
		if sizeLink.val == 1:
			sizeY.val = linkScale[1] * (boxY1.val - boxY0.val)
		updatePrevBox();

	elif evt==22:
		if boxX1.val < boxX0.val: boxX0.val = boxX1.val - 0.01
		if sizeLink.val == 1:
			sizeX.val = linkScale[0] * (boxX1.val - boxX0.val)
		updatePrevBox();

	elif evt==23:
		if boxY1.val < boxY0.val: boxY0.val = boxY1.val - 0.01
		if sizeLink.val == 1:
			sizeY.val = linkScale[1] * (boxY1.val - boxY0.val)
		updatePrevBox();

	# MaxN Slider
	elif evt==24:
		if sizeLink.val == 1:
			sizeZ.val = linkScale[2] * maxN.val
		updatePrevBox();

	Redraw()

Register(gui, event, bevent)