This is why I can't sleep at night.

This post in Peggy's journal: http://shatterstripes.livejournal.com/546349.html motivated me to get off my ass, and to do some more work on libfract [or maybe to-be-named GnuFract, but I'm thinking of making the core parts BSD license, or even public domain.] Something that I has been in the back of my head for the last few years, is wondering how difficult it would be to modify POV-RAY to do more than just quaternion Julia sets. Not difficult at all it turns out...

If my reasoning is correct, then literally, the only change to turn the current Julia set code into Mandelbrot set code is:
(Also with the added bonus of being able to permutate the initial value of Z Set parameters to <0,0,0,0> for vanilla M-set.)
(Actually... now that I think about this, this may not actually be correct.)

--- /tmp/povray-3.6.1/source/quatern.cpp        2004-08-02 16:11:37.000000000 -0700
+++ source/quatern.cpp  2006-03-25 07:55:38.000000000 -0800

@@ -137,10 +137,10 @@
 
     tmp = 3.0 * x2 - d;
 
-    Sx[i] = x = x * (x2 - 3.0 * d) + Julia->Julia_Parm[X];
-    Sy[i] = y = y * tmp + Julia->Julia_Parm[Y];
-    Sz[i] = z = z * tmp + Julia->Julia_Parm[Z];
-    Sw[i] = w = w * tmp + Julia->Julia_Parm[T];
+    Sx[i] = x = x * (x2 - 3.0 * d) + Sx[0] + Julia->Julia_Parm[X];
+    Sy[i] = y = y * tmp + Sy[0] + Julia->Julia_Parm[Y];
+    Sz[i] = z = z * tmp + Sz[0] + Julia->Julia_Parm[Z];
+    Sw[i] = w = w * tmp + Sw[0] + Julia->Julia_Parm[T];
   }
 
   return (true);

And also at "@@ -206,10 +202,10 @@", "@@ -305,10 +299,10 @@", and "@@ -401,10 +393,10 @@". Omitted here for brevity.
And a similar set of changes to source/hcmplx.cpp, if you like hypercomplex M-sets.

235     Sx[i] = xw - yz + Sx[0] + HCompl->Julia_Parm[X];
236     Sy[i] = 2.0 * (x * y - z * w) + Sy[0] + HCompl->Julia_Parm[Y];
237     Sz[i] = 2.0 * (x * z - w * y) + Sz[0] + HCompl->Julia_Parm[Z];
238     Sw[i] = 2.0 * (x * w + y * z) + Sw[0] + HCompl->Julia_Parm[T];

But... this may be wrong. The first several plots that I've made, didn't turn out looking like I had expected them to. I need to sit down at some point when I have the time, and work the math out to check.

Actually, now that I've had a full day's sleep; I realize what's wrong with this. For Z0 = 0+0i+0j+0k this is correct, but when Z0 ≠ 0+0i+0j+0k, that value gets added into the sum with every iteration, not just the first one. I'll fix it later when I have time.

I don't have time to write these words that I'm typing right now either, but quick explanation: In a Julia set, you start Z off as 0 + 0i and then add some constant of your choice, c, every iteration. In a Mandelbrot set, you also start Z off as zero, but c is your current point in space that you are plotting. (i.e. When you are at pixel 0+66i, you have c=0+66i, when you move to the next pixel, 0+67i, then c=0+67i. In a Julia set the c never changes.) For fun, you can start Z off with something other than zero, to get a distorted M-set, but you only use that value once. (Explains why the attractor basin shrank so much (orbits are escaping faster). I made this mistake before, a long long time ago in a 2-D {ℂ} program; I thought that this shape looked familiar)

Here's a few of the tests I made this weekend:
(Don't have time to make thumbnails right now.)

( WARNING: 12M of images )