Question: finding the paramterisation of the conic intersection of two spheres

Im using geom3d trying to find a way to paramterise the equation of a circle in three dimensions how ever there is no function in that ive found that can do this. Geom3d is great at drawing and giving the equations of the sphere that i require and can even plot easily using the 'intersection' command but this never gives the actual parameterisation of the circle itself.

The equations of the spheres are easy to find and using them ive been attempting to code a way to give a paramterisation of the circle from code ive found online, but this only gives me the paramterisation of the positive/negative semicircle of what I want. I can obtain both halves of the circle and plot them together but I need the paramterisation of the full circle as I intend to use this. Does anyone know a better way of achiving this? 

(also Im not sure how to extract the equation of the plane from the geom3d[intersection] comand, its there under detail but im not sure how to have it assigned to its own object.

restart;
with(geom3d);
with(plots);
geom3d[point](A_1, -8.5, 0, 0);
point(A_2, 8.5, 0, 0);
point(B_1, 0., -5.5, sqrt(166/4));
point(B_2, 0., 5.5, sqrt(166/4));
sphere(Sb1_10, [B_1, 10]);
sphere(Sb2_10, [B_2, 10]);
sphere(Sa2_10, [A_2, 10]);
sphere(Sa1_12, [A_1, 12]);
sphere(Sa2_12, [A_2, 12]);
intersection(circ, Sb1_10, Sb2_10);
detail(circ);
draw([B_1, B_2, A_1, A_2, Sb1_10, Sb2_10, circ], axes = boxed);
plane(plane1, Sb1_10, Sb2_10);
M1 := Plane(circ);
detail(M1);
detail(plane1);
M := plane(K, Plane(circ), [x, y, z]);
form(M);
Equation(M);

S1 := Equation(Sb1_10, [x, y, z]);
S2 := Equation(Sb2_10, [x, y, z]);
s := eliminate({S1, S2}, y);
S_1 := solve(s[1, 2], {x(t), z(t)});
                   

assign(s[1, 1]);
assign(S_1[1]);
Eq1 := [normal(x), normal(y), normal(z)];
space1 := spacecurve(Eq1, t = -500 .. 500, color = red, scaling = constrained, axes = boxed, numpoints = 20000);
   

unassign('x');
unassign('z');
unassign('y');
assign(s[1, 1]);
assign(S_1[2]);
Eq2 := [normal(x), normal(y), normal(z)];
space2 := spacecurve(Eq2, t = -500 .. 500, color = green, scaling = constrained, axes = boxed, numpoints = 20000, title = "2");
 

unassign('x');
unassign('z');
unassign('y');
assign(s[2, 1]);
assign(S_1[1]);
Eq3 := [normal(x), normal(y), normal(z)];
space3 := spacecurve(Eq3, t = -500 .. 500, color = pink, scaling = constrained, axes = boxed, numpoints = 20000, title = "3");
 

unassign('x');
unassign('z');
unassign('y');
assign(s[2, 1]);
assign(S_1[2]);
Eq4 := [normal(x), normal(y), normal(z)];
space4 := spacecurve(Eq4, t = -500 .. 500, color = orange, scaling = constrained, axes = boxed, numpoints = 20000, title = "4");
 

Eq := Eq1 + Eq2 + Eq3 + Eq4;
spaceC := spacecurve({Eq}, t = -500 .. 500, color = blue, scaling = constrained, axes = boxed, numpoints = 20000);
display({space1, space2, space3, space4, spaceC});
 

restart;
with(geometry);
circle(c1, [point(A, 0, 0), point(B, 2, 0), point(C, 1, 2)], 'centername' = O1);
f1 := (x - 8.5)^2 + (y + 0)^2 + (z + 0)^2 = 12^2;
f2 := (x + 8.5)^2 + (y + 0)^2 + (z + 0)^2 = 12^2;
                                      2    2    2      
                f1 := (x - 8.5)  + y  + z  = 144

                                       2    2    2      
                f2 := (x + 8.5)  + y  + z  = 144

restart;
with(plottools);
with(plots);
f1 := x^2 + y^2 + z^2 = 1;
f2 := x + y + z = 1;
Circle := intersectplot(f1, f2, x = -10 .. 10, y = -2 .. 2, z = -5 .. 15, color = red, thickness = 3, numpoints = 10000);
S1 := implicitplot3d(f1, x = -20 .. 20, y = -20 .. 20, z = -10 .. 20, style = patchnogrid, color = blue);
S2 := implicitplot3d(f2, x = -20 .. 20, y = -20 .. 10, z = -10 .. 20, style = patchnogrid, color = gold, transparency = 0.5);
display(S1, S2, Circle, scaling = constrained, axes = boxed);
whattype(Circle);
detail(Circle);
   

sol := eliminate({f1, f2}, z);
Sol := solve(sol[2, 1], {x(t), y(t)});
assign(sol[1, 1]);
assign(Sol);
Eq := [normal(x), normal(y), normal(z)];
spacecurve(Eq, t = -500 .. 500, color = red, scaling = constrained, axes = boxed, numpoints = 20000);
 

restart;
with(geom3d);
with(plots);
f1 := -28.25000000 + x^2 + y^2 + z^2 + 11.0*y - z*sqrt(166) = 0;
plane1 := y = 0;
Circle := intersectplot(f1, plane1, x = -10 .. 10, y = -2 .. 2, z = -5 .. 15, color = red, thickness = 3, numpoints = 10000);
sol := eliminate({f1, plane1}, y);
Sol := solve(sol[2], {x(t), z(t)});
 

assign(sol[1]);
assign(Sol[2]);
Eq1 := [normal(x), normal(y), normal(z)];
C2 := spacecurve(Eq1, t = -500 .. 500, color = red, scaling = constrained, axes = boxed, numpoints = 20000);
unassign('x');
unassign('z');
unassign('y');
assign(sol[1]);
assign(Sol[1]);
Eq2 := [normal(x), normal(y), normal(z)];
C1 := spacecurve(Eq2, t = -500 .. 500, color = green, scaling = constrained, axes = boxed, numpoints = 20000);
 

Eq3 := Eq1 + Eq2;
C3 := spacecurve(Eq3, t = -500 .. 500, color = green, scaling = constrained, axes = boxed, numpoints = 20000);
display(C1, C2, C3);