--- fnintern.cpp.orig	2009-02-26 00:09:01.000000000 +0100
+++ fnintern.cpp	2010-01-07 18:35:50.000000000 +0100
@@ -223,6 +223,7 @@
 DBL f_noise3d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 76
 DBL f_pattern(FPUContext *ctx, DBL *ptr, unsigned int fn); // 77
 DBL f_noise_generator(FPUContext *ctx, DBL *ptr, unsigned int fn); // 78
+DBL f_mandelbulb(FPUContext *ctx, DBL *ptr, unsigned int fn);
 
 void f_pigment(FPUContext *ctx, DBL *ptr, unsigned int fn, unsigned int sp); // 0
 void f_transform(FPUContext *ctx, DBL *ptr, unsigned int fn, unsigned int sp); // 1
@@ -314,6 +315,7 @@
 	{ f_noise3d,                 0 + 3 }, // 76
 	{ f_pattern,                 0 + 3 }, // 77
 	{ f_noise_generator,         1 + 3 }, // 78
+	{ f_mandelbulb     ,         2 + 3 }, // 79
 	{ NULL, 0 }
 };
 
@@ -325,7 +327,7 @@
 	{ NULL, 0 }
 };
 
-const unsigned int POVFPU_TrapTableSize = 79;
+const unsigned int POVFPU_TrapTableSize = 80;
 const unsigned int POVFPU_TrapSTableSize = 3;
 
 
@@ -1348,4 +1350,51 @@
 	}
 }
 
+DBL f_mandelbulb(FPUContext *ctx, DBL *ptr, unsigned int) // 79
+{
+	// PARAM_X  .. PARAM_Z  - x,y,z function parameters
+	// PARAM(1) - power of n
+	// PARAM(2) - iterations
+	DBL x,y,z,newx, newy, newz, r, theta, phi, rn, thetan, phin, sinthetan, cx,cy,cz, xkw_ykw;
+	int n,cnt,iterations;
+	n = (int)PARAM(0);
+	iterations = PARAM(1);
+	cnt = iterations;
+	cx = PARAM_X;
+	cy = PARAM_Y; 
+	cz = PARAM_Z;
+	x = 0.0;
+	y = 0.0;
+	z = 0.0;
+	for(;;)
+	{
+		if (cnt <= 0)
+			break;
+
+		xkw_ykw = x*x + y*y;
+		r = sqrt(xkw_ykw + z*z);
+
+		if (r >= 2.0)
+			break;
+
+		theta = atan2(sqrt(xkw_ykw), z);
+		phi = atan2(y,x);
+		rn = pow(r,n);
+		thetan = theta * n;
+		phin = phi * n;
+		sinthetan = sin(thetan);
+		newx = cx + rn * sinthetan * (cos(phin));
+		newy = cy + rn * sinthetan * (sin(phin));
+		newz = cz + rn * (cos(thetan));
+		// now we reassign them
+		x = newx;
+		y = newy;
+		z = newz;
+		cnt--;
+	}
+	if (cnt <= 0) { return -1.0; } else {return 0.0; }
+	/* return (1.0*cnt) / iterations; */
+}
+
+
 }