Spheres in Euclidean Spaces

by Rolf Sulanke

February 24, 1999, last revision November 1, 2000

Summary

This notebook contains  basic definitions of spheres as objects of euclidean and  Riemannian spherical geometry. It needs the packages euvec.m, spiralsf.m. As an application the construction and plotting of a sphere through four points in the euclidean 3-space is given. Furthermore, it contains a recursive definition of the generalized geographical parameter representations of n-spheres in the (n+1)-dimensional euclidean space. Some concepts needed in Möbius geometry,  the conformal geometry of the n-sphere, are introduced in euclidean terms. These concepts are: Stereographic projection and its inversion, reflections at hyperspheres (also called inversions),  spiral transformations, spiral surfaces, in particular the spiral cylinder, which is of interest in Möbius differential geometry: it is a homogeneous surface with constant negative curvature, see [3]. Spiral surfaces have been considered already by S. Lie  (1891), and E. Vessiot (1926), these and other references see  [3]. 
Changes Feb 2000:  A version of Erhard Schmidt's orthogonalization: esorthonorm, is introduced in section 4.1. The procedure normalized is renamed in dotnorme, normed in normalize, and the new concept normed is now an option for the orthogonalization esorthonorm. The procedures called sphere1, sphere3, spherrefl are renamed: sphere1 -> subsphere, sphere3 -> subspheremf, spherrefl -> sphericalreflection. The procedures sphere, sphereplot3D are introduced in euvec.m.

References

[1] Gray, Alfred. Modern Differential Geometry of Curves and Surfaces. CRC PRESS, 1993 (German translation: Spektrum Akademischer Verlag 1994)

[2] Kuska, Jens-Peer. http://phong.informatik.uni-leipzig.de/ ~kuska/mview3d.html.

[3] Sulanke, Rolf. Möbius geometry V: Homogeneous surfaces in the Möbius space (S^3). Colloquia Math. J. Bolyai 46, Topics in Differential Geometry, Debrecen 1984, 1141-1154.

Keywords

vector objects, random vectors, rank,  orthoframes,  unitvectors, norming, cross product (general), outzero, nullvector, standard base ,  hyperplanes, stereographical projection, invers stereographical projection, 3-sphere, spheres, spheres through four points, inversion at hyperspheres, flat torus, geodesics on the flat torus, torus, Erhard Schmidt's othogonalization, orthogonal complement in the euclidean 4-space.  Parameter representation for n-spheres, spiral transformations, spiral surfaces, MathGl3d, ThreeScript.

Address

sulanke@mathematik.hu-berlin.de

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

Initialization

The needed packages

If you want tu use all the notebooks of the item "Spheres", ensure that the packages
euvec.m
pseuvec.m
mspher.m
mcirc.m
spiralsf.m
and the master package init.m are in your working directory. These packages are contained in the item "Spheres".
 For the present notebook we only need the packages euvec.m,  spiralsf.m.

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.m only,  carry out the next command.

Needs["euvec`"]

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

For the next command:
DO NOT CONFIRM THE EVALUATION OF THE INITIALIZATION CELLS !!!

If you want to use the package euvec.m only, remove the file init.m  from the working directory and carry out the next command.

Needs["euvec`"]

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

1 Introduction

1.1 The 3-sphere

In this notebook we develop some  linear algebraic tools which can be applied for calculations in any dimension and for graphics in the threedimensional case. The commands appropriate for general use are collected in the packages euvec.m. The procedure for creating spiral surfaces is contained in the package spiralsf. m.  One has to fix the dimension dim of the underlying vector space, and its index ind; these global constants enter in almost all of the defined functions and modules. In our case their values are

dim  =  4; ind  =  0;

For an elementary approach to the Riemannian geometry of the standard threedimensional sphere  ("3-sphere") it is convenient to consider the 3-sphere as the submanifold of the unit vectors in the fourdimensional  euclidean space  ("4-space"):

vec[x]

vec[x] . vec[x]  ==  1

1.2 The usages

Concepts in euvec.m

Valid for arbitrary values of the dimension dim are

?nul

nul[]

nul[7]

?vec

vec[y]

?randomv

randomv[]

?stb

?cross

?dotnorme

The next command generates an error message:

dotnorme[nul[]]

?outzero

?rank

?sterproj

?invsterproj

?hyperplane

?stb

stb[i]

Since we defined dim = 4 the next command generates an error message:

stb[5]

stb[5, 7]

stb[4]

hyperplane[Table[stb[j], {j, 4}]]

?unit

?sphericalreflection

?esorthonorm

?innerprod

?normed

?print

?neglect

?renorm

?orthobase

orthobase[stb[5, 7]]

orthobase[nul[5]]

?smoothing

The following commands work in dimensions  4 or 3  only:

?ortho1

?orthoerg

?orthoframe

?tg4frame

?sphere

?sphereplot3D

?subsphere

?subspheremf

Concepts in spiralsf.m

If the master package init.m is not in use first import

Needs["spiralsf`"]

?spiral

?spiralsurface

?plotspiralsurface

Full definitions

For getting the usage together with the complete definition, enter e.g.

?? sterproj

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 in the original form in unevaluable cells.  It is recommended to look into the packages.

1.3 Some tools

Vector objects

In euvec.m a general vector is defined:

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

Example

vec[Sin]

N[%]

In euvec.m a  random vector with values in the unit n-cube is defined as

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

randomv[]

vec[randomv]

Reflections at arbitrary spheres (inversions)

A point pair x, y is in inversion with respect to a hypersphere with radius r and center a iff a, x, y are collinear and fulfil the equation

|y - a| |x - a| = (r^2).

Here |z| means the norm of the vector z.  We remember

?Infinity

Infinity

This leads to the definition (contained in euvec.m):

?sphericalreflection

sphericalreflection[{Infinity},  a_List, r_:  1]  :=  
        a;
\(sphericalreflection[x_List, \ a_List, \ r_:  1]\  := \ 
      If[Not[x == a], 

        a\  + \ \((\((x\  - \ a)\)\ *\ 
                r^2)\)/\((\((x - a)\) . \((x - a)\))\), {Infinity}, 

        a\  + \ \((\((x\  - \ a)\)\ *\ 
                r^2)\)/\((\((x - a)\) . \((x - a)\))\)]\)

sphericalreflection[{2, 2, 2},nul[3]]

sphericalreflection[vec[x],nul[]]

The reflection at the unit sphere with further examples will be considered in section 5.

sphericalreflection[vec[x],  vec[a],  r]

vy  =  %;

Test: The spherical reflections fulfill the following identities:

Simplify[(vy  - vec[a]) . (vec[x]  -  vec[a]) ==  
        r^2]

sphericalreflection[vec[a],  vec[a],  
        r] == {Infinity}

sphericalreflection[{∞}, vec[a], r]== vec[a]

Clear[vy]

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

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

Attributes[smoothing]

Example: Chop and smoothing

There exists a similar Mathematica operation

?? Chop

In spite of the fact that the attribute Listable is not reported, Chop seems to have this attribute. Therefore instead of smoothing[a, b] we could have used Chop[a, b]. The behavior of smoothing differs from that of Chop; for our purposes smoothing is sometimes more appropriate, see  the examples below.

Chop[N[Exp[-40]]]

Chop[{N[Exp[-40]],  N[Exp[-35]]}]

Chop[{N[Exp[-40]],  Exp[-35]}]

Chop[{N[I*Exp[-40]],  Exp[-35]}]

Chop[{N[Exp[-40]], I*Exp[-35]}]

In the next example Chop produces an unexpected result; since t may be complex, the argument should have been returned:

Chop[N[Pi*10^(-19)]*Cos[t]]

N[Pi*10^(-19)]*Cos[I*100. ]

In difference to Chop smoothing evaluates the inserted expressions as expected:

smoothing[{N[Exp[-40]],  Exp[-35]}]

smoothing[{N[I*Exp[-40]],  Exp[-35]}]

smoothing[{N[Exp[-40]], I*Exp[-35]}]

smoothing[N[Pi*10^(-19)]*Cos[t]]

smoothing[N[Pi*10^(-19)]]

dotnorme

dotnorme normes a vector with respect to the Dot product. It can be applied also for symbolic calculations:

?dotnorme

dotnorme[25*unit[4][vec[s]]]

Simplify[%]

Simplify[% . %]

Error message:

dotnorme[nul[]]

dot norme has been generalized to arbitrary, may be indefinite scalar products:

?normalize

normalize[nul[]]

2 Stereographic projection

2.1  Definition of the general stereographic projection

We want to visualize the geometry of the 3-sphere. By stereographic projection we go to the euclidean 3-space. Then we apply the 3D-graphics tools of Mathematica. Here is a definition of the general stereographic projection of the n-sphere \ north pole onto the n-dimensional Euclidean space. We put the north pole at the point nul[n+1] + stb[n+1] and project onto the equatorial hyperplane; the north pole itself will be projected to Infinity.

?sterproj

The definition contained in the package euvec.m is:

sterproj[v_List] := 
      Module[{m = Length[v], l, w}, 
l = 
          nul[m]; 
w = v[[m]]; 
If[
          smoothing[w - 1] == 0, 
          For[i = 1, i < m, i++, l[[i]] = Infinity], 
          l = Simplify[stb[m, m] + (v - stb[m, m])/(1 - w)], 
          l = Simplify[
              stb[m, m] + (v - stb[m, m])/(1 - 
                      w)]]; 
l = 
          Take[l, m - 1]; 
l]

Examples

stb[1, 3]

sterproj[stb[1, 3]]

?unit

unit[2][{s, t}]

sterproj[unit[2][{s, t}]]

testvec[s_, t_]  =  AppendTo[%,  0]

ParametricPlot3D[
      testvec[s, t], {s, 0, 2  Pi},  {t, - Pi/2 + 0.0000001, 
        0}]

sterproj[unit[3][{s, t, u}]]

WARNING! Without any message sterproj gives generally incorrect results for incorrect, not normalized arguments:

sterproj[stb[2]  + stb[4]]

sterproj[(stb[2]  + stb[4])/Sqrt[2. ]]

sterproj[nul[]]

The image of the north pole stb[n,n]  is the infinite point, e.g.

sterproj[stb[3, 3]]

The other basis vectors are "preserved" in the equatorial hyperplane, e.g.

stb[3, 7]

sterproj[stb[3, 7]]

The general formula is:

vec[n]

sterproj[vec[n]]

The elements of the equator remain fixed:

vec[n]  -  n[4]*stb[4]

sterproj[(vec[n]  -  n[4]*stb[4])]

2.2 Example: Geodesics on the flat torus

We consider the torus as a direct product of two circles in orthogonal planes, producing a surface  in the unit 3-sphere in the euclidean 4-space. Obviously, its inner metric is flat. By  stereographic projection onto the euclidean 3-space we obtain the usual torus as a surface of revolution. Its parameter representation is (with 0 < c < 1):

tor[c_][s_, t_] := sterproj[{Cos[t], Sin[t], 0, 0}*c +
{0, 0, Cos[s], Sin[s]}Sqrt[1 - c^2]]

tor[c][s, t]

tor[0][s, t]

tor[1][s, t]

ParametricPlot3D[tor[0.8][s, t], {s, 0, 2 π}, {t, 0, 
        2 π}, PlotRange -> All, 
      PlotPoints  ->  50]

gr1= %

Thϵ torus as submanifold of (S^3) is locally flat (local euclidean). Therefore its geodesics are the isogonal trajectories of its generating circles. Their images under the stereographic projection are:

kron[a_,  c_][t_] := 
        sterproj[{Cos[t],  Sin[t],  0,  0}*
              c  +  
{0,  0,  Cos[a t],  Sin[a t]} 
              Sqrt[1  -  c^2]]

kron[a,  c][t]

Since stereographic projection is not an isometry, the curves kron[a,c][t] are in general not geodesics on the image torus in the 3-space.

kron[ .5,  0.8][t]

ParametricPlot3D[kron[ .5,  .8][t], {t, 0, 10 π}, 
      PlotRange -> All, PlotPoints -> 100]

gr2 = %

Show[{gr1, gr2}]

gr3 = %

If you want to show or to animate this Graphics3D outside of Mathematica, use ThreeScript (see section 7.1), which has been imported with the Initialization.

?ThreeScript

ThreeScript["torexam.ts", gr3]

Unfortunately, ParametricPlot3D doesn't admit thickening the plot of the curve:

Options[ParametricPlot3D]

But with wire frame geometry the geodesic is well visible looking at the ts-file torexam.ts under MathGL3D or MathLive.  Alternatively, link MathGL3d (see section 7.2) and apply the OpenGLViewer with MVShow3D[gr3, MVLineTubeSize -> 0.02]

MVShow3D[gr3,  MVLineTubeSize  ->  0.02, 
        MVTextureMapType  ->  MVSphereMapping]

Now we show the torus by a family of Meridians

Table[
      ParametricPlot3D[tor[0.8][s, t], {s, 0, 2 π}], {t, 0, 
        2 π,  .1}]

Show[%]

tnet  =  %

Show[{tnet, gr2}]

gr4  = %

MVShow3D[gr4]

Here one can see that the intersection angle along the curve is constant. (Remark: Taking Appearance -> WireFrame needs not to represent the images of the generating circles, as gr4 does.)

PROOF of this statement.

The tangential vector of the curve is

tangkron[a_, c_][t_] := 
      D[kron[a,  c][t1], t1] /. t1 -> t

tangkron[a, c][t]

The tangential vectors of the parameter lines of the torus are

dttor[c_][s_, t_] :=  
      D[tor[c][s, t1], t1] /. t1 -> t

dttor[c][s, t]

dstor[c_][s_, t_] :=  
      D[tor[c][s1, t], s1] /. s1 -> s

dstor[c][s, t]

kron[a, c][u]

tor[c][s, t]

The parameter representation of the curve kron[a,c][u] on the torus is s=a u, t = u:

kron[a, c][u] == tor[c][a*u, u]

We show that the cosinus of the angle in which curve kron[a,c][u] intersects the parameter line is constant along the curve:

testf[t_] :=  
      dotnorme[tangkron[ .5,  .8][t]] . 
        dotnorme[dstor[ .8][  .5*t, t]]

Only numerical disturbations appear:

Plot[testf[t],  {t, 0, 2*Pi}]

Plot[Chop[testf[t] - testf[0]]  +  1,  {t, 0, 2*Pi}]

testf[0]

Clear[testf]

Test of the angle invariance

flattor[c_][s_,  t_] :=\({Cos[t], \ Sin[t], \ 0, \ 0}*
            c\  + 
{0, \ 0, \ Cos[s], \ Sin[s]} Sqrt[1\  - \ c^2]\)

dsflatt[c_][s_, t_]  :=  
      D[flattor[c][s1, t], s1] /. s1 -> s

dsflatt[c][s, t]

dtflatt[c_][s_, t_] :=  
      D[flattor[c][s, t1], t1] /. t1 -> t

dtflatt[c][s, t]

The tangential basis of the euclidean torus has been defined with the initialization:

dstor[c][s, t]

dttor[c][s, t]

dstor[c][s, t] . dstor[c][s, t]

Simplify[
      Together[%],  s > 0 && t > 0 && c > 0 && c < 1]

g11[c_][s_, 
        t_]  := ((-1 + (c^2)) (-3 + 
                (c^2) - (-1 + (c^2)) Cos[2 s] + 
                4 (Sqrt[1 - (c^2)]) Sin[s]))/(2 (1 - 
                4 (Sqrt[1 - (c^2)]) Sin[s] - 
                6 (-1 + (c^2)) Sin[s]^2 - 
                4 (((1 - (c^2)))^(3/2)) Sin[s]^3 + (((-1 + (c^2)))^2) Sin[s]^4))

dstor[c][s, t] . dttor[c][s, t]

dttor[c][s, t] . dttor[c][s, t]

g22[c_][s_, 
        t_]  :=  (((c^2) Cos[t]^2)/(((1 - (Sqrt[1 - (c^2)]) Sin[s]))^2)) + (((c^2) Sin[t]^2)/(((1 - (Sqrt[1 - (c^2)]) Sin[s]))^2))

Dividing g11,g22 by the norm squares (1-c^2), resp. c^2 of the tangent vectors of flattor[c][s,t], and substracting:

Simplify[g11[c][s, t]/(1 - c^2)  -  g22[c][s, 
            t]/(c^2)]

and this proves the (well known) conformality of the map sterproj in this case.

The isogonal trajectories are not geodesics

The last  result confirms the fact that the stereographical projection is a conformal map.  We show that the isogonal trajectories of the parameter lines, in general, are not  geodesics. To this aim we calulate the curvature vector of the curve kron[a,c][t]. The unit tangential vector is

tkron[a_, c_][t_]  :=  dotnorme[tangkron[a, c][t]]

The binormalvector of a regular curve cx[t] of the euclidean 3-space with an arbitrary parameter t is given by (see A. Gray [1]):

binormal[cx_][t_]  :=  
      Simplify[Cross[
D[cx[tt], tt],  D[cx[tt], {tt, 2}]]]/
          Simplify[
              Factor[Cross[D[cx[tt], tt], D[cx[tt], {tt, 2}]] . 
                  Cross[D[cx[tt], tt],  D[cx[tt], {tt, 2}]]]]^(1/2)  /. 
        tt -> t

Now we consider the isogonal trajectory as defined above

cx[t_]  :=  kron[0.5, 0.8][t]

binormal[cx][t]

bincx[t_]  =  Simplify[Chop[%],  t > 0]

If kron[a,c][t] would be a geodesic, the tangential component of this vector should vanish; but the following calculation shows that this is not the case.

ftest[t_] := dstor[ .8][ .5*t, t] . bincx[t]

Plot[ftest[t],  {t, 0, 2*Pi}]

Clear[ftest]

ftest[t_] := dttor[ .8][ .5*t, t] . bincx[t]

Plot[ftest[t],  {t, 0, 2*Pi}]

This shows that the isogonal trajectory cx is not a geodesic.

Clear[ftest]

2.3 The inverse stereographic projection

The following function maps an arbitrary point of the equatorial hyperplane to its preimage under the stereographic projection:

?invsterproj

\(invsterproj[y_List] := \ 
      Module[{n = Length[y] + 1, x = Append[y, 0], \ 
          u = y . y}, 
		x*2/\((1 + u)\)\  + 
          stb[n, n]*\((u - 1)\)/\((u + 1)\)]\)

invsterproj[vec[y, 3]]

Simplify[% . %]

Examples

Simple tests

invsterproj[nul[4]]

Lines are mapped onto circles

invsterproj[{1, t}]

ParametricPlot3D[{(2/(2 + (t^2))), ((2 t)/(2 + (t^2))), 
        ((t^2)/(2 + (t^2)))},  {t,  -100,  100}, 
      PlotRange  ->  All, PlotPoints  ->  1000]

invsterproj[stb[1] Cos[t] + stb[2] Sin[t] + stb[4]]

sterproj[
      invsterproj[stb[1] Cos[t] + stb[2] Sin[t] + stb[4]]]

sterproj[invsterproj[{a, b, c, d}]]

invsterproj[sterproj[%]]

Simplify[%]

Simplify[% /. (a^2) + (b^2) + (c^2)  ->  1 - (d^2)]

We remember the warning

?sterproj

3 Euclidean representations of spheres

3.1 A parameter representation of the unit n-sphere (S^n)

We define parameter representations for the unit vectors in subspaces of the d-dimensional euclidean space:

unit1[s_,  d_:  dim]  :=  
      stb[1, d]*Cos[s] +  stb[2, d]*Sin[s]

unit2[s_,  t_,  d_:  dim]  :=  
      unit1[s,  d]*Cos[t]  + stb[3, d]* Sin[t]

unit3[s_, t_, u_,  d_:  dim]  :=  
      unit2[s, t,  d]*Cos[u]  + stb[4, d]* Sin[u]

3.1.1 A parameter representation of (S^n), recursive Definition

The following recursive definition, contained in euvec.m,  yields the usual parameter representation of the unit n-sphere as a hypersurface in the (n+1)-dimensional euclidean space. For n =1,2,3 it almost coincides with the definitions above. The parameters appear as the first n elements of the List-argument.

\(\(unit[n_]\)[s_List]\  := \ 
      Module[{un, 
            i}, 
		un[1] = {Cos[s[\([1]\)]], \ 
              Sin[s[\([1]\)]]}; 
			For[i = 2, \ 
            i <= \ n, \ \(i++\), 
		un[i] = 
              AppendTo[un[i - 1], 0]*Cos[s[\([i]\)]] + 
                stb[1 + i, i + 1]*Sin[s[\([i]\)]]]; 
				Return[un[n]]] /; 
        Length[s] >= n\)

Examples

unit[2][{s, t}]

unit1,uni2,unit3 yield vectors of the fixed dimension d, completed with zeros, if necessary:

unit2[s, t]

unit[2][{s, t}]

TrueQ[unit2[s, t, 3]  == unit[2][{s, t}]]

unit[3][{s, t, u}]

unit3[s, t, u]

unit[4][vec[s]]

If Length[s_List] > n only the first n components are considered:

vec[s]

unit[2][vec[s]]

If Length[s_List] < n the evaluation failes:

unit[5][vec[s]]

3.2 The standard parameter representation of spheres in the euclidean 3-space

The general sphere with radius r and centre  m has the parameter representation with geographical coordinates u, v, defined in euvec.m:

?sphere

sphere[m_List, r_][u_, v_] := {m[[1]] + r*Cos[u]*Cos[v], 
    m[[2]] + r*Cos[v]*Sin[u], m[[3]] + r*Sin[v]} /; Length[m] == 3

Of course, the 2-spheres in the 3-space can directly be plotted; the sphere with centre m and radius r (optionally a part of it)  is

?sphereplot3D

sphereplot3D[m_List,  r_,  u0_:  0,  u1_:  2 Pi, v0_: - Pi/2, 
        v1_:  Pi/2,  opt___]  :=  
      ParametricPlot3D[sphere[m, r][u, v], {u, u0, u1}, {v, v0, v1}, 
        opt]

Using unit[n][s], one easily generalizes the definition of the parameter representation to the case of arbitrary n-spheres:

\(\(hypersphere[m_List, \ r_]\)[s_List]\  := \ 
      Module[{di\  = \ Length[m] - 1}, 
		m\  + \ r*\(unit[di]\)[s]] /; 
        Length[s] == Length[m] - 1\)

Examples

sphere[nul[3], 1][u, v]

sphereplot3D[nul[3], 1]

sphereplot3D[{0, 0, 5},  2]

sphereplot3D[{0, 0, 5},  2, 0,  Pi/2]

sphereplot3D[{0, 0, 5},  2, 0, Pi/2, 0, Pi/2,  Boxed  ->  False, 
      Axes  ->  False,  PlotPoints  ->  30]

hypersphere[{m1, m2}, r][{s}]

hypersphere[{m1, m2, m3}, r][{u, v}]

hypersphere[stb[1, 4],  r][{u, v, w}]

3.3 Parameter representations of spheres in the 3-sphere (S^3)

orthobase

The parameter representation unit3 of (S^3) is a vector function with values  in the 4-space; geometric objects in (S^( 3))can be plotted only after stereographic projection onto the 3-space. Each sphere (= hypersphere of (S^( 3))) is the intersection of (S^( 3)) with a hyperplane

vec[n] . vec[x]  ==  p

with  normal vector  vec[n] ϵ  (S^( 3)), and distance p, -1 < p < 1, from the center. The ambiguity in the signs of vec[n] and p can be used to distinguish the orientations of the spheres. The manifold of spheres in (S^( 3)) is 4-dimensional: three parameters of vec[n], and p ((n+1-dimensional in case of (S^( n))), or centre and radius. To get a parameter representation of the sphere corresponding to vec[n] and p we need three vectors orthogonal to vec[n]. Applying orthobase we shall get them.

?orthobase

orthobase[vn_]  :=  
      Module[{d  =  Length[vn],  tt , to,  
          v},  
tt  =  
          IdentityMatrix[d]; 
PrependTo[tt,  
          vn]; 
to  =  
          esorthonorm[tt]; 
v  =  
          to[[1]]; 
tt = 
          Take[to,  -(d - 1)]; 
AppendTo[tt,  
          v]; 
Return[tt]]

The Erhard Schmidt orthogonalization procedure esorthonorm is defined in the package euvec.m, see also section 4.

Simple cases

orthobase[stb[1]]

Det[%]

orthobase[stb[4]]

Det[%]

orthobase[stb[2]*a  + stb[4]*b]

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

Simplify[Det[%]]

orthobase[nul[]]

orthobase[stb[2, 7]]

Symbolic calculation

vec[v]

of  = orthobase[vec[v]]

Simplify[%, v[2]^2 + v[3]^2 + v[4]^2  == 1 - v[1]^2 && 
              v[1] > 0 && v[2] > 0 && v[3] > 0 && v[4] > 0]

Control of the scalar products:

Simplify[Table[of[[i]] . of[[j]], {i, 4}, {j, 4}]]

or

Simplify[of . Transpose[of]]

Det[of]

Simplify[%, v[1] > 0 && v[2] > 0 && v[3] > 0 && v[4] > 0]

Direct Definition

Using the formula for orthobase[vec[n]] as definition one could try defining orthobase directly, as a function of  vec[n].

orthobasis4[wx_List] := 
      Module[{c = wx . wx,  vx}, 
If[
          Length[wx]  ≠  4,  
          Print[wx,  " Incorrect argument!"]; 
          Return[]]; 
If[c == 0, 
          Return[IdentityMatrix[4]]]; 
If[c == 1, 
          vx  =  wx, vx  =  wx/Sqrt[c], 
          vx  =  wx/Sqrt[c]]; 
{{Sqrt[
              vx[[2]]^2 + vx[[3]]^2 + 
                vx[[4]]^2], -((vx[[1]]*vx[[2]])/
                  Sqrt[vx[[2]]^2 + vx[[3]]^2 + 
                      vx[[4]]^2]), -((vx[[1]]*
                      vx[[3]])/
                  Sqrt[vx[[2]]^2 + vx[[3]]^2 + 
                      vx[[4]]^2]), -((vx[[1]]*
                      vx[[4]])/
                  Sqrt[vx[[2]]^2 + vx[[3]]^2 + 
                      vx[[4]]^2])}, {0, 
            Sqrt[(vx[[3]]^2 + vx[[4]]^2)/(vx[[2]]^2 + 
                    vx[[3]]^2 + 
                    vx[[4]]^2)], -((vx[[2]]*vx[[3]]*
                      Sqrt[(vx[[3]]^2 + 
                              vx[[4]]^2)/(vx[[2]]^2 + 
                              vx[[3]]^2 + 
                              vx[[4]]^2)])/(vx[[3]]^2 + 
                      vx[[4]]^2)), -((vx[[2]]*vx[[4]]*
                      Sqrt[(vx[[3]]^2 + 
                              vx[[4]]^2)/(vx[[2]]^2 + 
                              vx[[3]]^2 + 
                              vx[[4]]^2)])/(vx[[3]]^2 + 
                      vx[[4]]^2))}, {0, 0, 
            Sqrt[vx[[4]]^2/(vx[[3]]^2 + 
                    vx[[4]]^2)], -((vx[[3]]*
                      Sqrt[vx[[4]]^2/(vx[[3]]^2 + 
                              vx[[4]]^2)])/
                  vx[[4]])}, {vx[[1]], vx[[2]], vx[[3]], 
            vx[[4]]}}]

WARNING! Doing so, for vectors with vx[[4]] == 0 we obtain indeterminate expressions:

orthobasis4[stb[1] +  stb[2]]

orthobasis4[stb[1] +  stb[2] +  stb[3]]

orthobasis4[stb[1] +  stb[2] +  stb[4]]

orthobasis4[{a, b, c, 0}]

Clear[a, b, c]

orthobasis4[{a, b, c, x}]

We can determine the basis as  a limit:

Limit[%,  x  ->  0]

Thus, we try to refine the definition as follows:

Clear[orthobasis4]

orthobasis4[wx_List] := 
      Module[{c = wx . wx,  vx,  ex  =  0,  ob,  main, sing, 
          t}, 
If[Length[wx]  ≠  4,  
          Print[wx,  " Incorrect argument!"]; 
          Return[]]; 
If[c == 0, 
          Return[IdentityMatrix[4]]]; 
If[c == 1, 
          vx  =  wx, vx  =  wx/Sqrt[c], 
          vx  =  wx/Sqrt[c]]; 
If[
          vx[[4]]  ==  0,  ex  =  1;  
          Goto[sing] ]; 
Label[
          main]; 
ob  = 
{{Sqrt[
                vx[[2]]^2 + vx[[3]]^2 + 
                  vx[[4]]^2], -((vx[[1]]*vx[[2]])/
                    Sqrt[vx[[2]]^2 + vx[[3]]^2 + 
                        vx[[4]]^2]), -((vx[[1]]*
                        vx[[3]])/
                    Sqrt[vx[[2]]^2 + vx[[3]]^2 + 
                        vx[[4]]^2]), -((vx[[1]]*
                        vx[[4]])/
                    Sqrt[vx[[2]]^2 + vx[[3]]^2 + 
                        vx[[4]]^2])}, {0, 
              Sqrt[(vx[[3]]^2 + vx[[4]]^2)/(vx[[2]]^2 + 
                      vx[[3]]^2 + 
                      vx[[4]]^2)], -((vx[[2]]*vx[[3]]*
                        Sqrt[(vx[[3]]^2 + 
                                vx[[4]]^2)/(vx[[2]]^2 + 
                                vx[[3]]^2 + 
                                vx[[4]]^2)])/(vx[[3]]^2 + 
                        vx[[4]]^2)), -((vx[[2]]*
                        vx[[4]]*
                        Sqrt[(vx[[3]]^2 + 
                                vx[[4]]^2)/(vx[[2]]^2 + 
                                vx[[3]]^2 + 
                                vx[[4]]^2)])/(vx[[3]]^2 + 
                        vx[[4]]^2))}, {0, 0, 
              Sqrt[vx[[4]]^2/(vx[[3]]^2 + 
                      vx[[4]]^2)], -((vx[[3]]*
                        Sqrt[vx[[4]]^2/(vx[[3]]^2 + 
                                vx[[4]]^2)])/
                    vx[[4]])}, {vx[[1]], vx[[2]], 
              vx[[3]], vx[[4]]}}; 
If[
          ex == 0,  Return[ob]]; 
If[
          ex == 1,  
          Return[Limit[ob,  t  ->  0]]]; 
Label[
          sing]; 
vx  = 
          vx /. vx[[4]] ->  t; 
Goto[
          main];
]

orthobasis4[stb[1]]

Table[%[[i]] . %[[j]],  {i, 4}, {j, 4}] // 
      MatrixForm

This is not the simplest solution:

orthobase[stb[1]]

With the new definition we do not get indeterminate expressions:

orthobasis4[stb[1] +  stb[2]]

orthobasis4[stb[1] +  stb[2] +  stb[3]]

orthobasis4[stb[1] +  stb[2] +  stb[4]]

orthobasis4[{a, b, c, 0}]

Clear[a, b, c]

orthobasis4[{a, b, c, x}]

Timings:

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

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

Very often orthobasis4 is much faster than orthobase:

Timing[orthobase[vec[vv]]]

Timing[orthobasis4[vec[vv]]]

Simplify[Together[%]]

Simplify[
      Table[%[[2]][[i]] . %[[2]][[j]], {i, 4}, {j, 
          4}]]

Remove[of]

Sphere given by vec[n], p

Let vec[n], p be as at the beginning of  3.3. Obviously, vec[n]*p is the center and Sqrt[1-p^2] is the radius of the sphere in the hyperplane vec[n].x=p.  A parameter representation of this sphere as a surface in the 4-space is given by subsphere[n,p][v,w].

?subsphere

subsphere[stb[1]*Cos[t]  + stb[2]*Sin[t],  1/6][v, w]

(% - (stb[1]*Cos[t]  + stb[2]*Sin[t])/
            6) . (% - (stb[1]*Cos[t]  + stb[2]*Sin[t])/
            6)

Simplify[Together[%]]

Error message :

subsphere[stb[5, 7], 0][v, w]

The next command gives the unit sphere with center in the origin in the hyperplane spanned by stb[2], stb[3], stb[4]; the same gives the next command.

subsphere[nul[], 0][v, w]

subsphere[stb[1], 0][v, w]

subsphere[stb[1], p][v, w]

Simplify[(% - stb[1]*p) . (% - stb[1]*p)]

subsphere[stb[1], 0.5][v, w]

subsphere[stb[1]*n[1] +  stb[2]*n[2], p][v, w]

subsphere[{1, 1, 0, 0}/Sqrt[2], 1/2][v, w]

subsphere[{1, 0, 0, 0}, 1/2][v, w]

subsphere[stb[4],   .9][v, w]

Timing[subsphere[stb[2]*Cos[t] + stb[3]*Sin[t], Random[]][v, 
        w]]

Timing[subsphere[stb[2]*Cos[t] + stb[3]*Sin[t], 0.8][v, 
        w]]

Symbolic calculations

subsphere[vec[n],  p][v, w]  /. 
      vec[n] . vec[n]  ->  1

Simplify[%,  n[1] > 0 && n[2] > 0 && n[3] > 0 && n[4] > 0]

Tangential framing

At first we only analyze the situation and find formulas, which  then will be used, in a simplified way, to define the frame explicitely as a function of the geographical coordinates s,t,u. The frame is a natural tangential frame at unit3[s,t,u]. We  derive  the parameter representation of (S^(3 )), and norm the tangential vectors. We use

dotnorme[b]

dotnorme[vec[n]]

unit3[s,  t, u]

tg1[s_, t_, u_]  :=  D[unit3[x,  t, u], x] /. x  ->  s

tg1[s, t, u]

% . %

ntg1[s_, t_, u_] :=  dotnorme[tg1[s, t, u]]

ntg1[s, t, u]

Simplify[%,  
      s ∈ Reals && t ∈ Reals && 
        u ∈ Reals]

Since t, u fulfil the inequality -Pi/2< t, u < Pi/2 we may set

Clear[ntg1];
ntg1[s_, t_, u_] :=  {-Sin[s],  Cos[s], 0, 0}

tg2[s_, t_, u_]  :=  D[unit3[s,  x, u], x] /. x  ->  t

ntg2[s_, t_, u_] :=  dotnorme[tg2[s, t, u]]

ntg2[s, t, u]

Simplify[%,  
      s ∈ Reals && t ∈ Reals && 
        u ∈ Reals]

Since u fulfils the inequality -Pi/2< u < Pi/2 we may set

Clear[ntg2];
ntg2[s_, t_, u_] :=  {-Cos[s]*Sin[t],  -Sin[s]*Sin[t],  
        Cos[t], 0}

tg3[s_, t_, u_]  :=  D[unit3[s,  t, x], x] /. x  ->  u

ntg3[s_, t_, u_] :=  dotnorme[tg3[s, t, u]]

Simplify[ntg3[s, t, u],  
      s ∈ Reals && t ∈ Reals && 
        u ∈ Reals]

We set

Clear[ntg3];
ntg3[s_, t_, 
        u_] :=  {-Cos[s] Cos[t] Sin[u], -Cos[t] Sin[s] Sin[
            u], -Sin[t] Sin[u], Cos[u]}

Clear[tgframe]

tgframe[s_, t_, u_]  :=  {ntg1[s, t, u], ntg2[s, t, u], 
        ntg3[s, t, u], unit3[s, t, u]}

tgframe[s, t, u] // MatrixForm

Having this result we can formulate the direct definition (remember 0<=s< 2Pi, -Pi/2<= t,u <= Pi/2).

This definition  contained in the package euvec.m is:

?tg4frame

MatrixForm[tg4frame[s, t, u]]

tg4frame yields a negatively oriented frame:

Simplify[Det[tg4frame[s, t, u]]]

Another parameter representation of the spheres, using the tangential frames

The tangential framing gives the possibility of the following parameter representation of arbitrary spheres, defined in euvec.m:

?subspheremf

We called this procedure subspheremf since the argumants s,t,u,p may serve as coordinates on the 4-dimensional manifold of all subspheres of (S^3).  The procedure is much faster, and more effectiv,  than the composition of subsphere with unit[3]:

Timing[subspheremf[s, t, u,  p][v, w]]

Timing[subsphere[unit[3][{s, t, u}], p][v, w]]

Examples

unit3[0, 0, 0]

For Abs[p] > 1 we get an imaginary sphere

subspheremf[0, -Pi/2, 0, 2][v, w]

and for Abs[p] = 1 we get a point only:

subspheremf[0, -Pi/2, 0, -1][v, w]

unit3[0, t, Pi/2]

unit3[s, t, Pi/2]

subspheremf[s, t, Pi/2, p][v, w]

We may plot the sphere in the intersecting hyperplane by e.g.

Timing[
      ParametricPlot3D[
        Take[subspheremf[0.4, 0.3, Pi/2, 0.9][v, w], 3],  {v,  0,  
          2 Pi}, {w,  -Pi/2, Pi/2}]]

a = Random[]; b = Random[]; c = Random[]; 
    Print[{a, b, c}];
ParametricPlot3D[
      Take[subspheremf[a, b, Pi/2, c][v, w], 3],  {v,  0,  
        2 Pi}, {w,  -Pi/2, Pi/2}]

Now we apply the stereographical projection to plot the corresponding sphere in the euclidean 3-space:

a = Random[]; b = Random[]; c = Random[]; 
    Print[{a, b, c}];
ParametricPlot3D[
      sterproj[subspheremf[a, b, c, 0.3][v, w]],  {v,  0,  
        2 Pi}, {w,  -Pi/2, Pi/2}, PlotRange  ->  All, 
      PlotPoints  -> 60]

Clear[a, b, c]

hyperplane[Table[stb[i],  {i, 4}]]

It follows that the north pole stb[4] of the standard unit hypersphere (S^3) belongs to the sphere

subsphere[{1, 1, 1, 1}/2, 1/2][v, w]

ParametricPlot3D[
      sterproj[subsphere[{1, 1, 1, 1}/2, 1/2][v, w]],  {v,  0,  
        2 Pi},  {w, -Pi/2,  Pi/2}]

Since the sphere containes the north pole of the hypersphere (S^3) we got by stereographic projection sterproj  a plane of the euclidean 3-space. The next example is the standard unit sphere in the 3-space (the equatorial sphere of (S^3).

subsphere[{0, 0, 0, 1}, 0][v, w]

sterproj[subsphere[{0, 0, 0, 1}, 0][v, w]]

ParametricPlot3D[
      sterproj[subsphere[{0, 0, 0, 1}, 0][v, w]],  {v,  0,  
        2 Pi},  {w, -Pi/2,  Pi/2}]

sterproj[subspheremf[s, t, u,  p][v, w] ]

Now we create a configuration of spheres:

\(Table[
      ParametricPlot3D[
        sterproj[\(subspheremf[i\ *2*Pi/5, 0, 0, \ 0.9]\)[v, 
            w]], 
			{v, \ 0, \ 2  Pi}, \ {w, \(-Pi\)/2, Pi/2}], \ {i, 
        5}]\)

Show[%]

sph5v  =  %

MVShow3D[sph5v]

ThreeScript["sph5.ts",  sph5v]

3.4 The sphere through four points

Hyperplanes through four points

Four points of the euclidean 4-space in general position define a hyperplane hp[n,p], n, p as in  3.3, namely

?hyperplane

\(hyperplane[pt_List]\  := \ 
      Module[{n, q, p, d = Length[pt]}, 
		q\  = \ 
          Table[pt[\([i]\)]\  - \ pt[\([d]\)], \ {i, 
              d\  - \ 1}]; 
			n\  = \ cross[q]; 
	
        If[Length[n]\  != \ d, \ Goto[fin]]; 
		If[n\  === \ nul[], 
          Return[\ Print["The points ", \ 
              pt, \ " are not in general position!"]]]; 
n\  = \ n/Sqrt[n . n]; 		
p\  = \ 
          n . pt[\([1]\)]; 
Return[{n, \ p}]; 
	Label[fin]; 

        Return[]]\)

Remark: The definition "hyperplane" is valid in any dimension.

Tests

For dim = 4, in case -1 < p < 1, the pair {n,p} defines the sphere with parameter representation subsphere[n, p][v,w].

?subsphere

hyperplane[{stb[1],  stb[2],  0.9*stb[3],  stb[4]}]

subsphere[%[[1]],  %[[2]]][v, w]

ParametricPlot3D[
      sterproj[%],  {v, 0, 2 Pi}, {w,  -Pi/2,  Pi/2}]

hyperplane[{stb[1],  stb[2],  stb[3],  -stb[3]}]

subsphere[%[[1]],  %[[2]]][v, w]

ParametricPlot3D[
      sterproj[%],  {v, 0, 2 Pi}, {w,  -Pi/2,  Pi/2}]

The next two examples return only messages...

hyperplane[{stb[1],  stb[2],  stb[2],  -stb[3]}]

hyperplane[{stb[1, 5],  stb[2, 5],  
        stb[3, 5],  -stb[3, 5]}]

hyperplane[{stb[1, 5],  stb[2, 5],  stb[3, 5],  -stb[3, 5], 
        stb[5, 5]}]

Sphere through four points via hyperplane

The following Module combines  stereographical projection, hyperplane and subsphere to get a parameter representation of the sphere through four points of the euclidean 3-space:

\(\(sphere4pts[pt_List]\)[v_, w_] := 
      Module[{a, 
          b}, 
			If[\((Length[pt]\  != 
                4)\)\  || \((Length[pt[\([1]\)]]\  != \ 3)\), \ 
          Print["Incorrect arguments!!!"]; 	Return[]]; 
			a\  = \ 
          Table[invsterproj[pt[\([i]\)]], \ {i, 4}]; 
			b\  = \ 
          hyperplane[a]; 
		If[b == \ Null, \ 	Return[]]; 
			Return[
          sterproj[\(subsphere[b[\([1]\)], \ b[\([2]\)]]\)[v, 
              w]]]
]\)

This gives, of course, the standard unit sphere:

sphere4pts[{stb[1,  3],  stb[2, 3], stb[3, 3],  -stb[3, 3]}][
      v, w]

\(ParametricPlot3D[%, {v, \ 0, \ 2\ Pi}, \ {w, \ \(-Pi\)/2, \ Pi/2}, \ 
      PlotRange\  -> \ All, \ 
PlotPoints\  -> \ 50]\)

ptslist =  {stb[1,  3],  stb[2, 3], stb[3, 3],  
        2  stb[3, 3]}

pts  = Map[Point, ptslist]

gr4pt  =   Graphics3D[Join[{PointSize[0.04]}, pts]]

Show[gr4pt]

grpts  =  %

sphere4pts[ptslist][v, w]

\(ParametricPlot3D[%, {v, \ 0, \ 2\ Pi}, \ {w, \ \(-Pi\)/2, \ Pi/2}, \ 
      PlotRange\  -> \ All, \ 
PlotPoints\  -> \ 40]\)

sph4pt  =  %

Show[{sph4pt, grpts}]

t4pt  =  %

MVShow3D[t4pt, MVPointSphereSize  ->  0.8]

Apparently, by numerical reasons, one point seems not to lay on the sphere.

The next three examples return only messages...

sphere4pts[{stb[1,  3],  stb[2, 3], stb[3, 3],  stb[3, 3]}][v, 
      w]

sphere4pts[{stb[1],  stb[3], stb[3],  -stb[3]}][v, 
      w]

sphere4pts[{stb[1],  stb[3], stb[3]}][v, w]

If wanted one can insert sphere4pts[] in a plot command and get the picture of the sphere through four points by

\(pt4sphere[pt_List, opt___]\  := 
      ParametricPlot3D[
        Evaluate[\(sphere4pts[pt]\)[v, w]], 
		{v, \ 0, \ 
          2\ Pi}, {w, \ \(-Pi\)/2, \ Pi/2}, \ opt]\)

\(pt4sphere[{{0.5, \ 1, \ 2. }, \ {\(-1\), \ 1.5, \ 0.7}, \ {2. , \ 
          1. , \ 0. }, \ stb[3, 3]}, 
PlotPoints\  -> \ 50, 
      PlotRange\  -> \ All]\)

ttgr  =  %

pts  = 
      Map[Point, {{0.5,  1,  2. },  {-1,  1.5,  0.7},  {2. ,  
            1. ,  0. },  stb[3, 3]}]

gr4ptt  =   Graphics3D[Join[{PointSize[0.04]}, pts]]

Show[gr4ptt]

grptts  =  %

Show[{ttgr, grptts}]

testgrtt  =  %

MVShow3D[testgrtt, MVPointSphereSize  ->  1.4]

A direct approach

Another, maybe simpler way to plot the sphere through four points in the 3-space is first to calculate its centre and radius and then to apply sphere.

\(centre[pt_List] := 
      Module[{p = Length[pt], d, i, j, k, ab, at, mat, mb, \ c}, 
		If[
          p\  != \ 4, \ Print["Incorrect argument!"]; \ 
          Goto[fin]]; 
			For[i = 1, \ i < 4, \ \(i++\), 
          at[i] = pt[\([i + 1]\)] - pt[\([i]\)]; 
				ab[
              i] = \((\((pt[\([i + 1]\)] . 
                      pt[\([i + 1]\)])\) - \((pt[\([i]\)] . pt[\([i]\)])\))\)/
              2]; 
				mat\  = \ Table[at[j], {j, 1, 3}]; 
		d\  = \ 
          Det[mat]; 
		If[d == \ 0. , \ 
          Print["Points lie on a plane!"]; 
Goto[
            fin]]; 		
		mb\  = \ 
          Table[ab[k], {k, 1, 3}]; 
				c\  = \ 
          LinearSolve[mat, mb]; 
	Return[c]; 
	Label[fin]; 
	
        Return[]]\)

radius[pt_List] :=  
      Module[{c  = centre[pt]} , 
If[c ==  Null,  
          Return[]]; 
Sqrt[(pt[[1]] - 
                c) . (pt[[1]] - c)]]

\(sphere4ptsplot3D[p_List, s0_:  0, \ s1_:  2\ Pi, t0_: - Pi/2, \ 
        t1_:  Pi/2, \ opt___]\  := 
      Module[{m, r}, 
m\  = \ 
          centre[p]; 
If[m\  == Null, \ 
          Return[]]; 
r\  = \ radius[p]; 
	
        sphereplot3D[m, r, s0, \ s1, t0, \ t1, \ opt]]\)

Examples, tests

sqr  = {{0.5,  1,  2. },  {-1,  1.5,  0.7},  {2. ,  1. ,  
          0. },  stb[3, 3]}

centre[sqr]

cpt  =  %

centre[sqr]

sqr[[1]]

Table[(cpt - sqr[[j]]) . (cpt - sqr[[j]]),  {j, 1, 
        4}]

radius[sqr]

radius[sqr]

%^2

sphere4ptsplot3D[sqr]

testgrtt  =  %

Show[{testgrtt, grptts}]

testgrtt1  =  %

MVShow3D[testgrtt1, MVPointSphereSize  ->  1.4]

Now the sphere fits much better to the points!

p2  = {nul[3], stb[1, 3], stb[2, 3],  stb[3, 3]}

centre[p2]

cp2  =  %

Table[(cp2 - p2[[i]]) . (cp2 - p2[[i]]),  {i, 1, 
        4}]

radius[p2]

sphere4ptsplot3D[p2,  0,  Pi,  -Pi/2, Pi/2,  
      PlotPoints ->  40]

Error messages

sphere4ptsplot3D[Take[p2, 2]]

p3  = {nul[3], stb[1, 3], stb[2, 3],  2*stb[2, 3]}

centre[p3]

radius[p3]

sphere4ptsplot3D[p3]

4 Orthogonality

4.1 Erhard Schmidt's orthogonalization

The procedure

esorthonorm follows Erhard Schmidt's orthogonalization, in difference to which the given vector sequence m_List needs not to be linearly independent. If a vector linearly depends on the foregoing vectors, esorthonorm generates the zero vector and eliminates it.  The module esorthonorm is contained in the package euvec.m.The numerical behavior of esorthonorm can be influenced by the option neglect, which sets the parameter of smoothing. 
WARNING: Because of the partly numerical character of the procedure, an inappropriate setting of the option neglect may cause false results. One always should control the result of esorthonorm checking the orthogonality and completeness of the produced sequence!

Options[esorthonorm]

?smoothing

?neglect

?print

?esorthonorm

\(\(\(esorthonorm[m_List, 
        opts___]\)\(\ \)\(:=\)\(\ \)\(
\)\(\ \ \)\(Module[
		{ma = \ m, 
          mb = {}, k\  = \ Length[m], rd = Length[m], \ j\  = \ 0, \ i, ne, 
          prod, \ pr, \ bn, \ sm, aus, \ h, \ 
          f}, 
		{prod, nor, pr, \ 
            bn}\  = \ 
			\({innerprod, normed, \ print, \ 
                neglect} /. {opts}\) /. 
            Options[esorthonorm]; 
If[
          Head[m[\([1]\)]\ ] == \ List, 
          rd\  = \ Length[m[\([1]\)]]]; 
	
        f[h_]\  := \ smoothing[h, 10^bn]; 
	
        sm[aus_]\  := \ MapAll[f, aus]; 
		Label[
          ne]; 
If[
          Length[mb] ≥ \ Min[rd, k], \ 
          Return[mb]]; 
		j\  = \ j + 1; 
If[j > k, 
          Return[mb]]; 
			ma[\([j]\)] = 
          sm[Simplify[
              ma[\([j]\)]\  - 
	
                Sum[mb[\([i]\)]*
                    prod[ma[\([j]\)], \ mb[\([i]\)]]/
                      prod[mb[\([i]\)], \ mb[\([i]\)]], \ 
	{i, 1, \ 
                    Length[mb]}]]]; \ 		

        If[pr\  == True, Print[ma[\([j]\)]]]; 

        If[prod[ma[\([j]\)], \ ma[\([j]\)]] == \ 0, 
          Goto[ne], 
\ \ \ \ \ \ \ \ If[nor\  == \ True, 
            AppendTo[mb, 
              sm[Simplify[
                  ma[\([j]\)]/
                    Sqrt[prod[ma[\([j]\)], 
                        ma[\([j]\)]]]]]], 
				AppendTo[mb, 
              Simplify[ma[\([j]\)]]]], 
				If[
            nor\  == True, 
				AppendTo[mb, 
              sm[Simplify[
                  ma[\([j]\)]/
                    Sqrt[prod[ma[\([j]\)], ma[\([j]\)]]]]]], \ 			

            AppendTo[mb, Simplify[ma[\([j]\)]]]]]; 		
		Goto[
          ne]]\)\(
\)\)\)

Simple tests

First we give some simple examples.

esorthonorm[{nul[]}]

esorthonorm[Table[stb[i]*i,  {i, 4}],  print  ->  True]

esorthonorm[Table[stb[i]*(i - 1),  {i, 4}], normed  ->  False, 
      print  ->  True]

Linearly dependent vectors are omitted:

at  = 
      AppendTo[Table[randomv[]*(j - 2), {j, 7}],  nul[]]

neglect -> 0 amd neglect -> default give the same results:

ot1 = esorthonorm[at, print  ->  True, neglect  ->  0]

Table[
        smoothing[ot1[[i]] . ot1[[j]],  1], {i, Length[ot1]},  {j, 
          Length[ot1]}] // MatrixForm

Chop[%]

ot2 = esorthonorm[at, print  ->  True]

ot2 - ot1

Also with normed -> False (and neglect -> 0) one gets the same results:

at  = 
      AppendTo[ Table[randomv[]*(j - 2), {j, 7}],  nul[]]

ot3 = 
      esorthonorm[at, normed  ->  False,  print  ->  True, 
        neglect  ->  0]

Table[
        smoothing[ot3[[i]] . ot3[[j]],  1], {i, Length[ot3]},  {j, 
          Length[ot3]}] // MatrixForm

Chop[%]

Now again with the default value of neglect:

ot4 = 
      esorthonorm[at, normed  ->  False,  print  ->  True]

ot3 - ot4

The Legendre polynomials with esorthonorm

The procedure esorthonorm  may be applied also for function spaces, for which a positive semidefinite scalar product must be defined in terms of Mathematica.

prodint[f_, g_]  :=  Integrate[f*g, {x,  -1, 1}]

prodint[1, x^2]

seq[h_]  :=  Table[x^i, {i, 0, h}]

seq[3]

esorthonorm[seq[1], innerprod  ->  prodint]

esorthonorm[seq[2],  innerprod  ->  prodint]

esorthonorm[seq[3],  innerprod  ->  prodint]

esorthonorm[seq[4], innerprod  ->  prodint]

esorthonorm[seq[5], innerprod  ->  prodint]

esorthonorm[seq[6],  innerprod  ->  prodint]

esorthonorm[seq[7],  innerprod  ->  prodint]

Timing[
      esorthonorm[seq[20],  innerprod  ->  prodint,  
        print  ->  True]]

Length[%[[2]]]

eslegpol20  =  %%[[2]];

Table[
        prodint[eslegpol20[[i]], eslegpol20[[j]]], {i, 21}, {j, 21}] // 
      MatrixForm

Comparison with the implemented LegendreP

?LegendreP

The functions LegendreP implemented in Mathematica are normed by the condition

LegendreP[_, 1]

Thus they are orthogonal but not normed with respect to prodint:

legpol20 = 
        Table[LegendreP[j, x],  { j, 0, 20}];

Table[prodint[legpol20[[i]], legpol20[[j]]], {i, 21}, {j, 21}] // 
      MatrixForm

LegendreP[20, 1]

eslegpol20[[21]] /. x -> 1

Up to the norming they coincide with the orthonormal polynomials obtained by esorthonorm:

Table[
      LegendreP[i, x] - 
        eslegpol20[[i + 1]]/(eslegpol20[[i + 1]] /. 
              x -> 1),  {i, 0, 20}]

Simplify[%]

A last example

Sometimes sequences containing numerical and functional vectors, cause difficulties. The next example runs normally

v1  =  randomv[3]

Timing[
      esorthonorm[
        Join[{{0, Cos[t], Sin[t]}, v1}, Table[stb[i, 3],  {i, 2}]],  
        print  ->  True,  normed  ->  False]]

ov  = Chop[%[[2]]]

Table[
      Simplify[ov[[i]] . ov[[j]]], {i, Length[ov]}, {j, 
        Length[ov]}]

Chop[%]

prodsov[t_] = %;

prodsov[Random[]]

Remember that we used the option normed -> False. Why we first have got small imaginary parts?

v1  =  randomv[3]

Timing[
      esorthonorm[
        Join[{{0, Cos[t], Sin[t]}, v1}, Table[stb[i, 3],  {i, 2}]],  
        print  ->  False]]

abc[t_]  =  %[[2]];

Chop[abc[t]]

Simplify[%]

Length[abc[t]]

Table[
      Dot[abc[t][[i]], abc[t][[j]]], {i,  3}, {j,  
        3}]

Chop[%]

prods[t_]  =  %;

prods[Random[]]

Chop[%]

Chop[prods[6*Random[]]]

Remove[ov, prods, abc]

The same with GramSchmidt, a procedure defined in

<< "LinearAlgebra/Orthogonalization.m"

?GramSchmidt

Options[GramSchmidt]

We recommend the user to load this package, and to compare GramSchmidt with esorthonorm, running the same tests including the timing for both procedures, e.g.

v1  =  randomv[3]

Timing[
      GramSchmidt[
        Join[{{0, Cos[t], Sin[t]}, v1}, 
          Table[stb[i, 3],  {i, 2}]]]]

Simplify[%[[2]]]

ov1[t_] =  %;

Chop[ov1[t]]

ov2[t_] = %;

Length[%]

The last vector is zero.

prods[t_] := 
      Table[Simplify[ov2[t][[i]] . ov2[t][[j]]], {i, 3}, {j, 
          3}]

prods[t]

Chop[%]

(% /. t -> Random[]) // MatrixForm

Chop[prods[Random[]]]

4.2 The Euclidean orthogonal complement in the 4-space

For later applications in e.g. Möbius geometry:
To construct the circle which corresponds to the intersection of two hyperplanes in the 4-space and the unit hypersphere we need an orthonormal solution of
v.x == 0, w.x == 0.
This is given by orthoerg defined in euvec.m:

?orthoerg

Examples

orthoerg[stb[3], stb[1]]

orthoerg[{0, Cos[t], Sin[t], 0}, {1, 0, 0, 0}]

Simplify[%,  t > 0]

Clear[v1, v2]

v1  = {Cos[t], 0, Sin[t], 0}; v2 = {1, 0, 0, 1};

orthoerg[v1, v2]

{v1 . %[[1]],  v2 . %[[2]], v1 . %[[2]],  
      v2 . %[[1]],  
      Table[%[[i]] . %[[j]],  {i, 2}, {j, 2}] // 
        MatrixForm}

Simplify[%]

Clear[v1, v2]

v1  =  randomv[]; v2  =  randomv[];

Timing[orthoerg[v1, v2]]

%[[2]]

{v1 . %[[1]],  v2 . %[[2]], v1 . %[[2]],  
      v2 . %[[1]],  
      Table[%[[i]] . %[[j]],  {i, 2}, {j, 2}] // 
        MatrixForm}

MapAll[smoothing, %]

Clear[v1]

v1  =  {0, Cos[t], Sin[t], 0}

v2

orthoerg[v1, v2]

vo  =  %;

Simplify[{v1 . vo[[1]],  v2 . vo[[2]], v1 . vo[[2]],  
        v2 . vo[[1]]}]

Chop[%]

Plot[%[[2]]  +  1, {t, 0, 2 Pi}]

Obviously, also the second term vanishes, but Mathematica does not recognize this.

Table[
      Chop[Simplify[vo[[i]] . vo[[j]]]], {i, 2}, {j, 2}]

Here all is ok. Sometimes the simplification is unsatisfactory, what depends on the random vector. The varry ou the next (or similar) commands

Plot[%[[2, 2]], {t, 0, 2 Pi}]

Plot[Chop[%%[[2, 2]] - 1] + 1, {t, 0, 2 Pi}]

OK, up to  small numerical errors.

5 Reflection at  spheres

Already in 1.3 we considered inversions at arbitrary hyperspheres. The reflection at the unit sphere can be defined much easier:

refl[x_List] := x/(x.x)

5.1 Examples

Tests

refl[{5, 5, 0, 0}]

refl[{2, 2, 2. }]

refl[vec[x]]

Simplify[Together[(% . %)*(vec[x] . vec[x])  - 1]]

Planes and spheres are transformed in planes or spheres

The unit square in the plane z = 1 is transformed into a spherical region:

ParametricPlot3D[
      refl[{u,  v,  1}],  {u,  0,  1},  {v,  0,  1}]

tgr1  =  %

Table[ParametricPlot3D[refl[{u, v, k}], {u, 0, 1}, {v, 0, 1}], 
{k, 1, 6}]

Show[%]

tgr2  = %

The infinite point is mapped in the origin:

Limit[refl[{0, 0, t}],  t ->  Infinity]

6 Spiral surfaces

6.1 The spiral  groups in the 3-Space

In what follows we need the package (automatically supplied by init.m)

Needs["spiralsf`"]

?spiral

The spiral groups are 1-parameter subgroups of the conformal group of the euclidean 3-space depending on a  parameter a:

MatrixForm[spiral[t, a]]

6.2 Fixed Points and Orbits

The only fixed point is the origin:

D[spiral[t, a], t]

m0 = % /. t -> 0

Solve[m0.vec[x,3] == nul[3], vec[x,3]]

The orbit of a point pt is the loxodromic spiral on the circular conus generated by the line through the origin and pt.

\(ParametricPlot3D[
      spiral[t, \ 0.1] . {1, \ 0, \ 1}, \ {t, \ \(-8\)\ Pi, \ 2\ Pi}, 

      PlotPoints\  -> \ 150]\)

conspir =  %

con[t_, u_] := spiral[t,  0.1] . {1,  0,  1} u

ax  =  con[t, u]

ax[[1]]^2 + ax[[2]]^2 - ax[[3]]^2

Simplify[%]

conb[t_, u_] := spiral[t,  0.1] . {1,  0, b} u

bx  =  conb[t, u]

bx[[1]]^2 + bx[[2]]^2 - (bx[[3]]^2)/b^2

Simplify[%]

Here is the corresponding cone:

\(ParametricPlot3D[
      con[t, u], \ {t, \ \(-6\)\ Pi, \ 0}, {u, \ 0, \ 2. }, 

      PlotPoints\  -> \ 30]\)

conus  =  %

Show[{conus, conspir}]

conspir2  =  %

The circlar conus can be created as usual, independent of the spiral:

ParametricPlot3D[{Cos[t],  Sin[t], 1}*u,  {t, 0, 2 Pi}, {u, 0, 
        2. }]

conus2  =  %

Show[{conus2, conspir}]

conspir3  =  %

Apply MathGL3d (or MathLive via ThreeScript) for animating these graphics! (See section 7)

MVShow3D[conspir2,  MVLineTubeSize  ->  0.01]

6.3 The spiral cylinder

A cylinder erected on the spiral is given by

spy[u_, v_, r_:10] := {Cos[u], Sin[u], 1/v} Exp[u/r]

Of course, that is a comparatively simple surface:

ParametricPlot3D[spy[u, v,15], {u, 0, 8 Pi}, {v, 0.5, 5}, 
PlotPoints -> 90]

spc  =  %

An apparently simpler (the usual) parameter representation of a spiral  cylinder  is given by

linspy[u_, v_, r_:10] := {Cos[u], Sin[u], 0} Exp[u/r]+ stb[3,3]*v

Of course, that is a comparatively simple surface:

ParametricPlot3D[linspy[u, v,15], {u, 0, 8 Pi}, {v, 0, 5}, 
PlotPoints -> {110, 20}]

The spiral cylinder is of special interest in Möbius geometry: it is a homogeneous surface. That means: the subgroup of the Möbius group preserving the cylinder acts transitively on this manifold. Moreover, it admits an invariant metric of constant negative curvature. 
For visualizing the behavior at the infinite point, we reflect the spiral cylinder at the unit sphere

reflspy[u_, v_, r_:10] := refl[spy[u, v, r]]

refllinspy[u_, v_, r_:10] := refl[linspy[u, 1/v, r]]

and obtain a conformally equivalent, interesting surface (zoom the created  Graphics3D!):

reflspy[u, v]

ParametricPlot3D[reflspy[u, v], {u, 0, 8 Pi}, {v, -.5, .5}, 
PlotPoints -> 90]

reflspcy  =  %

ParametricPlot3D[refllinspy[u, v,15], {u, 0, 8 Pi}, {v, -.5, .5}, 
PlotPoints -> 90]

refllinc  =  %

Apply MathGl3d (or MathLive via ThreeScript) for animating this graphics! (See section 7)

MVShow3D[reflspcy]

MVShow3D[refllinc]

6.4 Spiral Surfaces

The definition of the spiral surface generated by a curve xz[u] in the x, z-plane is given in spirals.m:

?spiralsurface

\(\(spiralsurface[func_, \ a_:  0.1]\)[u_, \ 
        t_]\  := \ 
\ \ \(\(spiral[t, \ a]\  . \ {\(func[u]\)[\([1]\)], \ 
          0, \ \(func[u]\)[\([2]\)]}\)\(\ \)\)\)

The corresponding plotting command is also imported from spirals.m:

\(\(\(\ \ \)\(\(\(plotspiralsurface[func_, \ a_:  0.1]\)[\ t0_, \ t1_, \ 
        u0_, \ u1_, \ 
        opt___]\)\(\ \)\(:=\)\(\ \)\(
\)\(	\)\(ParametricPlot3D[
        Evaluate[\(spiralsurface[func, \ a]\)[u, \ t]], \ 
	{t, \ t0, \ 
          t1}, {u, \ u0, \ u1}, \ opt]\)\(\ \)\)\)\)

Examples generated by these definitions

f1[s_] := 1  +  Sin[s]; f2[s_] :=  Cos[s]; 
    fct[s_] :=  {f1[s], f2[s]}

spiralsurface[fct, 2][u, t]

spiralsurface[fct][u, t]

plotspiralsurface[fct][-4 Pi, 2  Pi, 0,  2 Pi,  
      PlotPoints  ->  60]

ulitka  =  %

MVShow3D[ulitka]

ThreeScript["ulitka.ts",  ulitka]

xz[s_]  :=  {3  + s^2 , 1 + s^3}

spiralsurface[xz,   .5][u, t]

\(\(plotspiralsurface[xz,  .22]\)[\(-2\)*Pi, \ 2*Pi, \ \(-2\), \ 2, 

      PlotRange\  -> \ All, PlotPoints -> 50]\)

scubesf =  %

MVShow3D[scubesf]

fxz[u_]  :=  {1  +  0.2*Cos[u],  0.3*Sin[u]}

spiralsurface[fxz,   .3][u, t]

plotspiralsurface[fxz,   .17][-2 Pi,  4 Pi,  0,  2 Pi,  
      PlotPoints  ->  60,  PlotRange  ->  All]

spsf1 =  %

MVShow3D[spsf1]

ThreeScript["spsf1.ts", spsf1]

This plots a spiral section of a circular cone, 1 turn:

fxz[u_]  := {u, u};
plotspiralsurface[fxz, 0.2][0, 2*Pi,  0,  1,  
      PlotPoints  ->  45,  PlotRange  ->  All]

6.5 Spiral surfaces generated by Alfred Gray's curves collection

To generate other spiral surfaces we use parameter representations of certain plane curves. The parameter representations of these and other curves can be found in   A. Gray´s book [1].

\(\ \(\(astroid[a_]\)[t_] := \ {a*Cos[t]^3, 
          a*Sin[t]^3}\ ;\)\ \ \ \ \ \ \ \ \ \ \ \ \ 
\)
\(\ \ \ \ \(\(astroid[a_, b_]\)[t_] := \ {a*Cos[t]^3, 
          b*Sin[t]^3};\)
\)
\(\ \ \ \ \(\(astroid[n_, a_, b_]\)[t_] := \ {a*Cos[t]^n, 
          b*Sin[t]^n};\)
\)

ParametricPlot[astroid[1][u],  {u,  0,  2 Pi},  
      AspectRatio  ->  Automatic]

spiralsurface[astroid[1]][u, t]

plotspiralsurface[astroid[1], 0.06][ 0, 4 Pi,0, 2*Pi,
PlotPoints -> 60, PlotRange -> All]

ast  =  %

MVShow3D[ast]

fct[u_] :=  {2,  4}  +  0.5 astroid[1][u];

\(\(plotspiralsurface[fct, \ 0.06]\)[\ 0, \ 4\ Pi, 0, \ 2*Pi, 

      PlotPoints\  -> \ 40, \ PlotRange\  -> \ All]\)

spsf2 = %

MVShow3D[spsf2]

Clear[fct]

fct[u_] :=  {1.5  +  Sin[u], 
          1.5  +  Cos[u] Sin[u]};

\(\(plotspiralsurface[fct, \ 0.2]\)[\ \(-2\)\ Pi, \ 4\ Pi\ , \ 
      0, \ \ 2\ Pi, \ 
PlotRange\  -> \ All, \ 
      PlotPoints\  -> \ {25, \ 45}]\)

spsf3 = %

eight[
        u_] :=  {2, 2} + {Sin[u], Sin[u]*Cos[u]}

\(\(plotspiralsurface[eight, \ 0.3]\)[0, \ 4\ Pi, \ 0, \ 2\ Pi, \ 
\ 
      PlotRange\  -> \ All, \ PlotPoints\  -> \ {25, \ 45}]\)

spsf4 =  %

\(\(\(diamond[a_]\)[
          t_] := \ {2, \ 2}\  + 
\ \ \ \ \ \ \ \ {a*Cos[t]*Sqrt[Cos[t]^2], 
            a*Sin[t]*Sqrt[Sin[t]^2]};\)
\)
\(\ \ \ \ \(\(diamond[n_, a_, b_]\)[
          t_] := \ {2, 
            2}\  + 
\ \ \ \ \ \ \ \ {a*Sqrt[Cos[t]^2]^\((n\  - \ 1)\)*
              Cos[t], 
\ \ \ \ \ \ \ \ \ b*Sqrt[Sin[t]^2]^\((n\  - \ 1)\)*
              Sin[t]};\)\)

\(ParametricPlot[\(diamond[2, \ 1, \ 1]\)[u] // Evaluate, \ {u, \ 0.1, \ 
        2\ Pi\  + \ 0.1}, 
AspectRatio\  -> \ Automatic]\)

\(ParametricPlot[\(diamond[5, \ 1, \ 1]\)[u] // Evaluate, \ {u, \ 0.1, \ 
        2\ Pi\  + \ 0.1}, 
AspectRatio\  -> \ Automatic]\)

\(\(plotspiralsurface[diamond[4, \ 1, \ 1], \ 0.3]\)[\ 0, \ 4\ Pi, \ 
      0.1, \ 2\ Pi\  + \ 0.1, \ 
\ PlotRange\  -> \ All, \ 
      PlotPoints\  -> \ {25, \ 45}]\)

spsf5 = %

\(\(\(\ \)\(\(clothoidprime[n_, a_]\)[
        t_] := \ 
\ \ \ \ \ \ \ \ \ {a*Sin[t^\((n + 1)\)/\((n + 1)\)], 
        a*Cos[t^\((n + 1)\)/\((n + 1)\)]}\)\)\)

\(\(\(\ \)\(\(clothoid[n_, a_]\)[
        t_] := {1.8, \ 
          1.8}\  + 
\ \ \ \ \ \ \ \ \ Integrate[\(clothoidprime[n, a]\)[
            tt], {tt, 0, t}]\)\)\)

ParametricPlot[clothoid[1,  1][u] // Evaluate,  {u,  -5,  
        5},  AspectRatio  ->  Automatic]

\(\(plotspiralsurface[clothoid[\ 1, \ 1], \ 0.3]\)[\ \(-6\)\ Pi, \ \(-\ 
        Pi\), \(-6\), \ 6, \ 
\ PlotRange\  -> \ All, \ 
      PlotPoints\  -> \ {45, \ 45}]\)

spsf6  = %

MVShow3D[%]

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[conspir3,  MVLineTubeSize  ->  0.03]

MVShow3D[tgr3,  MVPointSphereSize  ->  1.8]

MVShow3D[ex123, MVPointSphereSize -> 1.8]

MVShow3D[tgr2]

Options[MVShow3D]

MVClose[]