(*
Package for Clifford Algebra calculus in Mathematica                       
Authors: Josep Manel Parra Serra,  email: jmparra@hermes.ffn.ub.es
         Lloren\c{c} Rosell\'{o} Saur\'{\i},  email:lrosello@goliat.upc.es 
         

Address: Departament de F\'{\i}sica Fonamental
         Universitat de Barcelona
         Diagonal 647  
         E-08028  Barcelona (Spain)
         FAX: 34 3 402 11 49
*)

BeginPackage["Clifford`"]

CliffordAlgebra::usage =
"CliffordAlgebra[a0^2,a1^2,...,an^2] defines the values of the
squares of the orthogonal basis vectors that generate the algebra.
Alias=CA."

form::usage = "gives the squares of the orthogonal basis vectors,
that is the coefficients of the quadratic form."

CliffordQ::usage = "CliffordQ[cliffor] returns True if cliffor is an
element of the Clifford algebra. False if not."

CP::usage = "CP[v1, v2,...,vN] is the Clifford product of N factors
in the Clifford algebra previously defined with
DefineCliffordAlgebra."

NPart::usage = "NPart[v, n] gives the n-degree component of v. For
instance, the bivector component of v is given by NPart[v, 2]. Alias
Pu"

ComplexConjugate::usage = "ComplexConjugate[v] gives the complex
conjugate of v. One assumes that all functions are real. Alias CC."

RPart::usage = "RPart[v] gives the real part of v. Alias ReP"

IPart::usage = "IPart[v] gives the imaginary part of v. Alias ImP"

e::usage = "e[i1,...,in] is a the simple multivector in the Clifford
algebra that results of the product of n distinct generators in the
specified order."

CTeX::usage = "TeX[v], where 'v' is in the Clifford Algebra, is a
modification of the function TeXForm in order to present in a good
TeX form the basis multivectors."

CommonFactor::usage = "CommonFactor[v, list] extracts as common
factors those e[i1,...,ik] with  i1,...,ik contained in the 'list'.
CommonFactor[vector] extracts as common factors all the multivector
units. Alias CF."

CPower::usage = "CPower[v, n] is the nth. Clifford power of v. n
should be a positive integer."

Rev::usage = "Rev[v] gives the reversion of v."

GradeInvol::usage = "GradeInvol[v] gives the grade involution of v.
Aliases GI, MI"

GradeRev::usage = "GradeRev[v] gives the Clifford conjugate of v; it
is the composition of reversion and the grade involution."

Changei::usage = "Changei[v, i] changes e[i] to -e[i] in v."

Numeric::usage = "Numeric[vector, n] gives the numeric value of the
coefficients with a nth-order precision."

ExteriorProduct::usage = "ExteriorProduct[v1,..., vN] evaluates the
exterior product of the N factors in the Clifford algebra. Alias EP."

InnerProduct::usage = "IP[v1,..., vN] evaluates the interior or
metric-contracted product of the N Clifford factors. Alias MP (metric
product)."

Contract::usage = "Contract[i1, i2,..., iN] contracts the indexes
i1,i2,..., iN accordingly the defined Clifford algebra. It is mainly
an auxiliary function."

SubindexOn::usage = "With this switch SubindexForm[] produces a
subindexed expression for the basis multivectors e[i1, ..., in] ."

SubindexOff::usage = "NoSubindexForm[] cancels SubindexForm[]."


PreReadOn::usage = "DefinePreRead[] defines a $PreRead in such a way
that PC[vector1,...,vectorN] is equivalent to < <vector1,...,vectorN
> >"

PreReadOff::usage = "QuitPreRead[] cancels the $PreRead."

RegularSession::usage = "RegularSession[] evaluates the group product
table of the geometric units. It is necessary to evaluate the regular
matrix representation of any Clifford algebra element. With this
function the algebra considered is the generated by e[1],..., e[n]."

RegularSession0::usage = "RegularSession0[] evaluates the group
product table of the geometric units. It is necessary to evaluate the
regular matrix representation of any Clifford algebra element. With
this function the algebra considered is the generated by
e[0],e[1],..., e[n]."

RegularMatrix::usage = "RegularMatrix[v] gives the regular matrix
representation of v. It requires having effected RegularSession[] or
RegularSession0[] accordingly the Clifford algebra we are using.
Alias RM."

GroupTable::usage = "Gives the product table of the set of 2^n unit
monomials in the standard mathematica ordering (1, (e[0]),e[1], ... ,
e[n],(e[0,1])...e[1,2] ... e[(0),1,...n]."


MetricTensor::usage = "MetricTensor[{h0, h1,..., hn}, {q0, q1,..., qn}]
is used to introduce a non-Cartesian metric in arbitrary orthogonal
coordinates. h0, h1,..., hn are Lam‚'s coefficients (the positive
norm of the holonomic coordinate tangent vectors, and q0,q1,..., qn
are the names given to the coordinates. The signature of the space
has been given before with DefineCliffordAlgebra. Note that, as it
occurs with DefineCliffordAlgebra, the zero component must be
defined, even when no use of it is intended. In this case it suffices
to put h0=1, q0=q0"

Ricci::usage = "Ricci[i, j, k] is Ricci's rotation coefficient,
accordingly to the derivation formula e[i](e[j])= Ricci[i,j,k] e[k]."

Nabla::usage = "Nabla[v] evaluates the Clifford-nabla or
Clifford-gradient of v. It does not evaluate ExtCod[v] as Nabla2[v]
does."

ExtCod::usage = "ExtCod[v] gives a list where the first element is
the exterior differential and the second is the codifferential of v."

ExtD::usage = "ExtD[v] gives the exterior differential of v. It
accepts the output of ExtCod as argument."

Codif::usage = "Codif[v] gives the codifferential of v. It accepts
the output of ExtCod as argument."

Nabla2::usage = "Nabla2[v] evaluates the Clifford-nabla or
Clifford-gradient of v through ExtCod by addition of the two parts.
It accepts the output of ExtCod as argument."



$PreRead=
StringReplace[#,{"CA"->"CliffordAlgebra","Pu"->"NPart",
"CC"->"ComplexConjugate", "ReP"->"RPart", "ImP"->"IPart",
"CF"->"CommonFactor","GI"->"GradeInvol","MI"->"GradeInvol",
"GR"->"GradeRev","EP"->"ExteriorProduct","MP"->"InnerProduct",
"RM"->"RegularMatrix","MT"->"MetricTensor"}]&

Begin["`Private`"]

Unprotect[IntegerQ, MatchQ, Expand]
SetAttributes[IntegerQ, Listable]
SetAttributes[MatchQ, Listable]
SetAttributes[Expand, Listable]
Protect[IntegerQ, MatchQ, ExpandQ]

CliffordQ[element_] := Module[{x = Expand[element]},
        If[x == 0 || NumberQ[x], Return[True]];
        If[x[[0]] == e || x[[0]] == Times, x = {x}];
        And @@ MatchQ[List @@ x, y_. e[p___]
                 /; ((And @@ IntegerQ[{p}]) && FreeQ[y, e[q___]])]]

Dif[S_, T_] := Complement[Union[S, T], Intersection[S, T]]

CliffordAlgebra[b__] := (form = List[b];)

M[a_. e[S___],b_. e[T___]] :=
  Times[a, b, Signature[{S}], Signature[{T}],
       Apply[Times, Flatten[Outer[If[#1 <= #2,1,-1]&, Sort[{S}], 
       Sort[{T}]]]],
       Apply[Times, Map[form[[1+#]]&, Intersection[{S}, {T}]]],
       Apply[e, Dif[{S}, {T}]]]

M[0, e[S___] b_.] := 0
M[e[S___] b_., 0] := 0
     (* The supression of the following three lines accelerates the
calculations *)
M[e[S___] b_., n_?NumberQ]:= n b Signature[{S}] e @@ Sort[{S}]
M[n_?NumberQ, e[S___] b_.]:= n b Signature[{S}] e @@ Sort[{S}]
M[n_?NumberQ, m_?NumberQ]:= n m

pc[x_, y_] := Distribute[M[x,y]]

CP[x_?CliffordQ, y__?CliffordQ] := Fold[pc, Expand[x], Expand[{y}]]

Attributes[CP] = {Flat}

Suprimir[a_. e[x___], n_] := a If[Length[{x}] == n, e[x], 0]
NPart[vector___, n_] := Distribute[Suprimir[vector, n]]

ComplexConjugate[v_] := v /. Complex[x_,y_] -> x - I y 
RPart[v_] := Expand[(v + ComplexConjugate[v])/2]     
IPart[v_] := Expand[(v - ComplexConjugate[v])/(2I)]

CTeX[x_] :=
        TeXForm[x /. {e[y__] -> Subscripted[e[SequenceForm[y]]], 
        e[] -> 1}] /. e -> "{\bf e}"

CommonFactor[x_, y_List] :=
Collect[x, Apply[e, Distribute[{{}, {#}}& /@ y, List, List, 
        List, Union],1]]

CommonFactor[x_] := CommonFactor[x, Range[Length[form]]-1]

pot[x_,n_]:= If[n == 0, 1, x]
Attributes[pot] = {Listable}

CPower[x_,n_Integer] :=
        Fold[CP[#1,#1, #2]&, 1, pot[x, IntegerDigits[n, 2]]]
CPower[x_,0]:=e[]

gira[x_. e[y___]]:= x e @@ Reverse[{y}]
Rev[x_] := CP[e[],Distribute[gira[x]]]

ap[x_. e[y___]] := (-1)^Length[{y}] x e[y]
GradeInvol[x_]:= CP[e[],Distribute[ap[x]]]

GradeRev[x_]:= GradeInvol[Rev[x]]

Unprotect[MemberQ]
MemberQ[{}, 0] := True
Protect[MemberQ]
can[e[x___] a_., n_] := If[MemberQ[{x}, n], -a e[x], a e[x]]
can[a_?NumberQ, n_]:= a
Changei[x_, n_] := Distribute[can[x, n]]
Changei[x_List, n_]:= Changei[#,n]& /@ x

Attributes[Rev] = Listable
Attributes[GradeInvol] = Listable
Attributes[GradeRev] = Listable

num[a_. e[x___], n_Integer]:= N[a, n] e[x]
Numeric[x_, n_Integer]:= Distribute[num[x, n]]

pext[a_. e[i__], b_. e[j__]] := a b Signature[{i,j}] e @@ Sort[{i,j}]
pext[a_. e[], b_. e[i__]]:=  a b e[i]
pext[a_. e[i__], b_. e[]]:=  a b e[i]
pext[a_. e[], b_. e[]] := a b e[]
pext[0, a_] := 0
pext[a_, 0] := 0
prodext[x_, y_] := Distribute[pext[x, y]]
ExteriorProduct[x_?CliffordQ, y__?CliffordQ] := 
        Fold[prodext, Expand[x], Expand[{y}]]

pint[a_. e[i___], b_. e[]] := 0
pint[b_. e[], a_. e[i___]] := 0
pint[0, a_] := 0
pint[a_, 0] := 0
pint[a_. e[i__], b_. e[j__]] :=
        If[Or @@ Map[MemberQ[{i}, #]&, {j}], a b CP[e[i], e[j]], 0]
prod[x_, y_] := Distribute[pint[x,y]]
InnerProduct[x_?CliffordQ, y__?CliffordQ] := Fold[prod, Expand[x], 
             Expand[{y}]]

Contract[x___] := CP @@ Apply[e, {#}& /@ {x}, 1]

SubindexOn[]:= (Format[e[x__]]:= Subscripted[e[SequenceForm[x]]])

SubindexOff[]:= (Format[e[x__]]=.)


PreReadOn[]:=
        ($PreRead = StringReplace[#, {"**" -> "~CP~"}]&);
        
(* Just another alternative to write the Clifford product: 
   < <...,...,...> >
PreReadOn[]:=
        ($PreRead = StringReplace[#, {"< <" -> "CP[", "> >" -> "]"}]&);
*)
PreReadOff[]:= ($PreRead =.)

RegularSession[]:= (
  GroupTable = Module[
  {lo = Apply[e, Sort[Distribute[{{}, {#}}& /@ Range[Length[form]-1], 
   List, List, List, Union]],1],
   signs, mat},
   signs = Inner[CP, lo, lo, List] /. e[] -> 1;
   mat = Transpose[Outer[CP, lo, lo]];
   mat = Inner[Times[#1,#2]&, signs, mat, List]];)

RegularSession0[]:= (
  GroupTable = Module[
  {lo = Apply[e, Sort[Distribute[{{}, {#}}& /@ Range[0,Length[form]-1], 
  List, List, List, Union]],1],
   signs, mat},
   signs = Inner[CP, lo, lo, List] /. e[] -> 1;
   mat = Transpose[Outer[CP, lo, lo]];
   mat = Inner[Times[#1,#2]&, signs, mat, List]];)

RegularRep[e[x___]]:= (GroupTable /. e[x] -> 1) /. e[y___] -> 0
RegularRep[a_. e[x___]]:= a RegularRep[e[x]]

RegularMatrix[x_]:= RegularRep[CP[e[], x]]
RegularMatrix[x_Plus]:= RegularRep /@ CP[e[], x]

MetricTensor[{h__},{v__}]:= Module[{n},
                n = Length[{v}]-1; variables := {v}; ax := {h};
                subst:=Table[variables[[i+1]]->x[i],{i,0,n}];
                substinv:=Table[x[i]->variables[[i+1]],{i,0,n}];
                dim = Length[{v}];
                ]

lame[i_Integer]:= ax[[i+1]] /. subst;

g[i_, j_, k_] := If[i != j, If[j != k, 0], g[i, i, k]]
g[i_, j_, j_] := 0
g[i_, i_, i_] := 0
g[i_,j_,i_] := D[lame[i],x[j]]/(lame[i] lame[j])
g[i_,i_,j_] := -form[[i+1]] form[[j+1]] D[lame[i],x[j]]/(lame[i] lame[j])

Ricci[i_Integer,j_Integer,k_Integer] := g[i,j,k] /. substinv

d[n_Integer]:= Sum[eD[i,n],{i,0,n-1}]

cder[eD[i_, n_], 0]:= 0

cder[eD[i_,n_], A_. e[]] := D[A,x[i]]/lame[i] e[]

cder[eD[i_,n_], A_. e[j_]] :=
 Sum[A g[i,j,k] e[k],{k,0,n-1}] + 
 (D[A,x[i]]/lame[i]) e[j]


cder[eD[i_,n_], A_. e[j_, k_]] :=
   CP[cder[eD[i,n], e[j]],A e[k]] +  CP[A e[j], cder[eD[i,n], e[k]]] +
  (D[A,x[i]] /lame[i]) e[j,k]

cder[eD[i_,n_], A_. e[u_, v_, w__]] :=
         CP[ cder[eD[i,n], e[u]],A e[v,w]] +
         CP[A e[u],cder[eD[i,n],e[v]],e[w]] +
         CP[A e[u],e[v],cder[eD[i,n],e[w]]] +
        (D[A,x[i]]/lame[i]) e[u,v,w]

nabla[eD[i_,n_],vector_]:=form[[i+1]] CP[e[i], cder[eD[i,n],vector]]

Nabla[vector_] := Module[{operator}, operator= d[dim];
        Distribute[nabla[operator, vector /. subst]]
        ] /. substinv

f[a_. e[x___], n_Integer]:= If[Length[{x}] > n, {a e[x],0}, {0, a e[x]}]

s[v_, n_Integer]:= Distribute[f[v, n]]

ss[a_. e[x___], d_]:= s[Nabla[a e[x]], Length[{x}]]

ExtCod[x_]:= Distribute[ss[x, dim]]

ExtD[x_] := ExtCod[x][[1]]
ExtD[x_List] := x[[1]]

Codif[x_]:= ExtCod[x][[2]]
Codif[x_List] := x[[2]]

Nabla2[x_]:= Plus @@ ExtCod[x]
Nabla2[x_List]:= Plus @@ x           

End[]

EndPackage[]