% draw a snowball in 3D
close all;
p=exp(2*pi*i*[1:3]/3);
x=real(p)';
y=imag(p)';
z=[0,0,0]';
tri1=[x,y,z];
figure; hold on; axis equal;
patch(x,y,z,[1,.7,.7],'LineWidth',2)
center=sum(tri1)/3;
r=sqrt(norm(tri1(1,:)-tri1(2,:))^2 -norm(tri1(1,:)-center)^2);
vert1=tri1(1,:);
vert2=tri1(2,:);
vert3=tri1(3,:);
side1=vert2-vert1;
side2=vert3-vert2;
side3=vert1-vert3;
r=sqrt(norm(side1)^2-norm(center-vert1)^2);
normal=cross(side1,side2);
normal=normal/norm(normal);
vert4=center+ normal*r;
tri2=[tri1(1,:);tri1(2,:);vert4];
tri3=[tri1(2,:);tri1(3,:);vert4];
tri4=[tri1(3,:);tri1(1,:);vert4];
patch(tri2(:,1),tri2(:,2),tri2(:,3),[.7,1,.7],'LineWidth',2)
patch(tri3(:,1),tri3(:,2),tri3(:,3),[1,1,.7],'LineWidth',2)
patch(tri4(:,1),tri4(:,2),tri4(:,3),[.7,1,1],'LineWidth',2)

% given three points in 3-space compute midpoints and then the 
% tetrahedron  with smaller triangle as its base

p1=(vert1+vert2)/2;
p2=(vert2+vert3)/2;
p3=(vert3+vert1)/2;
center=(p1+p2+p3)/3;
side1=p2-p1;
side2=p3-p2;
side3=p1-p3;
r=sqrt(norm(side1)^2-norm(center-p1)^2);
normal=cross(side2,side1);
normal=normal/norm(normal);
vert5=center+normal*r;
tri2=[p1;p2;vert5];
tri3=[p2;p3;vert5];
tri4=[p3;p1;vert5];
tri1=[p1;p2;p3;p1];
plot3(tri1(:,1),tri1(:,2),tri1(:,3),'k','LineWidth',2)
%patch(tri2(:,1),tri2(:,2),tri2(:,3),[.7,1,.7],'LineWidth',2)
%patch(tri3(:,1),tri3(:,2),tri3(:,3),[1,1,.7],'LineWidth',2)
%patch(tri4(:,1),tri4(:,2),tri4(:,3),[.7,.7,1],'LineWidth',2)

%return

% do another side

p1=(vert1+vert2)/2;
p2=(vert2+vert4)/2;
p3=(vert4+vert1)/2;
center=(p1+p2+p3)/3;
side1=p2-p1;
side2=p3-p2;
side3=p1-p3;
r=sqrt(norm(side1)^2-norm(center-p1)^2);
normal=cross(side2,side1);
normal=-normal/norm(normal);
vert8=center+normal*r;
tri2=[p1;p2;vert8];
tri3=[p2;p3;vert8];
tri4=[p3;p1;vert8];
tri1=[p1;p2;p3;p1];
plot3(tri1(:,1),tri1(:,2),tri1(:,3),'k','LineWidth',2)
%patch(tri2(:,1),tri2(:,2),tri2(:,3),[.7,1,.7],'LineWidth',2)
%patch(tri3(:,1),tri3(:,2),tri3(:,3),[1,.7,.7],'LineWidth',2)
%patch(tri4(:,1),tri4(:,2),tri4(:,3),[.7,.7,1],'LineWidth',2)

p1=(vert2+vert3)/2;
p2=(vert3+vert4)/2;
p3=(vert4+vert2)/2;
center=(p1+p2+p3)/3;
side1=p2-p1;
side2=p3-p2;
side3=p1-p3;
r=sqrt(norm(side1)^2-norm(center-p1)^2);
normal=cross(side2,side1);
normal=-normal/norm(normal);
vert6=center+normal*r;
tri2=[p1;p2;vert6];
tri3=[p2;p3;vert6];
tri4=[p3;p1;vert6];
tri1=[p1;p2;p3;p1];
plot3(tri1(:,1),tri1(:,2),tri1(:,3),'k','LineWidth',2)
%patch(tri2(:,1),tri2(:,2),tri2(:,3),[.7,1,.7],'LineWidth',2)
%patch(tri3(:,1),tri3(:,2),tri3(:,3),[1,.7,.7],'LineWidth',2)
%patch(tri4(:,1),tri4(:,2),tri4(:,3),[.7,.7,1],'LineWidth',2)

p1=(vert3+vert1)/2;
p2=(vert1+vert4)/2;
p3=(vert4+vert3)/2;
center=(p1+p2+p3)/3;
side1=p2-p1;
side2=p3-p2;
side3=p1-p3;
r=sqrt(norm(side1)^2-norm(center-p1)^2);
normal=cross(side2,side1);
normal=-normal/norm(normal);
vert7=center+normal*r;
tri2=[p1;p2;vert7];
tri3=[p2;p3;vert7];
tri4=[p3;p1;vert7];
tri1=[p1;p2;p3;p1];
plot3(tri1(:,1),tri1(:,2),tri1(:,3),'k','LineWidth',2)
%patch(tri2(:,1),tri2(:,2),tri2(:,3),[.7,1,.7],'LineWidth',2)
%patch(tri3(:,1),tri3(:,2),tri3(:,3),[1,.7,.7],'LineWidth',2)
%patch(tri4(:,1),tri4(:,2),tri4(:,3),[.7,.7,1],'LineWidth',2)

tri=[vert5;vert6;vert7]
patch(tri(:,1),tri(:,2),tri(:,3),[.7,1,.7],'LineWidth',2)
tri=[vert5;vert7;vert8]
patch(tri(:,1),tri(:,2),tri(:,3),[1,1,.7],'LineWidth',2)
tri=[vert5;vert8;vert6]
patch(tri(:,1),tri(:,2),tri(:,3),[.7,1,1],'LineWidth',2)
tri=[vert6;vert7;vert8]
patch(tri(:,1),tri(:,2),tri(:,3),[1,.7,.7],'LineWidth',2)

