(* ::Package:: *)

(* ::Title:: *)
(*Plucker coordinates in \[DoubleStruckCapitalP]^3*)


(* ::Chapter:: *)
(*0. Introduction*)


(* ::Subchapter:: *)
(*Preliminary general routines*)


(* ::Section:: *)
(*Matrix functions*)


(* ::Subsection:: *)
(*Display matrices using MatrixForm*)


MF[x_]:=Map[MatrixForm,x]
MF2[x_]:=Map[MatrixForm,x,2]


(* ::Subsection:: *)
(*Introduce zero vector and zero matrix*)


ZeroMatrix[n_]:=Table[0,{i,n},{j,n}]
ZeroVector[n_]:=Table[0,{i,n}]


(* ::Subsection:: *)
(*Direct sum of matrices*)


AppendZeroes[M_,n_]:=Map[Join[#,Table[0,{n}]]&,M]
AppendZeroes::usage = "AppendZeroes[M,n] takes a list of 
lists M and appends n zeros to each list in M.  If M is 
a matrix, this has the effect of adding n columns of zeros 
at the end";

PrependZeroes[M_,n_]:=Map[Join[Table[0,{n}],#]&,M]
AppendZeroes::usage = "PrependZeroes[M,n] takes a list of 
lists M and adjoins n zeros to the beginning of each list in M.  
If M is a matrix, this has the effect of adding n columns of 
zeros at the beginning";

DirectSum[A_,B_]:=Join[AppendZeroes[A,Length[B]],PrependZeroes[B,Length[A]]]
DirectSum::usage="DirectSum[A,B] returns the direct sum of the square matrices A and B.";


(* ::Section:: *)
(*Type checking and conversion routines*)


(* ::Subsection:: *)
(*Check types for skew-symmetric 4\[Cross]4 matrices and vectors in (\[DoubleStruckCapitalR]^6):*)


FourSkewSymQ[M_]:=(Dimensions[M]== {4,4}) && (M + Transpose[M]==ZeroMatrix[4])
SixVectorQ[v_]:=(Dimensions[v]=={6})


(* ::Subsection:: *)
(*Converting between skew-symmetric 4\[Cross]4 matrices and vectors in \[DoubleStruckCapitalR]^6:*)


ToFourSkewSymMatrix[v_]:={{0,v[[1]],v[[2]],v[[3]]},{-v[[1]],0,v[[4]],v[[5]]},{-v[[2]],-v[[4]],0,v[[6]]},{-v[[3]],-v[[5]],-v[[6]],0}}
ToSixVector[M_]:={M[[1,2]],M[[1,3]],M[[1,4]],M[[2,3]],M[[2,4]],M[[3,4]]}


(* ::Section:: *)
(*Generating symbolic input names*)


StringToSixVector[s_String]:=
Map[ToExpression[StringJoin[s,#]]&,{"12","13","14","23","24","34"}]


StringToFourVector[s_String]:=
Map[ToExpression[StringJoin[s,#]]&,{"1","2","3","4"}]


(* ::Section:: *)
(*Projective equivalence*)


EquivalentQ[a_List,b_List]:=
(ExtOuter[a,b]== Table[0,Length[a],Length[b]])/;(Length[a]== Length[b])
EquivalentQ::usage="EquivalentQ[a,b] returns True if the vectors 
a and b are projectively equivalent.";


(* ::Chapter:: *)
(*1. Coordinates on Subscript[\[DoubleStruckCapitalG], 2](\[DoubleStruckK]^4)*)


(* ::Subchapter:: *)
(*Data types for lines in \[DoubleStruckCapitalP]^3(\[DoubleStruckK]) *)


(* ::Section:: *)
(*Identify  \[CapitalLambda]^4(\[DoubleStruckK]^4)  with skew-symmetric 4x4 matrices*)


(* ::Subsection:: *)
(*The inner product on so(4) *)


(* ::Subsubsection:: *)
(*Here is the usual inner product on so(4):*)


InnerProduct[M_,N_]:= Simplify[-(Tr[M . N]/2)]


(* ::Subsection:: *)
(*Determinant and Pfaffian on so(4) *)


(* ::Subsubsection:: *)
(*The symmetric bilinear form corresponding to the Pfaffian*)


PF[M_,N_]:=M[[1,4]]N[[2,3]]+
M[[2,3]]N[[1,4]]+
M[[3,1]]N[[2,4]]+
M[[2,4]]N[[3,1]]+
M[[1,2]]N[[3,4]]+
M[[3,4]]N[[1,2]]


Pf[M_]:= Simplify[PF[M,M]/2]


(* ::Subsubsection:: *)
(*The square of the Pfaffian is the determinant*)


(* ::Subsection:: *)
(*The Plucker involution relates the nondegenerate symmetric bilinear forms given by the inner product and the Pfaffian bilinear form*)


(* ::Subsubsection:: *)
(*The Plucker involution on both vectors and matrices*)


PluckerInvolution[v_]:=
{v[[6]],-v[[5]],v[[4]],v[[3]],-v[[2]],v[[1]]}/;SixVectorQ[v]
PluckerInvolution[M_]:=
ToFourSkewSymMatrix[PluckerInvolution[ToSixVector[M]]]


(* ::Subsubsection:: *)
(*This linear mapping is an involution:*)


(* ::Subsubsection:: *)
(*This linear mapping is an isometry:*)


(* ::Subsubsection:: *)
(*This linear mapping relates the bilinear forms corresponding to the Pfaffian and the standard inner product on so(4):*)


(* ::Section:: *)
(*Products of vectors*)


(* ::Subsection:: *)
(*Given vectors v,w \[Element] \[DoubleStruckK]^4, their exterior outer product v\[And]w \[Element] \[CapitalLambda]^2\[DoubleStruckK]^4 \[LongLeftRightArrow] so(4)*)
(*represents the 2-plane \[LeftAngleBracket]v,w\[RightAngleBracket] \[Subset] \[DoubleStruckK]^4 they span.*)


ExtOuter[v_,w_]:=Outer[Times,v,w]-Outer[Times,w,v]


(* ::Subsubsection:: *)
(*The skew-adjoint endomorphism of \[DoubleStruckK]^4 defined by the exterior outer product v\[And]w has image the span \[LeftAngleBracket]v,w\[RightAngleBracket] \[Subset] \[DoubleStruckK]^4 and kernel the orthogonal complement \!\(\*SuperscriptBox[\(\[LeftAngleBracket]v, w\[RightAngleBracket]\), \(\[Perpendicular]\)]\). *)


(* ::Subsection:: *)
(*The projective class of the the exterior outer product v\[And]w depends only on  \[LeftAngleBracket]v,w\[RightAngleBracket] :*)


(* ::Subsection:: *)
(*The Pfaffian bilinear form applied to a pair of exterior products*)
(*is the determinant of the corresponding quadruple of vectors*)


(* ::Subsection:: *)
(*Plucker product of vectors*)


PluckerProduct[v_,w_]:=PluckerInvolution[ExtOuter[v,w]]


(* ::Subsection:: *)
(*The ternary Orth product*)


(* ::Subsubsection:: *)
(*This function is the "generalized cross product" which takes three 4-vectors*)
(*(for example, the homogeneous coordinates of a triple of points in \[DoubleStruckCapitalP]^3 in general position) and returns a fourth 4-vector orthogonal to the triple (for example, the homogeneous coordinates of the plane spanned by that triple).*)


Orth[a_,b_,c_]:= Det[{a,b,c,{i,j,k,l}}]/.{i->{1,0,0,0},j->{0,1,0,0},k->{0,0,1,0},l->{0,0,0,1}}


(* ::Chapter:: *)
(*2. Geometry of \[DoubleStruckCapitalP]^3(\[DoubleStruckK])*)


(* ::Text:: *)
(*We define datatypes for points, lines and planes in \[DoubleStruckCapitalP]^3(\[DoubleStruckK]). Points are projective equivalence classes of nonzero vectors in \[DoubleStruckK]^4and planes are projective equivalence classes of nonzero covectors in the dual vector space SuperStar[(\[DoubleStruckK]^4)].*)
(*Both are represented by lists of 4 scalars.  Lines are projective equivalence classes of rank 2 *)
(*skew-symmetric 4x4 matrices (nonzero singular elements of \[GothicS]\[GothicO](4,\[DoubleStruckK])).*)


(* ::Section:: *)
(*Points and Planes*)


PlaneThrough[Point[a_],Point[b_],Point[c_]]:=Plane[Orth[a,b,c]]
PlaneThrough::usage="PlaneThrough[A,B,C] is the plane containing the points A,B and C.";


IntersectionLine[Plane[a_],Plane[b_]]:=Line[PluckerInvolution[ExtOuter[a,b]]]
IntersectionLine::usage="IntersectionLine[A,B] is the line of intersection of the planes A and B.";


(* ::Section:: *)
(*Lines*)


LineThrough[Point[a_],Point[b_]]:=Line[ExtOuter[a,b]]
LineThrough::usage="LineThrough[A,B] is the line containing the points A and B.";


(* ::Section:: *)
(*Incidence routines*)


(* ::Subsection:: *)
(*Points and Planes*)


PointOnPlaneQ[Point[a_],Plane[b_]]:=(b.a==0)
PointOnPlaneQ::usage="PointOnPlaneQ[Point[a],Plane[b]] returnes True if the 
Point having homogeneous coordinates a(a 4-vector) lines on the plane whose
homogeneous coordinates form the covector b also a 4-vector).";


(* ::Subsection:: *)
(*Lines, Points, and Planes*)


PointOnLineQ[Point[a_],Line[b_]]:=(PluckerInvolution[b].a==ZeroVector[4])
PointOnLineQ::usage="PointOnLineQ[Point[a],Line[b]]returns True if
the Point having homogeneous coodinates a (4-vector) lies on the Line having
Plucker coordinates b (a 4\[Cross]4 skew-symmetric matrix).";


LineOnPlaneQ[Line[a_],Plane[b_]]:=(a.b==ZeroVector[4])
LineOnPlaneQ[Plane[b_],Line[a_]]:=LineOnPlaneQ[Line[a],Plane[b]]
LineOnPlaneQ::usage="LineOnPlaneQ[Line[a],Plane[b]]returns True if
the Line having Plucker coordinates a (a 4\[Cross]4 skew-symmetric matrix) 
lies on the Plane having homogeneous coordinate b a (4-vector).";


SkewQ[Line[a_],Line[b_]]:= (Chop[PF[a,b]] != 0)
SkewQ::usage="SkewQ[Line[a],Line[b]] returns True if the lines having
respective Plucker coordinates a, b (skew-symmetric 4\[Cross]4 matrices) are skew in \!\(\*SuperscriptBox[\(\[DoubleStruckCapitalP]\), \(3\)]\).";


(* ::Section:: *)
(*Intersecting Lines and Planes*)


IntersectionPoint[Plane[a_],Line[b_]]:=If[LineOnPlaneQ[Line[a],Plane[b]]==False,
Return[Point[b.a]],
Return[ ]]
IntersectionPoint::usage="IntersectionPoint[Plane[a],Line[Mb] returns the 
point of intersection of Plane[a] and Line[b]. This point exists only if 
LineOnPlaneQ[Line[a],Plane[b] tests False.";


(* ::Section:: *)
(*Preliminaries*)


PointList= {
Origin= Point[{0.,0.,0.,1.}], (* the point (0,0,0) *)
UnitX=Point[{1.,0.,0.,1.}],    (* the point (1,0,0) *)
UnitY=Point[{0.,1.,0.,1.}],    (* the point (0,1,0) *)
UnitZ=Point[{0.,0.,1.,1.}],    (* the point (0,0,1) *)
XIdealPoint=Point[{1.,0.,0.,0.}],
YIdealPoint=Point[{0.,1.,0.,0.}],
ZIdealPoint=Point[{0.,0.,1.,0.}]
};
PlaneList={
IdealPlane= Plane[{0.,0.,0.,1.}], (* the plane W = 0 (at infinity) *)
XYPlane=Plane[{0.,0.,1.,0.}],       (* the plane z = 0 *)
YZPlane=Plane[{1.,0.,0.,0.}],        (* the plane x = 0 *)
ZXPlane=Plane[{0.,1.,0.,0.}]             (* the plane y = 0 *)
};
CoordinateAxes={xAxis,yAxis,zAxis} =
Map[LineThrough[Origin,#]&,
{UnitX,UnitY,UnitZ}];
IdealCoordinateLines={
IdealXLine=LineThrough[YIdealPoint,ZIdealPoint],
IdealYLine=LineThrough[ZIdealPoint,XIdealPoint],
IdealZLine=LineThrough[XIdealPoint,YIdealPoint]};
IdealPlane=Plane[{0.,0.,0.,1.}];
XLevel[x_]:=Plane[{1.,0.,0.,-x}];
YLevel[y_]:=Plane[{0.,1.,0.,-y}]; 
ZLevel[z_]:=Plane[{0.,0.,1.,-z}]; 


(* ::Section:: *)
(*Projective Transformations*)


LineProjectiveTransform[Line[SkewSym_],A_]:=Line[A.SkewSym.Transpose[A]];
