Elementary Möbius Geometry II

Circles

by 
Rolf Sulanke

Starting date  February 9, 2011
Finished  March 11, 2011
Revised December 17,2011

Introduction

This notebook  continues  the notebook [emg1] about elementary Möbius geometry.  Read first the introduction of this notebook.  We consider now  the circles under the action of the Möbius group.   Again the three-dimensional  case is emphasized.   The calculations and examples illustrate and complete the presentation of the subject in our book [OS3], Section 2.7.  There one finds the foundations of elementary Möbius  geometry in terms of linear algebra.

The Möbius group acts transitively on the Möbius space: the sphere (S^n) as a point manifold, and on the manifold of the k-dimensional subspheres.   A circle is a 1-dimensional subsphere; thus it is modelled as  the set of isotropic 1-dimensional subspaces of a three-dimensional Minkowski space. The manifold of all circles in the n-dimensional Möbius space  is equivariant isomorphic to the manifold of all 3-dimensional pseudo-Euclidean subspaces of the (n+2)-dimensional Minkowski space.  Equivariantly one may  pass to the orthogonal complements being (n-1)-dimensional Euclidean subspaces of the (n+2)-dimensional Minkowski  space. Thus. in case n=3, the circle space appears as a 6-dimensional manifold, equivariantly modelled as an open orbit in the Grassmann manifold (Subscript[G, 5,3]) of the 3-dimensional subspaces of the 5-dimensional pseudo-Euclidean vector space of index 1.

In Chapter 2 the circles are considered as geometric objects in the 3-dimensional Euclidean space.  The Chapters 3 -6  may be considered as an extended update of the corresponding chapters in the notebook  mcircles.nb of the collection [Spheres].  The contents of each chapter is described in its first paragraph.

Keywords

Möbius group,  isotropic cone, isotropic-orthogonal coordinates, pseudo-orthonormal coordinates, circles,  circle pairs, interlaced circles, circle space, pseudo-Euclidean space,  pseudo-orthogonal Lie algebra,  generalized  Coxeter invariant, circle surfaces, tubes, geodesics of the circle manifold.

Copyright

This notebook and the accompanying packages are  public. Authors who intend to publish a changed or completed version of them should do it under their own names with the condition that they cite the original notebook with the internet address or other source where they got it. I am not responsible for errors or damages originated by the use of the procedures contained in my notebooks or packages; everybody who applies them should test carefully whether they are appropriate for his purposes.

Initialization

Before starting to work interactively with this notebook first time, read this section carefully and  make the necessary preparation. Later it suffices to call the menu item "Evaluation. Evaluate Initialization Cells".

The needed packages

For working with the notebook you need the packages
euvec.m, vectorcalc.m,  neuvec.m,  mspher.m, moeb.m, mcirc.m.
You may download these and other packages from my Homepage .  From there you find all the packages mentioned above together with the code of the notebook in the file emg.tgz.
Before initializing the notebook ensure that these packages and the current notebook are laying in your working directory, which should be on the $Path of Mathematica on your Computerm, or

Insert the path of your working directory into the  appropriate of the next two cells,  e.g.

For Windows:

SetDirectory["C:\\Users\\Rolf Sulanke\\mymath\\emg"];

For Linux:

SetDirectory["~/mymath/wdir"]

Now give the cell corresponding to your  operating system the properties "Cell Evaluatable" and "Initialization cell" (Menu Cell/Cell properties), and inactivate these properties for the  cells corresponding to the other operating system. If this is done, save the notebook. Next time you may start the notebook directly with the evaluation of the initialization, as follows:

The Initialization

Before starting to work interactively with the notebook

Activite Evaluation/Evaluate Initialization from the menu.

The first command loads the packages mentioned above:

Needs["mcirc`"]

In this notebook we  shall  consider the 3-dimensional  Möbius space as the elliptical quadric,  the 3-dimenssional sphere, in the 4-dimensional real projective space.  Therefore we set the dimension dim and the index ind  of the vector space as follows:

dim::usage="dim is the dimension of the vector space under consideration; it must be set within the Global Context.";

?dim

dim =5;

ind ::usage="ind is the index of the pseudo-Euclidean vector space under consideration; it must be set within the Global Context.";

ind =1;

Examples:

dim

?vec

vec[x]

?dual

dual[vec[x]]

ind

pssp[vec[x],vec[y]]

$Assumptions

We always assume that the global constants and variables appearing in the next command are real numbers:

$Assumptions=s∈Reals&&t∈Reals&&u∈Reals&&v∈Reals&&w∈Reals&&x∈Reals&&y∈Reals&&z∈Reals&&a∈Reals&&b∈Reals&&c∈Reals&&k∈Reals&&h∈Reals;

$Assumptions

1. List of Symbols and their Usages

In this Chapter one finds tables of the symbols introduced in the imported packages and in the Global Context.
To get the usages click on the name! If this does not work correctly, enable Dynamic Updating in the Evaluation Menu.

1.1.  Symbols in the Package vectorcalc.m

The package vectorcalc.m contains some functions specializing or supplementing the Mathematica built-in functions in the field of linear algebra, affine or projective geometry.

?vectorcalc`*

1.2. Symbols in the Package euvec.m

The package euvec.m contains functions for applications to the Euclidean analytic geometry.

?euvec`*

1.3. Symbols in the Package neuvec.m

The package neuvec.m contains enhancements of Mathematica  for applications to pseudo-Euclidean geometry.

?neuvec`*

Options[orthonorm]

ide[]

ide[]*IdentityMatrix[dim]

%//MatrixForm

Thiscoincides with the Gram matrix of a pseudo-orthonormal Basis:

Table[pssp[stb[i],stb[j]],{i,1,dim},{j,1,dim}]//MatrixForm

dual[vec[x,7]]

pssp[vec[x,7],vec[y,7],3]

Options[orthonorm]

??neglect

Context[neglect]

1.4. Symbols in the Package mspher.m

The package mspher.m contains functions relating Möbius and Euclidean geometry of spheres in the 3-dimensional Möbius space.

?mspher`*

1.5. Symbols in the Package moeb.m

The functions of the package moeb.m relate spseudo-orthonormal and isotropic-orthonormal coordinates of the pseudo-Euclidean vector space. Furthermore funcions describing Euclidean and hyperbolic geometry as subgeometries of the Möbius geometry.

?moeb`*

1.5.  Symbols in the Package mcirc.m

The package mcirc.m containes functions specific for the geometry of circles in the 3-dimensional Euclidean space and the Möbius space.

?mcirc`*

??control

1.6.  Symbols in the Global Context

Evaluate the next cell to get the actual list of the Global`*  symbols in the Global Context.

?Global`*

2.  Circles in the Euclidean 3-space

Mathematica has built-in  circles in dimension  two only:

?Circle

Thus,  first  we will  introduce circles as objects of the  three-dimensional Euclidean space. In 2.1 we consider two parmeter representations of  arbitrary circles; since a circle in the 3-space depends on 6 parameters the circle space is a 6-dimensional manifold. In 2.2 the corresponding plot definitions are presented.  In subsection 2.3 the circle (or line, or point) defined by a point tripel of the 3-space is calculated and plotted. As an application we construct tubes of space curves, especially torus knots, in the last subsection.

2.1 Parametrisation of the Circle Space. The General Circle

2.1.1. The Unit Sphere with Center o and its Tangential Unit Vectors

unit[2][{u,v}] is the usual parameter representation of the unit sphere in the 3-space with the geographical coordinates u, v: 
   0 <= u < 2Pi,  -Pi/2 < v < Pi/2.

unit[2][{u,v}]

Clear[sph1];sph1 =ParametricPlot3D[unit[2][{u,v}],{u, 0,2 Pi}, {v,-Pi/2,Pi/2}, Mesh -> False, PlotStyle -> Opacity[.3]]

For all values of u,v,t the following vector has norm 1:

?unitvec

Context[unitvec]

unitvec[u,v][t]

Simplify[unitvec[u,v][t].unitvec[u,v][t]]

unitvec[u,v][t] is orthogonal to unit[2][{u,v}]:

Simplify[unitvec[u,v][t].unit[2][{u,v}]]

2.1.2. Examples

For u,v fix, t variable, unitvec[u,v][t] is the parameter representation of the unit circle in the tangent vector space of the unit sphere at the point with parameters u,v.

tangcircle[u_,v_][t_] :=unit[2][{u,v}] + unitvec[u,v][t]

ParametricPlot3D[unitvec[0,Pi/3][t],{t,0,2 Pi}]

ParametricPlot3D[tangcircle[0,Pi/3][t],{t,0,2 Pi}]

Show[{%%,%},PlotRange -> All]

ParametricPlot3D[unitvec[0,Pi/2][t],{t,0,2 Pi}]

ParametricPlot3D[unitvec[Pi/2,Pi/4][t],{t,0,2 Pi},PlotStyle ->  {Thick,Red}]

ParametricPlot3D[tangcircle[Pi/2,Pi/4][t],{t,0,2 Pi},PlotStyle -> {Thick,Green}]

Show[{%%,%},PlotRange -> All]

cp1 = %;

pt = Graphics3D[{PointSize[.015],Point[unit[2][{Pi/2,Pi/4}]]}];

Show[{sph1,cp1,pt},PlotRange -> All]

2.1.3. Definition of circle3D

?circle3D

circle3D[x_,y_,z_,u_,v_,r_][t_]:={x, y,z} + unitvec[u,v][t]*r

circle3D[x,y,z,u,v,r][t]

A similar, but more general  definition of Circle3D as a special case of an ellipse in the 3-dimensional space is contained in the package Hyperbolic/Circle3D.m by Sylvio Levy  (mathsource 1992). Furthermore, it includes a plot definition. This package seems to run only under  Mathematica v. 2.

2.1.4. Examples

ParametricPlot3D[circle3D[2,0,0,Pi/2, Pi/4, 2][t],{t,0,2 Pi}]

u1=2 π RandomReal[];v1=((RandomReal[] π)/2);ParametricPlot3D[circle3D[0,0,0,u1,v1,1][t],{t,0,2 π}]

Table[u1=2 π RandomReal[];v1=((RandomReal[] π)/2);ParametricPlot3D[circle3D[0,0,0,u1,v1,1][t],{t,0,2 π},PlotStyle -> Red],{i,10}];

Show[%,PlotRange -> All]

Show[{sph1,%}]

Table[u1=2 π RandomReal[];v1=((RandomReal[] π)/2);ParametricPlot3D[circle3D[i/5,0,0,u1,v1,1][t],{t,0,2 π},PlotStyle -> Red],{i,10}];

Show[%,PlotRange -> All]

Table[u1=2 π RandomReal[];v1=((RandomReal[] π)/2);ParametricPlot3D[circle3D[i/5,0,0,u1,v1,i/5][t],{t,0,2 π},PlotStyle -> Red],{i,10}];

Show[%,PlotRange -> All]

Clear[u1,v1]

2.1.5. Another Definition: gencircle

In the next definition the radius of the circle is the norm of v:

?gencircle

gencircle[p_List,v_List][t_]:=Module[{a = v.v,b1, b2},
			If[Length[p]!=3||Length[v]!=3, Print["Incorrect arguments!!!"];Return[]];
	  	If[a== 0, Return[p]];
			If[v[[1]] == 0 && v[[2]] ==0, Return[p + {Cos[t], Sin[t], 0}*Sqrt[a]]];
			b1={-v[[2]], v[[1]], 0}/Sqrt[a - v[[3]]^2];
			b2 = Cross[v, b1];
			Return[p+b1*Sqrt[a]*Cos[t] +b2 Sin[t]]]

Context[gencircle]

2.1.6. Examples

ParametricPlot3D[gencircle[{2,0,0},2*unit[2][{Pi/2,Pi/4}]][t],{t,0,2 Pi}]

randomv[3]

(p1 = randomv[3]; v1 =-randomv[3];
	Print[p1,", ", v1,", ", Sqrt[v1.v1]]; 
ParametricPlot3D[gencircle[p1,v1][t],{t,0,2 Pi}])

Clear[p1,v1]

2.2 The Plot Definitions

?circle3D

plotcircle3D[x_,y_,z_,u_,v_,r_, opts___]:= ParametricPlot3D[
		Evaluate[circle3D[x,y,z,u,v,r][t]], {t, 0, 2 Pi}, opts]

?gencircle

plotgencircle[p_List,v_List, opts___]:= ParametricPlot3D[
		Evaluate[gencircle[p,v][t]], {t, 0, 2 Pi}, opts]

Examples

plotgencircle[null[3], stb[3,3]]

plotgencircle[null[3], 25.*stb[2,3]]

plotgencircle[ 2randomv[3],5 randomv[3]]

plotcircle3D[1,1,1,Pi/2,0, 2]

plotcircle3D[1,5 RandomReal[],1,(π/2),π RandomReal[],2]

2.3 Circles and Lines through Three Points

Definitions: Position Vector, Center, Radius, and Circle

The position vector of a circle is a unit vector orthogonal to the plane of the circle; it is uniquely defined up to sign, and it can be calculated from three points of the circle as follows:

posvec3pts[p1_List,p2_List,p3_List] :=
	Module[{ a},
		a = Cross[p2-p1, p3-p1];
		If[a ≠null[3],Return[Chop[a/Sqrt[a.a]]], Return[null[3]],Return[Chop[a/Sqrt[a.a]]]]]

Context[posvec3pts]

For the calculation of the center m of the circle through three points we start with  the  formula

m = p3 + (p1 - p3) x + (p2 - p3) y

Evaluating the conditions

(m-p3).(m-p3)== (m-p1).(m-p1);
(m-p3).(m-p3)== (m-p2).(m-p2);

we get a linear system for the unknowns x, y. By the Cauchy-Schwarz-inequality its determinant dn is zero if and only if the points p1, p2, p3 belong to a line. By definition,  then the center is {∞,∞,∞}. Otherwise, Cramer´s rule yields the solution.  These considerations lead to the following module:

center3pts[p1_List,p2_List,p3_List]:= 
  Module[{nn=Length[p1], q1= p1 - p3, q2 =p2 - p3,dn,dn1,dn2},
    dn = Det[{{q1.q1, q1.q2},{q1.q2,q2.q2}}];
    If[Chop[dn] == 0, Return[Table[Infinity,{nn}]]];
    dn1 = Det[{{q1.q1/2, q1.q2},{q2.q2/2,q2.q2}}];
    dn2 = Det[{{q1.q1, q1.q1/2},{q1.q2,q2.q2/2}}];
    Return[p3+(dn1*q1 + dn2*q2)/dn]]

Knowing the center it is easy to calculate the radius:

radius3pts[p1_List,p2_List,p3_List] :=
	Module[{a},
		a = center3pts[p1,p2,p3];
		If[a.a == Infinity, Return[Infinity]];
		Return[ Sqrt[(p1-a).(p1 - a)]]]

WARNING: Always, if there exist a line through the three points, e.g. if the points coincide, center3pts and radius3pts are infinite!

center3pts[stb[1,3],stb[3,3],stb[1,3]]

center3pts[stb[1,2],stb[2,2],-stb[1,2]]

radius3pts[stb[1,2],stb[2,2],-stb[1,2]]

radius3pts[stb[1,3],-stb[1,3],2*stb[1,3]]

center3pts and radius3pts work in arbitrary dimension, e,g. for dim = 5:

radius3pts[stb[1],stb[3],stb[2]]

center3pts[stb[1],stb[3],stb[2]]

Now we  give a definition of a 3D-circle through three points. The circle may eventually be a line (circle of infinite radius):

?circle3pts

circle3pts[p1_List,p2_List,p3_List][t_] :=
	Module[{q2 , q3 , a, m,r,nn =Length[p1]},
If[p1.p1 == Infinity,Return[p2 + t*(p3-p2)]];
If[p2.p2 == Infinity,Return[p1 + t*(p3-p1)]];
If[p3.p3 == Infinity,Return[p1 + t*(p2-p1)]];
m = center3pts[p1,p2,p3];
If[m.m ==Infinity,If[Chop[p2 -p1]==null[nn], Return[p1 + t*(p3 - p1)], Return[p1 + t*(p2 - p1)]]];
r = Sqrt[(p1-m).(p1 - m)];	
	q2 = p2-p1; q3 = p3-p1;
a = Orthogonalize[{q2,q3}];
		q2  = Normalize[a[[1]]];
		q3 =  Normalize[a[[2]]];
		m +(q2 *Cos[t] + q3*Sin[t])*r];

Here is the corresponding plot command:

plotcircle3pts[p1_List,p2_List,p3_List,opt___]:= ParametricPlot3D[Evaluate[circle3pts[p1,p2,p3][t]],{t, -Pi, Pi},opt]

Examples

circle3pts[{1,0,0},{0,1,0},{0,-1,0}][t]

plotcircle3pts[{1,0,0},{0,1,0},{0,-1,0}]

circle3pts[{1,0,0},{1,0,0},{0,-1,0}][t]

plotcircle3pts[{1,0,0},{1,0,0},{0,-1,0}]

vv=randomv[];ww=randomv[]; uu=randomv[];xx=circle3pts[uu,vv,ww][t]

mm = center3pts[uu,vv,ww]; rr =radius3pts[uu,vv,ww];

Chop[Simplify[rr^2-(xx-mm).(xx-mm)]]

Chop[N[%]]

smoothing[N[%%]]

FullSimplify[%]

Plot[%+1,{t,0,2*Pi}, PlotStyle->{Red,Thick}]

plotcircle3pts[{1,0,0},{0,1,0},{0,-1,1}]

circle3pts[{3,0,0},{3,1,0},{3,-1,0}][t]

plotcircle3pts[{3,0,0},{3,1,0},{3,-1,0}]

a=RandomReal[];b=RandomReal[];c=RandomReal[];
Print[{a,b,c}];
pa={a,0,0};pb={0,b,0};pc={0,0,c};
gr1=plotcircle3pts[pa,pb,pc];
gr2=Graphics3D[{PointSize[0.03`],Point[pa],Point[pb],Point[pc]}]
Show[{gr1,gr2}]

center3pts[pa,pb,pc]

radius3pts[pa,pb,pc]

posvec3pts[pa,pb,pc]

%.%

pa=randomv[3];pb=-randomv[3];pc=2*randomv[3];
Print[{pa,pb,pc}];
gr1=plotcircle3pts[pa,pb,pc];
gr2=Graphics3D[{PointSize[0.03`],Point[pa],Point[pb],Point[pc]}]
Show[{gr1,gr2}]

a=RandomReal[];b=RandomReal[];
Print[{a,b}];
pa={a,2 b,0};pb={0,-b,a};pc=((pa+pb)/2);
gr1=plotcircle3pts[pa,pb,pc];
gr2=Graphics3D[{PointSize[0.03`],Point[pa],Point[pb],Point[pc]}]
Show[{gr1,gr2}]

{pa,pb,pc}

If the three points coincide, circle3pts returns the point:

circle3pts[{1,0,0},{1,0,0},{1,0,0}][t]

plotcircle3pts[{1,0,0},{1,0,0},{1,0,0},PlotStyle-> {Red,Thick}]

If only two points coincide we get the line through these points:

circle3pts[{1,0,0},{1,0,0},{2,0,0}][t]

plotcircle3pts[{1,0,0},{1,0,0},{2,0,0}]

plotcircle3pts[{1,3,0},{1,3,0},{2,0,1}]

circle3pts[{1,3,0},{1,3,0},{2,0,1}][t]

circle3pts[{1,0,0},{2,0,0},{1,0,0}][t]

circle3pts[{1,0,0},{2,0,0},{2,0,0}][t]

posvec3pts[{1,0,0},{0,1,0},{0,-1,0}]

center3pts[{1,0,0},{0,1,0},{0,-1,0}]

radius3pts[{1,0,0},{0,1,0},{0,-1,0}]

For three collinear points the position vector is the null vector:

posvec3pts[{3,0,0},{3,1,0},{3,-1,0}]

center3pts[{3,0,0},{3,1,0},{3,-1,0}]

radius3pts[{3,0,0},{3,1,0},{3,-1,0}]

a=4 RandomReal[];b=4 RandomReal[];c=4 RandomReal[];
Print[{a,b,c}];
pa={a,0,0};pb={0,b,0};pc={0,0,c};
pv=posvec3pts[pa,pb,pc];
pm=center3pts[pa,pb,pc];
Print["Position vector = ",pv,", ","center = ",pm,", ","radius = ",radius3pts[pa,pb,pc]];
gr1=plotcircle3pts[pa,pb,pc];
gr2=Graphics3D[{PointSize[0.03],Point[pa],Point[pb],Point[pc]}];
gr3=Graphics3D[Arrow[{pm,pm+pv}]];
Show[{gr1,gr2,gr3},PlotRange->All]

a=RandomReal[];b=RandomReal[];
Print[{a,b}];
pa={a,b,0};pb={RandomReal[],b,RandomReal[] a};pc=((pa+pb)/2);
Print[center3pts[pa,pb,pc],", ",radius3pts[pa,pb,pc]];
gr1=plotcircle3pts[pa,pb,pc,PlotStyle -> {Red,Thick}];
gr2=Graphics3D[{PointSize[0.03],Point[pa],Point[pb],Point[pc]}];
Show[{gr1,gr2}]

pa=2 randomv[3];pb=-3 randomv[3];pc=((pa+pb)/2);
Print[center3pts[pa,pb,pc],", ",radius3pts[pa,pb,pc]];
gr1=plotcircle3pts[pa,pb,pc];
gr2=Graphics3D[{PointSize[0.03`],Point[pa],Point[pb],Point[pc]}];
Show[{gr1,gr2}]

Clear[a,b,c,pa,pb,pc,pv,pm,gr1,gr2,gr3]

2.4 An Application: Tubes of a Space Curve

We plot a tube around the helix.

?circle3D

c1 =ArcSin[0.5/Sqrt[1.25]]

ParametricPlot3D[{Cos[s],Sin[s],0.5` s},{s,-6,6}]

ParametricPlot3D[circle3D[Cos[s], Sin[s], 0.5 s,s+Pi/2,c1,0.1][t],{s,-6,6},{t, 0, 2 Pi},PlotPoints -> 60]

heltube= %;

The following example will be generalized in the next definition:

ParametricPlot3D[gencircle[{Cos[s], Sin[s], 0.5 s},{-Sin[s],Cos[s],0.5}*0.1][t],{s,-6,6},
	{t, 0, 2 Pi},PlotPoints -> 40]

heltube2=%;

The following definition is more general than that given in Alfred Gray's book [G1], because we do not need the normal vectors; regularity is not assumed, but applied implicitely. The parameter representation of the curve must be one times differentiable.

?tube

tube[x_, a_,b_,r_, opts___] := Module[{t,s, w,v, u},
					v[s_] := D[x[u],u]/.u -> s;
			    w[s_]:=If[v[s].v[s]=== 0.,Return[v[s]],
				  Return[ v[s]/Sqrt[v[s].v[s]]]];
				ParametricPlot3D[gencircle[x[s],r*w[s]][t], {s,a,b},{t,0,2 Pi},opts]]

Context[tube]

helix1[s_] := {Cos[s], Sin[s], 0.5*s}

Timing[tube[helix1,0, 6 Pi,0.3, PlotPoints -> 50]]

One finds Mathematica definitions of numerous examples of space curves  in Appendix A.3 of Alfred Gray's book [G1].  We construct tubes of some of these space curves.

Tubes of ast3d

Here is a curve ast3d with singularities, a space curve analog of the astroid:

curv[s_] := {((Cos[s])^3),1.5` ((Sin[s])^3),Cos[2 s]}

ParametricPlot3D[curv[s], {s,0,10}]

tube[curv,0, 2 Pi,0.08, PlotRange -> All]

Show[%, Boxed->False, Axes -> False]

ast3 = %;

Tubes of a bicilinder curve

The next curve bicylinder is the intersection of two cylinders:

curv[s_] := {Cos[s],Sin[s],(Sqrt[4-((Sin[s])^2)])}

ParametricPlot3D[curv[s], {s,0,10}]

tube[curv,0, 2 Pi,0.08, PlotRange -> All]

Tubes of a loxodrome

Now we shall consider a loxodrome, obtained by conformal embedding or inverse stereographical projection of a logar:thmic spirale, see the notebook [lox].

lox[a_][t_] :={((2 (Sqrt[2]) (E^(a t)) Cos[t])/(2+(E^(2 a t)))),((2 (Sqrt[2]) (E^(a t)) Sin[t])/(2+(E^(2 a t)))),-((-2+(E^(2 a t)))/(2+(E^(2 a t))))}

ParametricPlot3D[lox[.3][s], {s,-10,10},PlotRange -> All, PlotStyle -> {Red, Thick}]

loxodr1 = %;

Show[{sph1,loxodr1}, Boxed-> False, Axes -> None]

tube[lox[.3],-10, 10,0.03,PlotPoints ->50, PlotRange -> All]

Show[{sph1,loxodr1, %}, Boxed-> False, Axes -> None]

Example. Tubes of  Torus  Knots

?tube

The general circular torus is the surface of revolution

Clear[torussf];torussf[a_,b_][u_,v_]:= {(a+b*Cos[u])*Cos[v],(a+b*Cos[u])*Sin[v],b*Sin[u]}

ParametricPlot3D[torussf[1, 1/3][u,v], {u,0,2*Pi}, {v,0,2*Pi}, Mesh-> None, PlotStyle-> Opacity[.3]]

torus1=%;

A torus knot is a curve on the torus with the parameter representation

torkn[a_,b_][p_,q_][t_] :={(a+b*Cos[p*t])*Cos[q*t],(a+b*Cos[p*t])*Sin[q*t],b*Sin[p*t]}

A tube of radius r  of such a curve with tube radius p and circleparameter s can be defined as

torusknottube[a_,b_][p_,q_,r_,opts___]:= tube[torkn[a,b][p,q],0,2*Pi,r,opts]

Here a,b are the parameters defining the torus, p,q are integers defining the knot, and r is the tube radius

torusknottube[1,1/3][5,3,0.03,PlotPoints ->{80,20}]

Show[{%, torus1}, Axes-> None, Boxed -> False]

3.  Circles Defined by 2-Dimensional Euclidean Vector Subspaces

In Möbius geometry a circle in the n-sphere, or by stereographical projection, a circle or a line in the Euclidean n-space, is given by a 3-dimensional pseudo-Euclidean vector subspace: the generators of the isotropic cone of this 3-space are the (projective) points of the circle. Equivalently it is defined by the orthogonal complement, that is a (n-1)-dimensional Euclidean subspace of the (n+2)-dimensional pseudo-Euclidean vector space. The sets of these spaces are open orbits of the Möbius group in the corresponding  Grassmann manifolds. The Möbius group can be identified with two components of the pseudo-othogonal group O(1,4). In this section we construct  Mathematica  functions establishing  the equivariant bijection between the manifold of Euclidean 2-spaces of the 5-dimensional pseudo-Euclidean vector space of index 1 and the manifold of circles in the 3-dimensional Möbius space explicitely. In the next section we apply this bijection to calculate the Möbius invariants of pairs of circles. Elementary Möbius geometry is presented in [OS3], Section 2.7.

3.1 From the Subspace to the Circle

In this subsection we start with an arbitrary basis of an Euclidean 2-subspace, construct a pseudo-orthonormal basis of its orthogonal complement and obtain, using this basis and stereographical projection,  parameter representation and plots of the corresponding circle.

3.1.1. The Definition pscomplement

Let v, w be a pair of  vectors spanning a subspace of a pseudo-Euclidean vector space. Then one obtaines a pseudo-orthonormal basis of tits  orthogonal sub space as follows:

?dual

?pscomplement

pscomplement[v_List]:=chsort[orthonorm[NullSpace[Table[dual[v[[i]]],{i,1,Length[v]}]]]]

pscomplement[v_List,w_List]:=Module[{mm},
	mm= pscomplement[{v,w}];
	If[ch[mm[[1]]]<0, Return[mm], Print["Arguments don't span an Euclidean subspace!"];Return[],
 Return[mm]]]

pscomplement[{stb[1]+stb[5],stb[5]}]

pscomplement[{stb[1]+stb[5],stb[2]}]

pscomplement[stb[1]+stb[5],stb[2]]

pscomplement[stb[3],stb[2]]

Now we test pscomplement wit random matrices

?randommatrix

rm = randommatrix[3]

Table[Chop[pssp[%[[i]],%[[j]]]],{i,1,Length[%]},{j,1,Length[%]}]

nsp=pscomplement[rm]

Table[Chop[pssp[%[[i]],%[[j]]]],{i,1,Length[%]},{j,1,Length[%]}]

orthonorm[rm]

Table[Chop[pssp[%[[i]],%[[j]]]],{i,1,Length[%]},{j,1,Length[%]}]

Table[Chop[Simplify[pssp[rm[[i]],nsp[[j]]]]],{i,1,Length[rm]},{j,1,Length[nsp]}]

The functions work in pseudo-Euclidean spaces of higher index, too:

dim= 7; ind = 3;

rm = randommatrix[4]

orthonorm[rm]

Table[Chop[pssp[%[[i]],%[[j]]]],{i,1,Length[%]},{j,1,Length[%]}]

nsp=pscomplement[rm]

Table[Chop[pssp[%[[i]],%[[j]]]],{i,1,Length[%]},{j,1,Length[%]}]

Table[Chop[Simplify[pssp[rm[[i]],nsp[[j]]]]],{i,1,Length[rm]},{j,1,Length[nsp]}]

pscomplement[{{1,0,0,0,0,0,1},{0,1,0,0,0,1,0}}]

Table[Chop[pssp[%[[i]],%[[j]]]],{i,1,Length[%]},{j,1,Length[%]}]

Evaluate the next 5 cells several times :

rm1 = randomv[];rm2=randomv[];nn={rm1,rm2};orthonorm[nn]

Table[Chop[pssp[%[[i]],%[[j]]]],{i,1,Length[%]},{j,1,Length[%]}]

pscomplement[nn]

Table[Chop[pssp[%[[i]],%[[j]]]],{i,1,Length[%]},{j,1,Length[%]}]

pscomplement[rm1,rm2]

Clear[rm,rm1,rm2,nn,nsp]

dim = 5; ind = 1;

Remark.  For two equal vectors one obtains  correct results:

pscomplement[stb[2],stb[2]]

pscomplement[{stb[2],stb[2]}]

The arguments v,w of pscomplement, circle2spv, plotcircle2spv  need not to be orthogonal or normed:

pscomplement[stb[2]*c, stb[2]+stb[3]]

a= {t,1,0,0,0}; b =stb[2]*u + stb[5]*v;

pscomplement[a,b]

Simplify[Table[pssp[%[[i]],%[[j]]],{i,3},{j,3}]]//MatrixForm

WARNING! If one considers symbolic expressions for the given vectors v,w, one has to check whether they really span an Eucldean subspace, In our example this is not the case e.g. for t = 1, or some valuea of u,v in case t = .9:

pssp[a,a]*pssp[b,b] -pssp[a,b]^2

%/.{t->.9, u-> 2, v->.1}

For these arguments the vectors don't span an Euclidean subspace, what is hidden by the symbolic calculation:

pscomplement[{.9,1,0,0,0},{0,2,0,0,.1}]

3.1.2.  The Plot Command:  threepoints, circle2spv, plotcircle2spv

The following Modules yield the parameter representation of  the corresponding circle in the Euclidean 3-space, and the corresponding plot command; all these definition are protected and contained in mcirc.m.

?circle2spv

circle2spv[v_List,w_List][t_]:= Module[{psc},
    psc = pscomplement[v,w];
    If[psc == Null, Return[]];
      stepro[-psc[[2]]*Cos[t] -psc[[3]]*Sin[t] + psc[[1]]]];

?plotcircle2spv

plotcircle2spv[v_List,w_List,opt___]:= ParametricPlot3D[circle2spv[v,w][t], {t,0,2*Pi}, opt]

Another way: we take three  points of the circle corresponding to the subspace [v,w].

?threepoints

threepoints[v_,w_]:= Module[{b=pscomplement[v,w]},
If[b == Null, Return[];Break[]];
		{stepro[b[[3]]+  b[[1]]], stepro[b[[3]] - b[[1]]], stepro[b[[1]]+ b[[2]]]}]

Combining threepoints with circle3pts, plotcircle3pts of section 2 one obtaines the parameter representation circle3pts2spv[v,w][t], and the corresponding plotting procedure plotcircle3pts2spv:

circle3pts2spv[v_,w_][t_]  := circle3pts[threepoints[v,w][[1]],threepoints[v,w][[2]], threepoints[v,w][[3]]][t]

plotcircle3pts2spv[v_,w_,opt___]:=ParametricPlot3D[Evaluate[circle3pts2spv[v,w][t]],{t,-Pi,Pi},opt]

The possible options are those of ParametricPlot3D:

Options[ParametricPlot3D]

3.1.3. First Examples

Timing[circle2spv[stb[2], stb[3]][t]]

Timing[circle2spv[stb[2], stb[4]][t]]

Timing[circle2spv[stb[3], stb[4]][t]]

Timing[circle2spv[stb[2], stb[5]][t]]

ParametricPlot3D[{0,-Sin[t],-Cos[t]},{t, -Pi,Pi}]

circle2spv[stb[1], stb[3]][t]

?circle3pts2spv

?threepoints

threepoints[stb[2], stb[3]]

pscomplement[stb[2], stb[3]]

Timing[circle3pts2spv[stb[2], stb[5]][t]]

plotcircle3pts2spv[stb[2], stb[3]]

Timing[circle2spv[stb[3], stb[4]][t]]

Now we get the usual  parameter representation of the same line

Timing[circle3pts2spv[stb[3], stb[4]][t]]

Timing[plotcircle2spv[stb[3], stb[4]]]

Timing[plotcircle3pts2spv[stb[3], stb[4]]]

Timing[plotcircle3pts[stb[1,3], stb[2,3],-stb[1,3]]]

The next example shows random circles through the point pair {{-1,0,0},{1,0,0}}:

Table[plotcircle3pts[stb[1,3],- stb[1,3],randomv[3]], {20}];

Show[%,PlotRange->All]

Now we construct a uniform family of circles through the point pair {{-1,0,0},{1,0,0}}:

Quiet[Table[plotcircle3pts[stb[1,3],- stb[1,3],{0,Cos[t*Pi], Sin[t*Pi]}] , {t,0,1, Pi/50}]];

Show[%]

cirbdl = %;

3,1,4. The Frame Circles

The unit circle in the x,y-plane, centered at o:

grc1=plotcircle2spv[stb[4], stb[5],PlotStyle ->{Orange,Thick}]

This is the circle of radius 1 in the x, z-plane with center at the origin:

grc2=plotcircle2spv[stb[3],stb[5],PlotStyle->{Thick,Orange}]

This is the circle of radius 1 in the y, z-plane with center at the origin:

grc3=plotcircle2spv[stb[2],stb[5],PlotStyle->{Thick,Orange}]

circle2spv[stb[2],stb[5]][t]

This is the x-axis:

grc4=plotcircle2spv[stb[3],stb[4],PlotStyle->{Thick,Blue}]

This is the y-axis:

grc5=plotcircle2spv[stb[2],stb[4],PlotStyle->{Thick,Blue}]

This is the z-axis:

grc6=plotcircle2spv[stb[2],stb[3],PlotStyle->{Thick,Blue}]

Show[{grc1,grc2,grc3,grc4,grc5,grc6}, PlotRange -> {{-1.3,1.3},{-1.3,1.3},{-1.3,1.3}},Boxed -> False, Axes -> False,AspectRatio -> Automatic]

framecir = %;

3.1.5. Random Circle, Generated by a 2-Space Spanned by orthopair[random vector]

Clearly, the following vector must be spacelike with high probability:

bvec={2 *RandomReal[],-2 *RandomReal[],1.,3 *RandomReal[],RandomReal[]};ch[bvec]

?orthopair

Starting with this random vector we choose an orthopair, scetch the corresponding circle and the corresponding thrreepoints:

wr=orthopair[{2 RandomReal[],-2 RandomReal[],1.`,3 RandomReal[],RandomReal[]}];v=wr[[1]];w=wr[[2]];mwr=threepoints[v,w];pa=mwr[[1]];pb=mwr[[2]];pc=mwr[[3]];gr1=plotcircle3pts[pa,pb,pc];gr2=Graphics3D[{PointSize[0.03`],Point[pa],Point[pb],Point[pc]}];Show[{gr1,gr2}]

Here is the same procedure starting with a random vector randomv[]; if this vector is not spaselike error messages will be generated. Repeat the procedure severaltimes!

wr=orthopair[randomv[]];v=wr[[1]];w=wr[[2]];mwr=threepoints[v,w];pa=mwr[[1]];pb=mwr[[2]];pc=mwr[[3]];gr1=plotcircle3pts[pa,pb,pc];gr2=Graphics3D[{PointSize[0.03`],Point[pa],Point[pb],Point[pc]}];Show[{gr1,gr2}]

Remove[wr,v,w,mwr,pa,pb,pc,gr1,gr2, bvec]

3.2 From the Circle to the Subspace

For the geometrical interpretations of the concepts and results defined or calculated in the pseudo-Euclidean 5-space it is interesting to have expressions of the spacelike vectors as functions of the center {x,y,z}, the radius r, and the angles u,v which define the position vector of the circle in Euclidean 3-space. In this and the following sections we find and apply such functions.  We begin with a very simple, direct definition which is best appropriate for all applications. Defining a distinguished, from the Euclidean point of view, basis of the Euclidean 2-space corresponding to a given circle, the vectors of which depend on the mentioned  Euclidean invariants we get what we wanted.  The choice of such a basis means the choice of two spheres which intersect orthogonally in the given circle.  There are two such "spheres" canonical associated to the circle: The plane  containing the circle, and the sphere with the same radius and center containing the circle. Let us consider an arbitrary circle

Clear[a,b,c,x,y,z,u,v,r,t]

circle3D[x,y,z,u,v,r][t]

The spacelike vectors defining the associated spheres are orthogonal and span the spacelike 2-space of the circle. The vectors belonging to general planes and spheres are calculated in the notebook [emg1]. For a plane we have

?planevec

planevec[a,b,c,p]

?spt

spt[stb[5] + stb[1]]

stepro[stb[5] + stb[1]]

pssp[vec[nv],stb[5] + stb[1]]

pssp[planevec[a,b,c,p],stb[5] + stb[1]]

We know that the property nv[[1]] == nv[[5]] is characteristic for those spacelike vectors nv which represent planes. The reason is that all the subspheres through the north pole spt[stb[1] + stb[5]] of the 3-sphere, and only these, are projected onto planes by the stereographic projection.  Since we have the general identity

{a,b,c}.{a1,b1,c1} == pssp[{p,a,b,c,p},{p1,a1,b1,c1,p1}]

we may define the normed vector nv easily if we find Hesse's normal form

{a,b,c}.vec[x,3] == p

of the plane containing the circle, with

{a,b,c}.{a,b,c} == 1

Here vec[x,3] denotes the running point of the plane.  We write down the equation of the plane containing the circle with center {x, y, z} and position vector  unit[2][{u, v}] :

(vec[x,3]-{x,y,z}).unit[2][{u,v}] == 0

and recognize

{a,b,c} ==  unit[2][{u,v}]

p ==  unit[2][{u,v}].{x,y,z}

This leads to the definition

circleplane[x_,y_,z_,u_,v_] :={x Cos[u] Cos[v]+y Cos[v] Sin[u]+z Sin[v],Cos[u] Cos[v],Cos[v] Sin[u],Sin[v],x Cos[u] Cos[v]+y Cos[v] Sin[u]+z Sin[v]}

Simplify[pssp[circleplane[x,y,z,u,v],circleplane[x,y,z,u,v]]]

?circleplane

circleplane[0,0,0,0,0]

euklidsphereplot3D[circleplane[0,0,0,0,0]]

gp1 = %;

We emphasize that circleplane does not depend on the radius.  The second vector we know already:

?spherevec

Thus, a  vector corresponding to the sphere with the same radius and center as the circle does not depend on the position variables u,v; it is

spherevec[x,y,z,r]

circlespacevectors[x_,y_,z_,u_,v_,r_] :={circleplane[x,y,z,u,v],spherevec[x,y,z,r]}

circlespacevectors[0,0,0,0,0,1]

euklidsphereplot3D[%[[2]]]

Show[{%,gp1}]

plotcircle2spv[circlespacevectors[0,0,0,0,0,1]]

Clear[a,b,c,p]

circlespacevectors[x,y,z,u,v,r]

Simplify[Table[pssp[%[[i]],%[[j]]],{i,2},{j,2}]]

Obviously, the circlespacevectors are explicitly defined in terms of  Euclidean equivariants of the circle: center, position vector, radius. By its geometric definition as spacelike unit vectors such that the first corresponds to the plane of the circle they form  an orthonormal basis of the associated spacelike 2-space which is uniquely defined up to multiplications with  -1. If they have the mentioned property we call them adapted spacelike vectors of the circle. For a line any "sphere" containing it must be a plane. Therefore each vector v of the corresponding spacelike 2-space has the characteristic property (Subscript[v, 1]) == (Subscript[v, 5]), and any orthonormal basis of this 2-space will be considered as adapted to the line.

Now we show that the Euclidean plotcircle3D defined in section 2 and the plot generated by the spacelike vectors coincide:

?plotcircle3D

rv = randomv[6]

tv = circlespacevectors[rv[[1]],rv[[2]], rv[[3]], rv[[4]], rv[[5]], 10*rv[[6]]]

plotcircle3D[rv[[1]],rv[[2]], rv[[3]], rv[[4]], rv[[5]], 10*rv[[6]]]

plotcircle2spv[tv[[1]], tv[[2]]]

Remove[rv,tv]

Example: A circular cylinder

We consider the circles generating a cylinder of radius 1 around the z-axis:

zz[u_]:= circlespacevectors[0,0,u,0,Pi/2,1]

zz[u]

Table[pssp[zz[u][[i]], zz[u][[j]]], {i,2}, {j,2}]

Simplify[%]

circle2spv[zz[u][[1]],zz[u][[2]]][t]

Timing[Table[plotcircle2spv[zz[u][[1]],zz[u][[2]]],{u, 0,5,1/2}]]

Show[%[[2]],PlotRange ->All]

Plotting the cylinder:

Timing[ParametricPlot3D[Evaluate[circle2spv[zz[u][[1]],zz[u][[2]]][t]], {u, 0.,5},{t, 0.01, 2 Pi}, PlotPoints -> 50]]

Remove[zz]

3.3. Adapting the Frame

Often one meets the following situation: one has the Euclidean 2-space given by an arbitrary basis, and one wants to find the corresponding circle. One way to do that is the following: Transform the given frame in an adapted frame (defined in the last section), and then read of the Euclidean parameters from the adapted frame with the help of which one may get the parameter representation and the plot of the circle.  The adaption of the frame may be done by the following module:

?adaptsplframe

adaptsplframe[vx_List,vy_List] := Module[{aob,d1,d2,q = 0,v,ortho,nor},
    If[Simplify[pssp[vx, vx]*pssp[vy, vy]-pssp[vx, vy]^2<= 0],
      Print["The span is not 2-dimensional Euclidean!"];Return[],Goto[ortho],
	 Print["May be the span is not 2-dimensional Euclidean!"];Return[]];
	Label[ortho];
	aob = orthonorm[{vx,vy}, normed -> False];
    d1=Simplify[aob[[1]][[5]] - aob[[1]][[1]]];
    d2= Simplify[aob[[2]][[5]] - aob[[2]][[1]]];
    If[d1== 0, q=1;Goto[nor]];
    If[d2==0,v=aob[[1]]; aob[[1]] = aob[[2]];aob[[2]]= v;q=1];
    Label[nor];
    aob[[1]] = aob[[1]]/Sqrt[pssp[aob[[1]], aob[[1]]]];
    aob[[2]] = aob[[2]]/Sqrt[pssp[aob[[2]], aob[[2]]]];
    If[q==1,Return[Simplify[aob]]];
    d1=Simplify[aob[[1]][[5]] - aob[[1]][[1]]];
    d2= Simplify[aob[[2]][[5]] - aob[[2]][[1]]];
    q = Sqrt[d1^2 + d2^2];
    v = aob[[1]]; 
    aob[[1]] = Simplify[ -v*d2/q + aob[[2]]*d1/q];
    aob[[2]] =Simplify[ v*d1/q + aob[[2]]*d2/q];
    Return[aob]]

Tests

adaptsplframe[stb[1], stb[5]]

adaptsplframe[stb[5], stb[5]]

adaptsplframe[stb[5], a*stb[1]+ b*stb[5]]

adaptsplframe[stb[3], stb[4]]

adaptsplframe[{0,Sqrt[2]*Cos[t],Sqrt[2]*Sin[t],0,1},{0,0,0,0,1}]

adaptsplframe[stb[5],{0,Sqrt[2]*Cos[t],Sqrt[2]*Sin[t],0,1}]

adaptsplframe[stb[5],{0,Cos[t],Sin[t],0,1}]

adaptsplframe[stb[5],{a,Cos[t],Sin[t],0,a}]

adaptsplframe[stb[5],{1,Cos[t],Sin[t],0,1}]

adaptsplframe[{Sinh[s],0,0,0,Cosh[s]},{a,Cos[t],Sin[t],0,a}]

mm = {{Sinh[s],0,0,0,Cosh[s]},{a,Cos[t],Sin[t],0,a}};

dd=pssp[mm[[1]],mm[[1]]]*pssp[mm[[2]],mm[[2]]]-pssp[mm[[1]],mm[[2]]]^2

FullSimplify[%]

In general, the span needs not to be Euclidean, thus the Module yields Null.

In the next example adaptsplframe changes the order and normalizes only:

adaptsplframe[{0,0,0,Sqrt[1+s^2],s},{0,Cos[t],Sin[t],0,0}]

vx = randomv[]; vy = randomv[]; ot =adaptsplframe[vx,vy];ot//MatrixForm

ot1 =adaptsplframe[-vx,vy]

Table[pssp[%[[i]], %[[j]]], {i,2}, {j,2}]

plotcircle2spv[vx,vy,PlotRange -> All]

plotcircle2spv[ot[[1]],ot[[2]],PlotRange -> All]

euklidsphereplot3D[ot[[1]],-4,4,-4,4,Mesh-> None]

eb = %;

euklidsphereplot3D[ot[[2]],-1,1,-1,1,PlotStyle -> Opacity[.4], Mesh -> None]

sp = %;

?plotcircle2spv

plotcircle2spv[ot[[1]],ot[[2]],PlotStyle -> {Red,Thick},PlotRange-> All]

gc=%;

Show[{eb,sp,gc},PlotRange -> All]

Clear[vx,vy,ot,s1,t1,ep,eb,gc]

3.4.  Radius, Center, and Position Vector as Functions of the Subspace

In this section we calculate the basic Euclidean parameters radius, center, and position vector of the circle as functions of the  vectors v,w spanning an Euclidean 2-space in  the pseudo-Euclidean 5-space. There are two ways to calculate these invariants: one can use the adaption procedure adaptsplframe, or one can use threepoints, and the concepts of section 2.3. For numerical calculations the speed of both methods is nearly the same. We have chosen the first method with the corresponding definitions as standard concepts, and included them into the package mcirc.m, since it seems to be better appropriate for symbolic calculations.

3.4.1. The Euclidean invariants with Frame Adaption

Clear[x,y,z,u,v,r]

?circlespacevectors

circlespacevectors[x,y,z,u,v,r]

posvec[vx_List,vy_List] := Module[{ab},
ab=adaptsplframe[vx,vy];
    If[ab[[2]][[1]] == ab[[2]][[5]], Return[{0,0,0}]];
    Return[Take[ab[[1]],{2,4}]]]

center[vx_List,vy_List] := Module[{ab},
ab=adaptsplframe[vx,vy];
    If[ab[[2]][[1]] == ab[[2]][[5]], 
      Return[{Infinity,Infinity,Infinity}]];
    Return[Take[ab[[2]],{2,4}]/(ab[[2]][[1]] - ab[[2]][[5]])]]

radius[vx_List,vy_List] := Module[{ab},
ab=adaptsplframe[vx,vy];
    If[ab[[2]][[1]] == ab[[2]][[5]], Return[Infinity]];
    Return[1/(ab[[2]][[5]] - ab[[2]][[1]])]]

3.4.2. The Euclidean Invariants with threepoints

We take three points of circle2spv[v,w][t] with fixed parameter values, and calculate the invariants on the base of these three points applying the functions defined in section 2.3.

??threepoints

posvec1[v_List,w_List]:= Module[{p1,p2,p3,b},
		b=threepoints[v,w];
 If[b==Null,Return[];Break[]];
p1= b[[1]];
If[p1== {∞,∞,∞}, Return[{0,0,0}]];
		p2 = b[[2]];
If[p2== {∞,∞,∞}, Return[{0,0,0}]];
		p3 =  b[[3]];
If[p3== {∞,∞,∞}, Return[{0,0,0}]];
		posvec3pts[p1,p2,p3]]

center1[v_List,w_List]:= Module[{p1,p2,p3,b},
		b=threepoints[v,w];
 If[b==Null,Return[];Break[]];	
p1= b[[1]];
If[p1== {∞,∞,∞}, Return[p1]];
		p2 = b[[2]];
If[p2== {∞,∞,∞}, Return[p2]];
		p3 =  b[[3]];
If[p3== {∞,∞,∞}, Return[p3]];
	center3pts[p1,p2,p3] ]

radius1[v_List,w_List]:= Module[{p1,p2,p3,a, b},
			b=threepoints[v,w];
 If[b==Null,Return[];Break[]];	
		p1= b[[1]];p2 = b[[2]];p3 =b[[3]];
If[p1== {∞,∞,∞}||p2== {∞,∞,∞}||p3== {∞,∞,∞}, Return[Infinity]];
a = center3pts[p1,p2,p3];
		If[a == {Infinity,Infinity,Infinity},Return[Infinity],Sqrt[(a-p1).(a-p1)],Sqrt[(a-p1).(a-p1)]]]

3.4.3.  Standard Tests

posvec[stb[3],stb[4]]

posvec1[stb[3],stb[4]]

posvec[stb[5],stb[3]]

posvec1[stb[5],stb[3]]

center[stb[3],stb[4]]

center1[stb[3],stb[4]]

center[stb[2],stb[5]]

center1[stb[2],stb[5]]

radius[stb[3],stb[4]]

radius1 gives only non negative values, see the definition.

radius1[stb[3],stb[4]]

radius[stb[4],-stb[5]]

radius[stb[4],stb[5]]

radius1[stb[4],-stb[5]]

radius1[stb[4],stb[5]]

radius[stb[3],stb[5]]

radius[-stb[3],stb[5]]

radius[stb[3],-stb[5]]

radius1[stb[3],stb[5]]

3.4.4. Test with Random Space Vectors

wr=orthopair[{2 RandomReal[],-2 RandomReal[],1.`,3 RandomReal[],RandomReal[]}];v=wr[[1]];
w=wr[[2]];{v,w}

Timing[posvec[v,w]]

Timing[posvec1[v,w]]

%-%%

Timing[center[v,w]]

Timing[center1[v,w]]

Chop[%-%%]

Timing[radius[v,w]]

Timing[radius1[v,w]]

Chop[%-%%]

Clear[v,w]

Calculations via threepoint:

wr=orthopair[{2 RandomReal[],-2 RandomReal[],1.`,3 RandomReal[],RandomReal[]}];v=wr[[1]];w=wr[[2]];mwr=threepoints[v,w];
pa=mwr[[1]];pb=mwr[[2]];pc=mwr[[3]];
Timing[pv=posvec1[v,w];pm=center1[v,w];Print["Position vector = ",pv,", ",pv.pv,", ","center = ",N[pm],", ","radius = ",radius1[v,w]];gr1=plotcircle2spv[v,w];gr2=Graphics3D[{PointSize[0.05`],Point[pa],Point[pb],Point[pc],Point[pm]}];gr3=Graphics3D[Arrow[{pm,smoothing[pm+pv]}]];Show[{gr3,gr1,gr2}]]

Calculations via adaptsplframe, for the same vectors:

Timing[pv =  posvec[v,w];
pm =  center[v,w];
Print["Position vector = ", pv,", ",pv.pv,", ","center = ",N[pm],", ","radius = ",radius[v,w]];
gr1=plotcircle2spv[v,w];
gr2 = Graphics3D[{PointSize[0.05],Point[pa],Point[pb],Point[pc],Point[pm]}];
gr3 = Graphics3D[Arrow[{pm, smoothing[pm +pv]}]];
Show[{gr3,gr1,gr2}]]

The time difference is very small: 0.027 s.

radius[v,w]

radius[-v,-w]

Clear[x,y,z,u,v,r,wr,v,w,mwr,pa,pb,pc,pv,pm,gr1,gr2,gr3]

The next example shows a random circle with the canonical adapted plane and sphere:

?euklidsphereplot3D

s1=RandomReal[];t1=RandomReal[] π;
ot=adaptsplframe[{0,0,0,(Sqrt[1+(s1^2)]),s1},{0,Cos[t1],Sin[t1],0,0}];
gt1=plotcircle2spv[ot[[1]],ot[[2]],PlotStyle ->{Black, Thick}];
gt2=euklidsphereplot3D[ot[[1]],-Abs[radius1[ot[[1]],ot[[2]]]],3 Abs[radius1[ot[[1]],ot[[2]]]],-2 Abs[radius1[ot[[1]],ot[[2]]]],2 Abs[radius1[ot[[1]],ot[[2]]]] , Mesh -> None , PlotStyle-> {Blue,  Opacity[.4]}];
gt3=euklidsphereplot3D[ot[[2]], -1,1,-1,1, Mesh ->None,PlotStyle ->{Red, Opacity[.4]}];
Show[{gt1,gt2,gt3},PlotRange -> All]

gtest = %;

Turn the image with a left mouse click!

3.5.  From  Three Points to the Subspace of the Corresponding Circle

3.5.1. Definitions

Let p1, p2, p3 be three pairwise distinct points in the Euclidean 3-space. Their preimages under the stereographic projection define uniquely a circle on the 3-sphere, and the corresponding isotropic vectors span a pseudo-Euclidean 3-space, the orthogonal complement of which is the Euclidean subspace we ask for.

?vspace

3.5.2. Examples

vspace[stb[1,3], {0,0,0}, {2,0,0}]

plotcircle2spv[%[[1]],%[[2]]]

vspace[stb[1,3], {1,0,0}, {2,0,0}]

vspace[stb[1], {1,0,0}, {2,0,0}]

Timing[plotcircle3pts[{1,0,0}, {0,1,0},{1,-1,0}]]

vspace[{0,1,0},{1,0,0}, {1,-1,0}]

Timing[plotcircle2spv[%[[1]],%[[2]]]]

vs=%%;

radius[vs[[1]],vs[[2]]]

N[%]

N[center[vs[[1]],vs[[2]]]]

% -{1, 0, 0}

Sqrt[%.%]

N[posvec[vs[[1]],vs[[2]]]]

vspace[{1,0,0}, {0,0,0},{-1,0,0}]

vs = %

plotcircle2spv[vs[[1]],vs[[2]]]

radius[vs[[1]],vs[[2]]]

center[vs[[1]],vs[[2]]]

plotcircle3pts[{1,0,0}, {0,0,0},{-1,0,0}]

posvec[vs[[1]],vs[[2]]]

a=RandomReal[];b=2 RandomReal[];c=-RandomReal[];Print[TableForm[{a,b,c}]];Timing[cr =plotcircle3pts[{a,0,0},{0,b,0},{0,0,c}]]

Show[{cr,Graphics3D[{PointSize[.03],Map[Point,{{a,0,0},{0,b,0},{0,0,c}}]}]}]

Timing[vspace[{a,0,0}, {0,b,0},{0,0,c}]]

vs = %[[2]]

Timing[cr =plotcircle2spv[vs[[1]],vs[[2]]]]

Show[{cr,Graphics3D[{PointSize[.02],Map[Point,{{a,0,0},{0,b,0},{0,0,c}}]}]}]

It results what could have been expected: The direct Euclidean way constructing the circle through 3 points is much faster than constructing the circle by the help of vspace.

Clear[a,b,c,vs]

4.  Möbius Invariants for Pairs of Circles

In this section we apply the general theory, see the paper [S00],  or the book [OS3], Section 2.7,  to the case of circles in the 3-sphere and the Euclidean 3-space. For the Möbius geometry of spheres see the notebook  emg.nb. Invariants for the other possible pairings ((S^k), (S^m)), 0 ≤ k, m ≤ 3, are treated in the notebook emg3.nb. The expressions of the Möbius invariants of two circles by their Euclidean invariants obtained in subsection 4.8 are analogeous to the Coxeter invariant of two spheres, see the notebook emg.nb. They seems to be new (see also [OS3], subsection  2.7.4).

4.1 Case Discussion

Notations: Let C0, C1 be two distinct circles, W0, W1 the corresponding 3-dimensional pseudo-Euclidean subspaces, U0, U1 their orthogonal (Euclidean) complements.

Problem: Express the geometric properties and invariants of C0, C1 by algebraic properties and invariants of U0, U1.

Since 
                                 Dim(W0⋂W1) + Dim(W0 + W1) = Dim W0 + Dim W1= 6,		(1)
we distinguish (remember: orthogonal complement[W0⋂W1] = U0 + U1)
1.  Dim(W0⋂W1) = 1	<=> 	 Dim(U0 +U1) = 4,				
2. Dim(W0⋂W1) = 2	<=>	 Dim(U0 +U1) = 3.					  	(2)
The geometric meaning of these cases is, because of (1),
Case 1 <=> Dim(W0 + W1) = 5 <=> there don’t exist any 2-sphere containing C0⋃C1 
			        <=> Dim(U0⋂U1) = 0.					  	(3)
Case 2 <=> Dim(W0 + W1) = 4 <=> there exists a 2-sphere containing C0⋃C1  
    			         <=>  Dim(U0⋂U1) = 1.					  	(4)
In each of these cases we have three subcases with obvious geometric interpretations:
a) W0⋂W1 Euclidean <=> C0⋂C1 is empty and C0, C1 can be separated by spheres 
		          <=> U0 + U1 is pseudo-Euclidean, 4-dimensional.			(5)
Case 1 and b)
W0⋂W1 pseudo-Euclidean <=> C0⋂C1 is empty, C0, C1 are interlaced
  			      <=> U0 + U1 is Euclidean, 4-dimensional.			(6)
case 2 and b) 
W0⋂W1 pseudo-Euclidean <=> C0⋂C1 ={p0,p1} - two points 
			    <=> U0 + U1 is Euclidean, 3-dimensional.			   	(7)
c) W0⋂W1 isotropic <=> C0⋂C1 = {p} - one point <=> U0 + U1 is isotropic.		(8)

Here "separated by spheres" means that there exist disjoint spheres S0, S1 with C0⊂S0, C1⊂S1; "interlaced" means that C0⋂C1 is empty, and such spheres do not exist: any two spheres S0,S1 with C0⊂S0, C1⊂S1 must intersect.
In the following we discuss the Möbius invariants of the pairs in each of the cases and give examples.

4.2 Invariants for Pairs of  Euclidean 2-Subspaces

In Section 3 we established an equivariant bijection between  the space of the circles of the Möbius space and the space of  2-dimensional Euclidean subspaces of the corresponding pseudo-Euclidean vector  space. Thus pairs of circles correspond bijectively to pairs of such subspaces. In  [S00] or [OS3], Section 2.7, is proved that the Eigenvalues of a certain selfadjoint operator of one of the Euclidean 2-spaces are a complete invariant system for the pairs of subspaces and therefore also for the pairs of circles.  In Subsection 4.2.5 below we define this operator, named the double projection,  using the orthogonal decompositions of the pseudo-Euclidean vector space defined by the subspaces, and construct functions for the calculation of the invariants. First we consider the cross product of orthonormal basis vectors of the subspaces and use it for a case criterion, which is illustrated by examples.

4.2.1. Case Criterion

If U0 and U1 are the Euclidean vector subspaces spanned by orthonormal bases :U0={v1, v2}, U1={w1, w2}, then an important equivariant function  for the case characterization is the pseudo-orthogonal cross product

C0xC1 := pscross[{v1, v2, w1, w2}].

From the general properties of the cross product one easily concludes: 

Lemma. C0xC1 is independent of the choice of the oriented orthonormal bases spanning U0 and U1; it changes the sign if the bases correspond to different orientations of U0 or U1. ch[C0xC1] depends only of the (non-oriented) circles C0,C1; it does not change under arbitrary, not necessarily orthogonal transformations of the bases..

Proof. The e+quivariance of the vector product is well known.  We transform the basis of U0 and calculate

Clear[x,y,u,v]; ftest[t_] := pscross[{vec[x] Cos[t]+ vec[y] Sin[t],-vec[x]Sin[t] +vec[y]Cos[t], vec[u], vec[v]}]

Simplify[D[ftest[t],t]]

or

Simplify[ftest[t]== ftest[0]]

The equivariance of the vector products yields the following invariant of the two oriented circles C0= C0[x,y], C1 =C1[u,v], where x,y, resp. u,v are oriented orthonormal bases of vspace[C0) resp. vspace[C1].

Clear[x,y,u,v]; mgramdet[x_,y_,u_,v_] :=pssp[pscross[{x,y,u,v}],pscross[{x,y,u,v}]]

By  (3) we have case 1 iff  C0xC1 does not vanish; from (5) to (8) and the orthogonality of C0xC1 to U0 + U1 we see, that the following disjunction holds:
case 1.a) <==> ch[C0xC1] = 1;
case 1.b) <==> ch[C0xC1]  = -1;
case 1.c) <==> ch[C0xC1] = 0;
case 2  <==> C0xC1 = o.
Remember:

?ch

Case 2 leads to

ch[null[]]

4.2.2. Example: Random Pair, Case 1.a

We define  four vector functions such that {va1[t], va2[s]}, {wa1[u], wa2[v]} are orthonormal pairs for any values of s,t,u,v:

(va1[t_] := stb[3]*Cos[t] + stb[4]*Sin[t]; va2 [s_]:= stb[1]* Sinh[s] + stb[2]* Cosh[s];
	wa1[s_] := stb[3]*Cosh[s] + stb[1]*Sinh[s]; wa2[t_] := stb[2]*Sin[t] + stb[4]* Cos[t];
C0xC1 = pscross[{va1[t],va2[s],wa1[u],wa2[v]}]);

Simplify[pssp[va1[t],va2[s]]]

Simplify[pssp[wa1[t],wa2[s]]]

C0xC1

This vector  is always  proportional to stb[5]; it vanishes iff the factors depend linearly

Simplify[pssp[%,%]]

chf[s_,t_,u_,v_]=mgramdet[va1[t],va2[s],wa1[u],wa2[v]]

For some “singular” values, e.g. s = 0 and  t=Pi/2,  we have case 2. For all other values of t,s we have case 1. a : Non-intersecting, separated circles.  In general we have

rank[{va1[t],va2[s],wa1[u],wa2[v]}]

ch[C0xC1]

since generally for variable x the character of the vector yields

ch[{0,0,0,0,x}]

in spite of the possible case x=0:

ch[{0,0,0,0,0}]

Hints:
By numerical reasons it sometimes happens that the random circles have "corners"; try to avoid that increasing the number of plotpoints.
Evaluate the following cell several times and rotate the picture.
If the difference of the radii of the circles is too large it may happen that one of them is too small to be visible.

s= RandomReal[] *Pi;t= RandomReal[] *Pi;u= RandomReal[] *Pi;v= RandomReal[] *Pi;
Print["chf = ",chf[s,t,u,v]]
gr1=plotcircle2spv[va1[s],va2[t],PlotStyle ->{Green,Thick}];
gr2=plotcircle2spv[wa1[u],wa2[v], PlotStyle ->{Red,Thick}];
Show[{gr1,gr2},PlotRange->All]

Clear[s,t,u,v,gr1,gr2]

4.2.3. Example: Random Pair, Case 1.b

We define other four vector functions such that {vb1[t], vb2[s]}, {wb1[u], wb2[v]} are orthonormal pairs for any values of s,t,u,v:

Clear[s,t,u,v]

(vb1[t_] := stb[3]*Cos[t] + stb[4]*Sin[t]; vb2 [s_]:= -stb[5]* Sin[s] + stb[2]* Cos[s];
	wb1[s_] := stb[5]*Cos[s] + stb[2]*Sin[s]; wb2[t_] := -stb[3]*Sin[t] + stb[4]* Cos[t];
C0xC1 = pscross[{vb1[t],vb2[s],wb1[u],wb2[v]}])

Simplify[%]

chfb[s_,t_,u_,v_]:=mgramdet[vb1[t],vb2[s],wb1[u],wb2[v]]

chfb[s,t,u,v]

Simplify[%]

It follows: If s-u≠(2k+1)ℼ/2 and t-v≠(2k+1)ℼ/2 then any two circles generated by the orthonormal pairs {v1[t], v2[s]}, {w1[u], w2[v]} are  interlaced. One can test it evaluating the next cell. 
Hints:
By numerical reasons it sometimes happens that the random circles have "corners"; try to avoid that increasing the number of plotpoints.
Evaluate the following cell several times and use mouse  click left to rotate the picture.

s= RandomReal[] *Pi;t= RandomReal[] *2*Pi;u= RandomReal[] *2*Pi;v= RandomReal[] *Pi;
Print["chfb = ",chfb[s,t,u,v]]
gr1=plotcircle2spv[vb1[s],vb2[t],PlotStyle -> {Green,Thick}];
gr2=plotcircle2spv[wb1[u],wb2[v], PlotStyle ->{Red,Thick}];
Show[{gr1,gr2},PlotRange->All]

Clear[s,t,u,v,gr1,gr2]

4.2.4. Example: Random Pair, Case 1.c

We define random circles intersecting at the origin:

(pa=randomv[3];pb = randomv[3]; qa = 2*randomv[3]; qb = 2*randomv[3];
	mvp = vspace[null[3],pa,pb]; mvq =vspace[null[3],qa,qb];
	Print["C0xC1 = ",Chop[pscross[{mvp[[1]],mvp[[2]],mvq[[1]],mvq[[2]]}]]];
		Print["ch = ",ch[pscross[{mvp[[1]],mvp[[2]],mvq[[1]],mvq[[2]]}]]];
	gr1 = plotcircle3pts[null[3],pa,pb,PlotStyle -> {Green,Thick}];gr2 = plotcircle3pts[null[3],qa,qb,PlotStyle -> {Red,Thick}];
	Show[{gr1,gr2}, PlotRange -> All])

Note that the vector C0xC1 for this example is always isotropic.

Clear[pa,pb,qa,qb,mvp, mvq, gr1,gr2]

4.2.5.  The Double Projection

Now we are going to describe the basic invariants for pairs of circles under the Möbius group. In the next sections we then shall give further examples for which we calculate these invariants. As for stationary angles of subspaces of the Euclidean space the invariants are found by considering the projections defined by the corresponding spacelike subspaces, see [S00], or [OS3]. The generalized stationary angles correspond to the stationary values of the scalar products  pssp[v,w] subject to the constraints pssp[v,v] = 1, pssp[w,w] = 1, v ∈ U0, w∈ U1.

pp[v1, v2, w1, w2]  is the matrix of the double projection pp= p2∘ p1, the composition of the orthogonal projections
		p1: U0 -> U1, and p2: U1 -> U0, 
with respect to the base {v1, v2}. Have in mind that the bases are assumed to be orthonormal.

?pp

pp[v1_, v2_, w1_, w2_] := 
    {{pssp[v1, w1]^2 + pssp[v1, w2]^2, 
      pssp[v1, w1]*pssp[v2, w1] + pssp[v1, w2]*pssp[v2, w2]}, 
     {pssp[v1, w1]*pssp[v2, w1] + pssp[v1, w2]*pssp[v2, w2], 
      pssp[v2, w1]^2 + pssp[v2, w2]^2}}

pp is a symmetric matrix with non negative numbers on the main diagonal:

pp[v1, v2, w1, w2]//MatrixForm

The invariants of pp are the determinant and the trace, or equivalently, its eigenvalues:

ppdet[v1_, v2_, w1_, w2_] := Det[pp[v1, v2, w1, w2]]

pptr[v1_, v2_, w1_, w2_] := Sum[pp[v1, v2, w1, w2][[i,i]], {i, 1, 2}]

Simplify[ppdet[v1,v2,w1,w2]]

We see that  ppdet  is the square of the determinant of the scalar products pssp[xi, yj], i, j = 1, 2, and therefore always non negative. The same is valid for the trace:

pptr[v1,v2,w1,w2]

Simplify[Eigenvalues[pp[v1, v2, w1, w2]]]

Snce the eigenvalues of a sesfadjoint operator (with a symmetric matrix with respecttoan orthonormal basis) always are real. the radicand of the expressions for the eigenvalues must be non-negativ. Indeed, we have

FullSimplify[pptr[v1,v2,w1,w2]^2-4*ppdet[v1,v2,w1,w2]]

and it follows immediately

vec[w1]

and the fact that the Möbius group acts transitively on the circles, we can always assume that the first circle is the unit circle grc1. Since vec[w1], vec[w2] are not orthonormal, in general, the following calculations give invariants under the additional assumption of the orthonormality of these vectors only.

Lemma.  The eigenvalues α1, α2  of pp satisfy  α1 ≥ α2 ≥ 0.

4.2.6. The Orthogonal Case

Two circles are said to be orthogonal, iff the corresponding Euclidean subspaces  U0, U1 are mutually orthogonal. Geometrically this means that the intersection angle of any two spheres containing each one of the circles equals Pi/2. Obviously, then we have case 1.b: the circles are interlaced, and the eigenvalues fulfil α1 = α2 = 0. Two circles are called isogonal, if the eigenvalues satisfy α1 = α2. The question about the existence  of non-orthogonal isogonal circles will be answered positively in section 4.4 with the construction of examples. Obviously, orthogonal circles exist:

gr1 = plotcircle2spv[stb[2], stb[3],PlotStyle ->{Green, Thick} ];Show[{gr1,grc1},Boxed -> False, Axes -> None]

pp[stb[2],stb[3],stb[4],stb[5]]//MatrixForm

We use spherical reflection to give an example of two circles equivalent to the considered pair.

?sphericalreflection

src1[t_,r_] := sphericalreflection[{Cos[t],Sin[t],0},{1/2,0,0},r]

src2[t_,r_] := sphericalreflection[{0,0,t},{1/2,0,0},r]

ParametricPlot3D[src1[t,.5],{t,0,2*Pi},PlotStyle -> {Orange,Thick}]

ci1 =  %;

ParametricPlot3D[src2[t,.5],{t,-40,40},PlotRange -> All,PlotStyle -> {Green,Thick}]

ci2 = %;

Show[{ci1,ci2},PlotRange -> All, Boxed-> False,Axes -> None]

We use vspace to find the orthonormal bases for the vector subspaces corresponding to the circles; to this aim we choose three points on each circle and apply:

?vspace

vv1 =vspace[src1[0,.5],src1[Pi,.5],src1[Pi/2,.5]]

vv2 =vspace[src2[0,.5],src2[2,.5],src2[-2,.5]]

The first basis vectors vv1[[1]], vv2[[1]] correspond to the x,y- resp. x,z-plane intersecting orthogonal. Obviously, also the other pairings of basis vectors are orthogonal, e.g.

Chop[pssp[vv1[[2]],vv2[[2]]]]

euklidsphereplot3D[vv1[[2]],-1,1,-1,1, Mesh -> False, PlotStyle-> Opacity[.3]]

sph1= %;

euklidsphereplot3D[vv2[[2]],-1,1,-1,1, Mesh -> False, PlotStyle-> Opacity[.3]]

sph2= %;

Show[{ci1,ci2,sph1,sph2},PlotRange -> All, Boxed-> False,Axes -> None]

4.2.7. Example: Two Non Orthogonal Circles Intersecting Orthogonally

We construct an example of two non orthogonal circles intersecting orthogonally. The Euclidean 2-space defining the unit circle with centre (1,0,0) in the x,y-plane is spanned by the orthonormal vectors

c100=circlespacevectors[1,0,0,0,Pi/2,1]

plotcircle2spv[%[[1]],%[[2]]]

cc1= %;

We consider this circle together with the z-axis:

plotcircle2spv[stb[2],stb[3]]

Show[{%,cc1}, PlotRange-> {{-1,2},{-1,1},{-1,1}}]

the matrix of the double projection has eigenvalues 0, 1:

pp[stb[2],stb[3],stb[4],c100[[2]]]

rank[{stb[2],stb[3],stb[4],c100[[2]]}]

Spheres, containing each one of these circles, are tangent to each other or intersect in an arbitrary angle:

Clear[a,b]; pssp[stb[2]*Cos[a] + stb[3]*Sin[a],stb[4]*Cos[b]+ c100[[2]]*Sin[b]]

Now we replace the z-axis by the x-axis, and calculate the double projection

plotcircle2spv[stb[4],stb[3]]

Show[{%,cc1}, PlotRange-> {{-1,2},{-1,1},{-1,1}}]

pp[stb[4],stb[3],stb[4],c100[[2]]]

The same matrix, but, of course,  the rank-condition is hurt:

rank[{stb[3],stb[4],stb[4],c100[[2]]}]

4.3 Case 1.a: Disjoint,  Non-Interlaced Circles

Take two disjoint spheres and a circle in each. There are lot of examples.

v1 = {1, 0, 1, -1, 0}; v2 = stb[2]; 
  w1 = {3/Sqrt[2], 1/Sqrt[2], 0, 
    1/Sqrt[2], 3/Sqrt[2]}; w2 = stb[3]; 
  C0xC1 = pscross[{v1, v2, w1, w2}]

ch[%]

It follows: case 1.a; the circles are disjoint and separated.

mgramdet[v1,v2,w1,w2]

ppdet[v1,v2,w1,w2]

pptr[v1,v2,w1,w2]

Here is the circle C0

Clear[C0,C1]

plotcircle2spv[v1,v2,PlotStyle -> {Black,Thick},PlotRange->{{-10.,10.},{-4.,4},{-4,8}}]

C0=%;

and here the "xircle" C1, a line:

plotcircle2spv[w1,w2, PlotStyle-> {Green, Thick},PlotRange->{{-10.,10.},{-4.,4},{-4,8}}]

C1=%;

Show[{C0, C1}]

pp[v1,v2,w1, w2]//MatrixForm

Eigenvalues[pp[v1,v2,w1, w2]]

eval=N[%]

The stationary values of the Möbius invariant of two spheres  S0,S1 with C0⊂S0, C1⊂S1 are the square roots of the eigenvalues (see below)

Sqrt[%]

As it follows from the general theory [S00], for case 1.a is characteristic that the maximal eigenvalue must be greater than 1, and has multiplicity 1. Therefore, in case 1.a may not exist isogonal circle pairs.

Eigenvectors[pp[v1,v2,w1, w2]]

Therefore the "eigenspheres" (spheres corresponding to the eigenvectors) containing the circle C0 are

eigenv1 = Simplify[normalize[v1*((1/4)+(1/8) (-19-(Sqrt[353]))) +v2]]

Simplify[pssp[%,%]]

eigenv2 =Simplify[ normalize[v1 *((1/4)+(1/8) (-19+(Sqrt[353]))) + v2]]

Simplify[pssp[%,%]]

Simplify[pssp[%%,eigenv1]]

They intersect orthogonally in the circle C0. We construct the corresponding picture.

euklidsphereplot3D[eigenv1,-1,1,-1,1,PlotStyle -> {Green,Opacity[.25]},Mesh -> None]

sp1=%;

euklidsphereplot3D[eigenv2,-1,1,-1,1,PlotStyle -> {Opacity[.2,Red]},Mesh -> None]

sp2 = %;

Show[{sp1,sp2, C0,C1},PlotRange -> {{-10.,10.},{-6.,6},{-6,5}}]

sp12 = %;

The spheres containing the circles C0, C1 are the one-parameter families S0[t]⊃C0, S1[s]⊃ C1, respectively, defined by the spacelike unit vectors

uv0[t_]:= eigenv1 Cos[t] +  eigenv2 Sin[t]; 
uv1[s_] := w1 Cos[s] + w2 Sin[s];

For the generalized cosinus of the intersection angle (see the notebook emg.nb) we calculate

Simplify[pssp[uv0[t], uv1[s]]]

fctr[t_,s_] = %;

Plot3D[fctr[t,s],{s,0, 2 Pi},{t,0,2 Pi}]

FindMaximum[fctr[t,s],{s,3},{t,3}]

fctr[Pi,2.8184857533995356`]

Sqrt[eval]

As stated above, the square root of the first eigenvalue coincides with the maximal generalized  cosinus of the intersection angle. The generalized angle is maximal for the spheres corresponding to -eigenv1 and

maxs =uv1[2.8184857533995356`]

euklidsphereplot3D[maxs,-10,10,-10,10,PlotStyle -> {Red,Opacity[.3]},Mesh -> None]

sphmax = %;

Show[{sphmax,sp1,C0,C1},PlotRange ->  {{-10.,10.},{-6.,6},{-6,5}}]

figmax=%;

figmax shows a sphere pair with maximal generalized cosinus. FindMinimum yields the negative square root of the first eigenvalue

FindMinimum[fctr[t,s],{s,0},{t,3}]

For t = Pi/2, 3Pi/2 we get correspondingly for the second eigenvector:

Plot[fctr[Pi/2,s ], {s, 0, 2 Pi}]

lt follows that all spheres containing C1 (this are planes in the Euclidean picture) intersect the sphere eigenv2 (See figure sp12).

D[fctr[Pi/2,s],s]

Simplify[N[%]]

Solve[%==0,s]

{fctr[Pi/2,-1.893903227357868`],fctr[Pi/2,1.2476894262319251`]}

Sqrt[eval[[2]]]

The absolute value is the sqare root of the second Eigenvalue. The generalized angle is minimal for the spheres corresponding to  eigenv2  and

mins=uv1[1.2476894262319251`]

euklidsphereplot3D[mins,-8,8,-8,8,PlotStyle -> Opacity[.6],Mesh-> None]

emin= %;

Show[{emin,sp2,C0,C1},PlotRange ->{{-10.,10.},{-6.,6},{-6,5}}]

figmin =%;

figmin shows a paitr of spheres minimazing the Coxeter invariant.

4.4 Case 1.b: Interlaced Circles

Already in section 4.2, the orthogonal case,  we have seen:  the standard unit circle in the x,y-plane and the z-axis are interlaced "circles". We give now examples of isogonal and non-isogonal interlaced circle pairs.

4.4.1. Example: Isogonal Circles

A first example  of an isogonal pair:

(w1 ={0,1,1,0,1}/Sqrt[3]; w2 = {0,1,-1,1,0}/Sqrt[3];
griso = plotcircle2spv[w1,w2])

Show[{grc1,griso},PlotRange -> Automatic]

isosph=%;

psCross[stb[4],stb[5],w1,w2]

ch[%]

We again got an interlaced pair of circles.  The unique eigenvalue is 1/3:

pp[stb[5], stb[4],w1, w2]

The spheres containing the circles C0, C1 are the one-parameter families S0[t]⊃C0, S1[s]⊃ C1 defined by the spacelike unit vectors

uv0[t_]:= stb[4] Cos[t] +  stb[5] Sin[t]; 
uv1[s_] := w1 Cos[s] + w2 Sin[s];

Any pair (S0[t], S1[s]]) intersects in a certain circle, since (C0,C1) are interlaced. For the cosinus of the intersection angle (see the notebook emg.nb) we calculate

Simplify[pssp[uv0[t], uv1[s]]]

fctr[s_,t_] = %;

Plot3D[fctr[s,t],{t,0,2 Pi}, {s,0, 2 Pi}]

We do not have distinguished eigenspheres: since the eigenvalue has multiplicity 2, any sphere S0[t] ⊃ C0 generates the same angle function with the family (S1[s]) (up to a phase translation). Therefore we called the circle pair (C0, C1) isogonal. In concordance with the  general theory  the maximum is the square root of the eigenvalue:

a=RandomReal[];Print[a];Plot[fctr[s,a],{s,-Pi,3*Pi}]

Plot[fctr[s,Pi/2],{s,-Pi,3*Pi}]

fctr[0, Pi/2]

N[%]

The Amplitude of the Cos of the intersection angle is the square root of the intersection angle, It does not depend on t, in difference to the behavior of these functions in the case of interlaced,but non isogonal circles, see below.
Since pssp[-n, m] = - pssp[n, m], we get all relative minimal values as the negatives of the relative maximal values, changing the orientation of one of the spheres.

4.4.2. Family of Circles,  Isogonal to the Unit Circle grc1

Now we construct a family of isogonal pairs containing examples for any eigenvalue between 0 and 1. By the general theory it follows, that any pair of isogonal circles is Möbius equivalent to a uniquely defined member of this family.
Argumentation: By theorem 6 of [S00] the eigenvalues are a complete system of invariants for pairs of circles not contained in a 2-sphere.  From the case discussion in this paper, cases 1 and 4, there exists at most one eigenvalue ≥ 1, and then it has multiplicity 1. Since the eigenvalues are all non-negative, in the isogonal case  the unique eigenvalue α fulfils 0 ≤ α < 1. For  α = 0 the circles are orthogonal, and this is a special case of isogonal. We shall see that in the family we are going to construct exists for any  α, 0 < α  < 1,  exactly one member.

iso1[s_] :={ 0, 1/Cosh[s] , 1/Sinh[s],0,1}/Sqrt[1 + 1/Cosh[s]^2+1/Sinh[s]^2];
		iso2[s_] := {0, 1/Sinh[s] , -1/Cosh[s],1,0}/Sqrt[1 + 1/Cosh[s]^2+1/Sinh[s]^2];
grisofam[s_] := {iso1[s],iso2[s]}

Simplify[grisofam[s] , s∈ Reals]

Simplify[pssp[iso1[s],iso2[s]],s>0]

Simplify[pssp[iso1[s],iso1[s]],s>0]

Simplify[pssp[iso2[s],iso2[s]],s>0]

Simplify[pp[stb[4], stb[5], iso1[s], iso2[s]]]

This diagonal matrix show that for any value of s the eigenvalues coincide; as a function of s the only eigenvalue is

f[s_] := 1/(1 + Csch[s]^2 + Sech[s]^2)

f[s] - (Cosh[s]^2 * Sinh[s]^2)/(Cosh[s]^2 + Sinh[s]^2 +Cosh[s]^2 *Sinh[s]^2)

Together[%]

Therefore avoiding the poles will give a better definition for f:

f[s_] := (Cosh[s]^2 * Sinh[s]^2)/(Cosh[s]^2 + Sinh[s]^2 +Cosh[s]^2 *Sinh[s]^2)

Plot[f[s], {s,0,5}]

f[0]

Any value val = f[s], 1> val> 0, is reached by a uniquely defined circle pair of our family.  The case s= 0 , f(0) = 0, does not give a circle since s = 0 is a singularity of grisofam:

grisofam[0]

Limit[f[s], s-> Infinity]

Limit[grisofam[s], s-> Infinity]

Limit[grisofam[s], s-> 0]

plotcircle2spv[%[[1]],%[[2]]]

The circlefamily joins the z-axis with the standard unit circle :

plotcircle2spv[stb[4],stb[5]]

Sometimes it is useful to create a circle family instead of a circle surface. In mcirc.m is defined

circlefamily[v_,w_,t_,a_,b_, d_:20, opt___] := Table[plotcircle2spv[v[t],w[t],opt],{t, a, b, (b-a)/d}]

circlefamily[iso1,iso2,t,0.1,1.8,10, PlotPoints->80]

isofam = {grc1,%};

Show[isofam,PlotRange ->  Automatic]

isofamily= %;

The circle family isofam is contained in the circle surface (needed more than 3 min):

Timing[ParametricPlot3D[circle2spv[grisofam[s][[1]], grisofam[s][[2]]][t],
	{s, 0.1,2.5}, {t,0, 2  Pi}, PlotPoints -> 40, PlotRange -> {{-5.,5.},{-5,5},{-10,10}}]]

isosf2 = %[[2]]

Show[{isosf2,grc1}]

Any two circles of the circlefamily grisofam are isogonal (and interlaced, e.g.

plotcircle2spv[grisofam[.4][[1]], grisofam[.4][[2]]]

isoc1 = %;

plotcircle2spv[grisofam[1.7][[1]], grisofam[1.7][[2]]]

isoc2=%;

Show[{isoc1,isoc2}, PlotRange -> Automatic]

isocs = %;

pp[grisofam[1.7][[1]], grisofam[1.7][[2]],grisofam[.4][[1]], grisofam[.4][[2]]]

Indeed, this can be proved by  symbolic calculation:

Simplify[pp[grisofam[s][[1]], grisofam[s][[2]],grisofam[t][[1]], grisofam[t][[2]]]]//MatrixForm

%[[1,1]] == %[[2,2]]

4.4.3. Intersection Angle for Spheres Containing  Interlaced, But Non Isogonal Circles

Clear[a,b,br1,br2,sph1, sph2]

Here is the family of spheres containing the  unit circle with center o in the x,y,plane:

sph1[a_] := {0,0,0,Cos[a], Sin[a]}

sph1ta=Table[euklidsphereplot3D[sph1[i*Pi/10],-1,1,-1,1,Mesh-> None,PlotStyle-> Opacity[.4]],{i,-5,5}];Show[sph1ta,PlotRange -> All]

fam1 = %;

Here are two orthonormal vectors defined "at random":

vv1 = randomv[];vv2=randomv[];orthonorm[{vv1,vv2}]

br1=%[[1]]; br2 = %[[2]];

pssp[br1,br2]

pssp[br1,br1]

pssp[br2,br2]

pscross[{stb[4],stb[5],br1,br2}]

ch[%]

In  the case ch[%]>= 0 choose other random pairs, until ch[%]= -1; then the random circle is interlaced with the unit circle, case 1 b); continue!

cc1 =plotcircle2spv[br1,br2];Show[{grc1,cc1},PlotRange -> All]

pp[stb[4],stb[5],br1,br2]

Eigenvalues[%]

The eigenvalues differ; the circles are not isogonal.

Again we calculate the Möbius invariant (the cosinus of the intersection angle) of two spheres each containing one of the  circles:

Clear[a,b];sph2[b_] :=br1*Cos[b] +br2*Sin[b]

fct2[a_,b_]:= pssp[sph1[a], sph2[b]]

fct2[a,b]

Plot3D[fct2[a,b],{a,0,2*Pi},{b,0,2*Pi}]

ptgr = %;

For fixed a (or b) we again have cos-functions; now the amplitude depends on a:

Clear[a,b];a=RandomReal[];Print[a];Plot[fct2[a,b],{b,0,2 π}]

Clear[a]

The amplitude of fct2[a,.] as function of b is

amplitude[b_] :=Sqrt[(-0.15100011924327345` Cos[b]-0.7210844733072147` Sin[b])^2 +(-0.4799738432849672` Cos[b]+0.6294473570124225` Sin[b])^2]

Plot[amplitude[b], {b,0,2*Pi}]

In the non isogonal case the intersection angle of two spheres containing the circles generates angle functions with varying amplitude.

4.5 Case 1c: Circles Intersecting in Exactly One Point

4.5.1. First Example

See also the random pair, case 1.c, in section 4.2. Now we take circles intersecting at the origin in orthogonal planes:

circlespacevectors[0,1,0,0,Pi/2, 1]

v1 = %[[1]]

v2 = %%[[2]]

plotcircle2spv[v1, v2, PlotStyle -> {Green,Thick}]

c0 = %;

radius[v1,v2]

center[v1,v2]

posvec[v1,v2]

circlespacevectors[0,0,1,0,0, 1]

w1 = %[[1]]

w2 = %%[[2]]

radius[w1,w2]

center[w1,w2]

posvec[w1,w2]

plotcircle2spv[w1,w2,PlotStyle -> {Red,Thick}]

c1 = %;

pscross[{v1,v2,w1,w2}]

The isotropic vector corresponding to the intersection must be orthogonal to v1,v2,w1,w2; thus, the intersection is

stepro[%]

ch[%%]

rank[{v1,v2,w1,w2}]

gpto=Graphics3D[{PointSize[.02],Point[{0,0,0}]}]

It follows: case 1.c; the circles are intersecting in one point,namely the origin,  and are in "general position", as show the graphics:

Show[{c0,c1,gpto},PlotRange->All]

The matrix of the double projection

pp[v1,v2,w1, w2]

has eigenvalues 1 and 0.  The corresponding eigenvectors

Eigenvectors[pp[v1,v2,w1, w2]]

give the "eigenspheres" (spheres corresponding to the eigenvectors) containing the circle C0:

eigenv1 =v1

eigenv2 = v2

Graphics of the two  eigenspheres are given by

{euklidsphereplot3D[eigenv1, -1,1,0,2],euklidsphereplot3D[eigenv2]}

Show[%, PlotRange -> {{-1,1},{0,2},{-1,1}}]

eigensph ={%, c0}

Show[%]

eigensph = %;

They intersect orthogonally in the circle C0. The spheres containing the circles C0, C1 are the one-parameter families S0[t]⊃C0, S1[s]⊃ C1 defined by the spacelike unit vectors

?mcirc`*

uv0[t_]:= eigenv1 Cos[t] +  eigenv2 Sin[t]; 
uv1[s_] := w1 Cos[s] + w2 Sin[s];

The spheres corresponding to any pair{uv0[t],uv1[s]} intersect since both contain the intersection o = C0⋂C1.  The scalar product yields the cosinus of the intersection angle phi

fctr[s_,t_] =pssp[uv0[t], uv1[s]]

phi[s_,t_] := ArcCos[fctr[s,t]]*180/Pi

The next procedure shows two spheres, corresponding to a random pair {uv0[tt],uv1[ss]}, prints their intersection angle phi (in degrees), plots the spheres their intersection (the blue circle), and the circles C0, C1

tt=Random[]*2*Pi; ss =Random[]*2*Pi; Print["phi = ", phi[ss,tt]];gtt=euklidsphereplot3D[uv0[tt],-1,1,-1,1,PlotStyle -> {Red,Opacity[.06]},Mesh-> None];gss=euklidsphereplot3D[uv1[ss],-1,1,-1,1,PlotStyle -> {Blue,Opacity[.06]},Mesh-> None];gst=plotcircle2spv[uv0[tt],uv1[ss], PlotStyle->{Thick, Blue}];
Show[{gss,gtt,gst,c0,c1, gpto},PlotRange -> All,Boxed -> False, Axes-> False]

Clear[ss,tt]

For the cosinus of the intersection angle (see the notebook emg.nb) we obtain

Plot3D[fctr[s,t],{s,0, 2 Pi},{t,0,2 Pi}]

The first eigenvector  of pp[v1,v2,w1,w2] corresponds to t=0 or Pi. The spheres uv1[s] imtersect the sphere eigenv1 with the   cosinus fctr[s, 0]. Its maximum is the square root of the eigenvlue, 1 in our case:

Plot[fctr[s, 0], {s, 0, 2 Pi}]

fctr[Pi/2, 0]

uv1[Pi/2]

euklidsphereplot3D[%]

vcenter[uv1[Pi/2]]

vradius[uv1[Pi/2]]

Here are six spheres through the circle C1:

Table[euklidsphereplot3D[uv1[s * Pi],0,1,-1,1, PlotPoints -> 20, PlotStyle-> Opacity[.4],Mesh->None], {s, 0., 0.5, 0.1}]

testgr1= %;

Show[{testgr1,c1},PlotRange -> All]

testgr1 = %;

euklidsphereplot3D[eigenv1, -1,7,-4,4, Mesh -> None]

ev1 =%;

Show[{testgr1,ev1,c0,gpto}, PlotRange -> All]

testgr2 = %;

This picture also visualizes the intersection angle of the spheres through C1 with the eigensphere to the eigenvalue 1, being the x,y-plane.

pssp[uv1[s],eigenv1]

euklidsphereplot3D[uv1[0], -1,7,-4,4, Mesh -> None]

Remark: The value s = 0 gives the y,z-plane, and s= Pi/2 yields the maximum value 1: the sphere is tangent to eigenv1, the x,y-plane:

{euklidsphereplot3D[eigenv1],euklidsphereplot3D[uv1[Pi/2]]}

Show[{%,c0,c1},PlotRange -> All]

Now we shall visualize the intersection angle with the second eigensphere, which  equals Pi/2 for all values of s:

fctr[s,Pi/2]

pssp[eigenv2,uv1[s]]

euklidsphereplot3D[eigenv2,-1,1,-1,1, PlotStyle->{Yellow, Opacity[.4]},Mesh-> None]

ev2 = %;

Show[{ev2,c0}, PlotRange-> All]

Show[{ev2,c0,euklidsphereplot3D[uv1[.3* Pi],0,1,-1,1, PlotPoints -> 20, PlotStyle-> Opacity[.4],Mesh-> None]},PlotRange -> All]

Show[Flatten[{ev2,c0,testgr1,gpto}],PlotRange -> All]

testgr3 = %;

Clear[testgr2,v1,v2,w1,w2,fctr,testgr1,testgr,testgr2,testgr3, grsph1,uv1,uv0]

4.5.2. Another Example

bsp =vspace[{1,0,0},{2,0.5,1.5},{2.8,0,0}]

v1 =stb[4]

v2 = stb[5]

w1 = bsp[[1]]

w2 = bsp[[2]]

{plotcircle2spv[v1, v2,PlotStyle-> {Green,Thick}],plotcircle2spv[w1,w2,PlotStyle-> {Red,Thick},PlotRange -> All]}

Show[%,PlotRange-> All]

expc2 = %;

pscross[{v1,v2,w1,w2}]

ch[%]

This corresponds to case 1.c. The unique intersection point {1,0,0} is represented by the isotropic vector pscross[{v1, v2, w1, w2}].

Chop[stepro[%%]]

ptc = Graphics3D[{PointSize[.03],Point[{1,0,0}]}];

Show[{expc2,ptc}]

expc2=%;

pp[v1,v2,w1, w2]

Eigenvalues[%]

Eigenvectors[pp[v1,v2,w1, w2]]

Therefore the "eigenspheres" (spheres corresponding to the eigenvectors) containing the circle C0 are

eigenv1 =normalize[-v1*0.532813269686138823` - v2*0.846232840089750837`]

eigenv2 =normalize[v2*0.532813269686138823` - v1*0.846232840089750837`]

{euklidsphereplot3D[eigenv1,-1,1,-1,1,PlotStyle->{Red,Opacity[.3]},Mesh -> None],euklidsphereplot3D[eigenv2,-1,1,-1,1,PlotStyle->{Blue,Opacity[.5]}, Mesh -> None]}

Show[%, PlotRange ->All]

eigensph = %;

Show[{eigensph,expc2}, PlotRange -> All]

They intersect orthogonally in the circle C0. The spheres containing the circles C0, C1 are the one-parameter families S0[t]⊃C0, S1[s]⊃ C1 defined by the spacelike unit vectors

uv0[t_]:= eigenv1 Cos[t] +  eigenv2 Sin[t]; 
uv1[s_] := w1 Cos[s] + w2 Sin[s];

The spheres uv0[s],uv1[t] always intersect.  For the cosinus of the intersection angle (see the notebook emg.nb) we calculate

pssp[uv0[t], uv1[s]]

fctr[s_,t_] = Simplify[Chop[%]]

Plot3D[fctr[s,t],{s,0, 2 Pi},{t,0,2 Pi}]

The first eigenvector  of pp[v1,v2,w1,w2] corresponds to t=0 or Pi, where the spheres uv1[s] have the followwing cosinus with the sphere eigenv1 :

Plot[fctr[s, 0], {s, 0, 2 Pi}]

The maximum 1 equals the square root of the first eigenvalue.  For t = Pi/2, 3Pi/2 we get correspondingly for the second eigenvector:

Plot[fctr[s, Pi/2], {s, 0, 2 Pi}]

The amplitudes correspond to the square roots of the corresponding eigenvalues: 1 and

Sqrt[0.0737033666969973921`]

Indeed, this is the amplitude:

ff[s_]:= fctr[s,Pi/2]

FindMaximum[ff[s],{s,0}]

Remove[bsb,eigenv1,eigenv2,eigensph,explc2,v1,v2,w1,w2,fctr,uv1,uv0,ff]

4.6 The Cases 2

As remarked in Subsection 4.1 above,  in this case the circles lie in a 2-sphere. Thus it would be easy  to describe their properties diminishing the dimension, setting dim = 4. We leave this as an exercise to the user and  treat here the circles as objects in (S^3).

4.6.1. Definition of vspace 2d

A general circle in the x,y-plane with center {x, y} and radius r has the vspace spanne by

circlespacevectors[x,y,0,0,Pi/2,r]

This formula leads to the definition (we consider the circles in the x,y-plane as special 3D-objects):

vspace2d[x_, y_, r_] := 
  {{(1 - r^2 + x^2 + y^2)/(2*r), x/r, y/r, 0, 
    -(1 + r^2 - x^2 - y^2)/(2*r)}, 
   stb[4]}

vspace2d[x,y,r]

pssp[%[[1]],%[[1]]]

Together[%]

(w1 =vspace2d[0,1,1][[1]]; w2 =vspace2d[0,1,1][[2]];
plotcircle2spv[w1,w2])

4.6.2. a) Two Disjoint Circles:

(v1 = stb[4]; v2 = stb[5] ;
	w1 =vspace2d[0,1/2, 3][[1]];w2 =vspace2d[0,1/2, 3][[2]];
C0xC1 = pscross[{v1,v2,w1,w2}])

Indeed, since dimension(U0 + U1) = 3 the vectors must be linearly dependent.

{w1,w2}

ma={v1,v2,w1}

Det[Table[pssp[ma[[i]], ma[[j]]],{i,3},{j,3}]]

It follows: case 2.a; the circles are disjoint and separated; the sum U0+U1 = span[ma] is pseudo-Euclidean.

{plotcircle2spv[v1,v2],plotcircle2spv[w1,w2]}

Show[%,PlotRange-> All]

pp[v1,v2,w1, w2]

Eigenvalues[%]

N[%]

Eigenvectors[pp[v1,v2,w1, w2]]

Therefore the "eigenspheres" (spheres corresponding to the eigenvectors) containing the circle C0 are v1,v2

{euklidsphereplot3D[v1],euklidsphereplot3D[v2]}

Show[%,PlotRange -> All]

They intersect orthogonally in the circle C0. The spheres containing the circles C0, C1 are the one-parameter families S0[t]⊃C0, S1[s]⊃ C1 defined by the spacelike unit vectors

uv0[t_]:=v1 Cos[t] +  v2 Sin[t]; 
uv1[s_] := w1 Cos[s] + w2 Sin[s];

For the generalized cosinus of the intersection angle (see the notebook emg.nb) we calculate

pssp[uv0[t], uv1[s]]

fctr[s_,t_] = %;

Plot3D[fctr[s,t],{s,0, 2 Pi},{t,0,2 Pi}]

The first eigenvector  of pp[v1,v2,w1,w2] corresponds to t=0 or Pi, where the spheres uv1[s] have the following generalized cosinus with the sphere v1:

Plot[fctr[s, 0], {s, 0, 2 Pi}]

We remark that v1 = stb[3]  = w2 spans U0⋂U1 and defines the "sphere" containing both circles: the x,y-plane.  uv1[0] defines the sphere intersecting the x,y-plane orthogonally in the circle C1.

euklidsphereplot3D[uv1[0]]

Of course s= Pi/2 yields the x,y-plane:

euklidsphereplot3D[uv1[Pi/2]]

For t = Pi/2, 3Pi/2 we get correspondingly for the second eigenvector v2 = stb[4] the unit sphere through o; :

Plot[fctr[s, Pi/2], {s, 0, 2 Pi}]

For s = Pi/2 - the x,y-plane uv1[Pi/2] and the unitsphere intersect orthogonally. At the points with Abs[fctr[s, Pi/2]] = 1 the sphere uv1[s] and the unit sphere are tangent to each other, and for Abs[fctr[s, Pi/2]] > 1 both spheres are disjoint. The amplitudes correspond to the square roots of the corresponding eigenvalues.

Sqrt[2.640625`]

4.6.3. b) Two Intersecting Circles:

(v1 = stb[4]; v2 = stb[5] ;
	w1 =vspace2d[0,1/2, 0.7][[1]];w2 =vspace2d[0,1/2, 0.7][[2]];
C0xC1 = pscross[{v1,v2,w1,w2}])

Indeed, since dimension(U0 + U1) = 3 the vectors must be linearly dependent.

{w1,w2}

ma={v1,v2,w1}

Det[Table[pssp[ma[[i]], ma[[j]]],{i,3},{j,3}]]

It follows: case 2.b; the circles intersect; the sum U0+U1 = span[ma] is Euclidean.

{plotcircle2spv[v1,v2],plotcircle2spv[w1,w2]}

Show[%, PlotRange -> All]

pp[v1,v2,w1, w2]

Eigenvectors[pp[v1,v2,w1, w2]]

Therefore the "eigenspheres" (spheres corresponding to the eigenvectors) containing the circle C0 are v1,v2

{euklidsphereplot3D[v1],euklidsphereplot3D[v2]}

Show[%]

They intersect orthogonally in the circle C0. The spheres containing the circles C0, C1 are the one-parameter families S0[t]⊃C0, S1[s]⊃ C1 defined by the spacelike unit vectors

uv0[t_]:=v1 Cos[t] +  v2 Sin[t]; 
uv1[s_] := w1 Cos[s] + w2 Sin[s];

For the generalized cosinus of the intersection angle (see the notebook emg.nb) we calculate

pssp[uv0[t], uv1[s]]

fctr[s_,t_] = %

Plot3D[fctr[s,t],{s,0, 2 Pi},{t,0,2 Pi}]

The first eigenvector  of pp[v1,v2,w1,w2] corresponds to t=0 or Pi, where the spheres uv1[s] have the following generalized cosinus with the sphere v1:

Plot[fctr[s, 0], {s, 0, 2 Pi}]

We remark that v1 = stb[3]  = w2 spans U0⋂U1 and defines the "sphere" containing both circles: the x,y-plane.  uv1[0] defines the sphere intersecting the x,y-plane orthogonally in the circle C1.

euklidsphereplot3D[uv1[0]]

Of course s= Pi/2 yields the x,y-plane:

euklidsphereplot3D[uv1[Pi/2]]

For t = Pi/2, 3Pi/2 we get correspondingly for the second eigenvector v2 = stb[4] the unit sphere through o :

ff[s_]:= fctr[s, Pi/2]

Plot[ff[s], {s, 0, 2 Pi}]

Of course, any two spheres S0[t], S1[s] intersect: the "generalized" cosinus is the usual cosinus of the intersection angle. The amplitude is the square root of the smaller eigenvalue:

Sqrt[0.7844897959183673`]

FindMaximum[ff[s],s]

4.6.4. c) Two Tangent Circles

(v1 = stb[4]; v2 = stb[5] ;
	w1 =vspace2d[0,1/2, 1/2][[1]];w2 =vspace2d[0,1/2, 1/2][[2]];
C0xC1 = pscross[{v1,v2,w1,w2}])

Indeed, since dimension(U0 + U1) = 3 the vectors must be linearly dependent.

{w1,w2}

ma={v1,v2,w1}

Det[Table[pssp[ma[[i]], ma[[j]]],{i,3},{j,3}]]

It follows: case 2.c; the circles are tangent to each other; the sum U0+U1 = span[ma] is isotropic: the vector v2 + w1 is orthogonal to U0+U1.

{plotcircle2spv[v1,v2],plotcircle2spv[w1,w2]}

Show[%,PlotRange -> All]

pp[v1,v2,w1, w2]

Eigenvectors[pp[v1,v2,w1, w2]]

Therefore the "eigenspheres" (spheres corresponding to the eigenvectors) containing the circle C0 are v1,v2, the x,y-plane and the standard unit sphere.

They intersect orthogonally in the circle C0. The spheres containing the circles C0, C1 are the one-parameter families S0[t]⊃C0, S1[s]⊃ C1 defined by the spacelike unit vectors

uv0[t_]:=v1 Cos[t] +  v2 Sin[t]; 
uv1[s_] := w1 Cos[s] + w2 Sin[s];

pssp[uv0[t], uv1[s]]

fctr[s_,t_] = Simplify[%]

Plot3D[fctr[s,t],{s,0, 2 Pi},{t,0,2 Pi}]

The first eigenvector  of pp[v1,v2,w1,w2] corresponds to t=0 or Pi, where the spheres uv1[s] have the following generalized cosinus with the sphere v1:

Plot[fctr[s, 0], {s, 0, 2 Pi}]

We remark that v1 = stb[3]  = w2 spans U0⋂U1 and defines the "sphere" containing both circles: the x,y-plane.  uv1[0] defines the sphere intersecting the x,y-plane orthogonally in the circle C1.

euklidsphereplot3D[uv1[0]]

Of course s= Pi/2 yields the x,y-plane:

euklidsphereplot3D[uv1[Pi/2]]

For t = Pi/2, 3Pi/2 we get correspondingly for the second eigenvector v2 = stb[4] the unit sphere through o:

Plot[fctr[s, Pi/2], {s, 0, 2 Pi}]

Of course, any two spheres S0[t], S1[s] intersect: the "generalized" cosinus is the usual cosinus of the intersection angle. Only for s = k Pi  we have S1[s] tangent to the unitsphere:

{euklidsphereplot3D[v2,-1,1,-1,1,PlotStyle -> Opacity[.2]], euklidsphereplot3D[uv1[0]]}

Show[%,PlotRange -> All]

testgr = %;

4.7 Further Examples and Applications

4.7.1. Standard Unit Circle Compared with a Circle Family with Increasing Radius

cv[t_] := circlespacevectors[1,0,0,Pi/2,0,t]

cv1[t_] =cv[t][[1]]

cv2[t_]= cv[t][[2]]

circlefamily[cv1,cv2,t,0.1,3]

Show[{grc1,%},PlotRange -> All]

cvpp[t_] := pp[stb[4],stb[5], cv1[t], cv2[t]]

cvpp[t]

Of course we have vor all t≠ 0 the case 1:

pscross[{stb[4],stb[5], cv1[t], cv2[t]}]

For all t, s  we have case 2:

pscross[{cv1[s], cv2[s], cv1[t], cv2[t]}]

Since the maximal eigenvalue is ((t^2)/4), we have for 0<t<2 circles interlaced with grc1, for t = 2 a circle intersecting grc1 in one point, and for t>2 two disjoint circles which may be separated by spheres.

4.7.2. Osculating Circles of a Plane Curve

Plane curves could be treated as special cases of space curves, see the next subsection. To use 2-dimensional graphics tools of Mathematica it is better, and more naural, to consider them as plane objects. Therefore we have to diminish the dimension and to adapt some modules to this aim:

dim = 4;

An isotropic representative of the plane curve is defined by

curve2D[cu_][t_] := invstepro[{cu[t][[1]],cu[t][[2]]}]

cc[t_]:={x[t],y[t]}

curve2D[cc][t]

Simplify[pssp[%,%]]

The k-th osculating space of a curve in the n-dimensional Möbius space is defined as the linear hull of the derivatives of the representing vectors from 0 to k; it is shown: if k ≧ 2, then it is  a pseudo-Euclidean subspace for regular curves, see [MDG] . Thus we define the i-th derivatve by

der[i_][cu_][t_] := D[curve2D[cu][t],{t,i}]

der[0][cc][t]

Simplify[der[1][cc][t]]

Simplify[pssp[%,%]]

Clear[cu]

The spacelike normal vector of the second order osculating space of the curve cu[t] at the point t  is the cross product

osccirc[cu_][t_]:=psCross[der[0][cu][t],der[1][cu][t],der[2][cu][t]]

osccirc[cc][t]

Simplify[Together[%]]

Simplify[pssp[%,%]]

For a circle the osculating circle at any of its points coincides with the circle itself:

fc[t_] := {Cos[t], Sin[t]}

Simplify[osccirc[fc][t]]

Befor plotting the circles we must adapt the circle2spv and the corresponding plotting command to the plane case:

circlespv[u_List][t_]:= Module[{uu,uv,vv,s},
uu = u/Sqrt[pssp[u,u]];
uv = Append[Append[Take[uu,3],0],uu[[4]]];
dim = 5;
vv[s_]=Take[ circle2spv[uv,{0,0,0,1,0}][s],2];
dim = 4;
vv[t]]

circlespv[stb[4]][t]

plotcirclespv[u_List,opts___]:= ParametricPlot[circlespv[u][t],{t,0,2*Pi},opts]

Here appears the standard unit circle

osccirc[fc][t]/. t-> RandomReal[]

plotcirclespv[%]

Now we consider an ellipse

Clear[a]; fc[t_]:= {Cos[t],a*Sin[t]}

a=.5;fc[t]

ParametricPlot[{Cos[t],a*Sin[t]},{t,0,2*Pi},PlotStyle ->{ Red}]

ell= %;

osccirc[fc][t]

Chop[%]

FullSimplify[Together[%],t>0]

Simplify[Chop[%],t>0&&t>Pi/2]

vder[t_] = %;

vder[t]

osc1 = vder[Pi/4]

plotcirclespv[osc1, PlotStyle-> Green]

oscc1 =%;

Here one recognizes that the osculating circle crosses the curve at the common tangential point, in general:

Show[{ell,oscc1},PlotRange -> All]

Run the next cells several times!

vder[RandomReal[]*2*Pi]

plotcirclespv[%, PlotStyle -> Green]

Show[{ell,%},PlotRange -> All]

osc2 = vder[0]

plotcirclespv[osc2, PlotStyle-> Green]

oscc2 =%;

At the extremal pints of the Euclidean curvature, the vertices, the osculating circles remains locally at one side of the curve:

Show[{ell,oscc2},PlotRange -> All]

osc3 = vder[Pi/2]

plotcirclespv[osc3, PlotStyle-> Green]

oscc3 = %;

Show[{ell,oscc3},PlotRange -> All]

Now we generate a family of osculating circles

Table[plotcirclespv[vder[N[j*2*Pi/10]], PlotStyle-> Green],{j,1,10}]

Show[{%,ell},PlotRange -> All,Axes -> False]

picellosc = %;

4.7.3. Osculating Circles and Osculating Spheres of a Space Curve

We define the isotropic representative of a space curve:

curve3D[cu_][t_] := invstepro[{cu[t][[1]],cu[t][[2]],cu[t][[3]]}]

cc[t_]:={x[t],y[t],z[t]}

cc[t]

curve3D[cc][t]

Simplify[pssp[%,%]]

stepro[%%]

The corresponding derivatives are

der3D[i_][cu_][t_] := D[curve3D[cu][t],{t,i}]

der3D[0][cc][t]

der3D[1][cc][t]

Simplify[Together[%]]

Simplify[pssp[%,%]]

Clear[cu]

oscsphere[cu_][t_]:=psCross[der3D[0][cu][t],der3D[1][cu][t],der3D[2][cu][t],der3D[3][cu][t]]

oscsphere[cc][t];

Simplify[Together[%]]

Simplify[pssp[%,%]]

fc[t_] := {Cos[t], Sin[t],t}

Simplify[oscsphere[fc][t]]

%/.t-> 3.

This is the spacelike  unit vector  defining the osculating sphere of the screw line at the point t=3:

osph1 = %/Sqrt[pssp[%,%]]

euklidsphereplot3D[osph1,-1,1,-1,1,PlotStyle->Opacity[.2], Mesh-> None]

osphere1=%;

ParametricPlot3D[fc[t],{t,0,2*Pi}, PlotStyle ->{Red, Thickness[0.01]}]

screw=%;

ptos= Graphics3D[{PointSize[.025], Point[fc[3.]]}];

Show[{osphere1,screw,ptos}, PlotRange -> All]

Fig1=%;

osccircle3D[cu_][t_] :=psCross[der3D[0][cu][t],der3D[1][cu][t],der3D[2][cu][t],oscsphere[cu][t]]

Simplify[osccircle3D[fc][t]]

%/.t->3.

osph2 = %/Sqrt[pssp[%,%]]

The spacelike unit vectors osph1,osph2 are orthogonal to the curve and ist first two derivatives at the point t=3; being orthonormal; they define the osculating circle of the screwline at this point

Chop[pssp[osph1,%]]

Here is the osculating circle

plotcircle2spv[osph2,osph1,PlotStyle -> {Green, Thick}]

oscirc1= %;

Show[{Fig1,oscirc1},PlotRange -> All]

Fig2 = %;

Fig2 shows the common position of the osculating circle and sphere of the screw line  at the point t03. The next picture shows the screwline together with the osculating circle:

Show[{screw,oscirc1},PlotRange -> All]

4.7.4. The same for an Elliptically Distorted Helix

Here is  the parameter representation of the distorted helix:

fc[t_]:= {a*Cos[t],Sin[t],t/2}

a=.5;fc[t]

ParametricPlot3D[fc[t],{t,-2,7},PlotStyle ->{ Red,AbsoluteThickness[4]}]

cva= %;

Simplify[oscsphere[fc][t]]

Simplify[Chop[%]]

FullSimplify[%,t>0&&t<Pi/2]

Simplify[Chop[%],t>0&&t<Pi/2]

vder3D[t_] = %;

vder3D[t]/.t-> 3.

vo1 = %/Sqrt[pssp[%,%]]

This is the spacelike unit vector defining the osculating sphere of the curve at the point t=3. Here is this sphere

euklidsphereplot3D[vo1,-1,1,-1,1,PlotStyle -> Opacity[.1], Mesh-> None]

vo1sphere = %;

osccircle3D[fc][t];

N[%/.t->3.]

vo2 = Chop[%/Sqrt[pssp[%,%]]]

This spacelike uniz vector together with the orthonormal vector vo1 defines the osculating circle of the curve at the point t=3:

Chop[Simplify[pssp[%,vo1]]]

plotcircle2spv[vo1,vo2, PlotStyle-> {Green, Thick}]

oscirc=%;

Show[{cva,vo1sphere,oscirc},PlotRange -> All]

Fig4 = %;

4.8 The Möbius Invariants Expressed by Euclidean Invariants

A complete system of Euclidean invariants of  a pair (C0,C1) of circles is obtained by the following entities:
The radii r of C0, R of C1,
the distance d of the  centers,
the angle a between the position vectors nc0, nc1,
the angle b between the position vector nc0, and the line connecting the centers,
the angle c between the position vector nc1, and the line connecting the centers.

4.8.1. The calculation of the expressions

We assume that the circle  C0 belongs to the  x,y-plane, has the center in the origin, and Radius r:

csv0 =circlespacevectors[0,0,0,0,Pi/2, r]

We choose the coordinate x-axis  in the Euclidean 3-space such that the center of the second circle belongs to the x,z-plane:

csv1= circlespacevectors[x,0,z,u,v,R]

The matrix of the corresponding double projection is

matrpp = pp[csv0[[1]],csv0[[2]],csv1[[1]],csv1[[2]]]

FullSimplify[matrpp]

Eigenvalues[%]

FullSimplify[%,r>0&&R>0 &&u∈ Reals&&v∈ Reals]

After 50 minutes no result. Therefore it seems better to evaluate the determinant and the trace of pp:

FullSimplify[ppdet[csv0[[1]],csv0[[2]],csv1[[1]],csv1[[2]]]]

invdet = (2*x*z*Cos[u]*Cos[v] + (r^2 + R^2 - x^2 + z^2)*Sin[v])^2/(4*r^2*R^2);

FullSimplify[pptr[csv0[[1]],csv0[[2]],csv1[[1]],csv1[[2]]]]

invtr = (1/(4*r^2*R^2))*(r^4 + R^4 - 2*R^2*x^2 + 2*r^2*(2*R^2 - x^2 + z^2) + 
    (x^2 + z^2)^2 + 2*R^2*((-(r^2 + z^2))*Cos[2*v] + 2*x*Cos[u]*Cos[v]*
       (x*Cos[u]*Cos[v] + 2*z*Sin[v])));

Clear[r,R,d,a,b,c]

Now we express the coordinates x,z,u,v by  the Euclidean invariants described above, and use these expression as substitution rules.  An elementary calculation yields

subrules = {x-> d*Sin[b], z->d*Cos[b], Sin[v]-> Cos[a], Cos[v] -> Sin[a],Cos[u]->(Cos[c]-Cos[a]*Cos[b])/(Sin[a]*Sin[b])};

FullSimplify[invdet//.subrules]

euforminvdet[r_, R_, d_, a_, b_, c_] := 
  ((-d^2 + r^2 + R^2)*Cos[a] + 2*d^2*Cos[b]*Cos[c])^2/(4*r^2*R^2)

FullSimplify[invtr//.subrules]

%//.Cos[2*v]->Cos[v]^2-Sin[v]^2

%//.subrules

FullSimplify[%]

euforminvtr[r_, R_, d_, a_, b_, c_] := 
  (d^4 + r^4 + 4*r^2*R^2 + R^4 + 2*r^2*R^2*Cos[2*a] + 
    2*d^2*(r^2*Cos[2*b] + R^2*Cos[2*c]))/(4*r^2*R^2)

The functions euforminvdet, euforminvtr are the wanted expressions; as usually one may find the eigenvalues of pp using these functions.

euforminvdet[r,R,d,a,b,c]

euforminvtr[r,R,d,a,b,c]

Finally, we consider  the case d=0, in which the angles b and c are not defined. In this case we have x = z = 0, and

invdet

%//.{x->0, z->0}

Replacing Sin[v] by Cos[a], as before, we see that the function euforminvdet is valid in this case, too, since the terms containing b or c vanish. Analogeously:

invtr

%//.{x->0, z->0}

Replacing Cos[2v]= Cos[v]^2 -Sin[v]^2 by Sin[a]^2 -Cos[a]^2 = -Cos[2a] we conclude again the validity of the expression euforminvtr.

Simplify[(%/. Cos[2*v] -> -Cos[2*a])-euforminvtr[r,R,0,a,b,c]]

4.8.2. Test

We take two circles with in part random parameters, find their circlespacevectors, and define their invariants ppdet, pptr. Then we find their Euclidean invariants and calculate the corresponding values of euforminvdet, euforminvtr, and compare their values with those of ppdet, pptr; furthermore, we compare the values of the given random radii r00, r11 with the radii r0,r1 calculated by the functions radius of the generated vectors. All the differences are numerically small, see the last cells of this subsection. 
Also we calculate the eigenvalues of pp (see the last cell); use Graphics/Rendering/Animate  to study the mutual position of the two random circles.

The random circles:

r00=2 RandomReal[];vc0=circlespacevectors[RandomReal[],RandomReal[],RandomReal[],2 π RandomReal[],(1/2)RandomReal[] π,r00]

r11=2 RandomReal[];vc1=circlespacevectors[RandomReal[],0,0,2 π RandomReal[],(1/2)RandomReal[] π,r11]

plotcircle2spv[vc0[[1]],vc0[[2]],PlotStyle->{Red,Thick}]

gc0=%;

plotcircle2spv[vc1[[1]],vc1[[2]],PlotStyle->{Green,Thick}]

gc1 =%;

Show[{gc0,gc1},PlotRange -> All]

shgr= %;

adet = ppdet[vc0[[1]],vc0[[2]],vc1[[1]],vc1[[2]]]

atr = pptr[vc0[[1]],vc0[[2]],vc1[[1]],vc1[[2]]]

The Euclidean equivariants of the circle C0:

r0=radius[vc0[[1]],vc0[[2]]]

pos0=posvec[vc0[[1]],vc0[[2]]]

cent0 =center[vc0[[1]],vc0[[2]]]

The Euclidean equivariants of the circle C1:

r1=radius[vc1[[1]],vc1[[2]]]

pos1=posvec[vc1[[1]],vc1[[2]]]

cent1 =center[vc1[[1]],vc1[[2]]]

The distance of the centers

vcentdiff01 = cent1 -cent0

d01 = Sqrt[vcentdiff01.vcentdiff01]

The angle between the position vectors

a01=ArcCos[pos0.pos1]

The angles between the positionvectors and the line connecting the centers:

b01 =ArcCos[pos0.vcentdiff01/d01]

c01 =ArcCos[pos1.vcentdiff01/d01]

euforminvdet[r0,r1,d01,a01,b01,c01]

adet-%

euforminvtr[r0,r1,d01,a01,b01,c01]

atr-%

{Abs[r00]-Abs[r0],Abs[r11]-Abs[r1]}

Eigenvalues[pp[vc0[[1]],vc0[[2]],vc1[[1]],vc1[[2]]]]

%[[1]] + %[[2]]

%%[[1]] *%%[[2]]

Clear[r0,r00,r1,r11,d01,a01,b01,c01,vc0,vc1,vcentdiff01]

5. Normal Forms of Pairs of Circles

In the general classification of pairs {k-sphere, l-sphere} contained in an n-sphere  ( see [OS3] or  [S00]  ) we obtained normal forms for the bases of the generating vector spaces, where the coefficients of one base with respect to the other are defined by the invariant  eigenvalues of the matrix pp. In this section we describe these normal forms for each of the cases mentioned in the subsection  4.1 Case discussion, above.  The advantage of the normal form representation is that the relations between the invariants and the geometry of the circle pair can be studied directly. For the first circle we always take the unit circle grc1 = plotcircle2spv [stb[4],stb[5]] with center o in the x,y-plane, defined with the initialization; the other circle is given by the corrresponding normal form of the base. We remark that the case numbering and the notations differ from those used in [S00]. 
ATTENTION: For getting the normal forms the definitions in subsection 5.1 are needed. Evaluate 5.1 first!

Show [grc1]

5.1 Some Vector Functions Needed for the Definition of the Normal Forms

Clear[a,b]

v3[a_] := {0,0,Sin[a],0,Cos[a]}

v4[a_] := {Sinh[a],0,0,Cosh[a],0}

v5[a_] := {0,Sin[a],0,Cos[a],0}

5.2 Case 1a: Disjoint, Non-Interlaced Circles

This case corresponds to case 1, n =5, m = l =1, in [S00].

pscross[{stb[4], stb[5],v4[a],v3[b]}]

pssp[v4[a],v3[b]]

For b = kπ , or a =0 we have case 2; for other values of a, b this vector is spacelike; therefore U0 + U1 is pseudo-Euclidean:case 1.a.

MatrixForm[Evaluate[pp[stb[4], stb[5], v4[a], v3[b]]]]

The corresponding circle is plotted by

circ1a[a_, b_,opts___] := plotcircle2spv[v4[a], v3[b],opts]

circ1a[1.5,1.7,PlotStyle ->{Thick, Green}]

For a = b = 0 we get ste standard circle grc1 again.

circ1a[0,0]

Show[{grc1,%%},PlotRange -> All]

We calculate center and radius of the circles circ1a[a,b]

center[v4[a], v3[b]]

Simplify[%,{ a∈Reals,b∈Reals}]

centfct[a_,b_] := {0,-Sin[2 b]/(Cos[2 b]+Cosh[2 a]),Sinh[2 a]/(Cos[2 b]+Cosh[2 a])}

All the centers lie in the y,z-plane.

radius[v4[a], v3[b]]

Simplify[%,{ a∈Reals,b∈Reals}]

radfct[a_,b_] := ((Sqrt[2]))/(√(Cos[2 b]+Cosh[2 a]))

Obviously, this function has poles for a = 0, b = Pi/2 + k*Pi, k integer

Plot3D[radfct[a,b],{a, -1, 1},{b, -Pi, Pi},PlotPoints ->60, PlotRange ->{-1,5}]

pp[v4[a1],v3[b1],v4[a2],v3[b2]]

Simplify[%]

pscross[{v4[a1],v3[b1],v4[a2],v3[b2]}]

Simplify[%]

It follows that any pair of circles of the 2-parameter family under consideration fulfil the conditions of case 1 a) - i.e. are in general position and may be separated by spheres - wth the exception of the values a1 = a2 or b1 = b2 + k*Pi (k integer),  where the circles are not in general position (case 2, or they are equal).

a1 = RandomReal[];b1 = RandomReal[];a2 = RandomReal[];b2 = RandomReal[];

cc1 =circ1a[a1,b1,PlotStyle -> {Thick,Red}]

cc2 =circ1a[a2,b2,PlotStyle -> {Thick,Green}]

Show[{cc1,cc2},PlotRange -> All]

Example of a circlefamily and a circle surface

vv[t_] := v3[t+1]

pssp[v3[t+1],v4[t]]

circlefamily[v4,vv,t,-3,3]

Show[%,PlotRange -> All]

cfam1a = %;

Centers and radii of this circle family are:

cf1aradius[t_] := radius[v4[t], vv[t]]

cf1acenter[t_] := center[v4[t], vv[t]]

cflaposvec[t_] := posvec[v4[t], vv[t]]

Simplify[cf1aradius[t]]

Plot[cf1aradius[t],{t,-8,8},PlotRange -> All]

One may localize the maximum of the radius graphically:

Plot[cf1aradius[t],{t,.26,.28}]

FindMaximum[cf1aradius[t],t]

Limit[cf1aradius[t],t-> Infinity]

Limit[cf1aradius[t],t-> -Infinity]

Simplify[cf1acenter[t]]

The centers are lying within the y, z - plane:

ParametricPlot[Take[cf1acenter[t],-2], {t, -8,8}, PlotPoints -> 100,PlotRange -> All]

Limit[cf1acenter[t], t-> Infinity]

Limit[cf1acenter[t], t-> -Infinity]

cflaposvec[t]

Limit[%,t->-Infinity]

All directiones parallel to the y,z-plane may appear as position vectors of our family:

ParametricPlot[Take[cflaposvec[t],-2],{t,-10,10}]

We compare the series of circles with grc1

Table[Show[{grc1,circ1a[a, a+1]}],{a, 0.7,3.2, 0.5}]

Show[%, PlotRange ->All,Axes->False, Boxed->False]

5.3 Case 1b: Interlaced Circles

5.3.1. Introduction

This case corresponds to case 3, n =5, m = l =1, in [6]. Run first 5.1 !

Clear[a,b]

pssp[v5[a],v3[b]]

pscross[{stb[4], stb[5],v5[a],v3[b]}]

pssp[%,%]

For b = kπ , or a = kπ  we have case 2; for other values of a, b this vector is timelike; therefore U0 + U1 is Euclidean:case 1.b.  For two arbitrary circle of the family the same is valid; the exceptions are  b-b1 = kπ , or a-a1 = kπ:

Simplify[pscross[{v5[a],v3[b],v5[a1],v3[b1]}]]

MatrixForm[Evaluate[pp[stb[4], stb[5], v5[a], v3[b]]]]

circ1b[a_, b_,opts___] := plotcircle2spv[v5[a],v3[b],opts]

For a = b = 0 we get the standard circle grc1 again. For a = b > 0 the pairing with grc1 yields the family of isogonal circles.

circ1b[0,0,PlotStyle -> {Thick, Red}]

circ1b[Pi/4,Pi/2,PlotStyle -> {Thick, Green}]

Show[{%,%%}]

5.3.2. Radius, center, position vector

Center and radius do not depend on a, all the centers lie on the y-axis:

center[v5[a], v3[b]]

radius[v5[a], v3[b]]

Obviously, this function has poles for  b = Pi/2 + k*Pi, k integer

Plot[-Sec[b],{b, -2Pi, 2Pi}]

It follows: All circles circ1b[a, (Subscript[b, 0])], (Subscript[b, 0])= constant, lie on the sphere with center {0,-Tan[ (Subscript[b, 0])],0} and radius Abs[Sec[ (Subscript[b, 0])]]. This "sphere" is a plane in case  (Subscript[b, 0])= Pi/2 + k*Pi:

v3[b]

euklidsphereplot3D[v3[Pi/2]]

The position vector does not depend on b:

posvec[v5[a], v3[b]]

Clear[a1,a2,b1,b2];pp[v5[a1],v3[b1],v5[a2],v3[b2]]

Simplify[%]

pscross[{v5[a1],v3[b1],v5[a2],v3[b2]}]

Simplify[%]

It follows that any pair of circles of the 2-parameter family under consideration fulfil the conditions of case 1 b) - i.e. are in general position and interlaced - wth the exception of the values a1 = a2+k*Pi or b1 = b2 + k*Pi (k integer),  where the circles are not in general position (case 2, or they are equal).

5.3.3. The Möbius  invariant

The Möbius invariant (generalized cosinus) of two spheres containing the standard unit circle resp. the circle circ[1,b] is

Clear[a,b,s,t]

sphinv[a_,b_][s_,t_] := pssp[stb[4]*Cos[t] +stb[5]*Sin[t],v5[a]*Cos[s]+v3[b]*Sin[s]]

sphinv[a,b][s,t]

{Cos[.3], Cos[1.5]}

Plot3D[sphinv[0.3,1.5][s,t], {t,0,2*Pi},{s,0,2*Pi}, PlotPoints -> 50]

graph1 = %;

FindMaximum[sphinv[0.3,1.5][s,t],{s,0},{t,0}]

sphinv[0.3,1.5][0,0]

FindMinimum[sphinv[0.3,1.5][s,t],{s,1.5},{t,1.5}]

sphinv[0.3,1.5][Pi/2,Pi/2]

The relative extrema of the function sphinv[a,b] are in the points (s,t) = (0,0), (Pi/2, Pi/2); their values coincide with the roots of the eigenvalues, as expected.

pp[stb[4],stb[5],v5[.3],v3[1.5]]

Sqrt[%]

5.2.4. A circle family and a circle suface

A family of to grc1  non-isogonal, with grc1 interlaced circles:

Table[circ1b[a, 2*a],{a, 1,2, 0.25}]

Show[{% ,grc1}, PlotRange ->All]

nf1b =%;

A corresponding surface:

ParametricPlot3D[Evaluate[circle2spv[v5[a], v3[1.2*a]][t]],{a, .1, .9}, {t, 0, 2 Pi},PlotPoints -> 25,PlotStyle-> Opacity[.3]]

ngr1b = %;

Show[{grc1, ngr1b},PlotRange -> All]

ngr1b1 = %;

5.3.5. The circle family and the circle suface with b = 1.2 (case 1b)

In the following circlefamily any pair of circles belong to case 2; the circles intersectin two points and belong to a fixed sphere (case 2b):

pp[v5[a],v3[1.2],v5[a1],v3[1.2]]

psCross[v5[a],v3[1.2],v5[a1],v3[1.2]]

Chop[%]

ParametricPlot3D[Evaluate[circle2spv[v5[a], v3[1.2]][t]],{a, .1, .9}, {t, 0, 2 Pi},PlotPoints -> 25]

sf1b = %;

Show[{%, grc1},PlotRange -> All]

sf1b1 = %;

vv[a_] := v3[1.2]

circlefamily[v5,vv,a,.1,.9,10,PlotStyle -> Red]

Show[%,PlotRange -> All]

cirf1b =%;

The sphere containing all these circles is

euklidsphereplot3D[v3[1.2],-1,1,-1,1,PlotStyle -> Opacity[.2], Mesh -> None]

sph12 =%;

Show[{sph12, cirf1b},PlotRange -> All]

cirfsph1b = %;

Show[{cirfsph1b, grc1},PlotRange -> All]

In the sphere  defined by v3[1.2] our family coincides with the bundle of circles through two points which we calculate now:

Solve[{pssp[vec[x], v5[a]]==0, pssp[vec[x],v3[1.2]]== 0, pssp[vec[x],v5[1.]]== 0,pssp[vec[x],vec[x]]== 0}, vec[x]]

vec[x]/.%

pts12 =%

pssp[pts12[[1]],pts12[[1]]]

Chop[%]

pssp[pts12[[2]],pts12[[2]]]

Chop[%]

In the Euclidean space these points have the coordinates

eupts = Chop[{stepro[pts12[[1]]], stepro[pts12[[2]]]}]

Map[Point,eupts]

Show[Graphics3D[{PointSize[0.03],%}]]

grpts12 = %;

Show[{cirfsph1b, grpts12},PlotRange -> All]

togb12 = %;

5.3.6. The circle family and the circle suface with a = 1 (case 2b)

In the following circlefamily any pair of circles belong to case 2; the circles intersectin two points and belong to a fixed sphere (case 2b):

pp[v5[1.],v3[b],v5[1.],v3[b1]]

Simplify[%]

psCross[v5[1.],v3[b],v5[1.],v3[b1]]

Chop[%]

ParametricPlot3D[Evaluate[circle2spv[v5[1.], v3[b]][t]],{b, .1, .9}, {t, 0, 2 Pi},PlotPoints -> 25,PlotStyle -> Opacity[.4]]

cirsf1b2 = %;

Show[{%,grc1},PlotRange -> All]

cirsfb2c1=%

For any a v5[a] defines a plane containing the y-axis: the family {v5[a]} coincides with the bundle of planes through the y-axis:

v5[a]

Table[euklidsphereplot3D[v5[a*Pi]], {a,0,1,.1}]

Show[%]

Especially we have

euklidsphereplot3D[v5[1.],-4,2,-3,1,PlotStyle -> Opacity[.1]]

pl1b =%;

Show[{pl1b,cirsf1b2},PlotRange -> All]

Rotate the figure and recognize:  the circle surface is a subset of the plane pl1b!

Show[{pl1b,cirsf1b2,grc1},PlotRange -> All]

For the Euclidean invarints we have

radius[v5[1.], v3[b]]

center[v5[1.], v3[b]]

The position vector is constant:

posvec[v5[1.], v3[b]]

In the plane pl1b our circle family coincides with the bundle of circles through two points which we calculate now:

Solve[{pssp[vec[x], v5[1.]]==0, pssp[vec[x],v3[1.]]== 0, pssp[vec[x],v3[a]]== 0,pssp[vec[x],vec[x]]== 0}, vec[x]]

sol1 = vec[x]/.%[[1]]

sol2 = vec[x]/.%%[[2]]

pssp[sol1,sol1]

{stepro[sol1],stepro[sol2]}

Chop[%]

pts = %;

Graphics3D[Line

Graphics3D[{Thickness[.007],Line[{pts[[1]],pts[[2]]}]}]

l1b = %;

Show[Graphics3D[{PointSize[0.03],Point[{0.5403023058681398`,0,-0.8414709848078965`}],Point[{-0.5403023058681398`,0,0.8414709848078965`}]}]]

pts1b = %;

Show[{pts1b,l1b,cirsf1b2,pl1b,grc1},Boxed -> False]

tog1b=%;

5.3.7. The isogonal case

Setting a = b we get the family of isogonal circle pairs, we plot a part of this family:

?circlefamily

circlefamily[v5,v3,a,0.,1.,10]

Show[%,PlotRange -> All]

grison=%;

Examples of isogonal circle pairs have been considered already in subsection 4.4.

A corresponding surface:

ParametricPlot3D[Evaluate[circle2spv[v5[a], v3[a]][t]],{a, 0, 1.}, {t, 0, 2 Pi},PlotPoints ->40,PlotStyle -> Opacity[.3]]

grisosf =%;

Show[grisosf, Boxed-> False,Axes-> False]

Export["grisosf.eps",%]

Radius, center, position vector

All the centers lie on the y-axis:

center[v5[a], v3[a]]

radius[v5[a], v3[a]]

Remember

Plot[%, {a,-Pi,Pi}]

Thus the radius goes to infinity very fast for a-> Pi/2; for a = Pi/2 we obtain the z-axis:

circ1b[Pi/2,Pi/2]

posvec[v5[a], v3[a]]

pp[v5[a1],v3[a1],v5[a2],v3[a2]]

Simplify[%]

pscross[{v5[a1],v3[a1],v5[a2],v3[a2]}]

Simplify[%]

It follows that any pair of circles of the 1-parameter family circ1b[a,a]  fulfil the conditions of case 1 b) - i.e. are in general position and isogonal - wth the exception of the values a1 = a2+k*Pi  (k integer),  where the circles coincide (up to orientation).

For the spherical invariant function one has the simple form, as mentioned already in subsection 4.4, in a special case.

sphinv[a,a][s,t]

Simplify[%]

Plot3D[sphinv[0.6,0.6][s,t], {t,0,2*Pi},{s,0,2*Pi}, PlotPoints -> 50]

graph2 = %;

Export["graph2.eps", graph2]

5.4 Case 2a: Non Intersecting Circles on a Sphere

b4 = stb[5]

The eigenvalues are obtained from the normal form

pp[stb[4], stb[5], v4[a], stb[5]]

Both circles lie on the standard unit sphere. U0 +U1 is spanned by stb[4], stb[5], v4[a]. We set

mat[a_] := {stb[4], stb[5],v4[a]}

and calculate

Det[Table[pssp[mat[a][[i]],mat[a][[j]]], {i,3}, {j,3}]]

It follows that  for a ≠ 0 the subspace U0 + U1 is pseudo-Euclidean: case 2.a:

Table[plotcircle2spv[v4[a],stb[5],PlotStyle -> {Thick, Red}],{a, 0.2,3.2, 0.5}]

Show[%,PlotRange -> All]

nf2a= %;

ParametricPlot3D[Evaluate[circle2spv[v4[a], stb[5]][t]],{a, 0, Pi}, {t, 0, 2 Pi},PlotStyle->{Opacity[.5],Blue},Mesh->None]

Show[{%,nf2a}]

5.5 Case 2b: 2-Point Intersecting  Circles in the x,y-Plane

Two such circle families have already be considered in subsection 5.3. This case one gets also specialising case 1a to a = 0. Then the eigenvalues are 1, cos[b,]^2:

pp[stb[4], stb[5], stb[4], v3[b]]

mat[b_] := {stb[4], stb[5],v3[b]}

and calculate

Det[Table[pssp[mat[b][[i]],mat[b][[j]]], {i,3}, {j,3}]]

It follows that  for b ≠ 0 the subspace U0 + U1 is Euclidean: case 2.b

vb[b_]:=stb[4]

circlefamily[vb,v3,b,0.1,1.,10,PlotStyle ->{ Red}]

Show[%, PlotRange -> All]

bdl4 = %;+

Show[{bdl4, grc1}]

nf2a= %;

All the circles belong to the x,y-plane and intersect the standard unit circle grc1. The common intersection points are the unit points x = ±1, y = 0. The normal forms yield the bundle of all circles in the x,y-plane through these points.

Solve[{pssp[mat[b][[1]],vec[ x]]==0, pssp[mat[b][[2]],vec[ x]]==0,
		pssp[mat[b][[3]],vec[ x]]==0 ,pssp[vec[x],vec[ x]]==0},vec[x]]

vec[x]/.%

The solutions are, as stated:

{stepro[{-1,1,0,0,0}], stepro[{1, 1,0, 0, 0}]}

Obviously, we have case 2b), and the argument b is the intersection angle of the two circles.  For b = k π  we obtain the trivial pair {grc1,grc1}.

plotcircle2spv[stb[4],v3[Pi]]

For b = Pi/2 we get the x-axis:

plotcircle2spv[stb[4],v3[Pi/2]]

Show[%, grc1]

5.6 Cases 1c and 2c: Circles Intersecting in One Point

For the normal form in this case we use v5[b] and the vector

b4 := {1,0,1,0,-1}

Therefore all circles of the considered family  lie on the sphere

euklidsphereplot3D[b4,-1,1,-1,1,PlotStyle -> Opacity[.2], Mesh-> None]

sphc = %;

MatrixForm[Evaluate[pp[stb[4], stb[5], v5[b], b4]]]

pscross[{stb[4], stb[5],v5[b], b4}]

For all values of b ≠ kπ  this is an isotropic vector; therefore we have case 1.c: all the corresponding circles intersect the unit circle in a uniquely defined point.

circ1c[ b_] := plotcircle2spv[v5[b], b4 , PlotStyle -> Green]

Table[ plotcircle2spv[v5[b*Pi], b4, PlotStyle -> {Thick,Green}], {b, 0.1, 1, 0.3}]

Show[%,PlotRange -> All]

sertest1 =%;

Show[{sertest1,grc1},PlotRange -> All]

test1c=%;

matc[b_] := {stb[4], stb[5],v5[b], b4}

Solve[{pssp[matc[b][[1]],vec[x] ]== 0,pssp[matc[b][[2]],vec[x] ]== 0,pssp[matc[b][[3]],vec[x] ]== 0,pssp[matc[b][[4]],vec[x] ]== 0}, vec[x]]

vec[x]/.%

All the circles intersect the standard unit circle in the uniquely defined point

stepro[{1,0, 1,0,0}]

All the circles of the family defined by {b4, v5[b]} intersect in two points:

Solve[{pssp[matc[b][[3]],vec[x] ]== 0,pssp[matc[b][[4]],vec[x] ]== 0,pssp[matc[0][[3]],vec[x]]== 0, pssp[vec[x], vec[x]]== 0}, vec[x]]

vec[x]/.%

{stepro[{-1,0,0,0,1}], stepro[{1,0,1,0,0}]}

Thus we have a family of circles through the point {0,1,0}. Especially we get

plotcircle2spv[v5[Pi], b4]

Show[{grc1, %}]

This is a tangent circle to grc1 in the point {0,1,0}: case 2.c. The matrix pp has the diagonal form

MatrixForm[Evaluate[pp[stb[4], stb[5], v5[Pi], b4]]]

pscross[{stb[4], stb[5],v5[Pi], b4}]

If b4 is replaced by the vector {1,0,1,0,1} one gets  a family of lines lying in the plane y = 1 through the point {0,1,0}.

euklidsphereplot3D[{1,0,1,0,1},-1,1,-1,1,PlotStyle -> Opacity[.4], Mesh -> None]

ge1=%;

v5[t]

radius[v5[t], {1,0,1,0,1}]

center[v5[t], {1,0,1,0,1}]

Table[plotcircle2spv[v5[s], {1,0,1,0,1}, PlotStyle ->{ Thick, Green}],{s,0,5,0.5}]

Show[{Out[140],ge1}]

Especially we get the tangent:

plotcircle2spv[v5[Pi], {1,0,1,0,1}]

Show[{grc1, %}]

Of course, we have the same invariants:

MatrixForm[Evaluate[pp[stb[4], stb[5], v5[Pi], {1,0,1,0,1}]]]

pscross[{stb[4], stb[5],v5[Pi], {1,0,1,0,1}}]

Clear[x,y,z,u,v,r,ma]

6.  Circle Surfaces, Circle Orbits, Geodesics in the Circle Space

A circle surface is a surface in the 3-space through any point of which exists a circle contained in the surface. With other words: A circle surface can be described  as a curve in the circle space. The spheres (and planes) are the simplest circle surfaces. The tubes considered in section 2.4 are circle surfaces. More general, envelopes of one-parameter sphere families: the channel surfaces, are circle surfaces.  A special case are the Dupin cyclides, which are channel surfaces in a twofold way. That means, that there exist two different sphere curves whose envelopes are the same channel surface: tori, circular coni,  circular cylinders, and their images under Möbius transformations. Other special circle surfaces are the orbits of  1-parameter subgroups of the Möbius group in the circle space.

The circle space is pseudo-Riemannian symmetric of rank two. The geodesics of the circle space are orbits of certain 1-parameter subgroups of the Möbius group.   As in all reductive spaces, the geodesics through an element C are the orbits of  one-parameter subgroups exp(At) C, t ∈ R, where in our case C is the standard unit circle grc1 and A an element of the orthogonal complement of the isotropy algebra of C in the Lie algebra o(5,1); here orthogonality is defined by the Killing form  of the Lie algebra o(5,1), see my notebook liealgeb.nb, in MathSource [MS].  It would be interesting to investigate the relation between the pseudo-Riemannian geometry of the circle space and the circle manifolds more basically. But this may be the subject of another notebook. Some results about the Möbius geometry of circle surfaces are contained in the diploma paper  [UM] of Uwe May; especially the geodesics defined in the next subsections can be found there. We emphasize that we do not give here a classification of all the geodesics with respect to the Möbius group. We only treat three simple cases of geodesics lying in a 2-sphere.

6.1 Spacelike Geodesics

The following matrix function defines the 1-parameter subgroup exp[Uh*t], where Uh is an element of the Lie algebra o(n,1) orthogonal to the isotropy algebra of grc1 and spacelike with respect to the Killing form of o(n,1) (see  the notebook liealgeb.nb mentioned above):

Uh = {{0,0,0,0,1},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{1,0,0,0,0}}

Uh//MatrixForm

gUh[t_]  := MatrixExp[Uh t]

MatrixForm[gUh[t]]

pssp[gUh[t].vec[x],gUh[t].vec[y]]

Simplify[%]

vh1[t_]:= gUh[t].stb[5]

vh1[t]

Obviously, stb[4] is fixed under gUh[t]. Therefore the orbit of {stb[4], stb[5]} is given by

hypgeodesic[t_,u_]:=circle2spv[vh1[t],stb[4]][u]

hypgeodesic[t,u]

ParametricPlot3D[hypgeodesic[t,u],{t,-2,2},{u,0,2*Pi}]

hypgeo = %;

This is the family of concentric circles in the x,y-plane containing the standard unit circle grc1 (for t = 0).

The group is the group of all dilatations in the Euclidean picture, leaving the origin  fixed.

invstepro[{0,0,0}]

stepro[gUh[t].invstepro[{0,0,0}]]

The invariants along the hyperbolic geodesic

invpp[s_,t_]:= pp[vh1[t],stb[4],vh1[s],stb[4]]

invpp[s,t]

Simplify[%]

The Invariants are the elements of the diagonal; the first eigenvalue is (for s ≠t) always greater than 1, we show this:

ftt[s_,t_] := (1/4) (E^(-2 (s+t))) ((((E^(2 s))+(E^(2 t))))^2)

Plot3D[ftt[s,t],{s,-1,1}, {t,-1,1}]

FindMinimum[ftt[s,t],{s,t}]

ftt1[s_,t_]= D[ftt[s,t],s]

ftt2[s_,t_]= D[ftt[s,t],t]

Simplify[ftt1[s,t]]

Simplify[ftt2[s,t]]

The derivatives are zero for s = t only; obviously, for any fixed t, ftt[s,t] reaches its minimum for s = t, and this minimum is

ftt[s,s]

The radius of the circle as function of t is

radius[vh1[t],stb[4]]

Simplify[%]

For the varying invariant expressed as a function of the radius we obtain

invpp[t, 0][[1,1]]/.(E^(2*t))-> r^2

%/.(E^(-2*t))-> r^(-2)

Plot[((((1/(2 r))+(r/2)))^2),{r, 0.2,5}]

Equivalent description on the 2-sphere

U3h = {{0,0,1,0,0},{0,0,0,0,0},{1,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0}}

MatrixForm[U3h]

gU3h[t_]  := MatrixExp[U3h t]

MatrixForm[gU3h[t]]

Simplify[pssp[gU3h[t].vec[x],gU3h[t].vec[y]]]

v3h1[t_]:=gU3h[t].stb[3]

v3h1[t]

hyp3geodesic[t_,u_]:=circle2spv[v3h1[t],stb[5]][u]

ParametricPlot3D[hyp3geodesic[t,u],{t,0.,1},{u,0,2 π}]

hyp3geo =%;

?circlefamily

xy[t_]:= stb[5]

circlefamily[v3h1,xy,t,0, 1.,10, PlotPoints -> 100]

Show[%,PlotRange -> All]

hgeofam = %;

6.2 Timelike Geodesics

The following matrix function defines the 1-parameter subgroup exp[Ue*t], where Ue is an element of the Lie algebra o(n,1) orthogonal to the isotropy algebra of grc1 and timelike with respect to the Killing form of o(n,1),see [MS]:

Ue = {{0,0,0,0,0},{0,0,0,1,0},{0,0,0,0,0},{0,-1,0,0,0},{0,0,0,0,0}}

Ue//MatrixForm

gUe[t_]:=MatrixExp[Ue t]

gUe[t]//MatrixForm

ve1[t_] := gUe[t].stb[4]

ve1[t]

ve2[t_] := stb[5]

circlefamily[ve1,ve2,t,-N[Pi],N[Pi], 30]

Show[%,PlotRange -> All]

ellgeo = %;

The standard circle rotates around the y-axis. The surface creation  yields the unit sphere:

ellgeodesic[t_,u_]:=circle2spv[ve1[t],stb[5]][u]

ellgeodesic[t,u]

ParametricPlot3D[Evaluate[ellgeodesic[t,u]],{t,-Pi/2,Pi/2},{u,0,2 Pi}, PlotPoints -> 30,PlotStyle-> Opacity[.3]]

The invariants along the elliptic geodesic

invellpp[s_,t_] := pp[ve1[s],stb[5],ve1[t],stb[5]]

Simplify[invellpp[s,t]]

As expected, the only essential invariant is the cosinus square of the intersection angle.

6.3 Isotropic Geodesics

The following matrix function defines the 1-parameter subgroup exp[U*t], where U is an element of the Lie algebra o(n,1) orthogonal to the isotropy algebra of grc1 and isotropic with respect to the Killing form of o(n,1) (see [7], liealgeb.nb ):

U ={{0,0,0,1,0},{0,0,0,0,0},{0,0,0,1,0},{1,0,-1,0,0},{0,0,0,0,0}}

MatrixForm[U]

isotropugr[t_]:= MatrixExp[U*t]

isotropugr[t]//MatrixForm

gh[t_]:= isotropugr[t]

Simplify[pssp[gh[t].vec[x], gh[t].vec[y]]]

v1c[t_] := gh[t].stb[4]

v2c[t_] := gh[t].stb[5]

v1c[t]

v2c[t]

plotcircle2spv[v1c[1.],v2c[1.]]

circlefamily[v1c,v2c,t,0,10, 20,PlotRange -> All,PlotStyle-> Green]

Show[%,PlotRange -> All]

circlefamily[v1c,v2c,t,-10,0, 20,PlotRange -> All,PlotStyle-> Green]

Show[%,PlotRange -> All]

circlefamily[v1c,v2c,t,-3,3, 30,PlotRange -> All,PlotStyle-> Green]

Show[%,PlotRange -> All]

isotrop = %;

Show[Graphics3D[{PointSize[0.02], Point[{0,1,0}]}]]

isoge ={%, isotrop}

Show[%, PlotRange -> All]

isoge = %;

euklidsphereplot3D[stb[5],-1,1,-1,1,Mesh -> False, Axes -> False, PlotStyle-> Opacity[.2]]

sph1=%;

Show[{sph1,isoge,grc1},Boxed -> False]

This is a 1-parameter circle family on the unit sphere consisting of all circles tangent to the equator at the point {0,1,0}. We check this using Solve:

Solve[{pssp[stb[5],vec[x]]== 0, pssp[stb[4],vec[x]]== 0, pssp[v1c[1.],vec[x]]== 0, pssp[v1c[.5],vec[x]]== 0},vec[x]]

vec[x]/.%

The solution is

pz =Flatten[%]

pssp[pz,pz]

Since the points correspond to isotropic vectors, the uniquely defined solution is given by x[2]= 0:

pzz = {t,0,t,0,0};stepro[pzz]

The representation of the unit sphere as the corresponding circle surface can be obtained with

ParametricPlot3D[circle2spv[v1c[t],v2c[t]][u], {t,-Pi,Pi},{u, 0, 2 Pi}]

The invariants along the isotropic geodesic

invisopp[s_,t_] := pp[v1c[s],v2c[s],v1c[t],v2c[t]]

Simplify[invisopp[s,t]]

As can be proved easily by elementary considerations, any two pairs of two tangent circles are Möbius equivalent.

6.4 Character of the Geodesics

The character of a geodesic is the character of the generating element in the pseudo-orthogonal Lie algebra. We calculate it using the corresponding Modules of the package liealg.m. Since the definitions in this package are not compatible with those of the already imported packages,  first quit the kernel. Then do not evaluate the initialization cells,  and import the package

SetDirectory["C:\\Users\\Rolf Sulanke\\mymath\\emg"];Needs["liealg`"]

?liealg`*

Since in this package the timelike vectors of the ordered bases ly at the end, set the dimensions and the index as

dim = 5; ind = 4;

this does not influence the pseudo-orthonormal Lie algebra. Redefine the generating matrices:

Uh={{0,0,0,0,1},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{1,0,0,0,0}};

Ue={{0,0,0,0,0},{0,0,0,1,0},{0,0,0,0,0},{0,-1,0,0,0},{0,0,0,0,0}};

U={{0,0,0,1,0},{0,0,0,0,0},{0,0,0,1,0},{1,0,-1,0,0},{0,0,0,0,0}};

Uh//MatrixForm

psokilling[Uh,Uh]

Uh is spacelike.

Ue//MatrixForm

psokilling[Ue,Ue]

Ue is timelike

U//MatrixForm

psokilling[U,U]

U is isotropic, as stated.

References

[G1] Alfred Gray. Modern Differential Geometry of Curves and Surfaces. CRC PRESS, 1993 (German translation: Spektrum Akademischer Verlag 1994)

[UM] U. May. Der Raum der Kreise in der Sphäre. Diplomarbeit, Humboldt-Universität, Berlin, 1993.

[OS3] A. L. Oniscik, R. Sulanke.  Projective and Cayley-Klein Geometries. Springer-Verlag. Berlin,Heidelberg. 2006
(German version may be downloaded from my homepage). Projektive und Cayley-Kleinsche Geometrien.  2003.

[Spheres] Sulanke, Rolf. Spheres.  1999.  http://mathsource.com/Content/Enhancements/Geometry/0210-092 . Obsolete. Better use the new versions from my homepage.

[S00] R. Sulanke. Möbius invariants for pairs ((Subscript[S, 1]^m), (Subscript[S, 2]^l)) of spheres in the Möbius space (S^n). Beiträge zur Algebra und Geometrie.41 (2000), No.1, 233-246.  Preprint No. 373, SFB 288, Berlin 1999.  http://xxx.lanl.gov/abs/math.MG/9902131.

[MS] R. Sulanke. Lie Algebras. http://www.mathsource.com/cgi-bin/msitem?0210-845 , or on my homepage

[MDG] R. Sulanke. Möbius Differential Geometry. In preparation. On my homepage.

[emg1] R. Sulanke. Elementary Möbius Geometry I. Points and Spheres. Mathematica Notebook emg.nb. 2011.  On my homepage.

[lox] R. Sulanke.  loxodromes.nb

Homepage

http://www2.mathematik.hu-berlin.de/~sulanke