Welcome to Fractal Forums

// A triplex number is a '3D' number
triplex multiply(triplex a, triplex b) {
    triplex n;
    double pi=3.14159265;
    double r    = pow(a.x*b.x + a.y*b.y + a.z*b.z , 0.5);
    double yang = atan2( pow(a.x*b.x + a.y*b.y, 0.5) ,      a.z  )  ;
    double zang = atan2(a.y , b.x);
    n.x = (r*r) * sin( yang*2 + 0.5*pi ) * cos(zang*2 +pi);
    n.y = (r*r) * sin( yang*2 + 0.5*pi ) * sin(zang*2 +pi);
    n.z = (r*r) * cos( yang*2 + 0.5*pi );
	return n;
}

MMF-4DMinibrot {
::cRN/Qjn29l5WvtOuRc83DQ+OY4+SLOnjOSzYJZvLIQ3tBFoAnsAb3+cDolps1e0tKJHHnP9d
  ouYzLTO5hEyhD/zRD5wfBS5dysBZ5P94DrWNUMUqErf+5/5qvsazTPXUXsvrZY9qLFHGOJgw
  wVnUFHPNIiiDXVKvq66FRry7kVqeBmq1YsjA2tKrTdoYoXs+J5rFHW9s87F1/cyX3+VSmdrf
  8hxpPusZy2himax6fVm99jdNnrPseVTrMrY4qIiW0K1wpmDiqzlDFty++HfoS22WUfca6q6B
  VnI8rkjyj1i4wHfIvpj8WOOel8tCtDQ4qWVX2JV23FN55ryLKV16gddVVew58q1avJx6uOmD
  G+yf0UWcAf63PLp5XThY/6VtvkPlxEr/lStZ5QxrqV/v/L8psRBafp/is96bi6m52vb0+ip9
  rG21tvP9LvPNUeZR7Ur9nPKuq6pWZN1/5kthOZdvx8KGyOJW/2nv+53/8Fdo2epTgUeZsV0Y
  Kq9FKxkX01PcbiVy6DXFPT/WVq3v1rmsos58gYjeVe7q4LRBpJjze4dRYA9XKhr/dkAmE51C
  1FxmAtxcREMFdZ0eOdAhi5zDldTBdpKfg6qlMcaqd6jUaTRjTnaNsf+heoaxVtgNtWelv/2j
  QeXT9ghr7pl2y3+MJd0+WX9ZmOlOhMN/bmok0oPVNvq0CMGGUOR3eO9feg6P63t5SGg5JqXr
  2TKS8J370zMMQnEmsPv0zZ05tOjRMjgZlWCCzuZ0JDZZpeVm04gag20Eb0ztvS2tEv1N7Lql
  dzn1o2HUtUx84WU2JZ9R19DCT1Z/L6YtUXROtJ2+qsc5sju7yT/9Rh7JyZHADHGu2qoTXHrL
  GOfQRGP3rolVfTzvUWuasZ/cpDN0wpOV/pmyDznxv1NMIWHjF10KukR1Foi/ziLTaMKY0Y5+
  UbQsO8v/XC/23wE6vxpJhrXGCnjT6+Aalvvn+GVk1PMP4VzOvP3hae5ezsGqcpT8N9JZ6ep5
  Dl7rimvcpc/PR3o8CdjR+vWteesgXf5oqW1Vkpvl7+cI70kkkv7rmPyOauq5Alzeq41VRJ30
  oXpOQbm4m4EDFGvXeut++KK4TDTvZqSSHJqaFRGTR+KtrTHG+N6qTZ5i1ix7W/WRtS2tY7QR
  myaqHOXVdVYEndtnmv/S3X9W79o5osqSevud0/ilN3xulNXI/v3/ElUbJHSu9TqxkLVvqK7F
  bvlPO1choD9Grfep8Y49ClbJlQnsSXf494kSPm9qJFM3i0aGNXhdTwIPBjsEMyOhXTKAw9HU
  tmgrmgnmgtKVtgVcS9iCALRRXRxZRtOxgWxKaJKKQMI1S0NMZ0NeB7GLR3YJK1LJYbasTaN2
  Naj9UN2S1Y30asA3Fg2hbirqJeqm4mYTsC3EBsNIyOzmykES9EO1KcTtUNVEHkAu5gtuR7WP
  R3KADR1jHbnF2Ki2Fks1ONszV4dcHF2Zp9OzqTS4dCMJIOxOTExVhF5XiRmM1e0DnAnstNYj
  z+Xkf1WEXoTW3Yqv2JwM6JLQcAG5lzj8q9i8L+ITbsPkM6k9KQlg0/JkbCC5SQo/Kg3vMTLv
  2D7gVvEUBZYwO3c0GvngNs5Idln5ao9yMe1rwGBEslJL5VeG5XfSmQwJNp9KEjsXkYxGIIy8
  OA4HQMhPgYC8ETgjYCsETwgYC+ETghYCsETghYCuETwmYCOETwkYC+ETwmYCOET4HSMBHiJ4
  TMBGiJwQMBTiJYRMBXiJ4SMBGiJYSMBXiJ4SMBXiJwQMBHiJYRMBPiJ4SMBeiJYSMBLiJ4RM
  BGiJwQMBTiJYRMBWiJ4SMBGiJYSMBXiJ4TMBXiJwQMBHiJYRMBfiJwQMBGiJYSMBLiJwRMBX
  iJwQMBTiJ4TMBGiJ4SMBeiJYSMBHiJwQMBOiJwRMBLiJwQMBGiJ4RMhPgYCWETwjYCfAxE8I
  mAHxEcJmgHxE4ImAHxE4ImgFxEYJmAHxE8ImwHQMBbiJwQMBeiJ4RMBOiJ4RMBOiJ8BET8HQ
  MxPgYi8ETkjYisET0gYi+ETkhYisETkhYiuET0mYiOET0kYi+ET0mYiOET8HSMRHiJ6TMRGi
  JyQMRTiJaRMRXiJ6SMRGiJaSMRXiJ6SMRXiJyQMRHiJaRMRPiJ6SMReiJaSMRLiJ6RMRGiJy
  QMRTiJaRMRWiJ6SMRGiJaSMRXiJ6TMRXiJyQMRHiJaRMRfiJyQMRGiJaSMRLiJyRMRXiJyQM
  RTiJ6TMRGiJ6SMReiJaSMRHiJyQMROiJyRMRLiJyQMRGiJ6RMxPgYiWET0jYifAxE9ImIHxE
  dJmoHxE5ImIHxE5ImoFxEZJmIHxE9Im4HQMRbiJyQMReiJ6RMROiJ6RMROiJaTMJjF19FaA3
  f9tPf9vN/OnlV7LU1jvketLHKyzP3r0HZiHf18qszlyOqPO3tv4dFVRt86/L1Uup4tc8jj9P
  +0/+T/xJ5hmL9zvH/J7G+of90TTfsLoz0hL9QRy9Obmpz9jCOWW31Ms8C3zfJvOSQPS1DP+w
  0D30XyT/lhyJqa5EV15bdlVu24DdJ/SeTH+ElOaOP4LRdTt6xHO2JPoTTjj1X10Mca8BS/9E
  ej4y63AenYTSMMbKi2smMGFtNJc2KsdEsONQ8udxR7imHCTxZ7b3FmGGmuYf72Fhi3GvLNa8
  g28Y72tMWSaaKERHlm/yhmBKxpXizlvrIETu+/xgBG9y
}

 
Assume position on ray is "vp+alpha*d" where vp is the viewpoint and d is our normalised direction vector for the viewing ray.
We want to scan along the ray segment from a start point dictated by our front cutting plane giving us a starting value for alpha to our end point
dictated by the back cutting plane giving us an end value for alpha.
We use an array "maxd[maxiter]" for storing maximum distances to step at each iteration level.
Proceed as follows:
 
For each ray:
  alpha = starting alpha (as dictated by front clip)
  Set all maxd[n] to a maximum step distance value (usually a user parameter) - n from 0 to maxiter
  step = 0
  done = false
  binary search = false
  alpha found = end alpha
  repeat
     compute point p = vp + alpha*d
     oldstep = step
     iterate point p (set i to number of iterations done here and remember final values for computing smooth iteration)
     if not reached maxiter
       compute point p = vp + (alpha+abit)*d  ;----- where abit is a user parameter, say 1e-10
       Iterate point p
       if reached maxiter
          if oldstep==0
            alpha found = current alpha
            done = true
         else
           step = oldstep/2
           binary search = true
           alpha = alpha - step
         endif
       else
          compute smooth iteration value for first iteration s0
          compute smooth iteration value for second iteration s1
          compute directional DE = 1.0/(1.0+(1/abit)*abs(s1-s0))
          if binary search
            step = oldstep/2
            if DE>solid threshold ;--------- solid threshold is a user parameter, say 1e-4
               alpha = alpha + step
            else
               alpha = alpha - step
            endif
            if step<min accuracy ;--- min accuracy is a user parameter, say 1e-8
               alpha found = current alpha
               done = true
            endif
          elseif directional DE < solid threshold ;--------- solid threshold is a user parameter, say 1e-4
             step = oldstep/2
             binary search = true
             alpha = alpha - step
          else
             if maxd[i]<directional DE
               directional DE = maxd[i]
             else
                maxd[i] = directional DE
             endif
             if directional DE < maxd[i+1] ;***********ADDED
                maxd[i+1] = directional DE
             endif
             alpha = alpha + (step=directional DE/scale value) ;---- I use scale value = @scale*2.5 where @scale is a user parameter
;                                                                                                               default 1.0
           endif
       endif
     elseif oldstep==0
        alpha found = current alpha
        done = true
     else
       step = oldstep/2
       binary search = true
       alpha = alpha - step
     endif
  until  done || alpha>=end position

double r = sqrt(xx*xx+yy*yy+zz*zz);
xx = xx/r;
yy = yy/r;
zz = zz/r;
/* --- rotation - axis Z  */
//vector normalization:					        
r1 = sqrt(xx*xx+yy*yy);
xx = xx/r1;
yy = yy/r1;
zz = zz/r1;
//complex numbers calculation
newx1 = xx*xx - yy*yy;
newy1 = 2.0*xx*yy;
newz1 = zz;
						        
//---- rotation XY ------
r2 = sqrt(newx1*newx1+newy1*newy1);
//normalization:
r3 = sqrt(r2*r2+newz1*newz1);
newx1 = newx1/r3;
newy1 = newy1/r3;
newz1 = newz1/r3;
r2 = r2/r3;
//complex numbers calculation
newr2 = r2*r2 - newz1*newz1;
newz2 = 2.0*r2*newz1;
//spread result for 2 vectors (x and y)
rr = newr2/r2;
newx2 = newx1*rr;
newy2 = newy1*rr;        

//multiplying by normalization factors						        
m=r*r*r1*r3;
newx = newx2*m;
newy = newy2*m;
newz = newz2*m;
													        
xx = newx + aa;
yy = newy + bb;
zz = newz + cc;