Geometry of Circles

by

Rolf Sulanke

April 14, 1999, last revision November 30, 2000

Summary

This notebook describes the Moebius invariants for pairs of circles in the 3-sphere, or, by stereographical projection, in the euclidean 3-space. Section 2 contains the needed concepts for the euclidean geometry of circles in the 3-space (or in the 3-sphere). On the level of the euclidean 3-space we give  a parametrization of the 6-dimensional space of circles. We construct two plot commands which plot the circles corresponding to these parameters. Furthermore, we define functions giving the circle, and the euclidean invariants radius, center and position vector of the circle through three points on the euclidean level. As an application, a plot command for tubes of general space curves is developed. Section 3 gives the basic concepts for circles in the 3-dimensional Möbius space. The circles are represented by  2-dimensional euclidean subspaces of the 5-dimensional pseudo-euclidean vector space of index 1, which are defined by orthonormal pairs of vectors. By stereographical projection of the 3-sphere onto the euclidean 3-space the circles  become usual euclidean 3D-graphics (important for considering the mutual position of circles in the space). We construct procedures giving an adapted frame of the euclidean 2-space in terms of the euclidean circle invariants, and vice versa giving the circle as function of the euclidean 2-space.  In Section 4 we define a complete system of Möbius-geometric invariants for pairs of circles in the 3-sphere, and try to find out their geometric meaning. Section 5 gives the normal forms of the circle pairs in relation to this system of invariants. Geodesics  are studied as circle orbits and plotted in section 6.  The last section imports  ThreeScript and MathGL3d, which can be applied for animations of the created 3D-graphics. 
Many, but not all constructs developed in this notebook are collected in the package mcirc.m.
This notebook is related to the notebooks eusphere.nb, which contains the metric geometry of spheres, and the notebook mspheres.nb, devoted to the Möbius geometry of spheres. It uses tools of  pseudo-euclidean linear algebra, developed in the notebook pseuklid.nb. 
For differential geometry with Mathematica in the euclidean 3-space see the textbook A. Gray [1]. The concepts of Möbius geometry applied here are contained in the paper [4], see also my homepage [6]. A systematic description of the circle space  can be found  in the diploma paper of U. May [3].

References

[1] Alfred Gray. Modern Differential Geometry of Curves and Surfaces. CRC PRESS, 1993 (German translation: Spektrum Akademischer Verlag 1994)
[2] Jens-Peer Kuska. http://phong.informatik.uni-leipzig.de/ ~kuska/mview3d.html.
[3] U. May. Der Raum der Kreise in der Sphäre. Diplomarbeit, Humboldt-Universität, Berlin, 1993.
[4] 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. 
[5] R. Sulanke. Lie Algebras. http://www.mathsource.com/cgi-bin/msitem?0210-845
[6] R. Sulanke. Homepage. http://www-irm.mathematik.hu-berlin.de/~sulanke/

Keywords

pseudo-euclidean geometry,  Moebius geometry, 3D-circles, circles defined by subspaces, orthogonal circles, isospherical circles, stationary angles, eigenspheres, invariants of pairs of circles, normal forms of pairs of circles, geodesics in the circle space, MathGL3d.

Address

sulanke@mathematik.hu-berlin.de

The author will be very thankful for hints, corrections, comments.

Initialization

The needed packages

Before initializing the present notebook ensure that the packages
euvec.m
pseuvec.m
mspher.m
mcirc.m
and the master package init.m are in your working directory. These packages are contained in the item "Spheres".
 For the present notebook the package spiralsf.m is not needed.

3D graphics, animations

In the notebook command lines using ThreeScript, or MVShow3D appear. They may used for 3D-representations and animations of Mathematica Graphics3D objects. To apply them you must have the necessary software on your computer. Look in the last section.

For Windows

For the commands in this subsection:
DO NOT CONFIRM THE EVALUATION OF THE INITIALIZATION CELLS !!!

If the result of the following command

Directory[]

does not coincide with the working directory, e.g. "C:\Eigene Dateien\mymath\spheres.v4", containing the needed packages, insert the Path of your working directory in the next command line and activate it first, without calling the initialization:

SetDirectory["C:\\Eigene Dateien\\mymath\\spheres4"]

If you want to use the package euvec.mpseuvec.m,mspher.m only,  carry out the next command.

Needs["mspher`"]

If you want  that all commands contained in the packages of the item "Spheres" are available, activate

Get["init.m"]

Now: Activite Kernel/Evaluation/Evaluate Initialization from the menu!

For Linux (and other  systems: not tested!)

Activite Kernel/Evaluation/Evaluate Initialization from the menu!

1 Preliminaries

1.1 Pseudo-euclidean vector algebra and other prerequisites

Möbius geometry is treated by tools of euclidean and pseudo-euclidean linear algebra. The basic concepts are developed in  the notebooks eusphere.nb,  pseuklid.nb and contained in the packages euvec.m,  pseuvec.m. Basic concepts of Möbius geometry, in particular the Möbius geometry of spheres, can be found in the notebook mspheres.nb and the package mspher.m. All constructs can be called as part of the master package init.m of this directory. The symbols are listed below, see section 1.2.

Before using them one has to fix the dimension dim of the underlying vector space and the index  ind of the (non-degenerate, symmetric) scalar product. These basic constants are now defined by

dim  =  5;  ind  =  1;

because in this notebook the 3-dimensional sphere is treated as the elliptical hyperquadric in the four-dimensional projective  space, which belongs to the five-dimensional real vector space. You can change these constants only before starting the evaluation of the initialization cells  by changing this command line.. But remember that these constants enter in almost all of the treated concepts; therefore they should be defined as global constants before using the packages.

The matrix of the scalar products of any orthonormed basis of the pseudo-euclidean vector space is

MatrixForm[Table[pssp[stb[i],  stb[j]], {i, 5}, {j, 5}]]

1.2 The usages

After a package has been loaded you can get information about the  listed concepts in the usual way, e.g.

?orthonorm

Context[orthonorm]

For getting the usage and the complete definition, enter e.g.

?? spt

Since the names are prepended by the ContextNames the definition obtained in this way are hard to read. Therefore you will find some of them below in the original form in unevaluable cells.  It is recommended to look in the packages.

Symbols in mcirc.m

Here you can get information about the commands which are defined in mcirc.m:

\(TableForm[{"adaptsplframe", "center", "posvec", "radius", "center3pts", "circle2spv", "circle3D", "circle3pts", "circlefamily", "circleplane", "circlespacevectors", "control", "gencircle", "plotcircle2spv", "plotcircle3D", "plotcircle3pts", "plotgencircle", "posvec3pts", "pp", "ppdet", "pptr", "pscomplement", "radius3pts", "threepoints", "tube", "unitvec", "vspace"}]\)

To get the usages activate the next cell!

\(\({\(?"adaptsplframe"\), \(?"center"\), \(?"posvec"\), \(?"radius"\), \(?"center3pts"\), \ \(?"circle2spv"\), \ \(?"circle3D"\), \(?"circle3pts"\), \(?"circlefamily"\), \(?"circleplane"\), \(?"circlespacevectors"\), \ \(?"control"\), \ \(?"gencircle"\), \ \ \ \(?"plotcircle2spv"\), \(?"plotcircle3D"\), \ \(?"plotgencircle"\), \(?"plotcircle3pts"\), \(?"posvec3pts"\), \ \(?"pp"\), \(?"ppdet"\), \ \(?"pptr"\), \(?"pscomplement"\), \(?"radius3pts"\), \(?"threepoints"\), \ 
\ \(?"tube"\), \ \(?"unitvec"\), \ \(?"vspace"\)};\)\)

Symbols in euvec.m

\(TableForm[{"cross", "hyperplane", "invsterproj", "dotnorme", "esorthonorm", "innerprod", "normed", "print", "neglect", "nul", "ortho1", "orthobase", "orthoerg", "orthoframe", "outzero", "randomv", "rank", "renorm", "smoothing", "sphere", "sphereplot3D", \ "subsphere", "subspheremf", "sphereplot3D", "sphericalreflection", "stb", "sterproj", "tg4frame", "unit", "vec"}]\)

To get the usages activate the next cell!

\(\({\(?"cross"\), \(?"hyperplane"\), \(?"invsterproj"\), \(?"dotnorme"\), \(?"esorthonorm"\), \(?"innerprod"\), \(?"normed"\), \(?"print"\), \(?"neglect"\), \(?"nul"\), \(?"ortho1"\), \(?"orthobase"\), \(?"orthoerg"\), \(?"orthoframe"\), \(?"outzero"\), \(?"randomv"\), \(?"rank"\), \(?"renorm"\), \(?"smoothing"\), \(?"sphere"\), \(?"subsphere"\), \(?"subspheremf"\), \(?"sphereplot3D"\), \(?"sphericalreflection"\), \(?"stb"\), \(?"sterproj"\), \(?"tg4frame"\), \(?"unit"\), \(?"vec"\)};\)\)

?? invsterproj

Symbols in pseuvec.m

\(\(\(TableForm[{"ch", "dual", "ide", "normalize", "orthonorm", "orthopair", "pr", "psCross", "pscross", "psfilter", "pssp"}]\)\(\ \)\)\)

To get the usages activate the next cell!

\(\({\(?"ch"\), \(?"dual"\), \(?"ide"\), \(?"normalize"\), \(?"orthonorm"\), \(?"orthopair"\), \(?"pr"\), \(?"psCross"\), \(?"pscross"\), \(?"psfilter"\), \(?"pssp"\)};\)\)

?psfilter

?? pscross

?? psCross

Symbols in mspher.m

\(\(\(TableForm[{"confinv", \ "gencos", \ "invstepro", \ "planev", \ "plotsphere", 
"euklidsphere", \ "pstg5frame", \ "euklidsphereplot3D", \ "vradius", \ "vcenter", \ "spacevec", 
"sphvec", \ "spt", "spunit", \ "stepro"}]\)\(\ \)\)\)

To get the usages activate the next cell!

\(\({\(?"confinv"\), \ \(?"gencos"\), \ \(?"invstepro"\), \ \(?"planev"\), \ \(?"plotsphere"\), 
\(?"euklidsphere"\), \ \(?"pstg5frame"\), \ \(?"euklidsphereplot3D"\), \ \(?"vradius"\), \ \(?"vcenter"\), \ \(?"spacevec"\), 
\(?"sphvec"\), \ \(?"spt"\), \(?"spunit"\), \ \(?"stepro"\)};\)\)

?? spunit

1.3 Some tools

A general vector, random vectors

vec[x_,  n_:  dim]  :=  Table[x[i],  {i,  n}]

vec[Sin]

N[%]

randomv[n_:  dim]  :=  Table[Random[], {i,  n}]

randomv[]

Finding orthogonal pairs

The following definition is given in pseuklid.m:

?orthopair

Examples

orthopair[{0.5, 0.5, 1, 0, 0}]

orthopair[{0.5, 0.5, 1, 1, 1}]

ow = %

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

orthopair[stb[1]]

orthopair[stb[5]]

Context[orthopair]

Orthogonal complement in the euclidean 4-space

?orthoerg

Control of orthonormality

If necessary,  one can apply the following control:

?control

control[v_List, w_List,  prod_:  pssp] := 
      Simplify[(((prod[v, v] - 1))^2) + 
          prod[v, w]^2 + (((prod[w, w] - 1))^2)]

If all is correct the result must be 0, or numerically very small. If the result depends on a parameter, integrate!

EXAMPLE:

orthoerg[{Cos[t], Sin[t],  0, 0},  stb[2, 4]]

control[%[[1]], %[[2]],  Dot]

control[stb[1],  stb[5]]

(va  =  randomv[];  Print[va, ", ",  ch[va]]; 
      vb  = orthopair[va]; 
      Print[vb]; 
control[vb[[1]],  
        vb[[2]]])

Remove[va,  vb]

Smoothing

The following function makes an almost zero  number to be zero; the sensibility can be regulated by the optional second argument, whose default value is (10^(-10)); for the difference to Chop see the remark in the notebook eusphere.nb.

smoothing[a_, b_:  10^-10] := 
        Module[{c = a}, If[b == 1, Return[c], If[Abs[c] < b, c = 0]];  
          c];

Attributes[smoothing]

2 Euclidean geometry of circles in the 3-space

2.1 Parametrisation of the circle space. The general circle

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}]

ParametricPlot3D[unit[2][{u, v}], {u,  0, 2 Pi},  {v, -Pi/2, 
        Pi/2}]

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

?unitvec

Context[unitvec]

unitvec[u, v][t]

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

Simplify[%]

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

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

Simplify[%]

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[{%%, %}]

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

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

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

Show[{%%, %}]

cp1  =  %

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.

Example

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

\(\((u1\  = \ 2\ Pi\ Random[]; \ v1\  = \ Random[]\ Pi/2; 

      ParametricPlot3D[\(circle3D[0, 0, 0, u1, \ v1, \ 1]\)[t], {t, 0, 
          2\ Pi}])\)\)

\(Table[\((u1\  = \ 2\ Pi\ Random[]; \ v1\  = \ Random[]\ Pi/2; 

        ParametricPlot3D[\(circle3D[0, 0, 1, u1, \ v1, \ 1]\)[t], {t, 0, 
            2\ Pi}])\), {i, 10}]\)

Show[%]

Clear[u1, v1]

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, \ 
          Return[Print["Incorrect arguments!!!"]]]; 
	\ \ 	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]]]\)

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[nul[3],  stb[3, 3]]

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

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

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

plotcircle3D[1, 5 Random[], 1, Pi/2, Pi Random[],  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\  ≠ nul[3], Return[smoothing[a/Sqrt[a . a]]], \ 
          Return[nul[3]], Return[smoothing[a/Sqrt[a . a]]]]]\)

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[{q1 =  p1  -  p3,  q2  = p2  -  p3, dn, dn1, 
          dn2}, 
dn  =  
          Det[{{q1 . q1,  q1 . q2}, {q1 . q2, 
                q2 . q2}}]; 
If[dn  ==  0,  
          Return[{Infinity,  Infinity, 
              Infinity}]]; 
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\  == {Infinity, Infinity, 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[1, 3], stb[1, 3]]

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

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}, 
If[
          p1\  == \ {Infinity, Infinity, \ Infinity}, 
          Return[p2\  + \ t*\((p3 - p2)\)]]; 
If[
          p2\  == \ {Infinity, Infinity, \ Infinity}, 
          Return[p1\  + \ t*\((p3 - p1)\)]]; 
If[
          p3\  == \ {Infinity, Infinity, \ Infinity}, 
          Return[p1\  + \ t*\((p2 - p1)\)]]; 

        m\  = \ center3pts[p1, p2, p3]; 
If[
          m\  == \ {Infinity, Infinity, Infinity}, 
          If[smoothing[p2\  - p1]\  == nul[3], \ 
            Return[p1\  + \ t*\((p3\  - \ p1)\)], \ 
            Return[p1\  + \ 
                t*\((p2\  - \ p1)\)]]]; 
r\  = \ 
          Sqrt[\((p1 - m)\) . \((p1\  - \ m)\)]; 	
	
        q2\  = \ p2 - p1; \ q3\  = \ p3 - p1; 
        a\  = \ Cross[q2, q3]; 
		q3\  = \ Cross[a, q2]; 
		q2\ \  = \ 
          q2/Sqrt[q2 . q2]; 
		q3\  = \ \ q3/
            Sqrt[q3 . q3]; 
		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[{3, 0, 0}, {3, 1, 0}, {3, -1, 0}][t]

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

a = Random[];  b = Random[];  c = Random[];
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.05], Point[pa], Point[pb], 
            Point[pc]}];
Show[{gr1, gr2}]

center3pts[pa, pb, pc]

radius3pts[pa, pb, pc]

posvec3pts[pa, pb, pc]

% . %

a = Random[];  b = Random[];
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}]

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}]

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*Random[];  b = 4*Random[];  c = 4*Random[];
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.05], Point[pa], Point[pb], 
            Point[pc]}];
gr3  =  Graphics3D[Line[{pm, pm  + pv}]];
Show[{gr1, gr2, gr3}]

a = Random[];  b = Random[];
Print[{a, b}];
pa = {a, b, 0}; pb = {Random[], b, Random[]*a}; 
    pc = (pa  +  pb)/2;
Print[center3pts[pa, pb, pc], ", ", 
        radius3pts[pa, pb, pc]];
gr1 = plotcircle3pts[pa, pb, pc];
gr2  =  
        Graphics3D[{PointSize[0.05], 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.05], Point[pa], Point[pb], 
            Point[pc]}];
Show[{gr1, gr2}]

Remove[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[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 =  %

Animation: Before using ThreeScript first import the package (see section 7.1 below)!

ThreeScript["heltube.ts",  heltube]

MVShow3D[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.3]\)[t], {s, \(-8\), 
        8}, 
	{t, \ 0, \ 2\ Pi}, PlotPoints\  -> \ 40]\)

heltube2 = %

ThreeScript["heltube2.ts",  heltube2]

MVShow3D[heltube2]

The following definition is more general as that given in Alfred Gray's book [1], 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]]\)

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

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 [1].  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]

ast3  =  %

ThreeScript["ast3.ts",  %]

MVShow3D[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:

curv[
        s_]  :=  {((2 s Cos[Log[s]])/(1 + (s^2))), ((2 s Sin[Log[s]])/(1 + (s^2))), ((1 - (s^2))/(1 + (s^2)))}

Simplify[curv[s] . curv[s],  s > 0]

N[Log[30]]

ParametricPlot3D[curv[s],  {s, 0.1, 30}, 
      PlotRange  ->  All]

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

tube[curv, 0.0001,  20, 0.02, PlotPoints  ->  40,  
      PlotRange  ->  All]

The following parameter transformation gives simpler formulas for the evaluations:

Simplify[curv[s] /.  s -> Exp[t],  
      t ∈  Reals]

curve[t_]  =  %

ParametricPlot3D[curve[s],  {s, -10, 10}, 
      PlotRange  ->  All]

tube[curve, -3*Pi,  3 *Pi, 0.02,  PlotRange  ->  All, 
      PlotPoints  ->  {40, 20}]

loxotube  =  %

MVShow3D[loxotube]

Tubes of torus knots

Now we introduce a torus knot:

curv[
        s_]  :=  {Cos[
            2 s] (1 + 2 Cos[3 s]), (1 + 2 Cos[3 s]) Sin[2 s], 
        Sin[3 s]}

ParametricPlot3D[curv[s],  {s, 0. , 2 Pi}]

knot1  =  %

tube[curv, 0,  2 Pi, 0.02,  PlotRange  ->  All,  
      PlotPoints -> {50, 25}]

knot =  %

ThreeScript["knot.ts", knot]

Animation: Before using the next command first link MathGL3d (see section 7.2 below)!

MVShow3D[knot]

MVShow3D[knot1,  MVLineTubeSize  ->  0.005]

Clear[curv]

The general circular torus is the surface of revolution

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}]

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]}

?LCM

?GCD

If p,q, are integers torkn is a closed curve, i.e. it is periodic with the smallest period 2*Pi/GCD[p,q].

ParametricPlot3D[torkn[1, 1/3][2, 3][t],  {t, 0, 2*Pi},  
      PlotPoints  ->  100]

kn23 = %

ThreeScript["kn23.ts", kn23]

MVShow3D[kn23,  MVLineTubeSize  ->   .02]

curv[t_]  := torkn[1, 1/3][2, 3][t]

tube[curv, 0, 2*Pi, 0.04, PlotPoints  -> {80, 20}]

kntub23  =  %

MVShow3D[kntub23]

ThreeScript["kntub23.ts", kntub23]

ParametricPlot3D[torkn[1, 1/3][5, 2][t],  {t, 0, 2*Pi}, 
      PlotPoints  ->  100]

kn52  =  %

ThreeScript["kn52.ts", kn52]

Clear[curv]

curv[t_] := torkn[1, 1/3][5, 2][t]

tube[curv, 0, 2*Pi, 0.03, PlotPoints  -> {90, 20}]

kntub52 = %

ThreeScript["kntub52.ts", kntub52]

ParametricPlot3D[torkn[1, 1/3][7, 3][t],  {t, 0, 2*Pi}, 
      PlotPoints  ->  100]

kn73  =  %

ThreeScript["kn73.ts", kn73]

Clear[curv]

curv[t_] := torkn[1, 1/3][7, 3][t]

tube[curv, 0, 2*Pi, 0.03, PlotPoints  -> {90, 20}]

kntub73 = %

MVShow3D[kntub73]

ThreeScript["kntub73.ts", kntub73]

Clear[curv]

3 Circles defined by 2-dimensional euclidean vector subspaces (dim = 5, ind = 1)

In Möbius geometry a circle on the 3-sphere, or by stereographical projection, a circle or a line in the euclidean 3-space, is defined 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 can be defined by the orthogonal complement, that is a 2-dimensional euclidean subspace of the 5-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(5,1). It is well known that the orbits have a natural invariant structure of a pseudo-Riemannian symmetric space.

3.1 From the subspace to the circle

In this section we start with an arbitrary basis of an euclidean 2-space, 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.

The definitions pscomplement, threepoints, circle2spv, plotcircle2spv

Let v, w be a pair of  vectors spanning an euclidean 2-space in the pseudo-euclidean 5-space. Then one obtaines an pseudo-orthonormal basis of the wanted orthogonal complementary 3-space as follows:

?pscomplement

pscomplement[v_List, 
        w_List]:= Module[{j = 0, w1, w2, w3, 
          chj}, 
Label[chj]; 
If[
          j == dim,  
          Print["pscompl not found! Vectors do not span an euclidean \
2-space?"]; Return[]; Break[], 
          j = j + 1];               
 
        w1  =  pscross[{v, w, stb[j], stb[5]}]; 
If[
          w1  ==  nul[5],  Goto[chj],  
          w2  =  pscross[{v, w, w1, stb[5]}],  
          w2  =  pscross[{v, w, w1, stb[5]}]]; 
w3  =  
          pscross[{v, w, w1, w2}]; 
Return[{w1/
              Sqrt[pssp[w1, w1]], w2/Sqrt[pssp[w2, w2]], 
            w3/Sqrt[-pssp[w3, 
                    w3]]}]]

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[v_List, w_List][t_] :=  
      Module[{psc}, 
psc  =  
          pscomplement[v, w]; 
If[psc  ==  Null,  
          Return[]]; 
stepro[-psc[[1]]*Cos[t]  - 
            psc[[2]]*Sin[t]  +  psc[[3]]]]

?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[v_, w_] := \ 
      Module[{b = pscomplement[v, w]}, 
	
        If[b == Null, Return[]; 
          Break[]]; 		
{stepro[
            b[\([3]\)] + \ \ b[\([1]\)]], \ 
          stepro[b[\([3]\)]\  - \ b[\([1]\)]], \ 
          stepro[b[\([3]\)] + \ b[\([2]\)]]}]\)

Combining threepoints with circle3pts, plotcircle3pts of section 2 one obtaines a parameter representation and the plotting procedure plotcircle:

?circle3pts

OBSOLET:

\(\(oldcircle2spv[v_, \ w_]\)[t_]\  := \ 
      Module[{p1, p2, p3, b}, 
		b\  = 
          threepoints[v, w]; 
	\ \ \ \ If[b == Null, 
          Return[]; Break[]]; 
		p1 = \ \ b[\([1]\)]; 
		p2\  = 
          b[\([2]\)]; 
		p3\  = \ \ b[\([3]\)]; 
		\(circle3pts[p1, p2, 
            p3]\)[t]]\)

?plotcircle3pts

\(oldplotcircle2spv[\ v_, w_, \ opt___] := \ 
      Module[{p1, p2, p3, b}, 
		b\  = 
          threepoints[v, w]; 
\ 
        If[b == Null, Return[]; 
          Break[]]; 
		p1 = \ \ b[\([1]\)]; 
		p2\  = 
          b[\([2]\)]; 
		p3\  = \ \ b[\([3]\)]; 
		plotcircle3pts[p1, 
          p2, p3, opt]]\)

The possible options are those of ParametricPlot3D:

Options[ParametricPlot3D]

3.2 Examples, tests

First examples

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

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

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

pscomplement[stb[2],  stb[2] + stb[4]]

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

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

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

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

The Frame Circles

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

grc1 = plotcircle2spv[stb[3],  stb[4]]

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

grc2 = plotcircle2spv[stb[2], stb[4]]

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

grc3 = plotcircle2spv[stb[1], stb[4]]

This is the x-axis:

grc4 = plotcircle2spv[stb[2], stb[3]]

This is the y-axis:

grc5 = plotcircle2spv[stb[1], stb[3]]

This is the z-axis:

grc6 = plotcircle2spv[stb[1], stb[2]]

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

framecir  =  %

Further examples

In this example circle2spv, plotcircle2spv are shown to be stable against linear combination of the arguments; it would give an incorrect result in case a=0, or b = 0:

pscomplement[stb[1]*b + stb[4]*c, stb[4]*a]

Timing[circle2spv[stb[1]*b + stb[4]*c, stb[4]*a][t]]

plotcircle2spv[stb[4]*c, stb[4]*a]

plotcircle2spv[stb[1]+stb[4]*c, stb[4]*0]

Of course, with unspecified parameters plotcircle may not work, in general:

plotcircle2spv[stb[1]*b + stb[4]*c, stb[4]*a]

Here the parmeter c is unessential:

plotcircle2spv[stb[1] + stb[4]*c, stb[4]]

Moreover, the last example shows  that circle2spv does not require neither orthogonality nor norming; generally, in  all  modules using circle2spv or the 2-space spanned by v, w, these vectors may be an arbitrary basis of the 2-space.

threepoints[stb[3],  stb[4]]

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

threepoints[stb[1],  stb[2]]

Timing[plotcircle2spv[ stb[1], stb[2]]]

\(\((wr\  = 
        orthopair[{2*Random[], \(-2\)*Random[], 1. , 3*Random[], \ 
            Random[]}]; \ 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.05], Point[pa], Point[pb], 
            Point[pc]}]; 
Show[{gr1, gr2}])\)\)

Timing[plotcircle2spv[ v, w]]

v = v +  0.1*stb[1];  w  = w +  0.01*stb[5];

Timing[plotcircle2spv[ v, w]]

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

3.3 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.  We have two such "spheres" canonical associated to the circle: The plane which containes 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 mspheres.nb. For a plane we have

?planev

planev[a, b, c, p]

spt[stb[4]  +  stb[5]]

stepro[stb[4]  +  stb[5]]

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

pssp[planev[a, b, c, p], stb[4]  +  stb[5]]

We know that the property nv[[4]] == 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[4] + stb[5]] of the 3-sphere, and only these, are projected onto planes by the stereographic projection.  Since

{a, b, c} . {a, b, c}  ==  
      pssp[{a, b, c, p, p}, {a, b, c, p, p}]

for all values of p we may define the normed vector nv easily if we found 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:

(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_]  := {Cos[u] Cos[v], Cos[v] Sin[u], 
        Sin[v], x Cos[u] Cos[v] + y Cos[v] Sin[u] + z Sin[v], 
        x Cos[u] Cos[v] + y Cos[v] Sin[u] + z Sin[v]}

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

?spacevec

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

spacevec[x, y, z, r]

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

Clear[a, b, c, p]

Obviously, the circlespacevectors are explicitly defined in terms of the euclidean invariants of the circle. 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, 4]) == (Subscript[v, 5]), and any orthonormal basis of this 2-space will be considered as adapted to the line.

Tests

Let us check our considerations:

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

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

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

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

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]

Simplify[%,  u > 0]

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

Show[%[[2]]]

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.4 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[vx_List, vy_List]\  := \ 
      Module[{aob, d1, d2, q\  = \ 0, v, nor}, 
aob\  = \ 
          orthonorm[{vx, vy}, \ 
            normed\  -> \ False]; 
If[
          Length[aob] < 
            2, 
Print["The span is not 2-dimensional and spacelike!"]; Return[]]; 
If[
          Simplify[pssp[aob[\([1]\)], \ aob[\([1]\)]]] ≤ \ 0\  || 
            Simplify[pssp[aob[\([2]\)], \ aob[\([2]\)]]] ≤ \ 
              0\ , 
Print["The span is not 2-dimensional and spacelike!"]; Return[]]; 
d1 = 
          Simplify[\(aob[\([1]\)]\)[\([5]\)]\  - \ \(aob[\([1]\)]\)[\([4]\)]]; 
d2 = \ 
          Simplify[\(aob[\([2]\)]\)[\([5]\)]\  - \ \(aob[\([2]\)]\)[\([4]\)]]; 
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]\)]\)[\([4]\)]]; 
d2 = \ 
          Simplify[\(aob[\([2]\)]\)[\([5]\)]\  - \ \(aob[\([2]\)]\)[\([4]\)]]; 
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[1],  stb[1]]

adaptsplframe[stb[1],  stb[4]]

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

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

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

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

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

Simplify[%,  s > 0 && a > 0]

This is not a direct result, but it coincides with the direct one:

Simplify[% /.  a  ->  0]

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

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

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

plotcircle2spv[vx, vy]

plotcircle2spv[ot[[1]], ot[[2]]]

Remove[vx, vy, ot, s1, t1]

3.5  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.

The euclidean invariants with  frame adaption

?circlespacevectors

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

\(posvec[vx_List, vy_List]\  := \ 
      Module[{ab}, 

        ab = adaptsplframe[vx, 
            vy]; 
If[\(ab[\([2]\)]\)[\([4]\)]\  == \ \(ab[\([2]\)]\)[\([5]\)], \ Return[{0, 0, 0}]]; 
Return[
          Take[ab[\([1]\)], 3]]]\)

\(center[vx_List, vy_List]\  := \ 
      Module[{ab}, 

        ab = adaptsplframe[vx, 
            vy]; 
If[\(ab[\([2]\)]\)[\([4]\)]\  == \ \(ab[\([2]\)]\)[\([5]\)], \ 
          Return[{Infinity, Infinity, Infinity}]]; 
Return[
          Take[ab[\([2]\)], 
              3]/\((\(ab[\([2]\)]\)[\([5]\)]\  - \ \(ab[\([2]\)]\)[\([4]\)])\)]]\)

\(radius[vx_List, vy_List]\  := \ 
      Module[{ab}, 

        ab = adaptsplframe[vx, 
            vy]; 
If[\(ab[\([2]\)]\)[\([4]\)]\  == \ \(ab[\([2]\)]\)[\([5]\)], \ Return[Infinity]]; 
Return[
          1/\((\(ab[\([2]\)]\)[\([5]\)]\  - \ \(ab[\([2]\)]\)[\([4]\)])\)]]\)

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)]]]

Examples

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

posvec[stb[1], stb[3]]

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

center[stb[2], stb[1]]

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

radius[-stb[4], stb[1]]

radius[stb[3], stb[1]]

(wr  = 
        orthopair[{2*Random[], -2*Random[], 1. , 3*Random[],  
            Random[]}];  v  =  wr[[1]]; 
      w  =  wr[[2]])

Timing[posvec[v, w]]

Timing[posvec1[v, w]]

Timing[center[v, w]]

Timing[center1[v, w]]

Timing[radius[v, w]]

Timing[radius1[v, w]]

Remove[v, w]

\(Timing[\((wr\  = 
          orthopair[{2*Random[], \(-2\)*Random[], 1. , 3*Random[], \ 
              Random[]}]; \ v\  = \ wr[\([1]\)]; w\  = \ wr[\([2]\)]; 
        mwr = threepoints[v, w]; 
pa = mwr[\([1]\)]; pb = mwr[\([2]\)]; 
        pc = mwr[\([3]\)]; 
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]}]; 

        gr3\  = \ Graphics3D[Line[{pm, \ smoothing[pm\  + pv]}]]; 

        Show[{gr3, gr1, gr2}])\)]\)

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]}];
gr3  =  Graphics3D[Line[{pm,  smoothing[pm  + pv]}]];
Show[{gr3, gr1, gr2}]

Remove[wr, v, w, mwr, pa, pb, pc, pv, pm, gr1, gr2, gr3]

\(Timing[\((wr\  = 
          orthopair[{2*Random[], \(-2\)*Random[], 1. , 3*Random[], \ 
              Random[]}]; \ v\  = \ wr[\([1]\)]; w\  = \ wr[\([2]\)]; 
        mwr = threepoints[v, w]; 
pa = mwr[\([1]\)]; pb = mwr[\([2]\)]; 
        pc = mwr[\([3]\)]; 
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]}]; 

        gr3\  = \ Graphics3D[Line[{pm, \ smoothing[pm\  + pv]}]]; 

        Show[{gr3, gr1, gr2}])\)]\)

\(Timing[\((wr\  = 
          orthopair[{2*Random[], \(-2\)*Random[], 1. , 3*Random[], \ 
              Random[]}]; \ v\  = \ wr[\([1]\)]; w\  = \ wr[\([2]\)]; 
        mwr = threepoints[v, w]; 
pa = mwr[\([1]\)]; pb = mwr[\([2]\)]; 
        pc = mwr[\([3]\)]; 
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]}]; 

        gr3\  = \ Graphics3D[Line[{pm, \ smoothing[pm\  + pv]}]]; 

        Show[{gr3, gr1, gr2}])\)]\)

Remove[wr, v, w, mwr, pa, pb, pc, pv, pm, gr1, gr2, gr3]

s1  =  Random[];  t1  =  Random[]*Pi;
ot  =  
        adaptsplframe[{Sqrt[1 + s1^2], 0, 0, 0, s1}, {0, Cos[t1], Sin[t1], 0, 
            0}];
gt1  =  
        plotcircle2spv[ot[[1]],  
          ot[[2]]];
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]]]]];
gt3 = euklidsphereplot3D[ot[[2]]];
Show[{gt1, gt2, gt3}]

gtest  =  %

MVShow3D[gtest]

Remove[ot, s1, t1]

3.6 From  three points to the subspace of the corresponding circle

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

\(vspace[p1_, p2_, p3_] := 
      Module[{q1, q2, q3, j = 1, w1, w2, chj}, 
        q1 = Simplify[invstepro[p1]]; 
q2 = 
          Simplify[invstepro[p2]]; 
q3 = 
          Simplify[invstepro[p3]]; 
Label[
          chj]; 
If[j == dim, 
          Print["vspace not found! Points not in general position?"]; 
          Return[], j = j + 1]; 
w1 = 
          Simplify[pscross[{q1, q2, q3, stb[j]}]]; 
If[
          smoothing[w1] == nul[dim], Goto[chj], 
          w1 = Simplify[w1/Sqrt[pssp[w1, w1]]], 
          w1 = Simplify[w1/Sqrt[pssp[w1, w1]]]]; \ \ \ \ \ \ w2 = 
          Simplify[pscross[{q1, q2, q3, w1}]]; 
If[
          w2 == nul[dim], Goto[chj]]; 
w2 = 
          Simplify[w2/Sqrt[pssp[w2, w2]]]; 		Return[
          smoothing[{w1, w2}]]]\)

Examples

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

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[{1, 0, 0},  {0, 1, 0}, {1, -1, 0}]

vs = %

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

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\  = \ Random[]; 
b\  = \ 2*Random[]; 
c\  = \ \(-\ Random[]\); \ 
      Print[TableForm[{a, b, c}]]; 

      Timing[plotcircle3pts[{a, 0, 0}, \ {0, b, 0}, {0, 0, 
            c}]]
)\)\)

\(\((a\  = \ Random[]; 
b\  = \ 2*Random[]; 
c\  = \ \(-\ Random[]\); \ 
      Print[TableForm[{a, b, c}]]; 

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

\(\(\(	\)\(vs\  = \ %[\([2]\)]\)\)\)

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

\(\(\(		\)\(Timing[
      plotcircle3pts[{a, 0, 0}, \ {0, b, 0}, {0, 0, c}]]\)\)\)

Remove[a, b, c, vs]

4 Position and invariants of two circles

In this section we apply the general theory, see the paper [4], to the case of circles in the 3-sphere and, of course, the euclidean 3-space. For the Möbius geometry of spheres see the notebook  mspheres.nb. Invariants for the other possible pairings ((S^k), (S^m)), 0 ≤ k, m ≤ 3, are treated in the notebook pairs.nb.

4.1 Case discussion

Notations: C0, C1 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 does not 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.		(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

Case criterion

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

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

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[nul[]]

From the general properties of the cross product one easily concludes: 
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 U1 or U1. ch[C0xC1] depends only of the (non-oriented) circles C0,C1.

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

D[ftest[t],  t]

Remove[ftest]

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[5]*\ Sinh[s]\  + \ stb[2]*\ Cosh[s]; 
	
      wa1[s_]\  := \ stb[3]*Cosh[s]\  + \ stb[5]*Sinh[s]; \ 
      wa2[t_]\  := \ stb[2]*Sin[t]\  + \ stb[4]*\ Cos[t]; 

      C0xC1\  = \ pscross[{va1[t], va2[s], wa1[s], wa2[t]}])\)\)

Simplify[%]

Simplify[ch[%]]

Remark: For s = 0 or t=(2k+1]Pi/4 we have case 2. For all other values of t,s we have case 1. a : Non-intersecting, separated circles.

gr1 = plotcircle2spv[va1[Random[]*Pi],  va2[Random[]*Pi]]
gr2  =  
      plotcircle2spv[wa1[Random[]*Pi],  wa2[Random[]*Pi]]
Show[{gr1, gr2},  PlotRange  ->  All]

separ1  =  %

Apply MathGL3D to visualize the mutual position!

MVShow3D[separ1]

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:

\(\((vb1[t_]\  := \ stb[3]*Cos[t]\  + \ stb[4]*Sin[t]; \ 
      vb2\ [s_] := \ \(-stb[1]\)*\ Sin[s]\  + \ stb[2]*\ Cos[s]; 
	
      wb1[s_]\  := \ stb[1]*Cos[s]\  + \ stb[2]*Sin[s]; \ 
      wb2[t_]\  := \ \(-stb[3]\)*Sin[t]\  + \ stb[4]*\ Cos[t]; 

      C0xC1\  = \ pscross[{vb1[t], vb2[s], wb1[s], wb2[t]}])\)\)

Simplify[%]

ch[%]

It follows: 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 four cells and applying MathGL3d, see section 7.2 below. (By numerical reasons it sometimes happens that the random circles have "corners"; try to avoid that increasing the number of plotpoints.)

gr1 = plotcircle2spv[vb1[Random[]*Pi],  vb2[Random[]*Pi]]
gr2  =  
      plotcircle2spv[wb1[Random[]*Pi],  wb2[Random[]*Pi]]
Show[{gr1, gr2},  PlotRange  ->  All]

inter1  =  %

MVShow3D[inter1]

Example: Random pair, case 1.c

We shall define random circles intersecting at the origin:

\(\((pa = randomv[3]; pb\  = \ randomv[3]; \ qa\  = \ 2*randomv[3]; \ 
      qb\  = \ 2*randomv[3]; 
	mvp\  = \ vspace[nul[3], pa, pb]; \ 
      mvq\  = vspace[nul[3], qa, qb]; 
	
      Print[smoothing[
          pscross[{mvp[\([1]\)], mvp[\([2]\)], mvq[\([1]\)], 
              mvq[\([2]\)]}]]]; 
		Print[
        ch[pscross[{mvp[\([1]\)], mvp[\([2]\)], mvq[\([1]\)], 
              mvq[\([2]\)]}]]]; 
	gr1\  = \ plotcircle3pts[nul[3], pa, pb]; 
      gr2\  = \ plotcircle3pts[nul[3], qa, qb]; 
	
      Show[{gr1, gr2}, \ PlotRange\  -> \ All])\)\)

onept  =  %

ParametricPlot3D[{x, y, 0}, {x, -1, 2},  {y,  -1, 
        2}]

Show[{onept, %}]

onept1  =  %

MVShow3D[onept1,  MVLineTubeSize  ->  0.01]

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

Stationary angles

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 [4].

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}.

\(\(\ 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}}\)\)

The invariants of pp are the determinant and the trace:

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}]

pp[v1, v2, w1, w2]

ppdet[v1, v2, w1, w2]

Simplify[%]

pptr[v1, v2, w1, w2]

We show that  ppdet  is the square of the determinant of the scalar products pssp[xi, yj], i, j = 1, 2, and therefore always non negative:

ppdet[x1, x2, y1, 
        y2]  - ( 
          pssp[x1,  y1]*pssp[x2,  y2]  -  
            pssp[x1,  y2]*pssp[x2,  y1])^2

Simplify[%]

It follows ppdet[v1_,v2_,w1_,w2_]>= 0,  but Mathematica does not now that our square is a real number:

TrueQ[ppdet[v1_, v2_, w1_, w2_] >=  0]

TrueQ[Simplify[ppdet[v1_, v2_, w1_, w2_]] >=  0]

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

Using the vector object

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

Eigenvalues[pp[stb[3],  stb[4],  vec[w1], vec[w2]]]

Simplify[%]

For shortness we write a,b,c,d instead of w1[3]=pssp[v1,w1], w2[4]= pssp[v2,w2],  w2[3]=pssp[v1,w2],  w1[4]=pssp[v2,w1], and consider the argument of Sqrt:

poly[a_, b_, c_, 
        d_]  :=  -4*(c*d  -  a*b)^2  +  (a^2  + b^2  + 
              c^2  + d^2)^2

Factor[poly[a, b, c, d]]

Since the square of the first terms in Simplify[Eigenvalues] above coincides with ((((a^2) + (b^2) + (c^2) + (d^2)))^2) and the first term is non negative we conclude for the eigenvalues α1, α2 of pp:

α1 ≥ α2 ≥ 0.

Simplify[((((a^2) + (b^2) + (c^2) + (d^2)))^2)  -  ((a^2) + 2 a b + 
              (b^2) + (c^2) - 2 c d + (d^2)) ((a^2) - 2 a b + (b^2) + 
              (c^2) + 2 c d + (d^2))]

If U0, U1 are mutually orthogonal, we say that the corresponding circles are orthogonal. 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[1],  stb[2]]; 
    gr2  =  plotcircle2spv[stb[3], stb[4]]; Show[{gr1, gr2}]

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

4.3 Case 1.a: disjoint,  non-interlaced circles

Example

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

\(\((v1\  = \ {1, \ 0, \ 1, \ \(-1\), \ 0}/Sqrt[3]; \ 
      v2\  = \ stb[2]\ ; 
	w1\  = \ {3\/\@2, 1\/\@2, 0, 1\/\@2, 3\/\@2}; \ 
      w2\  = \ stb[3]\ ; 

      C0xC1\  = \ pscross[{v1, v2, w1, w2}])\)\)

ch[%]

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

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

Show[%, PlotRange  -> All]

testgr  =  %

pp[v1, v2, w1,  w2]

Simplify[%]

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

N[%]

Sqrt[%]

As it follows from the general theory [4], 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  =  
      normalize[-v1*((-3 + (Sqrt[57]))/(4 (Sqrt[3])))  +  v2]

eigenv2  =  
      normalize[v1 *((3 + (Sqrt[57]))/(4 (Sqrt[3])))  +  v2]

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, 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 mspheres.nb) we calculate

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

Simplify[Together[%]]

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 eigenv1 .

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

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.

The two distinguished spheres are

{euklidsphereplot3D[eigenv1], euklidsphereplot3D[eigenv2]}

Show[%]

eigensph  =  %

MVShow3D[eigensph]

pssp[eigenv1, eigenv2]

Simplify[Together[%]]

The intersection is the circle C0.

plotcircle2spv[v1, v2]

Show[{%,  eigensph}]

testgr0  =  %

MVShow3D[testgr0,  MVLineTubeSize  ->  0.02]

Table[
      euklidsphereplot3D[uv1[s *2 Pi], -1, 1, -1, 1,  
        PlotPoints  ->  30],  {s,  0.3,  0.5,  0.1}]

testgr1 =  %

Show[testgr1]

testgr1  =  %

euklidsphereplot3D[eigenv2]

testgr2 =  Flatten[{testgr1, %}]

Show[%]

testgr2  =  %

radius[w1, w2]

N[%]

N[center[v1, v2]]

N[posvec[v1, v2]]

There are many possibilities for experimental graphics...

The projection of the eigenvalue eigenv2 onto the 2-space U1 is

preigenv2  =  
      w1*pssp[w1, eigenv2]  +  w2*pssp[w2, eigenv2]

Simplify[%]

npreigenv2  =  normalize[preigenv2]

Simplify[%]

N[pssp[npreigenv2, eigenv2]]

This is the square root of the eigenvalue. Since it is greater 1, the spheres defined by these vectors are disjoint and separate the circles

plotcircle2spv[v1, v2]

c1  =  %

euklidsphereplot3D[eigenv2]

s1  =  %

Show[{c1, s1}]

sc1  =  %

plotcircle2spv[w1, w2]

c2  =  %

euklidsphereplot3D[npreigenv2]

s2  =  %

Show[{c2, s2}]

sc2  =  %

Show[{sc1, sc2}]

test  =  %

Now link MathGL3d (see section 7.2), apply

MVShow3D[test, MVLineTubeSize  ->  0.01]

and look at the picture with Appearance/WireFrame, spin it around the axis. You will see the eigensphere eigenv2, containing the circle C0, and the sphere corresponding to the projection npreigenv2, containing C1. The two spheres separate the circles.

Another example

nv = N[circlespacevectors[2, 0, 0, Pi/4, Pi/4,  1]]

\(\((v1\  = \ stb[2]\ ; \ v2\  = \ stb[4]\ ; 
	w1\  = nv[\([1]\)]; \ 
      w2\  = \ \ nv[\([2]\)]; 

      C0xC1\  = \ pscross[{v1, v2, w1, w2}])\)\)

ch[%]

It follows: case 1.a.

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

Show[ %]

pp[v1,  v2, w1,  w2]

Eigenvalues[%]

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

eigenv =  %

eigenv1  =  
      normalize[-v1*0.25666793515702424` -  
          v2*0.9664996487646695`]

eigenv2  =  
      normalize[
        v2*0.25666793515702424` -  v1*0.9664996487646695`]

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];

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

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

fctr[s_, t_]  =  Chop[%]

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

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

Maximum and Minimum are defined by the square root of the eigenvalue:

Sqrt[2.1327822185373186`]

Table[
      euklidsphereplot3D[uv1[s *2 Pi], -1, 1, -1, 1,  
        PlotPoints  ->  30],  {s,  0.3,  0.5,  0.1}]

testgr1 =  %

euklidsphereplot3D[eigenv1]

testgr  =  Flatten[{testgr1,  %}]

Show[%]

testgr  =  %

The other eigenvector of pp[v1,v2,w1,w2]   corresponds to t=Pi/2 or 3 Pi/2, where the spheres uv1[s] intersect the eigensphere   in different angles, whose cosinus are defined by:

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

Maximum and Minimum are defined by the square root of the eigenvalue:

Sqrt[0.11721778146268136`]

euklidsphereplot3D[eigenv2]

There are many possibilities for experimental graphics...

testgr  =  Flatten[{testgr1,  %}]

Table[Show[{testgr[[i]],   testgr[[4]]}], {i, 3}]

Show[%]

testgr2  =  %

Remove[v1, v2, w1, w2, fctr, testgr,  grsph1, uv1, uv0]

4.4 Case 1.b: Interlaced circles

Already at the end of section 4.2 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.

Example: Isogonal circles

A first example  of an isogonal pair:

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

      griso\  = \ plotcircle2spv[w1, w2])\)\)

Show[{grc1, griso}]

isosph = %

We again got an interlaced pair of circles; for seeing this clearer apply MathGL3d, section 7.2.  The unique eigenvalue is

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

Eigenvalues[%]

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[3] Cos[t]  +   stb[4] 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 mspheres.nb) we calculate

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

Simplify[%]

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:

fctr[0,  Pi/4]

N[%]

fctr[0,  5 Pi/4]

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.

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.
Arguments: By theorem 6 of [4] 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. The value α  = 1 corresponds to a pair of tangent circles (case 2.c); all these pairs are Möbius equivalent, see also section 5.6.

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

grisofam[s]

iso2[s]

Simplify[pssp[iso2[s], iso2[s]]]

pp[stb[3],  stb[4],  iso1[s],  iso2[s]]

Eigenvalues[%]

Together[%]

Simplify[%]

f[s_] := Csch[s]^2/(1 + Csch[s]^2 + Sech[s]^2) + 
   Sech[s]^2/(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}]

The case s= 0 , f(0) = 1, does not give a circle since s = 0 is a singularity of grisofam:

grisofam[0]

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___]\  := \ 
      ParametricPlot3D[
		Evaluate[
		Table[\(circle2spv[v[t], \ 
                w[t]]\)[u], \ {t, \ a, \ b, \ \((b - a)\)/d}]], {u, 0, 
          2\ Pi}, \ opt]\)

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

isofam  =  {grc1, %}

Show[%]

isofam =  %

ThreeScript["isofam.ts", isofam]

The circle family isofam is contained in the circle surface

\(ParametricPlot3D[\(circle2spv[\(grisofam[s]\)[\([1]\)], \ \(grisofam[
              s]\)[\([2]\)]]\)[t], 
	{s, \ 0.1, 2.1}, \ {t, 0, \ 
        2\ \ Pi}, \ PlotPoints\  -> \ 50]\)

isosf  =  %

Show[{grc1, isosf}]

isosf =  %

Constructing smaller parts of the circle surface, e.g.  for 0.5 < s < 1.5, is recommended.

\(ParametricPlot3D[\(circle2spv[\(grisofam[s]\)[\([1]\)], \ \(grisofam[
              s]\)[\([2]\)]]\)[t], 
	{s, \ 0.4, 1.7}, \ {t, 0, \ 
        2\ \ Pi}, \ PlotPoints\  -> \ {30, 50}]\)

isosfb  =  %

Show[{grc1, isosfb}]

isosfb = %

ThreeScript["isosfb.ts", isosfb]

The meaning of  'isogonal'

Spheres, containing the circles:

Clear[fct, s, b, a,  sph2]

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

sph2[s_, b_]  :=  iso1[s]*Cos[b]  +  iso2[s]*Sin[b]

The cosinus of the intersection angle of two of such spheres is:

fct[s_][a_, b_] :=  pssp[sph1[a],  sph2[s, b]]

fct[s][a, b]

Simplify[%]

Clearly, for fixed a (or b) the resulting functions of  b (resp. a) are cos--functions with the same amplitude, depending on s only, and distinct phases.  For non isogonal circle pairs the amplitude is not constant (see the next subsection).

Plot3D[fct[1. ][a, b], {a, 0, 2*Pi}, {b, 0, 2*Pi}]

tgri =  %

MVShow3D[tgri]

Clear[a, b];  a  = 2*Pi* Random[];  
    Plot[fct[1. ][a, b], {b, 0, 2*Pi}]; a

Table[ (a  = 2*Pi*Random[];  
        Plot[fct[1. ][a, b], {b, 0, 2*Pi}]), {5}]

Show[%]

Intersection angle for spheres containing non isogonal circles

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

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

Here are two orthonormal vectors defined "at random":

br1  =  {-0.72439280526, 0, 0.6893874554, 0, 0}

br2 =  {0, -0.33744663644, 0, 0.9413446593, 0}

pssp[br1, br2]

pssp[br1, br1]

pssp[br2, br2]

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  =  %

MVShow3D[ptgr]

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

Clear[a, b];  a  =  Random[];  Plot[fct2[a, b], {b, 0, 2*Pi}]; 
    a

Clear[a]

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

amplitude[a_]  := 
      Sqrt[(0.689387455 Cos[a] )^2 + (0.9413446593 *
                Sin[a])^2]

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

Interlaced, non isogonal circles

Of course, the random circle pair, case 1.b, in section 4.2 are non-isogonal with probability 1:

\(\((vb1[t_]\  := \ stb[3]*Cos[t]\  + \ stb[4]*Sin[t]; \ 
      vb2\ [s_] := \ \(-stb[1]\)*\ Sin[s]\  + \ stb[2]*\ Cos[s]; 
	
      wb1[s_]\  := \ stb[1]*Cos[s]\  + \ stb[2]*Sin[s]; \ 
      wb2[t_]\  := \ \(-stb[3]\)*Sin[t]\  + \ stb[4]*\ Cos[t];)\)\)

\(Eigenvalues[
      pp[vb1[Random[]*Pi], \ vb2[Random[]*Pi], 
wb1[Random[]*Pi], \ 
        wb2[Random[]*Pi]]]\)

With probability 1 one gets two distinct eigenvalues, 0 < (Subscript[α, 1])<(Subscript[α, 2])<1.
We consider the following example:

\(\((v1\  = \ stb[3]\ ; \ v2\  = \ stb[4]\ ; 
	
      w1\  = \ {\(-\(3\/5\)\), 0, 0, 4\/5, 0}; \ w2\  = \ stb[2]\ ; 

      C0xC1\  = \ pscross[{v1, v2, w1, w2}])\)\)

ch[%]

It follows: case 1.b; the circles are interlaced.

plotcircle2spv[w1, w2]

Show[{grc1,  %}]

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

Eigenvalues[%]

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[3] Cos[t]  +   stb[4] 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 mspheres.nb) we calculate

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

fctr[s_, t_]  =  %

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

We have the eigenvector stb[3] of pp[v1,v2,w1,w2], which correspond to t=0 or Pi, where the spheres uv1[s] intersect the sphere stb[3] (that is the x,y-plane) orthogonally:

Table[
      euklidsphereplot3D[uv1[s *2 Pi], -1, 1, -1, 1 , 
        PlotPoints  ->  30],  {s,  0.3,  0.5,  0.1}]

testgr1 =  %

euklidsphereplot3D[stb[3]]

testgr  =  Flatten[{testgr1,  %}]

Show[%]

testgr  =  %

The other eigenvector of pp[v1,v2,w1,w2]  is stb[4], which correspond to t=Pi/2 or 3 Pi/2, where the spheres uv1[s] intersect the sphere stb[4] (that is the unit sphere with center o)  in different angles, whose cosinus are defined by:

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

Maximum and Minimum are defined by the square root of the eigenvalue, here 4/5.  This is the unit sphere

euklidsphereplot3D[stb[4]]

grsph1 = %

There are many possibilities for experimental graphics...

testgr  =  Flatten[{testgr1,  grsph1}]

Table[Show[{testgr[[i]],   testgr[[4]]}], {i, 3}]

Show[%]

testgr2  =  %

4.5 Case 1c: Circles intersecting in exactly one point

First example

See also the random pair, case 1.c, in section 4.2. Circles in orthogonal planes intersecting at the origin:

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

v1  =  %[[1]]

v2  =  %%[[2]]

plotcircle2spv[v1,  v2]

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

w1  =  %[[1]]

w2  =  %%[[2]]

plotcircle2spv[w1, w2]

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

ch[%]

It follows: case 1.c; the circles are intersecting in one point and are in "general position".

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

Show[%]

testgr  =  %

pp[v1, v2, w1,  w2]

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

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

eigenv1  = v1

eigenv2  =  v2

The two  eigenspheres are

{euklidsphereplot3D[eigenv1], euklidsphereplot3D[eigenv2]}

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

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

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 mspheres.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 generalized cosinus with the sphere eigenv1 .

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

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

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

Clearly, that means zero; we have only numerical errors,  remember:

fctr[s, t]

fctr[s, Pi/2]

It is a surprise by plotting zero to get a  -sin curve. Try better

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

pssp[eigenv1, uv1[s]]

uv1[s]

eigenv1

The amplitudes correspond to the square roots of the corresponding eigenvalues.

Table[
      euklidsphereplot3D[uv1[s *2 Pi], -1, 1, -1, 1,  
        PlotPoints  ->  40],  {s,  0.4,  0.5,  0.1}]

testgr1 =  %

Show[%]

euklidsphereplot3D[eigenv1]

testgr2 =  Flatten[{testgr1, %}]

Show[%]

testgr2  =  %

There are many possibilities for experimental graphics...

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

Another example

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

bsp  =  smoothing[%]

v1  = stb[3]

v2  =  stb[4]

w1  =  bsp[[1]]

w2  =  bsp[[2]]

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

Show[%]

exp1c2  =  %

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

ch[%]

This corresponds to case 1.c. The unique intersection point is {1,0,0}.

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], euklidsphereplot3D[eigenv2]}

Show[%,  PlotRange  -> All]

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

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 mspheres.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 followwing generalized cosinus with the sphere eigenv1 :

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

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

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

pssp[eigenv1, uv1[s]]

uv1[s]

eigenv1

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

Sqrt[0.0737033666969973921`]

Table[
      euklidsphereplot3D[uv1[s *2 Pi], -1, 1, -1, 1,  
        PlotPoints  ->  40],  {s,  0.1,  0.3,  0.1}]

testgr1 =  %

Show[%]

euklidsphereplot3D[eigenv1]

testgr2 =  Flatten[{testgr1, %}]

Show[%]

testgr2  =  %

There are many possibilities for experimental graphics...

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

4.6 The cases 2

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_] := 
  {{x/r, y/r, 0, 
    -(1 + r^2 - x^2 - y^2)/(2*r), (1 - r^2 + x^2 + y^2)/(2*r)}, 
   stb[3]}

vspace2d[x, y, r]

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

Together[%]

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

      plotcircle2spv[w1, w2])\)\)

a) Two disjoint circles:

\(\((v1\  = \ stb[3]; \ v2\  = \ stb[4]\ ; 
	
      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 linearily dependent.

ch[%]

{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[%]

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[%]

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 mspheres.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`]

b) Two intersecting circles:

\(\((v1\  = \ stb[3]; \ v2\  = \ stb[4]\ ; 
	
      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 linearily dependant.

{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[%]

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 mspheres.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}]

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

c) Two tangent circles:

\(\((v1\  = \ stb[3]; \ v2\  = \ stb[4]\ ; 
	
      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 linearily dependant.

{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[%]

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];

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}]

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],  euklidsphereplot3D[uv1[0]]}

Show[%]

testgr  =  %

Call MathGL3d (section 7.2) with testgr and use wire frame!

MVShow3D[testgr]

4.7 Further examples

The x-axis and the standard unit circle compared with  circle families

The x-axis

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

gro  =  %

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

cv[t]

cv1[t_]  := cv[t][[1]]; 
    cv2[t_] :=  cv[t][[2]]

circlefamily[cv1, cv2, t, -1, 1]

Show[{gro, %}]

cvpp[t_]  :=  pp[stb[2], stb[3],  cv1[t],  cv2[t]]

cvpp[t]

The circles for t ≠0 are disjoint to the x-axis; the maximal eigenvalue is 1+(t^2)>1corresponding to case 1a.

pscross[{stb[2], stb[3],  cv1[t],  cv2[t]}]

Replacing the x-axis by the standard unit circle we obtain

Clear[cvpp]

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

cvpp[t]

Eigenvalues[%]

R e m a r k. In this example the circle pair is not in general position:

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

The sphere

vn[t_]  := spacevec[0, 0, t/2,  Sqrt[1  +  (t/2)^2]]

vn[t]

contains both circles of the pair belonging to the parameter t. For the standard unit circle this is obvious. On the other hand, considering

compl[t_] :=  pscomplement[cv1[t], cv2[t]]

compl[t]

we obtain

Table[
      Simplify[pssp[vn[t],  compl[t][[i]]]],  {i, 3}]

what proves the statement. Therefore, the fact that there exist two eigenvalues greater or equal 1 does not contradict to the statement after formula (50) of the paper [4], since assumption B is not fulfilled, which is essential for the proof of the eigenvalue properties.

Eigenvectors[cvpp[t]]

Clear[eigenv1]

eigenv1[t_]  :=  normalize[{0, 0, -t/2,  1, 0}]

eigenv1[t]

This again is the sphere containing both circles of the pair (up to orientation):

Simplify[vn[t]  +  eigenv1[t],  t >  0]

Clear[eigenv2]

eigenv2[t_]  :=  normalize[{0, 0, 2/t, 1, 0}]

eigenv2[t]

For t= 1 we get the corresponding eigenspheres

euklidsphereplot3D[eigenv1[1. ]]

es1  =  %

euklidsphereplot3D[eigenv2[1. ]]

es2  =  %

Show[{es1, es2}]

vradius[eigenv1[t]]

Simplify[%]

Simplify[vcenter[eigenv1[t]]]

Table[euklidsphereplot3D[eigenv1[t]], {t,  0, 5}]

Show[%]

Simplify[vradius[eigenv2[t]]]

Simplify[vcenter[eigenv2[t]]]

The normed projections of eigenv1,eigenv2 on the vspace spanned by cv1,cv2 are

projeigenv1[t_]  :=  
      normalize[
        cv1[t]*pssp[eigenv1[t], cv1[t]] + 
          cv2[t]*pssp[eigenv1[t], cv2[t]]]

projeigenv2[t_]  :=  
      normalize[
        cv1[t]*pssp[eigenv2[t], cv1[t]] + 
          cv2[t]*pssp[eigenv2[t], cv2[t]]]

Simplify[pssp[projeigenv1[t], projeigenv2[t]]]

euklidsphereplot3D[projeigenv1[1. ]]

pes1  =  %

euklidsphereplot3D[projeigenv2[1. ]]

pes2  =  %

Show[{pes1, pes2}]

vradius[projeigenv1[t]]

Simplify[%]

Simplify[vcenter[projeigenv1[t]]]

Simplify[vcenter[projeigenv2[t]]]

Simplify[vradius[projeigenv2[t]]]

Clear[cv, cvpp, eigenv1, vn, eigenv2]

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, %}]

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

cvpp[t]

For all t >0 we have case 2:

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

fct[t_]  = pssp[%, %]

Together[Simplify[fct[t]]]

Plot[fct[t], {t, 0, 3}]

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.

circlefamily[cv1, cv2, t, 0.2, 3.2, 10]

Show[{%, grc1}]

cfv1  =  %

MVShow3D[cfv1]

5 Normal forms of pairs of circles

In the general classification of pairs {k-sphere, l-sphere} contained in an n-sphere  ( see [4]) we obtained normal forms for the bases of the generating vector spaces, where the coefficients of one base with respect to the other are the invariant  eigenvalues of the matrix pp. In this section we describe these normal forms for each of the cases mentioned in the subsection "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 = plotcircle [stb[3],stb[4]] 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, the numbering of the basis vectors, and the notations differ from those used in [4]. 
ATTENTION: For getting the normal forms the definitions in subsection 5.1 are needed. Run 5.1 first!

Show [grc1]

5.1 Some vector functions needed for the definition of the normal forms

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

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

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

5.2 Case 1a: disjoint, non-interlaced circles

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

pscross[{stb[3],  stb[4], 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[3],  stb[4],  v4[a],  v3[b]]]]

circ1a[a_,  b_]  :=  plotcircle2spv[v4[a],  v3[b]]

For a = b = 0 we get ste standard circle grc1 again.

circ1a[0, 0]

A first series

We compare a series of circles with grc1

Table[
      Show[{grc1, circ1a[a,  a + 1]}], {a,  0.2, 3.2,  0.5}]

Show[%,  PlotRange  -> All]

nf1a  = %

MVShow3D[nf1a,  MVLineTubeSize  ->  0.05]

A surface generated by the circles circ1a[a,a+1]:

ParametricPlot3D[
      Evaluate[circle2spv[v4[a],  v3[a + 1]][t]], {a,  0,  5},  {t,  
        0,  2 Pi}, PlotPoints  ->  90,  
      PlotRange  ->  All]

ngr1a  =  %

MVShow3D[ngr1a]

Limit[circle2spv[v4[a],  v3[a + 1]][0],  
      a  ->  Infinity]

Limit[circle2spv[v4[-a],  v3[-a + 1]][1. ],  
      a  ->  Infinity]

ThreeScript["nf1a.ts", nf1a]

ThreeScript["ngr1a.ts", ngr1a]

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

ngr1a1  =  %

MVShow3D[ngr1a1,  MVLineTubeSize  ->  0.03]

ThreeScript["ngr1a1.ts", ngr1a1]

Another series

circ1a[0, 1]

?center

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])}

?radius

radius[v4[a],  v3[b]]

Simplify[%, { a ∈ Reals, b ∈ Reals}]

radfct[a_, b_]  :=  ((Sqrt[2])/(Sqrt[Cos[2 b] + Cosh[2 a]]))

centfct[0, 1. ]

radfct[0, 1. ]

We compare a series of circles with grc1

Table[
      Show[{grc1, circ1a[a,  a + 1]}], {a,  0.2, 3.2,  0.5}]

Show[%,  PlotRange  -> All]

nf1a  = %

MVShow3D[nf1a,  MVLineTubeSize  ->  0.03]

A surface generated by the circles circ1a[a,a+1]:

ParametricPlot3D[
      Evaluate[circle2spv[v4[a],  v3[a + 1]][t]], {a,  0,  5},  {t,  
        0,  2 Pi}, PlotPoints  ->  30,  
      PlotRange  ->  All]

ngr1a  =  %

MVShow3D[ngr1a]

Limit[circle2spv[v4[a],  v3[a + 1]][0],  
      a  ->  Infinity]

Limit[circle2spv[v4[-a],  v3[-a + 1]][1. ],  
      a  ->  Infinity]

ThreeScript["nf1a.ts", nf1a]

ThreeScript["ngr1a.ts", ngr1a]

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

ngr1a1  =  %

MVShow3D[ngr1a1,  MVLineTubeSize  ->  0.03]

ThreeScript["ngr1a1.ts", ngr1a1]

5.3 Case 1b: interlaced circles

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

pscross[{stb[3],  stb[4], v5[a], v3[b]}]

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.

MatrixForm[Evaluate[pp[stb[3],  stb[4],  v5[a],  v3[b]]]]

circ1b[a_,  b_]  :=  
      plotcircle2spv[Evaluate[v5[a]], Evaluate[ v3[b]]]

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]

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  = %

MVShow3D[nf1b]

A corresponding surface:

ParametricPlot3D[
      Evaluate[circle2spv[v5[a],  v3[2*a]][t]], {a,  1,  2},  {t,  
        0,  2 Pi}, PlotPoints  ->  70]

ngr1b  =  %

Remark: If MathGL3d is already  loaded appears here a "MeshGraphics3D" object. This package changes the definition of Parametricplot3D. Of course, I could not find neither in the Help nor in the Mathematica book a description of this type of objects. MeshGraphics3D cannot be joined with a Graphics3D object:

Show[{grc1,  ngr1b}]

ngr1b1  =  %

MVShow3D[ngr1b1, MVLineTubeSize -> 0.03]

5.4 Case 2a

b4  =  stb[4]

The eigenvalues are obtained from the normal form

pp[stb[3],  stb[4],  v4[a],  stb[4]]

Both circles lie on the standard unit sphere. U0 +U1 is spanned by stb[3], stb[4], v4[a]. We set

mat[a_]  :=  {stb[3],  stb[4], 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[4]], {a,  0.2, 3.2,  0.5}]

Show[%]

nf2a =  %

ParametricPlot3D[
      Evaluate[circle2spv[v4[a],  stb[4]][t]], {a,  0,  Pi},  {t,  
        0,  2 Pi}]

ngr2a  =  %

ThreeScript["nf2a.ts", nf2a]

ThreeScript["ngr2a.ts", ngr2a]

5.5 Case 2b: 2-point intersecting circles in the x,y-plane

This case one gets specialising case 1a to a = 0. Then the eigenvalues are 1, cos^2[b]:

pp[stb[3],  stb[4],  stb[3],  v3[b]]

mat[b_]  :=  {stb[3],  stb[4], v3[b]}

and calculate

Det[
      Table[pssp[mat[b][[i]], mat[b][[j]]],  {i, 3},  {j, 
          3}]]

It follows that  for a ≠ 0 the subspace U0 + U1 is euclidean: case 2.b

Table[plotcircle2spv[stb[3], v3[b]], {b,  0.2, 3.2,  0.5}]

Show[%]

Show[{%,  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.

vec[x]

\(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]]\)

The solutions are, as stated:

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

Obviously, we have the 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[3], 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  :=  {0, 1, 0, -1, 1}

Therefore all the considered circles lie on the sphere

euklidsphereplot3D[b4]

MatrixForm[Evaluate[pp[stb[3],  stb[4],  v5[b],  b4]]]

pscross[{stb[3],  stb[4], 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 ]

Table[ 
      plotcircle2spv[v5[b*Pi],  b4],  {b,  0.1,  1,  0.3}]

Show[%]

sertest1  = %

Show[{sertest1, grc1}]

test1c = %

matc[b_]  :=  {stb[3],  stb[4], 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]]

All the circles intersect the standard unit circle in the uniquely defined point

stepro[{0,  1, 0, 0, 1}]

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]]

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

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[3],  stb[4],  v5[Pi],  b4]]]

pscross[{stb[3],  stb[4], v5[Pi],  b4}]

If b4 is replaced by the vector {0,1,0,1,1} one gets  a family of lines lying in the plane y = 1 through the point {0,1,0}.

euklidsphereplot3D[{0, 1, 0, 1, 1}]

v5[t]

radius[v5[t],  {0, 1, 0, 1, 1}]

center[v5[t],  {0, 1, 0, 1, 1}]

Table[
      plotcircle2spv[v5[s],  {0, 1, 0, 1, 1}], {s, 0, 5, 0.5}]

Show[%]

Especially we get the tangent:

plotcircle2spv[v5[Pi],  {0, 1, 0, 1, 1}]

Show[{grc1,  %}]

Of course, we have the same invariants:

MatrixForm[
      Evaluate[pp[stb[3],  stb[4],  v5[Pi],  {0, 1, 0, 1, 1}]]]

pscross[{stb[3],  stb[4], v5[Pi],  {0, 1, 0, 1, 1}}]

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

6 Circle surfaces, circle orbits, geodesics in the circle space

6.1 Introduction

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 canal surfaces, are circle surfaces.  A special case are the Dupin cyclides, which are canal surfaces in a twofold way. That means, that there exists two different sphere curves whose envelopes is the same canal surface: tori, circular coni,  circular cylinders, and their images under Moebius 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 [5].  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  [3] 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.2 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, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 
          0, 1}, {0, 0, 0, 1, 0}}

gUh[t_]   :=  MatrixExp[Uh t]

MatrixForm[gUh[t]]

vh1[t_] :=  gUh[t] . stb[4]

vh1[t]

Obviously, stb[3] is fixed under gUh[t]. Therefore the orbit of {stb[3], stb[4]} is given by

hypgeodesic[t_, u_] := circle2spv[vh1[t], stb[3]][u]

ParametricPlot3D[
      hypgeodesic[t, u], {t, -0.5, 1}, {u, 0, 2 π}]

hypgeo  =  %

This is the family of concentric circles in the x,y-plane containing the standard unit circle grc1 (for t = 0).

ThreeScript["hypgeo.ts", hypgeo]

The group is the group of all dilatations in the euclidean picture, leaving the "north pole" of (S^3) fixed.

The invariants along the hyperbolic geodesic

invpp[s_, t_] :=  pp[vh1[t], stb[3], vh1[s], stb[3]]

invpp[s, t]

Simplify[%]

The Invariants are the elements of the diagonal; the variable one 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}]

ftt1[s_, t_] :=  D[ftt[s1, t], s1] /. s1 -> s

ftt2[s_, t_] :=  D[ftt[s, t1], t1] /. t1 -> 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[3]]

Simplify[%]

For the varying invariant expressed as a function of the radius we obtain

invpp[t,  0][[1, 1]] /. (E^t) ->  
        r

% /. (E^(-t)) ->  r^(-1)

Plot[((((1/(2 r)) + (r/2)))^2), {r,  0.2, 5}]

Equivalent description on the 2-sphere

U3h  =  {{0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 1}, {0, 0, 0, 
          0, 0}, {0, 0, 1, 0, 0}}

MatrixForm[U3h]

gU3h[t_]   :=  MatrixExp[U3h t]

MatrixForm[gU3h[t]]

v3h1[t_] := gU3h[t] . stb[3]

v3h1[t]

hyp3geodesic[t_, u_] := circle2spv[v3h1[t], stb[4]][u]

ParametricPlot3D[
      hyp3geodesic[t, u], {t, -1. , 1}, {u, 0, 2 π}]

hyp3geo  = %

?circlefamily

xy[t_] :=  stb[4]

circlefamily[v3h1, xy, t, 0,  1. , 10,  
      PlotPoints  ->  100]

hgeofam  =  %

ThreeScript["hgeofam.ts", hgeofam]

6.3 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 liealgeb.nb ):

Ue  =  {{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}}

gUe[t_] := MatrixExp[Ue t]

gUe[t] // MatrixForm

ve1[t_]  :=  gUe[t] . stb[3]

ve1[t]

ve2[t_]  :=  stb[4]

circlefamily[ve1, ve2, t, -N[Pi], N[Pi],  30]

ellgeo  =  %

ThreeScript["ellgeo.ts", ellgeo]

The standard circle rotates around the y-axis. The surface creation  yields the unit sphere, with two omitted sections:

ellgeodesic[t_, u_] := circle2spv[ve1[t], stb[4]][u]

ParametricPlot3D[
      ellgeodesic[t, u], {t, -Pi/2, Pi/2}, {u, 0, 2 Pi},  
      PlotPoints  ->  30]

MVShow3D[%]

The strange behaviour of the picture at the borders t = Pi/2, -Pi/2 comes from the careful calculation of the expression ellgeodesic[t,u] by Mathematica :

ellgeodesic[t, u]

Simplification with the condition  t>-Pi/2&&t<Pi/2 yields a usual parameter representation of the unit sphere:

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

Clear[ellgeodesic]

ellgeodesic[t_, u_] :=  {Cos[t] Sin[u],  
        Cos[u], - Sin[t] Sin[u]}

ParametricPlot3D[
      ellgeodesic[t, u], {t, -Pi/2, Pi/2}, {u, 0, 2 Pi},  
      PlotPoints  ->  30]

The invariants along the elliptic geodesic

invellpp[s_, t_]  :=  pp[ve1[s], stb[4], ve1[t], stb[4]]

Simplify[invellpp[s, t]]

As expected, the only essential invariant is the cosinus square of the intersection angle.

6.4 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 [5], liealgeb.nb ):

U  = {{0, 0, 1, 0, 0}, {0, 0, 0, 0, 0}, {-1, 0, 0, 0, 1}, {0, 0, 
          0, 0, 0}, {0, 0, 1, 0, 0}}

MatrixForm[U]

isotropugr[t_] :=  MatrixExp[U*t]

isotropugr[t]

gh[t_] :=  isotropugr[t]

Simplify[
      Table[pssp[gh[t][[i]],  gh[t][[j]]], {i, 5}, {j, 
          5}]]

v1c[t_]  :=  gh[t] . stb[3]

v2c[t_]  :=  gh[t] . stb[4]

v1c[t]

v2c[t]

plotcircle2spv[v1c[1. ], v2c[1. ]]

circlefamily[v1c, v2c, t, 0, 10,  20, PlotRange  ->  All]

circlefamily[v1c, v2c, t, -10, 0,  20, 
      PlotRange  ->  All]

circlefamily[v1c, v2c, t, -3, 3,  30, 
      PlotRange  ->  All]

isotrop  =  %

Show[Graphics3D[{PointSize[0.05],  Point[{1, 0, 0}]}]]

isoge  = {%,  isotrop}

Show[%]

isoge  =  %

ThreeScript["isoge.ts", isoge]

MVShow3D[isoge, MVPointSphereSize -> 0.5, 
 
      MVLineTubeSize  ->  0.01]

This is a 1-parameter circle family on the unit sphere consisting of all circles tangent to the equator at the point {1,0,0}. We check this using Solve:

Solve[{pssp[stb[4], vec[x]] ==  0,  pssp[stb[3], vec[x]] ==  0,  
        pssp[v1c[1. ], vec[x]] ==  0,  pssp[v1c[ .5], vec[x]] ==  0}, 
      vec[x]]

The solution is

pz  =  {t, s, 0, 0, t}

pssp[pz, pz]

Since the points correspond to isotropic vectors, the uniquely defined solution is given by s = 0:

pzz  =  {t, 0, 0, 0, t}; stepro[pzz]

The representation as circle surface can be obtained with

ParametricPlot3D[circle2spv[v1c[t], v2c[t]][u],  {t, -Pi, 
        Pi}, {u,  0,  2 Pi}]

MVShow3D[%]

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 tw pairs of two tangent circles are Möbius equivalent.

7 ThreeScript, MathGL3d

If you want to generate *.ts files for animations outside of Mathematica, activate the next section. To use the animation program MathGl3d in this notebook, activate section 7.2.

7.1 Threescript

The ThreeScript commands generate *.ts files appropriate  for certain animation programs. The good freeware program  MathGL3d  created by Jens-Peer Kuska (for Windows95, Linux and other platforms) can be downloaded from [2].

Under Windows the commercial 3D-Animation programs MathLive and DynamicVisualizer are available. To produce .ts-files one

Needs["Graphics`ThreeScript`"]

?ThreeScript

The standard extension of these files is ts: "filename.ts" or "Path/filename.ts"

ThreeScript["isofam.ts",  isofam]

7.2 MathGL3d

The free animation program MathGL3d by Jens-Peer Kuska [2] can be used from a command line and within a Mathematica session. A new version of this program is in preparation. Jens-Peer Kuska has sent me a new startup file OpenGLViewer.m, and friendly permitted to add it to the item "Spheres". Replace the old file OpenGLViewer.m by this new OpenGLViewer.m. If you did this, changes of Mathematica's ParametricPlot3D, especially the sometimes disturbing generating of MeshGraphics3D objects, will not take place.
Establish the link to the current Mathematica session with the command:

Get["MathGL3d`OpenGLViewer`"]

If you want to close the link, enter

MVClose[]

The following command clears the Viewer window; this and other MathGL3d commands you may enter also with  palettes.

MVClear[];

MVShow3D[\[Placeholder], MVPolygonShading -> MVSmooth, 
      MVGrayBackground -> False, MVTexture -> None, 
      MVTextureMapType -> MVMeshUVMapping, 
      MVScaleTexture -> {1, 1}, 
      MVTranslateTextureTarget -> Automatic, 
      MVRotateTextureTarget -> Automatic, 
      MVScaleTextureTarget -> Automatic]

MVShow3D[gre1
, MVPointSphereSize -> 0.8,  
      MVLineTubeSize  ->  0.03]

MVShow3D[isoge, MVPointSphereSize -> 0.5, 
 
      MVLineTubeSize  ->  0.01]

MVShow3D[tgr3,  MVPointSphereSize  ->  1.8]

MVShow3D[ex123, MVPointSphereSize -> 1.8]

MVShow3D[hyp3geo]

Options[MVShow3D]

MVClose[]