POV Ray 3.7 code
//==========================================================
// 3D quaternion Julia-Mandelbrot-Mandelbrot (C = -0.7454294)
//  set generation by a cellular automaton (POV-Ray 3.7)
//==========================================================

#include "colors.inc"  // Standard Color definitions
#include "glass.inc"   //  Glass pigment effect
#include "metals.inc"   // Metal pigment effect
global_settings {assumed_gamma 1.0}
background { color rgb <00,0.0,0.0> }

camera {
    location <-20, 20, -300>
    look_at  <-5, 0,  0>
}

light_source {
	< -20, 20, -10>
	rgb <1.000000, 1.000000, 1.000000> * 2.0
}

light_source {
	< 30, 20, -10>
	rgb <1.000000, 1.000000, 1.000000> * 3.0
}

#declare R =  .5;      // set radius value
#declare L = 2;   // set square area side
#declare Cx = -0.7454294;
#declare Cz = 0;
#declare n = 150;     // set number of pixels per area side
#declare st = 1;   // set increment step
#declare N = 40;    // set number of cycles
#declare Th = -45;   // set rotation angles
#declare Ph = 30;

#declare K = array[2*n+1][2*n+1][2*n+1];   // escape rate 3D matrix
 #for (p, -n, n, st)
     #declare IncX = p*L/n;   // Julia X increment

    #for (q, -n, n, st)
       #declare IncY = q*L/n;    // Mand Y increment

    #for (r, -n, n, st)
    #declare IncZ = r*L/n;  // Julia Y increment

      #declare X = IncX;   // start JuliaX
      #declare Y = 0;   // start MandY
      #declare Z = 0;   // start MandZ
      #declare K[p+n][q+n][r+n] = 0;

     #for (k,1,N)
       #declare XX = X*X - Y*Y - Z*Z + Cx;  // cycle JuliaX
       #declare YY = 2*X*Y + IncY;  // cycle MandY
       #declare ZZ = 2*X*Z + IncZ;   // cycle JuliaZ
       #declare X = XX;
       #declare Y = YY;
       #declare Z = ZZ;
       #declare W = X*X +Y*Y + Z*Z;

     	#if (W < R)     // escape if
     	#declare K[p+n][q+n][r+n] = k;
  #end  // end if
  #end  // end for k
  #end  // end for q
  #end  // end for p
  #end  // end for r

  #declare Rnd_1 = seed (1153);
 #declare Rnd_2 = seed (553);
 #declare Rnd_3 = seed (876);

 #declare pp = 0;
 #declare qq = 0;
 #declare rr = 0;

// repalce   n*n*n/8   by    n*n*n/8*clock   for animation
  union{#for(j,0,n*n*n/8) // partial values n*n/4, n*n/2, n*n, n*n*n/64,
   n*n*n/8, n*n*n

  #declare p = pp + st*pow(-1,int(n*rand(Rnd_1)));
   #declare q = qq + st*pow(-1,int(n*rand(Rnd_2)));
   #declare r = rr + st*pow(-1,int(n*rand(Rnd_3)));

   #if ( abs(p) < n)
   	#if( abs(q) < n)
   		#if (abs(r) < n)
 #if(K[p+n][q+n][r+n] > 0)
  #declare pp = p;
  #declare qq = q;
  #declare rr = r;
 	 	sphere {
  		  < p, q, r >, 1   //  adding 3d axis
     			 texture {
    			  pigment { color Col_Glass_Yellow }
      		 }
      		 finish { ambient rgb <0.3,0.1,0.1>
     			  diffuse .3
     			 reflection .3
     			 specular 1 }   // plot sphere

    translate < 30, 15,10 >
  rotate < Th, Th, Ph > }
 #end // end if
  #end // end if
  #end // end if
  #end // end if
    #end}  // end for j