kill(all);

van=facebook.com/LearningCalculus$

      The normal to the plane is the vector (1,1,1)
      The point  (2,2,2) is on the plane in the center of the circle
      The radius of the circle is the distance between (2,3,1) and (2,2,2).

P(x,y,z):=[[x],[y],[z],[1]];

P(x,y,z)[1];

p:P(1,0,0);

p[1];

dist(P1, P2):= sqrt((P1[1]-P2[1])^2 +(P1[2]-P2[2])^2 +(P1[3]-P2[3])^2);

dist(P(2,3,1),P(2,2,2));

      So the circle of the intersection of x^2+y^2+z^2 = 14 and x + y + z = 6 has radius _/2
      So now all we have to do is rotate the sphere/plane pair whose normal is (1,1,1,0) to be
      coincident with the z axis.
      The plane is a distance 3.46... away from the origin.

dist(P(2,2,2),P(0,0,0));

float(%);

      Then we scale the coordinates to make the radius of the circle unity:
      The scale factor is 1/_/2 in x, y, and z
      then we check that the new radius of the circle to be projecteds is 1

      Then we translate  in z to bring the circle to the x, y plane and we have the correct transform
      We then transform equation 3 into this new coordinate system and find the point of intersection.
      We then invert the transform, apply it to this result and get the point 2, 3, 1.

      Tools for Rotation

RotateX(a):=matrix(
                    [1,      0,      0, 0],
                    [0, cos(a),-sin(a), 0],
                    [0, sin(a), cos(a), 0],
                    [0,      0,      0, 1]
                  );

RotateY(b):=matrix(
                    [cos(b), 0,-sin(b), 0],
                    [     0, 1,      0, 0],
                    [sin(b), 0, cos(b), 0],
                    [     0, 0,      0, 1]
               );

RotateZ(c):=matrix(
                 [cos(c),-sin(c), 0, 0],
                 [sin(c), cos(c), 0, 0],
                 [     0,      0, 1, 0],
                 [     0,      0, 0, 1]
               );

P(x,y,z):=[[x],[y],[z],[1]];

RotateZ(%pi).P(1,0,0);

P(2,3,1);

foo:RotateZ(atan(2/3)).P(2,3,1);

bar:RotateX(atan(sqrt(13))).foo;

ScaleXYZ(sx,sy,sz):=matrix(
                        [sx,  0,  0, 0],
                        [ 0, sy,  0, 0],
                        [ 0,  0, sz, 0],
                        [ 0,  0,  0, 1]
                       );

ScaleXYZ(1, 1, 1/sqrt(14)).bar;

ScaleXYZ(1, 1, 1/sqrt(14)).RotateX(atan(sqrt(13))).RotateZ(atan(2/3)).P(2,3,1);

ScaleXYZ(1, 1, 1/sqrt(14)).RotateX(atan(sqrt(13))).RotateZ(atan(2/3)).P(x,y,z);

float(sqrt(14));

3*x = 6;

ScaleXYZ(1, 1, 1/sqrt(14)).RotateX(atan(sqrt(13))).RotateZ(atan(2/3)).P(1,1,1);

float(%);

solve(%,(a,b));

f[i,j]:= if i = j then x^i else 0$
g[i,j]:= if i = j then 1 else 0$

A:genmatrix(f,4,4);

determinant(A);

eigenvalues(A);

eigenvectors(A);

gramschmidt(A);

invert(A).genmatrix(f,4,4);

A^^A;

A ^^ -1;

transpose(A);

load("nchrpl")$
ncharpoly(A,x);

solve(%,x);

load("eigen")$
similaritytransform(A);

A^2 + 2*A;