(* ::Package:: *)
(* :Context: L3toHorosphere` *)
(* :Title: Central Projection of Hyperbolic Space Onto a Horosphere *)
(* :Author: Marijana Babic and Srdjan Vukmirovic *)
(* :Version: Mathematica 6.0 *)
(* :Package Version: 1.0 *)
(* :Keywords:
Hyperbolic space, horosphere
*)
(* :History: *)
(* :Requirements: Mathematica 6.0 or higher *)
(* :Warnings: *)
(* :Summary:
This package implements visualization of Hyperbolic space in Half space model, in
Klein projective model, and by means of central projection onto horosphere. Horosphere is isometric to Euclidean plane. Therefore,
we can say that if we embedd computer screen into Hyperbolic space it may become a part
of horosphere. This package also contains implementation of isometries in hyperbolic space.
*)
(* :Sources: *)
BeginPackage["L3toHorosphere`"]
PSphere2Klein::usage = "PSphere2Klein[u] converts point u from Poincare sphere model (the one where lines are circles orthogonal to the apsolute) into Klein projective model.
In Klein model, the point has 4 coordinates "
Klein2PSphere::usage = "Klein2PSphere[u] converts point u from Klein projective model (4 coordinates) "
Klein2HS::usage = " Klein2HS[u] converts point u from Klein projective model (4 coordinates) into half space model "
HS2Klein::usage = "HS2Klein[u] converts point u from half space model into Klein projective model (4 coordinates) "
(*
PSphere2HS::usage = " PSphere2HS "
HS2PSphere::usage = " HS2PSphere "
*)
drawProjection::usage = "drawProjection[\[Omega], k, vertK, faces, options] projects polyhedra vith vertex coordinates vertK, and vertex connectivity given by faces onto horsphere z=k,
point (0,0,\[Omega]) of half space model. Options:
horizon->False: do not draw the circle that separates front and back of the observer
visibility->True: draw non-visible lines dashed. Visibility works only for convex, closed polyhedra.
color->Black: by default draw projection in black color. Color specification is RGBColor[r,g,b]
thickness->Thin: by default draw lines thin. Thickness specification is of form Thickness[r]
stereo->e: Project from point (e,0,\[Omega]) - useful for stereo viewing experiment. Default is 0."
modelHS::usage = "modelHS[\[Omega], k, vertK, faces, options] draws, in the half space model, a polyhedra and its central projection onto horosphere z=k, from point
(0,0,\[Omega]). The polyhedara has vertices coordinates vertK in Klein model and vertx connectivity given by faces. Options:
range->{{-5,5}, {-5,5}, {0,5}} is basically the plot range
drawSphere->False: whether to draw half sphere, representing plane that divides front and back of the observer.
drawAbsolute->False: whether to draw the absolute z = 0
drawHorosphere->False: whether to draw the horosphere z=k
drawRays->True: whether to draw the rays of projection"
edges::usage = "edges[faces] finds all edges of polyhedra with given faces. Purely combinatorial function. "
normalVector::usage = "normalVector[{pt1, pt2, pt3}] returns Minkowski normal vector of the plane determinde by three points. The returned normal vector is given in
homogenous coordinates, i.e. of the length 4"
m::usage = "m[a, b] returns hyperbolic midpoint of segment ab. The argument and result are 4 homogenous coordinates"
translate::usage = "translate[a, b] returns 4x4 matrix that represents translation from point a to point b. Points a and b are given in Klein coordinates.
The resulted matrix is applied to any point of Klein model by means of matrix multiplication."
rotate::usage = "rotate[a,b][\[Theta]] returns 4x4 matrix that represent the rotation around line ab by angle \[Theta]. Points a and b are given in Klein coordinates.
The resulted matrix is applied to any point of Klein model by means of matrix multiplication."
reflect::usage = "reflect[a] returns 4x4 matrix that represent the reflection with respect to point a. Points a is in Klein coordinates.
The resulted matrix is applied to any point of Klein model by means of matrix multiplication."
limitRotate::usage = "limitRotate[P][p][fi] returns 4x4 matrix that represent the limit (horocyclic) rotation. It is composition of two paralel planes
\[Alpha] and \[Beta] intersecting along line with 3D vector p and touching unit sphere at 3D point P. The fi is 'angle' of horocyclic rotation and belongs to
open interval (-Pi/2, Pi/2). If fi = Pi/2 (mod Pi), the error message will appear. Also p must be orthogonal to P.
The resulted matrix is applied to any point of Klein model by means of matrix multiplication."
Options[drawProjection] = {horizon->False, visibility->True, color->Black, thickness->Thin, stereo->0};
Options[modelHS] = {range->{{-5,5}, {-5,5}, {0,5}}, drawSphere->False, drawAbsolute->False, drawHorosphere->False, drawRays->True}
Begin["`Private`"]
(* PSphere (Poincare sphere) model is interior of unit sphere (called Absolute) embedded in R3. Lines and planes are parts of circles (or lines)
and spheres (or planes) orthogonal to the Absolute. Isometries are generated by inversion (reflections) with respect to planes. Only conversion of coordinates
from this model to another two models are supported *)
(* Klein projective model is hyperquadric x^2 + y^2 +z^2 -w^2=0 in projective 3 space with homogenous coordinates (x:y:z:w). Isometries in this model
are projective transformations that preserve the hyperquadric i.e. the bilinear form given by matrix diag(1:1:1:-1). The important is that these isometries are linear in
homogenous coordinates. In affine coordinates (x,y,z), when w = 1, Klein model becomes interor of unit sphere in R3 with lines and planes being parts of euclidean lines and planes. *)
(* The Klein model is used for perfoming isometries, visibility calculations..., since it is linear *)
(* HS (half space) model is upper half plane z>0 of R3, where the plane z=0 is the Absolute. Lines are semicircles orthogonal to the Absolute and planes are
hemispheres or halfplanes orthogonal to the Apsolute. In this model isometries are generated by inversions (reflections) with respect to the planes. This model is
used for projection onto horosphere because the plane z=k represents horosphere in this model *)
(************************************* Conversions beetween models ***********************************)
(* Poincare sphere model to Klein model *)
PSphere2Klein[u_]:=Append[(2u)/(1+u.u),1]
(* Klein model to Poincare sphere model *)
Klein2PSphere[s_]:=Module[{sn},
sn = Drop[s/s[[4]], -1];
sn/(1+Sqrt[1-sn.sn])
] (* Module *)
(* To get half space from unit sphere, we translate the sphere by -1 along z-axis. Then we perform inversion with respect
sphere with center in origin and radius 2 plus Oxy plane reflection. Finally we translate by 2 along z-axis and perfom additional
reflection with respect to Oxy plane.*)
inversion[{x_, y_, z_}]:=4/(x^2 + y^2 + z^2) {x, y, z}
(*Poincare sphere to half space model*)
PSphere2HS[{x_,y_,z_}]:=
(inversion[{x,y, z}-{0,0,1}]+{0,0,2}){1,1,-1}
HS2PSphere[{x_,y_,z_}]:= inversion[{x, y, z}{1,1,-1}-{0,0,2}]+{0,0,1}
Klein2HS[s_]:=PSphere2HS[Klein2PSphere[s]];
HS2Klein[u_]:=PSphere2Klein[HS2PSphere[u]];
(************************************* Isometries of Hyperbolic space ***********************************)
(* All isometries result in 4x4 matrix. Arguments are points Klein model *)
(* matrix of (Minkowski) scalar product in Klein projective model. Isometries of hyperbolic space (Klein projectiva model) preserve this scalar product *)
j=({
{1, 0, 0, 0},
{0, 1, 0, 0},
{0, 0, 1, 0},
{0, 0, 0, -1}
})
(* hyperbolic midpoint in Klein model *)
m[a_, b_]:=a Sqrt[(b.j.b)*(a.j.b)]+ b Sqrt[(a.j.a)*(a.j.b)];
(* normal vector to the plane determined by three points *)
normalVector[{pt1_, pt2_, pt3_}] := {1,-1,1,1}(Minors[{pt1, pt2, pt3}, 3][[1]]//Reverse)
(* 4x4 matrix that represents reflection with respect to the point p (given by 4 homogenous coordinates Klein model) *)
reflect[p_]:=Module[{pm},
pm = List /@ p;
IdentityMatrix[4] -(2 pm.Transpose[pm]. j)/(p.j.p)//N
]
(* Translation from point a to point b *)
translate[a_, b_]:=reflect[m[a,b]].reflect[a]
(* rotation around axes ab (a and b are given points) by angle \[Theta]:
First, we translate ab so it contains the origin, then we rotate around origin to align ab with z-axis, then we rotate by angle \[Theta] around z-axis
Then we rotate back and translate back. All rotations around origina are euclidean since euclidena and hyperbolic angles coincide in the origin. *)
rotate[a_,b_][\[Theta]_]:=Module[{an, bn, bt, u, v, w, m, r},
bt=translate[a,{0,0,0,1}]. b ;
bn=Drop[bt/bt[[4]], -1];
(* rotation matrix by angle \[Theta] around z-axis*)
r=({
{Cos[\[Theta]], Sin[\[Theta]], 0, 0},
{-Sin[\[Theta]], Cos[\[Theta]], 0, 0},
{0, 0, 1, 0},
{0, 0, 0, 1}
});
If[ bn[[1]]^2+bn[[2]]^2<0.00000000001,
(* If after translation ab is aligned to z-axis, we perform the rotation and translate back *)
translate[{0,0,0,1},a].r.translate[a,{0,0,0,1}],
(* else - make new orthonormal basis with z axes aligned to bn *)
w=bn/Norm[bn];
v=Cross[w,{0,0,1}]/Norm[Cross[w,{0,0,1}]] ;
u=Cross[v,w];
m=({ (* matrix of rotation that makes translated ab to be new z-axes *)
{u[[1]], u[[2]], u[[3]], 0},
{v[[1]], v[[2]], v[[3]], 0},
{w[[1]], w[[2]], w[[3]], 0},
{0, 0, 0, 1}
});
(* m is matrix of rotation so m^-1=m^T. We rotate by r around new z-axis and translate back *)
translate[{0,0,0,1},a].Transpose[m].r.m.translate[a,{0,0,0,1}]
](*If*)
](*Module*)
(* horocyclic rotation *)
limitRotate[P_][p_][fia_]:=Module[{na,nb, a1, a2, a3, b1, b2,b3},
{a1, a2, a3}=Cos[fia] Cross[P,p] + Sin[fia] P;
{b1, b2, b3}=Cos[0] Cross[P,p]+ Sin[0] P;
na={a1, a2, a3, a1 P[[1]]+a2 P[[2]]+a3 P[[3]]};
nb={b1, b2, b3, b1 P[[1]]+b2 P[[2]]+b3 P[[3]]};
reflect[nb].reflect[na]
]/; (N[P.p] < 0.00000001)
(* central projection from eyepoint (0,0,\[Omega]) onto a horosphere z=k of half space model.
The point is given in the HS model.
We assume that \[Omega] > k, that is THE EYEPOINT INSIDE THE HOROSPHERE *)
projection[\[Omega]_, k_][{x_,y_,z_}]:=Module[{m, n},
If[ x^2+y^2<=0.0000000001, (* if point belong to z-axis *)
If[z<\[Omega],
{x, y}, (* if M is below the eye-point its projection is the normal projection *)
{\[Infinity],\[Infinity]} (* if it is above it projects to infinity of the plane (must be a pair of numbers to be comparable - used in function arcAB *)
],
(* else, if point doesn't belong to z-axis, we project it along arc, i.e. hyperbolic line *)
m = (\[Omega]^2 - x^2 - y^2 - z^2)/(2 Sqrt[x^2 + y^2]);
n = -m + Sqrt[m^2 - k^2 + \[Omega]^2] ;
{(n x)/Sqrt[x^2 + y^2], (n y)/Sqrt[x^2 + y^2]}
] (* if *)
] (* projection *)
(* The projection of segment ab into horosphere is circular arc. That circular arc is intersection of horosphere z=k and
half-sphere with center on Absolute z=0 containing a, b and eyepoint O. The center of that half sphere is intersection of
medial planes of segments ab, Ob and the plane z=0. *)
centerArcAB[\[Omega]_, k_][{a1_, a2_}, {b1_, b2_}]:=Module[{t} ,
t=(a1 b1+a2 b2-k^2+\[Omega]^2)/(-2 a2 b1+2 a1 b2);
{t (a2 - b2)+ (a1 + b1)/2, t (b1 - a1) + (a2 +b2)/2} (* coordinates of center AB in the horosphere z=k *)
]
angle[t_]:= If[t<0, 2 Pi + t, t] ;(* function that converts angle from [-Pi, Pi] to [0, 2 Pi] *)
maxx= 1000;
Ray[a_,b_]:= Line[{a, a + maxx (b-a)}]
(* The projection of segment is circular arc. There are two arcs with given endpoints, therefore we must specify one point on the arc *)
(* Function for drawing circular arc ab containing point s in projection. Points a, b, and s are projections, i.e. belong to the horosphere z=k *)
(* argument vis tells if the line visible or not (in the latter case it will be drawn dashed *)
arcAB[\[Omega]_, k_][a_, b_, s_, vis_:False]:= Module[{c, v1, v2, vs, r, v1u, v2u, vsu} ,
(*We first check if some of endpints is in the infinity *)
(*1*)
If[a=={\[Infinity],\[Infinity]}, (* If 1 *)
If[b=={\[Infinity],\[Infinity]},
(* (1) If both endpoints are infinite we draw nothing *)
{},
(* (2) If a is infinite, b finite, we draw ray from b in the direction of s *)
Prepend[#, If[vis,{}, Dashed]]& [{Ray[b,s]}]
], (* If *)
If[b=={\[Infinity],\[Infinity]}, (* If 2 *)
(* (3) If a is finite, b infinite, we draw ray from a in the direction of s *)
Prepend[#, If[vis,{}, Dashed]]& [{Ray[a, s]}],
(* (4) If both endpoints are finite *)
If[ Abs[a[[1]] b[[2]] - a[[2]] b[[1]] ]<0.0000000000001,
(* we check whether eyepoint O, a, b belong to the plane orthogonal to Absolute *)
(* If yes, if s is beetween a and b we draw a segmente,
If no, we draw two concurent rays with endpoints a and b *)
If[Norm[s-a]0,
Prepend[#, If[vis,{}, Dashed]]& [{Line[{a, b}]}],
{Prepend[#, If[vis,{}, Dashed]]& [{Ray[a,2a-b]}],
Prepend[#, If[vis,{}, Dashed]]& [{Ray[b, 2b-a]}]}
], (* If *)
(* If no, projection of the segment is circular arc that we find in the sequel: *)
c = centerArcAB[\[Omega], k][a, b];
v1 = a - c; (* vectors from center of arc to points a and b *)
v2 = b - c;
vs =s -c; (* vector of point s on the arc *)
v1u = angle[N[Arg[{1,I}.v1]]]; (* convert the angles into range [-Pi, Pi], i.e.[0, 2 Pi] *)
v2u = angle[N[Arg[{1,I}.v2]]];
vsu = angle[N[Arg[{1,I}.vs]]];
{v1u, v2u}=Sort[{v1u, v2u}];
r = N[Sqrt[v1.v1]];
Prepend[#, If[vis,{}, Dashed]]& [{Circle[c,r,If[ vsuFalse, all edges are visible *)
]; (* If *)
Graphics[Append[ (* the projection of all edges (arcs) is returned *)
{OptionValue[color], OptionValue[thickness],
Table[arcAB[\[Omega],k][vertHS[[edgeList[[n]][[1]]]],
vertHS[[edgeList[[n]][[2]]]],
{e,0} + projection[\[Omega], k][Klein2HS[Append[midEdges[[n]], 1]]+ {-e,0,0}],
MemberQ[visible, n]
], (* arcAB *)
{n,1, Length[edgeList]}
]}, (* Table *)
If[ OptionValue[horizon],{Black, Thin, visibilityCircle}, {}]] (* add visibilyCircle to the final image if neccesary *)
] (* Graphics *)
] (* drawProjection *)
(* Center of arc ab normal to Absolute (that is segment ab in HS model) is obtained in the intersection
1. medial plane of ab, with equation:
(b1-a1)(x-(a1+b1)/2)+(b2-a2)(y-(a2+b2)/2)+(b3-a3)(z-(a3+b3)/2)=0
2. plane orthogonal to the Absolute z=0, containing a(a1,a2,a3) and b(b1,b2,b3)
with parametric equation: x=a1+s(b1-a1), y=a2+s(b2-a2), z=a3+s(b3-a3)+t
3. Absolute z=0.
*)
centerArcHS[{a1_, a2_, a3_}, {b1_, b2_, b3_}]:=Module[{s} ,
s=(-a3^2+(a1-b1)^2+(a2-b2)^2+b3^2)/(2 ((a1-b1)^2+(a2-b2)^2));
{a1+ s (b1-a1),a2+s (b2-a2) , 0}
];
(* Function for hyperboic segment ab in HS model
If ab is circular arc:
The plain containing arc is determined by vectors ab and vector e={0,0,1} which are not orthogonal.
Therefore, in that plane we find unit vector u orthogonal to e: w=
\!\(\*OverscriptBox["AB", "\[LongRightArrow]"]\)*e , u=(e*w)/(||e*w||).
u is of the form u(u1, u2, 0).
Now, the circle is parameterized by X(t)=c+ r Cos[t] u + r Sin[t] e
*)
arcHS[{a1_, a2_, a3_}, {b1_, b2_, b3_}]:=
Module[{ v1, r,c, u,e, tA, tB} ,
If[
(* if AB is orthogonal to Absolute that is a1=b1 i a2=b2,*)
(a1 -b1)^2+(a2 -b2)^2<0.000000000001,
Graphics3D[Line[{{a1,a2, a3}, {b1, b2, b3}}]], (* then hyperbolic segment ab is euclidean segment *)
(* else: *)
c =centerArcHS[{a1, a2, a3}, {b1, b2, b3}];
v1 = {a1,a2, a3} - c; (* vector from the center to point a*)
r = N[Sqrt[v1.v1]];
u=1/Sqrt[(b1-a1)^2+(b2-a2)^2] {b1-a1, b2-a2, 0};
e={0,0,1};
(*
We look for parameter tA and tB along the circle which correspond to endpoints a and b of the arc i.e.
c1 + r costA u1 = a1, c2 + r costA u2 = a2, c3 + r sintA = a3
Points a and b satisfy z>0, therefore tA, tB belong to {0, Pi}. Therefore tA is obtained from the first equation: tA=ArcCos[(a1-c1)/(r u1)].
This might cause problem only if u1=0, that is a and b belong to the plane orthogonal to x-axes, that is a1=b1. Then tA is determined from the second equation tA=ArcCos[(a1-c1)/(r u1)]. If u1 and u2 ar both 0, then ab is segment, not an arc, and this is the case considered above.
*)
If[a1==b1,
tA=ArcCos[(a2-c[[2]])/(r u[[2]])]; tB=ArcCos[(b2-c[[2]])/(r u[[2]])], (*else*)
tA=ArcCos[(a1-c[[1]])/(r u[[1]])]; tB=ArcCos[(b1-c[[1]])/(r u[[1]])]
];
Sort[{tA,tB}];
ParametricPlot3D[ c + r Cos[t]u + r Sin[t]e,{t,tA,tB}, PlotStyle->{Thickness[0.001]}]
] (* If *)
] (* Module *)
(* Function that draws ray from eyepoint a through point b, all to the Apsolute *)
rayHS[{a1_, a2_, a3_}, {b1_, b2_, b3_}]:=
Module[{c, v1, r, u,e, tA, tB, t0} ,
If[
(* If ab i orthogonal to the Absolute, i.e. if a1=b1 and a2=b2, it is Euclidean segment *)
(a1 -b1)^2+(a2 -b2)^2<0.0000000000001,
If [a3>b3,
(* then it is Euclidean segment from a to Absolute *)
Graphics3D[{ Opacity[0.3], Line[{{a1,a2, a3}, {b1, b2, 0}}] }],
(* else it is Euclidean segment from a up to infinity *)
Graphics3D[{ Opacity[0.3], Line[{{a1,a2, a3}, {b1, b2, b3+1}}] }]
] ,
(* otherwise it is a circular arc *)
c =centerArcHS[{a1, a2, a3}, {b1, b2, b3}];
v1 = {a1,a2, a3} - c; (* vector from center of the arc to the point a *)
r = N[Sqrt[v1.v1]];
u=1/Sqrt[(b1-a1)^2+(b2-a2)^2] {b1-a1, b2-a2, 0};
e={0,0,1};
If[a1==b1,
tA=ArcCos[(a2-c[[2]])/(r u[[2]])]; tB=ArcCos[(b2-c[[2]])/(r u[[2]])], (*else*)
tA=ArcCos[(a1-c[[1]])/(r u[[1]])]; tB=ArcCos[(b1-c[[1]])/(r u[[1]])]
];
t0=If[tAOpacity[0.3]]
] (* If *)
] (* Module *)
(* Function that draws picture in HS model *)
modelHS[\[Omega]_, k_, vertK_, faces_, OptionsPattern[]]:=Module [{sphere,horosphere, r, poly,absolute, rays, arcs, vertHS, edgeList,x1, x2, y1, y2},
r = OptionValue[range];
poly = {{r[[1,1]], r[[2,1]],0}, {r[[1,2]], r[[2,1]],0}, {r[[1,2]], r[[2,2]],0}, {r[[1,1]], r[[2,2]],0}};
horosphere = Graphics3D[{Opacity->0.5,Polygon[(# + {0,0,k})& /@ poly ]}];
absolute = Graphics3D[Polygon[poly]];
sphere =Graphics3D[{Opacity[0.3],Sphere[{0,0,0},\[Omega]]}];
edgeList = Union[ Sort /@ Flatten[edges /@ faces , 1]];
vertHS= Klein2HS[#]&/@vertK;
rays=rayHS[{0, 0, \[Omega]}, #]&/@ vertHS;
arcs=arcHS[vertHS[[#[[1]]]], vertHS[[#[[2]]]]
]&/@edgeList;
Show[
Graphics3D[{Point[#],Red, Point[{0, 0, \[Omega]}]}& /@ vertHS ],
projPolyhedraHS[\[Omega], k, vertK, faces],
Sequence @@ arcs,
If[OptionValue[drawRays], rays,{}],
If[OptionValue[drawSphere], sphere,{}],
If[OptionValue[drawHorosphere], horosphere,{}],
If[OptionValue[drawAbsolute], absolute,{}],
Boxed->False,PlotRange->OptionValue[range]
] (* Show *)
] (*Module*)
projPolyhedraHS[\[Omega]_, k_,vertK_, faces_]:= Module[{edgeList, projVert, arcs},
edgeList = Union[ Sort /@ Flatten[edges /@ faces , 1]];
projVert = N[projection[\[Omega], k][Klein2HS[#]]]& /@ vertK;
(* for each edge we draw arc *)
arcs=projectionArcHS[\[Omega],k][
projVert[[#[[1]]]],
projVert[[#[[2]]]],
projection[\[Omega], k][Klein2HS[Mean[{vertK[[#[[1]]]], vertK[[#[[2]]]]}]]]]
&/@edgeList;
Show[arcs]
]
rayHS[a_,b_,k_]:=Module[{ahs,bhs},
ahs=Append[a,k];
bhs=Append[b,k];
Graphics3D[Line[{ahs, ahs + maxx (bhs-ahs)}]]]
(* Function that draws circular arc ab, containing point S *)
projectionArcHS[\[Omega]_, k_][a_, b_, s_]:= Module[{c, v1, v2, vs, r, v1u, v2u, vsu, u1, u2} ,
(* We first check if some (or both) of a or b is point in infinity *)
If[a=={\[Infinity],\[Infinity]},
If[b=={\[Infinity],\[Infinity]},
(* If both are infinite we draw nothing *)
{},
(* else if a is infinite and b finite we draw an arc from b containing s *)
rayHS[b,s,k]
],(* If *)
If[b=={\[Infinity],\[Infinity]},
(* If a is finite and b infinite we draw ray from a in the direction s *)
rayHS[a, s,k],
(* if both a and b are finite *)
If[ Abs[a[[1]] b[[2]] - a[[2]] b[[1]] ]<0.000000000001, (* if O, a, b, belong to the plane orthogonal to Absolute *)
(* If they belong and projection of center is beetween a and b, we draw segment
otherwise we draw two rays *)
If[Norm[s-a]0,
Graphics3D[Line[{Append[a,k], Append[b,k]}]],
{rayHS[a,2a-b,k],rayHS[b, 2b-a,k]}
],
(* else, when the do not belong to the plane orthogonal to Absolute the projection is circular arc: *)
c = centerArcAB[\[Omega], k][a, b];
v1 = a - c;
v2 = b - c;
vs =s -c;
v1u = angle[N[Arg[{1,I}.v1]]];
v2u = angle[N[Arg[{1,I}.v2]]];
vsu = angle[N[Arg[{1,I}.vs]]];
{v1u, v2u}=Sort[{v1u, v2u}];
r = N[Sqrt[v1.v1]];
{u1, u2}=If[ vsu