--- fnintern.cpp	Fri Dec 18 10:02:13 2009
+++ /usr/src/povray-3.7.0.beta.34/source/backend/vm/fnintern.cpp	Fri Dec 
18 10:01:10 2009
@@ -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); // 79

  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,              5 + 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;


@@ -1270,6 +1272,43 @@
  	return Noise(Vec, ngen);
  }

+DBL f_mandelbulb(FPUContext *ctx, DBL *ptr, unsigned int) // 79
+{  // Coded by David Wagner,
+	// based on Daniel White's formula at 
http://www.skytopia.com/project/fractal/2mandelbulb.html#formula
+	// and Abram Hindle's patch at 
http://softwareprocess.us/index.cgi/Mandelbulb
+	DBL x = PARAM_X, y = PARAM_Y, z = PARAM_Z;
+	DBL p_r       = PARAM(0);// Radial power.
+	DBL p_theta	  = PARAM(1);// Power around theta (usually the same as p_r).
+	DBL p_phi     = PARAM(2);// Power around phi (usually the same as p_r).
+	int i_bailout = PARAM(3);// Maximum number of iterations.
+	DBL r_bailout = PARAM(4);// Assumed divergence radius, usually 3.
+	int i=1;
+	DBL r2_bailout, r2, rp, ptheta, pphi, rpsinptheta;
+	r2_bailout = pow(r_bailout,2);
+	r2 = pow(x,2) + pow(y,2) + pow(z,2) ;
+	while(i < i_bailout && r2 < r2_bailout) {
+		 ++i;
+		 rp          = pow(r2,0.5*p_r);
+		 ptheta      = p_theta * atan2(sqrt(pow(x,2) + pow(y,2)), z);
+		 pphi	       = p_phi   * atan2(y,x);
+		 rpsinptheta = rp * sin(ptheta);
+		 x = PARAM_X + rpsinptheta * cos(pphi);
+		 y = PARAM_Y + rpsinptheta * sin(pphi);
+		 z = PARAM_Z + rp * cos(ptheta);
+		 r2 = pow(x,2) + pow(y,2) + pow(z,2);
+   }
+   // The return value corresponds to a smoothed inverted iteration count,
+   // for a range from near zero for many iterations,
+   // through 1/3 at three iterations, 1/2 at two, and unity at one 
iteration;
+   // and greater than one if the starting point is already outside the 
assumed divergence radius.
+   if (i < i_bailout) {
+       return 1/( i + log(log(r2_bailout) / log(  r2)) / log(p_r) );
+   }
+   else {
+       return 1/( i + log(log(r2_bailout) / log(1+r2)) / log(p_r) );
+   }
+}
+
  void f_pigment(FPUContext *ctx, DBL *ptr, unsigned int fn, unsigned 
int sp) // 0
  {
  	VECTOR Vec = { PARAM_N_X(5), PARAM_N_Y(5), PARAM_N_Z(5) };

Post a reply to this message