% help function for sierp3

% "sierpiter.m" written by Albert Dvornik <bert@mit.edu>
% Redistribute freely, but it would be nice to let me know if you use it.


function [] = sierpiter(x,y,z,n,c)

%%% are we done recursing?
if (n <= 0)

%%% if we are, just draw the 4 sides of the tetrahedron
  side1 = [1,2,3];
  side2 = [1,2,4];
  side3 = [1,3,4];
  side4 = [2,3,4];
  patch(x(side1),y(side1),z(side1),c(1,:));  % color 1 in the cmap
  patch(x(side2),y(side2),z(side2),c(2,:));  % color 2
  patch(x(side3),y(side3),z(side3),c(3,:));  % color 3
  patch(x(side4),y(side4),z(side4),c(4,:));  % color 4

else

%%% find the midpoints of all eight edges...
  x12 = (x(1)+x(2))/2;  y12 = (y(1)+y(2))/2;  z12 = (z(1)+z(2))/2;
  x13 = (x(1)+x(3))/2;  y13 = (y(1)+y(3))/2;  z13 = (z(1)+z(3))/2;
  x14 = (x(1)+x(4))/2;  y14 = (y(1)+y(4))/2;  z14 = (z(1)+z(4))/2;
  x23 = (x(2)+x(3))/2;  y23 = (y(2)+y(3))/2;  z23 = (z(2)+z(3))/2;
  x24 = (x(2)+x(4))/2;  y24 = (y(2)+y(4))/2;  z24 = (z(2)+z(4))/2;
  x34 = (x(3)+x(4))/2;  y34 = (y(3)+y(4))/2;  z34 = (z(3)+z(4))/2;

%%% ...and recurse
  sierpiter([x(1),x12,x13,x14],[y(1),y12,y13,y14],[z(1),z12,z13,z14],n-1,c);
  sierpiter([x12,x(2),x23,x24],[y12,y(2),y23,y24],[z12,z(2),z23,z24],n-1,c);
  sierpiter([x13,x23,x(3),x34],[y13,y23,y(3),y34],[z13,z23,z(3),z34],n-1,c);
  sierpiter([x14,x24,x34,x(4)],[y14,y24,y34,y(4)],[z14,z24,z34,z(4)],n-1,c);

end
