###################################################################
#
# This is a solution to problem 5 of Lab 8.
# The program draws a tetrahedron and its
# circumscribed sphere.
#
# Input: 4 distinct points in space
#
# The auxiliary function "inter" (for computing the
# intersection of two 3D lines) had to be slightly
# modified.
#
###################################################################

inter := proc (L1, L2)

  local eq1, eq2, eq3, sol;
        eq1 := L1[1] - subs(t=s, L2[1]);
        eq2 := L1[2] - subs(t=s, L2[2]);
        eq3 := L1[3] - subs(t=s, L2[3]);
        sol := solve({eq1, eq2, eq3}, {s, t});
        RETURN ([subs(sol,L1[1]),subs(sol,L1[2]),subs(sol,L1[3])]);
end;

###################################################################

tetrasphere := proc (x1,y1,z1,  # coordinates of 1st point
                     x2,y2,z2,  # coordinates of 2nd point
                     x3,y3,z3,  # coordinates of 3rd point
                     x4,y4,z4)  # coordinates of 4th point

  local cs, r,          # center and radius of
                        # circumscribed sphere
        n, s1, s2,      # auxiliary vectors
        c1, c2, b1, b2, #
        tetra, sphere,  # the two plots
        points,         # plot of corner points
        m,              # intersection of bisectors
        d1, d2          # lines intersecting in cs
        ;
  s1 := [x2-x1,y2-y1,z2-z1];
  s2 := [x3-x1,y3-y1,z3-z1];
  c1 := [(x2+x1)/2,(y2+y1)/2,(z2+z1)/2];
  c2 := [(x3+x1)/2,(y3+y1)/2,(z3+z1)/2];
  n  := linalg[crossprod](s1,s2);
  b1 := evalm (c1 + t*linalg[crossprod](s1,n));
  b2 := evalm (c2 + t*linalg[crossprod](s2,n));
  m  := inter(b1,b2);
  d1 := evalm (m + t*n);

  s2 := [x4-x1,y4-y1,z4-z1];
  c2 := [(x4+x1)/2,(y4+y1)/2,(z4+z1)/2];
  n  := linalg[crossprod](s1,s2);
  b1 := evalm (c1 + t*linalg[crossprod](s1,n));
  b2 := evalm (c2 + t*linalg[crossprod](s2,n));
  m  := inter(b1,b2);
  d2 := evalm (m + t*n);

  cs := inter(d1,d2);
  r  := evalf(sqrt((cs[1]-x1)^2 + (cs[2]-y1)^2 +(cs[3]-z1)^2));

  points := plots[pointplot]({[x1,y1,z1],[x2,y2,z2],[x3,y3,z3],[x4,y4,z4]},
            color=red);
  tetra  := plots[polygonplot3d]([[x1,y1,z1],[x2,y2,z2],[x3,y3,z3],
            [x4,y4,z4],[x1,y1,z1],[x3,y3,z3],[x4,y4,z4],[x2,y2,z2]],
            color=green);
  sphere := plot3d([cs[1]+r*sin(x)*cos(y),cs[2]+r*cos(x)*cos(y),
            cs[3]+r*sin(y)],x=0..2*Pi,y=0..2*Pi);
  plots[display]({points,tetra,sphere},style=WIREFRAME);

end;

###################################################################