$fn=8;

ofaces = [
  [0, 1, 2], [0, 2, 4], [0, 4, 5], [0, 5, 1], 
  [3, 5, 4], [3, 1, 5], [3, 2, 1], [3, 4, 2],
];	

edges = [[0, 1], [0, 2], [0, 4], [0, 5], [3, 1], [3, 2], [3, 4], [3, 5], [1, 2], [2, 4], [4, 5], [5, 1]];


module octahedron(center, size) {
  points = [
    center + [l / 2, 0, 0], center + [0, l / 2, 0], center + [0, 0, l / 2],
    center - [l / 2, 0, 0], center - [0, l / 2, 0], center - [0, 0, l / 2],
  ];
  polyhedron(points=points, faces=ofaces, convexity=1);
}

// return a point on the (p1, p2) axis, next to p2, but away from p1 by alpha.
function extend(p1, p2, alpha) = p2 + alpha*(p2 - p1);

module tetrahedron(points) {
  assign(pts=[points[0],
          extend(points[0], points[1], .1),
          extend(points[0], points[2], .1),
          extend(points[0], points[3], .1)]) {
    echo(points);
    echo(pts);
    polyhedron(points=pts, faces=[[1, 2, 3],[0, 1, 3],[0, 2, 1],[0, 3, 2]], convexity=1);
  }
}

function sumv(v, i=0) = (i + 1 == len(v)? v[i]: v[i] + sumv(v, i + 1));
//function sumsqv(v, i=0) = (i + 1 == len(v)? v[i] * v[i]: v[i] * v[i] + sumsqv(v, i + 1));
function length(v) = sqrt( pow(v[0], 2) + pow(v[1], 2) + pow(v[2], 2));

magnify=4;
radius=pow(2, magnify);
l=pow(2,magnify+7);

module segments(center, points) {
        for (segment = edges) {
          assign(edge = (points[segment[0]] + points[segment[1]]) / 2 - center) {
            assign(b = acos(edge[2] / length(edge)),
                   c = (edge[0] == 0)? sign(edge[1]) * 90: (atan(edge[1] / edge[0]) + ((edge[0] > 0)? 0: 180))) {
              //echo(length(edge), b, c);
              translate(center) rotate([0, b, c]) rotate([0, 0, 22.5]) cylinder(h=length(edge)*1.1, r=radius);
            }
          }
        }
}
module volumes(center, points) {
for (face = ofaces) {
        tetrahedron([center, (points[face[0]] + points[face[1]]) / 2,
                               (points[face[1]] + points[face[2]]) / 2,
                               (points[face[2]] + points[face[0]]) / 2]);
        }
}

module carveout(points, n, type) {
  if (0 != n) {
    assign(center = sumv(points)/len(points)) {
      if (type==0)
        difference() {
          segments(center, points);
          %volumes(center, points);
        }
      else
        difference() {
          volumes(center, points);
          %segments(center, points);
        }
      for (point = points) {
        carveout([(point + points[0]) / 2, (point + points[1]) / 2, 
                  (point + points[2]) / 2, (point + points[3]) / 2,
                  (point + points[4]) / 2, (point + points[5]) / 2], n - 1, type);
      }
    }
  }
}

module foo(center, size, n, type) {
  points = [
    center + [l / 2, 0, 0], center + [0, l / 2, 0], center + [0, 0, l / 2],
    center - [l / 2, 0, 0], center - [0, l / 2, 0], center - [0, 0, l / 2],
  ];
  difference() {
    polyhedron(points=points, faces=ofaces, convexity=1);
    carveout(points, n, type);
  }
}

scale(.5/radius)
union () {
  translate([-1000, 0, 0])
    foo([0, 0, 0], l, 3, 0); // this I don't know if it is problem, but renders faster to give an idea about model
  translate([ 1000, 0, 0])
    foo([0, 0, 0], l, 3, 1); // this I don't know if it is problem, but renders faster to give an idea about model
}