Welcome to Fractal Forums

// mi is the iterations number
// imax is the iterations number for the current point
// rxy is the z axis trap
// expo is the exponent (that is, "p" in Zn+1=Zn^p+C)
//First get the DE of (cx,cy,cz)  
      imax=0;rxy=0;//declarations
      d2M=mandelbulb(cx,cy,cz,rxy,imax);
      if(imax==mi) d2M=0;//(cx,cy,cz) is inside the mandelbulb
      is=1-1/expo;//a magic number.
      d2M=d2M^is;//a magic operation
//////////////////////////////////////////////////////////////////////////////////////////////////
mandelbulb(x,y,z,&kz,&imax){// returns: the DE, kz:correction to the "stars", imax:the number of iterations actually used
    x0=x;y0=y;z0=z;                     
    dr=1;                                 
    i=0;
    rt=10000;// for (0,0,0) trap
    rz=10000;// for z axis trap
    while(i<mi) {
      r2=x*x+y*y+z*z;
      r =sqrt(r2);
      
      rt=min(rt,r);
      rz=min(rz,x*x+y*y);
      
      if(r>2^10) break;

      zo0 = asin( z/r );
      zi0 = atan2( y,x );

      zr = pow( r, expo-1 );
      zo = zo0 * expo;
      zi = zi0 * expo;
    
      dr = zr*dr*expo + 1;
      zr *= r;
      x = x0 + zr*cos(zo)*cos(zi);
      y = y0 + zr*cos(zo)*sin(zi);
      z = z0 + zr*sin(zo);

      i++;
    }
    imax=i;
    kz=1-(1/(16*rz+1))*(1-1/((2*rt)^expo+1));//this will be used in the ray marcher
    return 0.5*log(r)*r/dr;
}
//////////////////////////////////////////////////////////////////////////////////////////////////
juliabulb(x,y,z,&kz,&imax){
    x0=cx;y0=cy;z0=cz;
    dr=1;                                
    i=0;
    rt=10000;
    rz=10000;
    while(i<mi) {
      r2=x*x+y*y+z*z;
      r =sqrt(r2);
      rt=min(rt,r);
      rz=min(rz,x*x+y*y);
      
      if(r>2^10) break;

      zo0 = asin( z/r );
      zi0 = atan2( y,x );

      zr = pow( r, expo-1 );
      zo = zo0 * expo;
      zi = zi0 * expo;
    
      dr = zr*dr*expo ;
      zr *= r;
      x = x0 + zr*cos(zo)*cos(zi);
      y = y0 + zr*cos(zo)*sin(zi);
      z = z0 + zr*sin(zo);

      i++;
    }
    imax=i;
    dist=0.5*log(r)*r/(dr+2^-53);//Estimated Distance
    k=1-d2M/(rt^(expo-1)+d2M);//Correction to the estimated distance due to (0,0,0) trap
    kz=1-(1/(16*rz+1))*(1-1/((2*rt)^expo+1));//Correction to the estimated distance due to the "z" axis trap. used in the ray marcher.
    return dist*k;//This scaling can be done here 
}
/////////////////////////////////////////////////////////////////////////////////////
// The ray marcher: returns the point on the ray ((ox,oy,oz),(vx,vy,vz)) where the intersection occures 
nextd0(ox,oy,oz,vx,vy,vz,tmin,tmax)
{   
   ldist=min(0.05,max(errmax,param3*2^-8));//Don't care. ldist is "targeted" distance to the fractal.
   rxy=0;imax=0;
   t=tmin;
      for(i=0;i<maxiter;i++){
         d=juliabulb(ox+t*vx,oy+t*vy,oz+t*vz,rxy,imax);//or mandelbulb(x,y,z,rxy,imax)
         if (imax==mi) d=0;
         rxy-=1;rxy*=exp(-20*d^2);rxy+=1;// scale the "star" correction term by a function of d.
         dscale=(0.95*rxy+0.05);//the final scale. here: 0.05 is the minimal step distance. 0.95+0.05=1. 
                                          //One can use ((1-min_step_distance)*rxy+min_step_distance).  
         if(d-t*ldist<errmax*t) {// Reached an intersection: Refine
            d0=ad;d1=d;
            t0=t-ad;t1=t;
            for(j=0;j<5;j++){
               t=0.5*(t0+t1);
               d=juliabulb(ox+t*vx,oy+t*vy,oz+t*vz,rxy,imax);//or mandelbulb(x,y,z,rxy,imax)
               if((d0-t0*ldist)*(d-t*ldist)<0){t1=t;d1=d;}
               else{t0=t;d0=d;}
            }
            t=0.5*(t0+t1);//Why not? but t may be inside!
            break;
         }
         if(t>tmax) {break;}
         t+=d*dscale;//We step forward by the scaled distance
      }
   return t;
}

rxy-=1;rxy*=exp(-20*d^2);rxy+=1;// scale the "star" correction term by a function of d.