Cacher is the code snippet organizer for pro developers

vector rotateFractal(vector P; float angle)
{
  vector axis = normalize(chv("rotAxis"));
  matrix rotMat = ident();
  rotate(rotMat, radians(angle), axis);
  vector pos =P*rotMat;
  return pos;
}

float sdSphere( vector p)
{
  return length(p)-1.0;
}

float maxcomp(vector p )
{
  return max(p.x,max(p.y,p.z));
}

float sdBox( vector p; vector b )
{
  vector di = set(abs(p) - b);
  float mc = maxcomp(di);
  return min(mc,length(max(di,0.0)));
}
vector swap(float a; float b; float c)
{
  return set(a,b,c);
}
vector2 rotateF(float r; float u; float v)
{
    float sine, cosine;
    sine = sin(r);
    cosine = cos(r);
    
    float x = u*cosine - u*sine;
    float y = v*cosine + u*sine;
    
    return (vector2)set(x, y);
}
float apollonian( vector pos; int ptnum)
{
        vector p = pos;
        float scale = 1.0;        
        for( int i=0; i<8;i++ )
        {
                p = -1.0 + 2.0*((0.5*p+0.5)-floor(0.5*p+0.5));

                float r2 = dot(p,p);
                float k = 1.0/r2;
                p *=k;
                scale *=k;
        }
        
        return 0.25*abs(p.y)/scale;
}

float menger(vector pos; int ptnum; float angle)
{
int iterations = point(0, "frIt", ptnum);
int modCount = point(0, "frMc", ptnum);
float modSize = point(0, "frMs", ptnum);

float scale = point(0, "msSc", ptnum);
vector offset = set(point(0, "msOF00", ptnum),point(0, "msOF01", ptnum),point(0, "msOF02", ptnum));
vector boxScale = set(point(0, "msBO00", ptnum),point(0, "msBO01", ptnum),point(0, "msBO02", ptnum));
vector rotateV = set(point(0, "msRo00", ptnum),point(0, "msRo01", ptnum),point(0, "msRo02", ptnum));



vector P = rotateFractal(pos, angle);
if(modCount > 0)
{
  vector modulo;
  float modSizeOverTwo = modSize * 0.5;
  float limit = modSizeOverTwo + modSize * modCount;
 
  //X axis---------------------------------------------------------
  if(abs(P.x) < limit)
  {
    modulo.x = ((P.x + modSizeOverTwo)%modSize) - modSizeOverTwo;
  }
  else
  {
    if(P.x < 0.0) {modulo.x = (P.x + modSize * modCount);}
    else {modulo.x = (P.x - modSize * modCount);}
  }
  //Y axis----------------------------------------------------------
  if(abs(P.y) < limit)
  {
    modulo.y = ((P.y + modSizeOverTwo)%modSize) - modSizeOverTwo;
  }
  else
  {
    if(P.y < 0.0) {modulo.y = (P.y + modSize * modCount);}
    else {modulo.y = (P.y - modSize * modCount);}
  }
  //Z axis----------------------------------------------------------
  if(abs(P.z) < limit)
  {
    modulo.z = ((P.z + modSizeOverTwo)%modSize) - modSizeOverTwo;
  }
  else
  {
    if(P.z < 0.0) {modulo.z = (P.z + modSize * modCount);}
    else {modulo.z = (P.z - modSize * modCount);}
  }
 
  P = modulo;
}

vector z = P; // position

int n = 0;
float distance = 0.0;
float trap =length2(z);
//trap =-1.0;
while(n<iterations)
{
  z = abs(z);

  if(z.x < z.y) {z = set(swap(z.y, z.x, z.z));}
  if(z.x < z.z) {z = set(swap(z.z, z.y, z.x));}
  if(z.y < z.z) {z = set(swap(z.x, z.z, z.y));}
    
  z = scale * z - offset * (scale-1);
  
  //rotation
  if(rotateV.x != 0.0) 
  {
    vector2 zRotX = rotateF(rotateV.x, z.y, z.z);
    z.y = zRotX.x;
    z.z = zRotX.y;

  }
  if(rotateV.y != 0.0) 
  {
    vector2 zRotY = rotateF(rotateV.y, z.x, z.z);
    z.x = zRotY.x;
    z.z = zRotY.y;
  }
  if(rotateV.z != 0.0)
  {
    vector2 zRotZ = rotateF(rotateV.z, z.x, z.y);
    z.x = zRotZ.x;
    z.y = zRotZ.y;
  }
  
  if(z.z < -0.5 * offset.z * (scale-1))
  {
    z.z += offset.z * (scale-1);
  }

  trap = min(trap, length2(z));
  n++;
}
 return length(z)*pow(scale, float(-n)-1.0);
}

float DE(vector P; int ptnum; float angle)
{
  //float dist = menger(P, ptnum, angle);
  float dist = apollonian(P, ptnum);
  //float dist = sdSphere(P); 
  //float dist = sdBox(P, 2.0);
  return dist;
}

// Finite difference normal
vector getNormal(vector pos; int ptnum; float angle) {
    float normalDistance = 0.0002;
    vector e = set(0.0,normalDistance,0.0);
    
    return normalize(set(
                    DE(pos+swap(e.y, e.x, e.x),ptnum,angle)-DE(pos-swap(e.y, e.x, e.x),ptnum,angle),
                    DE(pos+swap(e.x, e.y, e.x),ptnum,angle)-DE(pos-swap(e.x, e.y, e.x),ptnum,angle),
                    DE(pos+swap(e.x, e.x, e.y),ptnum,angle)-DE(pos-swap(e.x, e.x, e.y),ptnum,angle)));
}

// ray marching
void rayMarch(vector from; vector dir; int ptnum; float angle; float list[]) 
{
    float t = 0.0;
    float h = 1.0;
    float maxSteps = 256*2;
    float maxDist = 10;
    float minDist = 0.0002;
    int steps = 0;
    for( int i=0; i<maxSteps; i++ )
    {
      if( h<minDist || t>maxDist ) break;
      h = DE(from + dir*t, ptnum, angle);
      t += h;
      steps = i;
    }
    float smoothStep = float(steps)+ (h/minDist);
    float ao = 1.1-smoothStep/(float(maxSteps)*0.4);
    vector normal = getNormal(from + dir*t, ptnum, angle);
    if(t > maxDist) t = -1;
    //return set(t, ao);
    list = array(t, ao, normal.x, normal.y, normal.z);
}
//float dist = rayMarch(@P, @N, @ptnum);
float ray [];
rayMarch(@P, @N, @ptnum, chf("angleDegrees"), ray);
float dist = ray[0];
@Cd = ray[1];
if(dist==-1) removepoint(0,@ptnum);
@P+=dist*@N;  
@N = set(ray[2], ray[3], ray[4]);