(* symplecticgeo.m *)
(* Version 1. *)
(* by Rolf Sulanke *)
(* Date: January 11, 2006 *)

(* This package contains definitions and procedures for 2n-dimensional 
symplectic vector algebra and the corresponding projective symplectic \
\
\
geometry. 
The functions and modules defined here are developed in 
the notebook symplectic.nb.*)

(* Keywords:  Gram matrix, index of  a subspace, isotropic subspaces, lines, \
\
\
\
null system, orthosymplectic basis, polar, pole,  projective symplectic \
\
\
space, 
 spherical symplectic geometry, sym, symplectic basis, symplectic \
\
\
orthogonalization, 
 symplectic scalar product, symplectic transvection, symplectic vector \
\
\
space,*)

(* 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. For graphics functions the dimension dim must be set to 4. *)

BeginPackage["symplecticgeo`", {"vectorcalc`"}]

(*The usages *)

ssp::usage ="ssp[vx,vy] is the symplectic scalar product of two d-dimensional \
\
\
\
\
\
vectors; the dimension d= Length[vx] must be even."

grammatrix::usage = "grammatrix[b_List] is the matrix of the scalar products \
\
\
\
\
\
of the vectors b[[i]] of List b with respect to the scalar product \
\
\
ssp."

gram::usage = "gram[d:dim] is Gram's matrix of the scalar products of the \
\
\
\
\
\
standard basis in dimension d."

ortho::usage = 
    "ortho[b] is a renumbering of the vectors of a symplectic basis b such \
\
\
\
\
\
that the result is an orthosymplectic basis. WARNING! ortho[b] yields \
\
\
false \
\
\
results if b is not a symplectic basis, without any error \
\
\
message!"

orthogram::usage="orthogram[d], d even, is the Gram matrix of orthosymplectic \
\
\
\
\
\
bases in the d-dimensional symplectic vector space."

symplectic::usage = 
    "symplectic[b] is a renumbering of the vectors of an orthosymplectic \
\
\
\
\
\
basis b such that the result is a symplectic basis. WARNING! \
\
\
\
symplectic[b] \
\
yields false results if b is not an orthosymplectic \
\
\
basis, \
without any \
error \
message!"

orthosymplectic::usage =
    "orthosymplectic is an option for the Module symortho; its default value \
\
\
\
\
\
is False. If set to True, symortho yields an orthosymplectic vector \
\
\
\
\
sequence, \
and a symplectic one else."

symortho::usage= "symortho[b]: for a given vector sequence b of a symplectic \
\
\
\
\
\
vector space with scalar product ssp the result symortho[b] is a \
\
\
\
symplectic \
\
vector sequence with the same span as that of b; setting \
\
the \
\
option \
\
orthosymplectic -> True, the resulting sequence is \
\
\
\
orthosymplectic. The \
rank \
and the defect of the span of b is \
\
calculated \
\
and printed, too."

polarspace::usage = "Given a List m of vectors, polarspace[m] is a basis of \
\
\
\
\
\
the vector subspace orthogonal to the span of m."

transvection::usage = "transvection[b,c] is the symplectic transvection \
\
\
\
\
\
corresponding to the vector b and the scalar c."

transv::usage = "transv[b,c] is the  matrix of transvection[b,c]."

sym::usage = "sym[x, y, u, v] yields a symplectic invariant of two projective \
\
\
\
\
\
symplectic lines spanned by the vectors (x,y) and (u,v) \
respectively; \
\
for \
\
\
isotropic lines it is Infinity or Indeterminate."

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

Begin["`Private`"]

(* dimm:= Global`dim; *)

ssp[vx_,vy_] := Module[{ n=Length[vx],d},
    If[Head[vx] != List|| Head[vy] != List||OddQ[n]||
        Length[vy]!= n, Print["Incorrect arguments ",vx, vy];
      Return[], Null,Print["Arguments correct?"];Return[]];
    d=n/2;
    Sum[vx[[i]]*vy[[d+i]]-vx[[d+i]]*vy[[i]],{i,1,d}]]

grammatrix[b_List] :=Table[ssp[b[[i]], b[[j]]],{i,Length[b]},{j,Length[b]}]

gram[d_:dim]:= 
  If[EvenQ[d]&& d > 0,grammatrix[IdentityMatrix[d]],
    Print["Incorrect argument!"]; Return[], Print["Argument?"]; Return[]]
    
    ortho[c_] := Module[{d,b},
    If[MatrixQ[c],d=Length[c],Print["The argument is not a \
\
\
\
\
matrix!"];Return[],
      Print["The argument is not a matrix!"];Return[]];
    If[OddQ[d],Print["Incorrect argument ", c];Return[]];
    b= nullmatrix[d];
    For[i= 1,i<= d/2,i++, b[[2*i]]= c[[d/2 + i]];
      b[[2*i-1]]=c[[i]]];
    b]
    
 ortho[c_] := Module[{d,b},
    If[MatrixQ[c],d=Length[c],Print["The argument is not a \
\
\
\
\
matrix!"];Return[],
      Print["The argument is not a matrix!"];Return[]];
    If[OddQ[d],Print["Incorrect argument ", c];Return[]];
    b= nullmatrix[d];
    For[i= 1,i<= d/2,i++, b[[2*i]]= c[[d/2 + i]];
      b[[2*i-1]]=c[[i]]];
    b]

orthogram[d_:dim] := 
  Module[{c = ortho[IdentityMatrix[d]]}, 
    Table[ssp[c[[i]],c[[j]]],{i,d},{j,d}]]
    
symplectic[c_] := Module[{d,b},
    If[MatrixQ[c],d=Length[c],Print["The argument is not a \
\
\
\
\
matrix!"];Return[],
      Print["The argument is not a matrix!"];Return[]];
    If[OddQ[d],Print["Incorrect argument ", c];Return[]];
    b= nullmatrix[d];
    For[i= 1,i<= d/2,i++, b[[i]]= c[[2*i-1]];b[[d/2+i]]=c[[2*i]]];
    b]  
    
Options[symortho]={orthosymplectic -> False};

symortho[a_,opts___]:= Module[{b,bb,m, ea ,iso={}, c={},or},
    or = orthosymplectic /.{opts}/.Options[symortho];
    b=noprops[a];
    If[b==  Null, Return[]];
    
    m= Length[b];
    If[m==0, Return[b]];
    If[OddQ[Length[b[[1]]]], Print["Incorrect dimension!"];Return[]];
    
    Label[start];
    ea = 0;
    b = outzero[b,10^(-20)];
    m= Length[b];
    
    If[m<= 1, iso=Join[iso,b];Goto[fin]];
    
    For[i=2,i<= m,i++,
      If[ssp[b[[1]],b[[i]]]!= 0,ea = i;Goto[nnn],Null,
        Print["Scalar products may not be evaluated!"];Return[]]];
    
    If[ea== 0,
      For[i=2,i<=m,i++, 
        If[wedge[b[[1]],b[[i]]]== nullmatrix[], b[[i]]=null[]]];
      AppendTo[iso,b[[1]]];b=Table[b[[i]],{i,2,m}];Goto[start]];
    
    Label[nnn];
    AppendTo[c,b[[1]]];AppendTo[c,b[[ea]]/ssp[b[[1]],b[[ea]]]];
    bb=Join [Table[b[[i]],{i,2,ea-1}],Table[b[[l]],{l,ea+1,m}]];
    m=Length[bb];
    For[i=1,i<= m,i++,
      b[[i]] = bb[[i]] -c[[-2]]*ssp[bb[[i]],c[[-1]]] +
          c[[-1]]*ssp[bb[[i]],c[[-2]]]];
    bb  = Table[Chop[b[[i]]],{i,m}];
    b = noprops[bb];
    Goto[start];
    Label[fin];
    m=Length[c]/2;Print["The rank is ", 2*m];
    ea=Sum[ssp[c[[2*i-1]],c[[2*i]]],{i,m}];
    If[Chop[m-ea] != 0,
      Print["Scalar products may not be evaluated!"];Return[]];
    If[Length[iso]!= 0,iso = outzero[RowReduce[iso]]];
    Print["The defect is ", Length[iso]];
    ea= Sum[Chop[ssp[c[[i]],iso[[j]]]],{i,2*m},{j,Length[iso]}]+
        Sum[Chop[ssp[iso[[i]],iso[[j]]]],{i,Length[iso]},{j,Length[iso]}];
    If[Chop[ea] != 0,Print["Scalar products may not be evaluated!"];
      Return[]];
    If[or==False&&  c!= {},c=symplectic[c]];
    
    c=Join[c,iso];
    Print["The dimension is ", Length[c]];
    Return[c]]
    
 polarspace[m_List]:= Module[{b},b = Map[(gram[].#)&,m]; Chop[NullSpace[b]]]
 
 transvection[b_,c_][x_] := Module[{bb =Length[b], xx = Length[x]},
    If[bb != xx|| OddQ[bb], Print["Incorrect arguments!"];Return[]];
    x + b*c*ssp[b,x]]

transv[b_,c_] := Transpose[Table[transvection[b,c][stb[i]],{i,4}]]

sym[x_,y_,u_,v_] :=(ssp[x,u]*ssp[y,v]-ssp[x,v]*ssp[y,u])/(ssp[x,y]*ssp[u,v])

End[]

Protect[gram,grammatrix,ortho,orthogram,orthosymplectic,polarspace,ssp,sym,\
symortho, \
symplectic,transv, transvection]

EndPackage[]