(* vectorcalc.m *)
(* Version 1. *)
(* by Rolf Sulanke *)
(* Date: January 6, 2006 *)

(* This package contains definitions and procedures for n-dimensional 
vector algebra. The functions and modules defined here can be found in 
the notebook symplectic.nb; they are useful for calculations in euclidean 
or pseudo-euclidean, hermitean or pseudo-hermitean geometries, too. Older 
versions are contained in the package euvec.m, and have been applied e.g. 
in the notebook Euklid.nb.*)

(* Keywords: vector objects, null vector, basis vectors, random vector,
stereographic projection, inverse stereographic
projection, matrix, random matrix, rank, transposition matrix, wedge product. \
\
\
*)

(*Most of the concepts depend on the dimension dim, which must be defined in
the  Global Context! *)

(* Free software: Feel free to apply, distribute or change this package. \
\
\
Please,
cite the source and notice the changes, or new concepts,  which you
 introduced.*)

(*Warning. The package vectorcalc.m can not substitute the package euvec.m in \
\
\
the item 
Spheres, MathSource, Enhancements/Geometry 0210-092. *)

BeginPackage["vectorcalc`"]

(*The usages *)
 
 (* dim::usage="dim is the dimension of the vector space under consideration; \
\
\

 it must be set within the Global Context." *)

stb::usage="stb[i, d:dim] is the i-th standard basis vector, 1 <= i <= d."

vec::usage="vec[x, n:dim] is the n-vector with components x[i]."

null::usage="null[d:dim] is the zero vector of length d."

dotnorme::usage="dotnorme[b] yields the unit vector with respect to the Dot
  product in direction of b, b!=o."

sterproj::usage="sterproj is the stereographical projection of the unit 
  n-sphere from the north pole stb[n+1] onto the euclidean (equatorial) \
\
\
\
\
n-space. 
  WARNING: The result is correct only for arguments of norm 1 !"

invsterproj::usage="invsterproj is the inverse of the general euclidean 
  stereographic projection sterproj."
  
projpt::usage = "projpt[vx] is the point in the affine model of the \
\
\
projective space corresponding to the vector vx; if the last component of \
\
vx \
equals 0, projpt returns ComplexInfinity."

reprprojpt::usage = "vx = reprprojpt[px] is a vector with last component 1 \
\
\
representing the projective point px = projpt[vx]."

projpt3D::usage = "projpt3D[a] is the Graphics primitive Point of the \
\
\
corresponding projective point corresponding to the 4-dimensional vector \
\
a."

spherpt3D::usage = 
    "If a is a vector of dimension 4, then spherpt3D[a] is the graphics3D \
\
\
primitive Point of the stereographical projection sterproj of the normed \
\
\
vector a.";

randomv::usage="randomv[n:dim] is a random n-vector with values in the  \
\
positive unit n-cube."

smoothing::usage = "smoothing[a, b] is a listable function returning zero if 
Abs[a]< b: 10^-10, and a else; for b = 1 a is returned: smoothing[a, 1] = a."

wedge::usage = "wedge[a,b] is the alternate product of vectors a,b given by \
\
the matrix of its coordinates."

outzero::usage="outzero[b_List,dd] is the sublist of the vectors b[[i]] with \
\
\
\
\
Abs[b[[i]].b[[i]]] > dd; the default value of dd is 10^(-20)."

basis::usage="basis[v] yields a sublist of v being a basis of the span of 
  the vectors v[[i]] of v; v must be a matrix."

noprops::usage = "noprops[c] deletes rows c[[j]] of the matrix c which are \
\
proportional to an element c[[i]] with i<j, and returns the matrix \
\
originated that way."

matrix::usage = "matrix[b,n:dim] is a square matrix of order n with elements \
\
\
\
\
b[i,j]."

nullmatrix::usage = "nullmatrix[d] is the square nullmatrix of order d; the \
\
\
\
\
default value of d is dim."

blf::usage = "blf[a, d:dim] is a symmetric matrix of order d with elements \
\
a[i,j]."

randommatrix::usage = "randommatrix[k,l]is a matrix with k rows and l \
\
\
\
\
columns, with random elements out of the unit intervall."

rank::usage="rank[m] calculates the rank of the rectangular matrix m of a not
   too large order (< 10000)."

trp::usage = "trp[i,j,n:dim] is the unit matrix of order n where the i-th and \
\
\
\
\
the j-th row are transposed." 

renorm::usage = 
  "renorm[b] for a given matrix b  is the matrix with  normed vectors 
b[[i]]/|b[[i]]|, where |x| denotes the usual euclidean norm of the vector \
\
\
\
x. \
Before norming the zero vectors of b are deleted."


(*The dimension must be set within the Global Context.*)


Begin["`Private`"]

dimm:= Global`dim;

(*A function which makes small entities to become zero.*)
  
  smoothing[a_, b_:10^-10]:=
  Module[{ c= a}, If[b == 1, Return[c], If[Abs[c]< b, c=0]]; c];
Attributes[smoothing]=Listable  

(*The standard basis.*)

stb[i_Integer,d_:dimm]:=IdentityMatrix[d][[i]]

(* A general vector *) 

vec[x_,n_:dimm]:=Array[x,n]

(* The null vector *)

null[d_:dimm]:=Table[0,{d}]

(* Norming a vector *)

dotnorme[b_]:=b/Sqrt[b.b]

(* Stereographical  projection *)

sterproj[v_List]:=Module[{m=Length[v],l,w},l=null[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]

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

projpt[vx_] := Module[{dd = Length[vx], v = vx.vx},
    If[Chop[v] == 0, Return["nopt"]];
    If[Chop[vx[[dd]]]== 0, Return[ComplexInfinity]];
    Take[vx/vx[[dd]],dd-1]]
    
reprprojpt[a_List] :=Join[a,{1}]

projpt3D[a_] := Point[projpt[a]]

spherpt3D[a_] := Point[sterproj[dotnorme[a]]]

(* Random vectors *)

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

(* Wedge product *)

wedge[a_,b_]:= 
  Module[{d=Length[a]},Table[a[[i]]*b[[j]]-a[[j]]*b[[i]],{i,1,d},{j,1,d}]]

(* Vector sequences *)

outzero[b_List, dd_:10^(-20)]:=
  Module[{v={},l=Length[b]},
    For[i=1,i<=l,i++,
      If[Not[Chop[b[[i]].b[[i]],dd]==0],AppendTo[v,b[[i]]],Null[],
        AppendTo[v,b[[i]]]]];Return[v]]

basis[v_]:=
  Module[{r,l,w,i,fdrop},
    If[Not[MatrixQ[v]],Print["Argument is not a matrix, no span!"];Return[]];
    r=rank[v];
    If[r==0,Print["The span is the null space!"];Return[]];
    w=outzero[v, 10^(-20)];
    fdrop[i_]:=Drop[w,{i}];
    Label[new];
    l=Length[w];
    If[r==l,Return[w]];
    For[i=1,i<=l,i++,
      If[rank[fdrop[i]]==r,w=fdrop[i];Goto[new]]]]

noprops[c_]:= Module[{m=1,cc, d},
    If[Length[c]<=1,Return[c]];
    If[Not[MatrixQ[c]], Print["incorrect argument!"];Return[]];
    d=Length[c[[1]]];
    cc=outzero[c];
    If[Length[cc]<=1,Return[cc]];
    Label[ttt];
    For[j=m+1,j<= Length[cc],j++, 
      If[Chop[wedge[cc[[m]],cc[[j]]],10^(-20)]== nullmatrix[d],
        cc[[j]]=null[d]]];
    cc=outzero[cc];m = m+1;
    If[Length[cc] <= m, Return[cc],Goto[ttt]]]

(* Matrices *)

matrix[b_,n_:dimm]:= Array[b,{n,n}]

nullmatrix[d_:dimm]:= Table[0,{i,1,d},{j,1,d}]

rank[m_] := Module[{ di , nsp},
    If[Not[MatrixQ[m]], Print["incorrect argument!"];Return[]];
    di = Length[m[[1]]];
                  nsp = NullSpace[m];
               di - Length[nsp]]
               
blf[a_,d_:dimm]:= Module[{b = matrix[a]},
    For[i = 1,i<d,i++,For[j=i+1,j<=d,j++,a[j,i]=a[i,j]]];
    b]

randommatrix[k_,l_] := Table[randomv[l], {k}]

trp[i_,j_,n_:dimm]:=Module[{m=IdentityMatrix[n]}, m[[i]]=stb[j,n];
    m[[j]]=stb[i,n];
    m]


(* The renorming procedure (useful for badly conditioned matrices) *)

renorm[b_?MatrixQ]:= Module[{bo = outzero[b],i },
	Table[dotnorme[bo[[i]]], {i, Length[bo]}]]

End[]

Protect[basis, blf, dotnorme, invsterproj, matrix,noprops,null, nullmatrix, \
\
\
outzero, projpt, reprojpt, projpt3D,spherpt3D,randommatrix, randomv, \
rank, \
renorm, smoothing, stb,sterproj,trp, \
vec, wedge]

EndPackage[]