import bpy
import mathutils
import bpy_extras
import bmesh
import random
import math

bpy.ops.object.select_all(action='SELECT')
bpy.ops.object.delete(use_global=False)

# Mandelbox parameters
scale = -1.5
min_radius = 0.5
fixed_radius = 1.0
max_iter = 40
bailout = 8.0

# Grid resolution and bounds
grid_size = 300
height = 300
extent = 2.5  # How far from center the grid reaches
step = (extent * 2) / grid_size
sphere_trim = False

def random_ccmp():
	return random.randrange(0, 1000) / 1000.0

def random_color():
	return (random_ccmp(), random_ccmp(), random_ccmp(), 1.0)

def vector_length(V):
	return math.sqrt(V[0] ** 2 + V[1] ** 2 + V[2] ** 2)

def mandelbox(z):
	C = z.copy()
	for i in range(max_iter):
		# Box fold
		for j in range(3):
			if z[j] > 1.0:
				z[j] = 2.0 - z[j]
			elif z[j] < -1.0:
				z[j] = -2.0 - z[j]
		
		# Sphere fold
		r2 = z.length_squared
		if r2 < min_radius ** 2:
			z *= (fixed_radius ** 2) / (min_radius ** 2)
		elif r2 < fixed_radius ** 2:
			z *= (fixed_radius ** 2) / r2

		# Scale and offset
		z = z * scale + C

		if z.length_squared > bailout:
			return False  # Escapes
	return True  # Bounded

verts = []
faces = []

def add_rp(A, B, C, D, E, F, G, H):
	voff = len(verts)
	verts.extend([A, B, C, D, E, F, G, H])
	faces_z = [(0, 1, 3, 2), (0, 1, 5, 4), (0, 2, 6, 4), (2, 3, 7, 6), (3, 1, 5, 7), (4, 5, 7, 6)]
	faces.extend([[n + voff for n in f] for f in faces_z])

def add_cube(P):
	cube_size = step/2
	A = (P[0] + cube_size, P[1] + cube_size, P[2] + cube_size)
	B = (P[0] - cube_size, P[1] + cube_size, P[2] + cube_size)
	C = (P[0] + cube_size, P[1] + cube_size, P[2] - cube_size)
	D = (P[0] - cube_size, P[1] + cube_size, P[2] - cube_size)
	E = (P[0] + cube_size, P[1] - cube_size, P[2] + cube_size)
	F = (P[0] - cube_size, P[1] - cube_size, P[2] + cube_size)
	G = (P[0] + cube_size, P[1] - cube_size, P[2] - cube_size)
	H = (P[0] - cube_size, P[1] - cube_size, P[2] - cube_size)
	add_rp(A, B, C, D, E, F, G, H)

for i in range(grid_size):
	for j in range(grid_size):
		for l in range(grid_size - height, grid_size):
			x = -extent + i * step
			y = -extent + j * step
			z = -extent + l * step

			is_inside = mandelbox(mathutils.Vector((x, y, z)))
		
			if is_inside and (not sphere_trim or (vector_length((x, y, z)) < 2.0)):
				add_cube((x, y, z))

mat = bpy.data.materials.new(name="VertexColorMaterial")
mat.use_nodes = True
bsdf = mat.node_tree.nodes.get("Principled BSDF")
color_node = mat.node_tree.nodes.new(type='ShaderNodeVertexColor')
color_node.layer_name = "Col"
mat.node_tree.links.new(color_node.outputs['Color'], bsdf.inputs['Base Color'])

mb_mesh = bpy.data.meshes.new(name="Mandelbox Mesh")
mb_mesh.from_pydata(verts, [], faces)
mb_mesh.validate(verbose=True)
mb_object = bpy_extras.object_utils.object_data_add(bpy.context, mb_mesh)
mb_object.data.materials.append(mat)

bpy.ops.object.mode_set(mode='EDIT')
mb_bmesh = bmesh.from_edit_mesh(mb_object.data)
color_layer = mb_bmesh.loops.layers.color.new("Col")
for i, face in enumerate(mb_bmesh.faces):
	if i % 6 == 0:
		d_center = vector_length((face.loops[0].vert.co.x, face.loops[0].vert.co.y, face.loops[0].vert.co.z))
		face_color = (0.5 + math.sin(math.pi/2 + 4.0 * d_center)/2, 0.5 + math.sin(4.0 * d_center)/2, 0.5 + math.sin(1.0 * d_center)/2, 1.0)
	for loop in face.loops:
		loop[color_layer] = face_color
bmesh.update_edit_mesh(mb_object.data)
bpy.ops.object.mode_set(mode='OBJECT')