Welcome to Fractal Forums

 double rad1 = 1.0;
 double rad2 = 0.2;
 double nSect = 2*pi/9;
 double fact = 3;

 norm = x*x+y*y+z*z;
 while ((norm<16)&&(iter<imax)) {

   (*nfunc)++;

   r1 = x*x + y*y;
   flag = 0;
   if (r1 < (rad1+rad2)*(rad1+rad2)) {
    r = sqrt(r1);

    c1 = rad1 * x/r;
    s1 = rad1 * y/r;

    x1 = (x - c1);
    y1 = (y - s1);
    z1 = z;

    r2 = x1*x1+y1*y1+z1*z1;
    if (r2 < rad2*rad2) {
     //maxnorm = (iter % 5); maxnorm /= 4.0;iter=imax+1;
     maxnorm = (theta2 + pi) / (2.0*pi) ;iter=imax+1;
     flag = 1;
    }
   }

   if (!flag) {

    theta = atan2(y,x);
    theta2 = nSect * round (theta / nSect);

    c1 = cos(theta2);
    s1 = sin(theta2);

    x1 = c1 * x + s1 * y;
    y1 = -s1 * x + c1 * y;
    z1 = z;

    x1 = x1 - rad1;

    x=fact*x1;y=fact*z1;z=fact*y1;

   }
   norm = x*x+y*y + z*z;
}

        double rad1 = 1.0;
	double rad2 = 0.1;
	double nSect = 2 * M_PI / 4;
	double fact = 3;
	double theta2 = 0.0;

	double r = sqrt(z.x * z.x + z.y * z.y);

	double c1 = rad1 * z.x / r;
	double s1 = rad1 * z.y / r;

	CVector3 z1 = CVector3(z.x - c1, z.y - s1, z.z);

	double r2 = z1.Length();

	if (r2 > rad2)
	{

		double theta = atan2(z.y, z.x);
		theta2 = nSect * round(theta / nSect);

		double c1 = cos(theta2);
		double s1 = sin(theta2);

		double x1 = c1 * z.x + s1 * z.y;
		double y1 = -s1 * z.x + c1 * z.y;
		double z1 = z.z;

		x1 = x1 - rad1;

		z.x = fact * x1;
		z.y = fact * z1;
		z.z = fact * y1;
	}
	else
	{
		?????  WHAT SHOULD BE HERE ?????
	}

	double R = sqrt(z.x * z.x + z.y * z.y);
	double R2 = rad1 - R;
	double t = R2 * R2 + z.z * z.z - rad2 * rad2;

double R = sqrt(z.x * z.x + z.y * z.y);
double R2 = rad1 - R;
double t = R2 * R2 + z.z * z.z - rad2 * rad2;
if (t < 0.03)
{
   z /= t;
}

	
        double rad1 = 1.0;
	double rad2 = 0.1;
	double nSect = 2 * M_PI / 9;
	double fact = 3;
	double theta2 = 0.0;

	double R = sqrt(z.x * z.x + z.y * z.y);
	double R2 = rad1 - R;
	double t = R2 * R2 + z.z * z.z - rad2 * rad2;

	if (t > 0.03)
	{
		double theta = atan2(z.y, z.x);
		theta2 = nSect * round(theta / nSect);

		double c1 = cos(theta2);
		double s1 = sin(theta2);

		double x1 = c1 * z.x + s1 * z.y;
		double y1 = -s1 * z.x + c1 * z.y;
		double z1 = z.z;

		x1 = x1 - rad1;

		z.x = fact * x1;
		z.y = fact * z1;
		z.z = fact * y1;
	}
	else
	{
		z /= t;
	}

uniform float BigRadius; slider[-10,1.,10]
float rad1=BigRadius;

uniform float SmallRadius; slider[-10,.1,10]
float rad2=SmallRadius*SmallRadius;    //only use rad2*rad2 soo.....

uniform float DonutFactor; slider[-10,3.,10]

uniform int Donuts;slider[1,9,20]
float nSect=6.283/float(Donuts);


void donut(inout vec3 z, inout float dr) {
	

	float R = sqrt(z.x * z.x + z.y * z.y);
	float R2 = rad1 - R;
	float t = R2 * R2 + z.z * z.z - rad2;  //rad2= SmallRadius^2

	if (t > 0.)   //.03 isn't as nice looking???  and it will usually be greater than zero, right??  
	{
		float theta = atan(z.y, z.x)/nSect;
		if (theta>-.5) {theta= floor (theta+.5);
		} else { theta= ceil(theta-.5);}
		theta = nSect * theta;

		float c1 = cos(theta);
		float s1 = sin(theta);

		float x1 = c1 * z.x + s1 * z.y;
		float y1 = -s1 * z.x + c1 * z.y;
		float z1 = z.z;

		x1 = x1 - rad1;

		z.x = DonutFactor * x1;
		z.y = DonutFactor * z1;
		z.z = DonutFactor * y1;
		dr=dr*(DonutFactor);
	}
	else if (t!=0.) 
	{
		z /= t;   
	}
	
}

  if (fractal->transformCommon.rotationEnabled)
  {
    z = fractal->transformCommon.rotationMatrix.RotateVector(z);
  }

  if (fractal->transformCommon.addCpixelEnabledFalse)
  {
    CVector3 tempFAB = c;
    if (fractal->transformCommon.functionEnabledx)
    {
            tempFAB.x = fabs(tempFAB.x);
    }
    if (fractal->transformCommon.functionEnabledy)
    {
            tempFAB.y = fabs(tempFAB.y);
    }
    if (fractal->transformCommon.functionEnabledz)
    {
            tempFAB.z = fabs(tempFAB.z);
    }
    tempFAB *= fractal->transformCommon.constantMultiplier000;
    if (z.x > 0) z.x += tempFAB.x;
    else z.x -= tempFAB.x;
    if (z.y > 0) z.y += tempFAB.y;
    else z.y -= tempFAB.y;
    if (z.z > 0) z.z += tempFAB.z;
    else z.z -= tempFAB.z;
  }
}

if (t > 0.03) 

if (t > 0.00) 

if (t > 0.03) 

		
double rad2 = fractal->donut.ringThickness;
double R = sqrt(z.x * z.x + z.y * z.y);
double R2 = fractal->donut.ringRadius - R;
double t = R2 * R2 + z.z * z.z - rad2 * rad2;
if(t > 0.03)
{
	z *= 2.0;
}
else
{
	z /= t;
}

if (t > 0.0)

void donut(inout vec3 z, inout float dr, inout bool donutcheck) {
	
	float c1; float s1;float t;
	float R = sqrt(z.x * z.x + z.y * z.y);

	if (test1) {  // radius check variable t=  msltoe's original radius check formula

		c1 = z.x- rad1 * z.x / R;		//rad1= Large Radius (of alleged tori)
		s1 = z.y- rad1 * z.y / R;	
		c1*=c1;
		s1*=s1;
		t=sqrt(c1+s1+z.z*z.z);

	} else {  // t = Buddhi's radius check 

		float R2 = rad1 - R;	
		t = abs(R2 * R2 + z.z * z.z - rad2);  //rad2= SmallRadius^2
	}

// I used the smaller radius to check t 
// instead of using SmallerRadius^2 * DonutFactor (.03 for original formula)

		if (t > SmallRadius)  { 			
			float theta = atan(z.y, z.x)/nSect;     //check if we are in specific section of donut

			if (theta>-.5) { theta= floor (theta+.5);}     //same as rounding theta
				   else   { theta=  ceil  (theta-.5);}    // same as rounding theta

			theta = nSect * theta;
	
			float c1 = cos(theta);
			float s1 = sin(theta);

			float x1 = c1 * z.x + s1 * z.y;
			float y1 = -s1 * z.x + c1 * z.y;
			float z1 = z.z;

			x1 = x1 - rad1;
	
			z.x = DonutFactor * x1;
			z.y = DonutFactor * z1;
			z.z = DonutFactor * y1;
			dr=dr*(DonutFactor);
		}
		else if (t!=0.) 
		{
			
			z *= t;      // yes... *=  not /=     I don't get a lot of hits for /=???
			
	   }
	
}

// different DEs!:

if (returnmode==1) {
		return 0.5*log(r)*r/dr;    // used Mandelbulb frag as a template... lazy :D
		
		
	} else if (returnmode==4) { 
		return .5*sqrt(z.x*z.x+z.y*z.y)/dr/DonutFactor;   // works nice with smaller radius = 0
	} else if (returnmode==3) { 
		return sqrt(z.x*z.x+z.z*z.z)/dr/DonutFactor*.5;   // works nice with smaller radius = 0
	} else if (returnmode==2) { 
		return sqrt(z.y*z.y+z.z*z.z)/dr/DonutFactor*.5;   // works nice with smaller radius = 0
	} else if (returnmode==5) { 
		return abs(z.x)/dr/DonutFactor*.5;    
	} else if (returnmode==6) { 
		return abs(z.y)/dr/DonutFactor*.5;
	 } else if (returnmode==7) { 
		 return abs(z.z)/dr/DonutFactor*.5;
	} else if (returnmode==8) { 
		return .5*r/dr/DonutFactor;
	}

vect3 z0 = ResultOfFractalIterations(z)
float r = length(z0)

vect3 z1 = ResultOfFractalIterations(z + vect3(delta, 0,0))
float r1 = length(z1)
float dr1 = fabs(r1 - r) / delta

vect3 z2 = ResultOfFractalIterations(z + vect3(0, delta, 0))
float r2 = length(z2)
float dr2 = fabs(r2 - r) / delta

vect3 z1 = ResultOfFractalIterations(z + vect3(0, 0, delta))
float r3 = length(z3)
float dr3 = fabs(r3 - r) / delta

double dr = sqrt(dr1*dr1 + dr2*dr2 + dr3*dr3)

distance = 0.5 * r / dr  (or distance = 0.5 * r * log(r) / dr for mandelbulbs)