(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



BeginPackage["GREATER2`"]

(* 

Written by Tom Hartman. Uses code from GREAT.m by Tristan Hubsch, EinsteinTensor.m by Pekka Janhunen, and diffgeo.m by Matt Headrick. 

Last updated June 2012.

*)







S = Simplify[#,TimeConstraint->1]&;
SMF = Simplify[MatrixForm[#], TimeConstraint->1]&;
MF = MatrixForm;
GSimplify = (#1&);


datapoints = {};
Sto[datapoint_]:= datapoints=Join[datapoints, {datapoint}];



Schwarzian[w_, z_]:= D[w,{z,3}]/D[w,z]-3/2 (D[w,{z,2}]/D[w,z])^2;
(* 3-pt function in 2d cft *)
cft3[h_, z_]:=(z[[1]]-z[[2]])^(-h[[1]]-h[[2]]+h[[3]]) (z[[2]]-z[[3]])^(-h[[2]]-h[[3]]+h[[1]]) (z[[3]]-z[[1]])^(-h[[3]]-h[[1]]+h[[2]])


IMetric::usage = "IMetric[g], with g an n.n-matrix (two lower indices),
  returns the inverse metric (two upper indices)."

Christoffel::usage = "Christoffel[g,x], with g a n.n-matrix and x
  n-vector of coordinates, gives the Christoffel symbol of the 2nd
  kind (1st upper, two lower indices)." 

Riemann::usage = "Riemann[g,x], with g a n.n-matrix and x n-vector of
  coordinates, gives the Riemann tensor (1st upper, three lower
  indices)."

Ricci::usage = "Ricci[g,x], with g a n.n-matrix and x n-vector of
  coordinates, gives the Ricci tensor (two lower symmetric indices)."

SCurvature::usage = "SCurvature[g,x], with g a n.n-matrix and x
  n-vector of coordinates, gives the Scalar Curvature."

EinsteinTensor::usage = "EinsteinTensor[g,x] with g a nxn-matrix
  (the metric with lower indices) and x n-vector (the coordinates)
  gives the Einstein tensor (a nxn-matrix) with lower indices."

SqRicci::usage = "SqRicci[g,x], with g a n.n-matrix and x
  n-vector of coordinates, gives the norm-square of the Ricci tensor."

SqRiemann::usage = "SqRiemann[g,x], with g a n.n-matrix and x
  n-vector of coordinates, gives the norm-square of the Riemann tensor."

helpGREAT::usage = "GREAT functions are: IMetric, Christoffel,
  Riemann, Ricci, SCurvature, EinsteinTensor, SqRicci, SqRiemann."


Simplif1 = Simplify[#1, TimeConstraint->1] &;
IMetric[metric_] := Simplif1[Inverse[metric]];

Christoffel[metric_,x_]:= Christoffel[metric,x] = 
  Module[{Dim, iMet, PreChristoffel, Christoffel, i, j, k},
         Dim = Length[x];
         iMet = IMetric[metric];
           (* Metric with upper indices *)
         PreChristoffel =
           Table[ D[metric[[k,i]],x[[j]]]
                + D[metric[[j,k]],x[[i]]]
                - D[metric[[i,j]],x[[k]]],
          	   {k,Dim}, {i,Dim}, {j,Dim} ];
          	 (* The \{k,ij\} Christoffel symbols *)
         PreChristoffel = Simplif1[PreChristoffel];
            (* The \Gamma^k_{ij} Christoffel symbols *)
         Christoffel = (1/2) iMet . PreChristoffel;
         (* Return the Christoffel symbol: *)
         Simplif1[Christoffel] ]

Riemann[metric_,x_]:=
  Module[ {Dim, iMet, ChrisSymbol, Riemann, PreRiemann,
          a, b, c, i, j, k},
          Dim = Length[x];
          iMet = IMetric[metric];
          ChrisSymbol = Christoffel[metric,x];
          PreRiemann = 
             Table[ D[ChrisSymbol[[a,i,c]],x[[b]]]
                  + Sum[ChrisSymbol[[k,i,c]]
                        * ChrisSymbol[[a,k,b]],
                        {k,Dim} ],
                    {a,Dim}, {i,Dim}, {b,Dim}, {c,Dim} ];
           	(* Riemann tensor (1st upper, three lower indices)
           	   is antisymmetrized PreRiemann: *)
          Riemann = Table[ PreRiemann[[a,i,b,c]]
                         - PreRiemann[[a,i,c,b]],
                           {a,Dim}, {i,Dim}, {b,Dim}, {c,Dim} ];
          (* Return the Riemann tensor: *)
          Simplif1[Riemann] ]

Ricci[metric_,x_]:=
  Module[ {Dim, Riem, Ricci, a, i, j},
          Dim = Length[x];
          Riem = Riemann[metric,x];
          Ricci = Table[ Sum[Riem[[a,i,a,j]], {a,Dim}],
                         {i,Dim}, {j,Dim} ];
          (* Return the Ricci tensor (two lower indices): *)
          Simplif1[Ricci] ]

SCurvature[metric_,x_]:=
  Module[ {Dim,iMet,CurvatureScalar,i,j},
          Dim=Length[x];
          iMet=IMetric[metric];
          CurvatureScalar=Tr[iMet.Ricci[metric,x]];
    (*Return Scalar Curvature:*)Simplif1[CurvatureScalar]]

EinsteinTensor[metric_,x_]:=
  Simplif1[Ricci[metric,x] - (1/2) SCurvature[metric,x] metric];

SqRicci[metric_,x_]:=
  Module[ {Dim, iMet, Ric, RRicci, i, j, k, l},
          Dim = Length[x];
          iMet = IMetric[metric];
          Ric = iMet.Ricci[metric, x];
          RRicci = Tr[Ric.Ric];
          (*Return norm - square of Ricci tensor :*)Simplif1[RRicci]]

SqRiemann[metric_, x_] := Module[{Dim, iMet, Riem,
   RRiem, i, j, k, l, m, n, p, q}, Dim = Length[x];
    iMet = IMetric[metric];
    Riem = Riemann[metric, x].iMet;
    RRiem = Sum[Riem[[i, j, l, k]]*
    Riem[[l, k, i, j]], {i, Dim}, {j, Dim}, {k, Dim}, {l, Dim}];
    (*Return norm - square of Riemann tensor :*)Simplif1[RRiem]]

helpGREAT:= Print["GREAT functions are: IMetric, Christoffel,
  Riemann, Ricci, SCurvature, EinsteinTensor, SqRicci, SqRiemann."]







ds2met::usage = "ds2met[ds2,dx] convert \!\(\*SuperscriptBox[\(ds\), \(2\)]\) expression to \!\(\*SubscriptBox[\(g\), \(\[Mu]\[Nu]\)]\) matrix."
ds2met[ds2_,dx_]:=
Module[{dim,ds2metdouble,icounter,jcounter, xds2, res},
dim = Length[dx];
xds2 = Expand[ds2];
res =
Table[Coefficient[xds2,dx[[icounter]] dx[[jcounter]]]/(2-KroneckerDelta[icounter,jcounter]),
{icounter,dim},{jcounter,dim}];
S[res]
];

Metric::usage = ToString[
"Metric[\!\(\*
StyleBox[\"ds2\",\nFontSlant->\"Italic\"]\), {\!\(\*SubscriptBox[
StyleBox[\"x\",\nFontSlant->\"Italic\"], \(1\)]\), \!\(\*SubscriptBox[
StyleBox[\"x\",\nFontSlant->\"Italic\"], \(2\)]\), ...}] returns the metric tensor specified by the line element \!\(\*
StyleBox[\"ds2\",\nFontSlant->\"Italic\"]\) with coordinates \!\(\*SubscriptBox[
StyleBox[\"x\",\nFontSlant->\"Italic\"], \(i\)]\). (Just like ds2met but does not require you to define differentials.)", StandardForm];
Metric[ds2_, X_]:=Module[{dX},
dX = Differential[X];
ds2met[ds2,dX]
];

Differential[X_List]:=Module[{},
	If[!MatchQ[Union[Head/@X], {Symbol}], Message[GREATER::invalidCoordinate,X]];
	Symbol["d"<>SymbolName[#]]&/@X
];


(* Forms to tensors; use \[Vee] = Vee as a abstract wedge product, Wedge as a component-wise wedge product*)
TFormQ[f_, dX_]:=Not[NotTFormQ[f,dX]];
NotTFormQ[f_, dX_]:=And@@(FreeQ[f, #]&/@dX);
TFormExRules[dX_] := {
HoldPattern[Vee[a_,b_+c_]]:> Vee[a,b] + Vee[a,c],
HoldPattern[Vee[a_+b_,c_]]:> Vee[a,c] + Vee[b,c],
HoldPattern[Vee[n_?(NotTFormQ[#,dX]&) * b_, c_]]:> n Vee[b,c],
HoldPattern[Vee[c_, n_?(NotTFormQ[#,dX]&) * b_]]:> n Vee[c,b]
};
formToTensor[form_, dX_, n_]:=Module[{dim,xform, res},
dim = Length[dX];
xform = Expand[form]//.TFormExRules[dX];
res = Table[Coefficient[xform, Vee[dX[[ia]],dX[[ib]]]] - Coefficient[xform, Vee[dX[[ib]],dX[[ia]]]], {ia,dim}, {ib,dim}];
res
]/;n == 2;


Unprotect[TensorProduct];
TensorProductT[a__]:= Outer[Times, a];
Protect[TensorProduct];





antisymmetrize[X_]:= Module[{perms, n, numperm, res,ia},
n = ArrayDepth[X];
perms = Permutations[Range[n]];
numperm = Length[perms];
res = Sum[
Transpose[X,perms[[ia]]]Signature[perms[[ia]]],
{ia,numperm}
]/numperm;
res
];
symmetrize[X_]:= Module[{perms, n, numperm, res,ia},
n = ArrayDepth[X];
perms = Permutations[Range[n]];
numperm = Length[perms];
res = Sum[
Transpose[X,perms[[ia]]],
{ia,numperm}
]/numperm;
res
];




IndexAtStart[T_, ind_]:=
Transpose[T,
(*Join[{ind},Range[1, ind-1], Range[ind+1, ArrayDepth[T]]]*)
Join[Range[2, ind],{1}, Range[ind+1, ArrayDepth[T]]
]
];
Contract[T1_, T2_, ind1_, ind2_]:=
Inner[Times, T1, IndexAtStart[T2, ind2], Plus, ind1];
SelfContract[T1_, ind1_, ind2_]:=
Tr[IndexAtStart[IndexAtStart[T1,ind1],ind2], Plus, 2]/;ind1<ind2;

Raise[T_, i_, Guu_]:=IndexMove[Contract[T, Guu, i, 1],Length[Dimensions[T]], i] ;


SetAttributes[CircleTimes,HoldAll];
CircleTimes[Part[T1_, ind1__], Part[T2_, ind2__]]:=
 Module[{sumind, res, ii, ilist1, ilist2},
ilist1 = List[ind1];
ilist2 = List[ind2];
sumind = Intersection[ilist1,ilist2];
res= Outer[Times, T1, T2];
Do[res = SelfContract[res,Position[ilist1, ii][[1,1]], Length[ilist1] + Position[ilist2,ii][[1,1]]],
{ii, sumind}];
res];
EpsilonTensor[4]:= Table[Signature[{ia,ib,ic,id}], {ia,4}, {ib,4}, {ic,4}, {id,4}];
EpsilonTensor[5]:= Table[Signature[{ia,ib,ic,id,ie}], {ia,5}, {ib,5}, {ic,5}, {id,5}, {ie,5}];
Epsilon[args_]:= Signature[args];
TraceFirst[T1_]:= SelfContract[T1, 1, 2];
(*IndexShift[T_, istart_, ishift_]:=
Transpose[T,
Join[Range[1, istart-1], Range[istart+1, istart+ishift], {istart}, Range[istart+ishift+1, ArrayDepth[T]]]];*)
IndexShift[T_, istart_, ishift_]:=
Transpose[T,
Join[Range[1, istart-1], {istart+ishift}, Range[istart, istart+ishift-1], Range[istart+ishift+1, ArrayDepth[T]]]];
IndexMove[components_List, from_, to_]:= Transpose[components, Join[Range[1,from-1], {to}, Range[from, to-1], Range[to+1,ArrayDepth[components]]]]/;to>from;
IndexMove[components_List, from_, to_]:= Transpose[components, Join[Range[1,to-1], Range[to+1,from], {to}, Range[from+1, ArrayDepth[components]]]]/;to<from;
IndexMove[components_List, from_, to_]:= components/;to==from;

ShowComponents[T_, X_, ind_, name_:"Y"]:=Module[{Tflat, irep, istr},
Tflat = Flatten[T];
For[ii = 1, ii<= Length[Tflat], 
If[UnsameQ [Tflat[[ii]],  0],
irep = IntegerDigits[ii-1, Length[X]];
irep = PadLeft[irep, Length[ind]]+1;
irep = (Part[X, #1]&)/@ irep;
istr = name;
For[ij = 1, ij <= Length[ind], 
whichscript = If[ind[[ij]]==-1, Subscript, Superscript];
istr = whichscript[istr, irep[[ij]]];
ij = ij + 1
];
istr = istr/.{Subscript[Subscript[a1_,b1_],c1_]:> Subscript[a1,ToString[b1]<>ToString[c1]]};
istr = istr/.{Superscript[Superscript[a1_,b1_],c1_]:> Superscript[a1,ToString[b1]<>ToString[c1]]};
Print[istr , "=", Tflat[[ii]]]

];
ii = ii + 1
];
];


(* Converts, say, dt + r dr to a component p-form {1, r} *)
pForm[f_, dX_]:= Module[{ia},
Table[Coefficient[f, dX[[ia]]], {ia, Length[dX]}]];

(*HodgeStar[wdd_, Gdd_, Guu_]:=Module[{n, detg, dims, wuu},
n = Length[Dimensions[w]];
detg = Sqrt[-Det[Gdd]];
dims = Length[Gdd];
wuu = Guu.wdd.Guu//S;
res = 1/Factorial[n]detg Table[Sum[Signature[{ia,ib,ic,id,ie}]wuu[[id,ie]], {id,dims}, {ie,dims}], {ia,dims}, {ib,dims}, {ic,dims}];
res]/;Length[Dimensions[wdd]]\[Equal]2;
HodgeStar::usage = "Provide upper-indexed form, get back lower-indexed dual form";
HodgeStar[wUPPER_, Gdd_, Guu_]:=Module[{n, detg, dims, ep},
n = Length[Dimensions[wUPPER]];
detg = Sqrt[-Det[Gdd]];
dims = Length[Gdd];
ep = EpsilonTensor[dims];
res =1/Factorial[n]detg Contract[ep, wUPPER, Length[wUPPER 
res =  Table[Sum[Signature[{ia,ib,ic,id,ie}]wuu[[id,ie]], {id,dims}, {ie,dims}], {ia,dims}, {ib,dims}, {ic,dims}];
res]/;Length[Dimensions[wUPPER]]\[Equal]2;
*)

(* Make epsilon, without any Sqrt[g]*)
MakeEpsilon[n_]:=Module[{ilist, tableIters, ii,  compoTilde},
	ilist = Array[ii, {n}];
	tableIters = Table[{jj, n}, {jj,ilist}];
	compoTilde = Table[Signature[ilist], ##]&@@tableIters;
compoTilde
];





SimpleDeriv[v_, x_] :=Module[{i},
GSimplify[Table[D[v, x[[i]]], {i, Length[x]}]]];
SimpleDerivt[v_, x_] :=Module[{i},
Simplify[Table[Dt[v, x[[i]]], {i, Length[x]}]]];

PD[v_,x_]:=SimpleDeriv[v,x];


LieD[v_, T_, x_, ind_] :=
Module[{dv },
dv = SimpleDeriv[v,x];
 GSimplify[v.SimpleDeriv[T, x]
+ 
 Sum[
If[ind[[iind]]== 1,
-IndexShift[Contract[dv, T, 1, iind], 1, iind-1],
IndexShift[Contract[dv, T, 2, iind], 1, iind-1]
],
{iind, 1, Length[ind]}
]]
];
LieCom[V1_, V2_, X_]:= LieD[V1, V2, X, {1}];
CoD[T_, met_, x_, ind_, \[Omega]dab_:ERROR] :=
Module[{\[CapitalGamma] , p1},
(*Print["Cod0"];*)
\[CapitalGamma] = Christoffel[met, x];
p1 = SimpleDeriv[T, x];
(*Print["CoD2"];*)
GSimplify[p1+
Sum[
Switch[ind[[iind]],
1, IndexShift[Contract[\[CapitalGamma], T, 3, iind], 1, iind],
-1, -IndexShift[Contract[\[CapitalGamma], T, 1, iind], 2, iind-1],
2, IndexShift[Contract[\[Omega]dab, T, 3, iind], 2, iind-1],
-2, -IndexShift[Contract[\[Omega]dab, T, 2, iind], 2, iind-1]
],
{iind, 1, Length[ind]}
]]
];

SpinorCoD[\[Psi]_, Gdd_, X_, \[Omega]dbb_, \[Gamma]a_]:= Module[{\[Gamma]aa},
\[Gamma]aa = 1/2 Table[ \[Gamma]a[[ia]].\[Gamma]a[[ib]] - \[Gamma]a[[ib]]. \[Gamma]a[[ia]], {ia,Length[X]}, {ib,Length[X]}];
 SimpleDeriv[\[Psi], X] +spinsign 1/4 Table[Sum[ \[Omega]dbb[[iu,ia,ib]]*\[Gamma]aa[[ia,ib]], {ia,Length[X]}, {ib,Length[X]}], {iu,Length[X]}].\[Psi]
];
SpinorCoDBar[\[Psi]_, Gdd_, X_, \[Omega]dbb_, \[Gamma]a_]:= Module[{\[Gamma]aa},
\[Gamma]aa = 1/2 Table[ \[Gamma]a[[ia]].\[Gamma]a[[ib]] - \[Gamma]a[[ib]]. \[Gamma]a[[ia]], {ia,Length[X]}, {ib,Length[X]}];
 SimpleDeriv[\[Psi], X] -1/4 Table[Sum[ \[Omega]dbb[[iu,ia,ib]]*(\[Psi].\[Gamma]aa[[ia,ib]]), {ia,Length[X]}, {ib,Length[X]}], {iu,Length[X]}]
];

KillingEquations[k_,met_,x_]:=
Module[{delk},  
delk = CoD[k,met,x, {-1}];
Simplify[delk + Transpose[delk]]];



BasisFormToTensor[dx_List, dX_]:=Length[dx]*antisymmetrize[Outer[Times, ##]&@@ (OneFormToBasisVector[#, dX]&/@ dx)];
OneFormToBasisVector[dx_, dX_]:= Table[If[i==Position[dX, dx][[1,1]],1,0], {i, Length[dX]}];
Form[expr_, X_List]:=Module[{rank, p, dX, compo, sig},
	dX = Differential[X];
	p = Expand[expr];
	If[FreeQ[expr,Wedge], p=Wedge[p]];
	rank = Union[Length/@Cases[{p}, Wedge[___], Infinity,Infinity, Heads->True]];

	compo = p//.{
		Wedge[d___, c_+a_, b___]:> Wedge[d,c,b] + Wedge[d,a,b],
		Wedge[d___, c_ * a_, b___] :> c Wedge[d,a,b]/;And@@(FreeQ[c,#]&/@dX),
		Wedge[d___, c_ , b___] :> c Wedge[d,b]/;And@@(FreeQ[c,#]&/@dX)
		
	}/.{Wedge[f__?(MemberQ[dX,#]&)] :> BasisFormToTensor[List[f], dX]};
	sig = Table[-1, {ArrayDepth[compo]}];
	compo
];


Options[ChangeCoords] = {Inverse->False};

ChangeCoords[V_, X1_, X2_,trans_, ind_] := Module[{jac},
jac = TJacobian[X2, X1, trans];
Simplify[(jac.V)/.Trans2rules[X1,X2,trans]]]/;ind=={1};
ChangeCoords[V_, X1_, X2_,trans_, ind_] := Module[{jac},
jac = TJacobian[X1, X2, trans];
Simplify[(V.jac)/.Trans2rules[X1,X2,trans]]]/;ind=={-1};
ChangeCoords[g_, X1_, X2_, trans_, ind_,OptionsPattern[]] :=Module[{jac},
jac = TJacobian[X1, X2, trans];
S[(Transpose[jac].g.jac)/.Trans2rules[X1,X2,trans]]]/;ind=={-1,-1};
ChangeCoords[g_, X1_, X2_, trans_, ind_] :=Module[{jac},
jac = TJacobian[X1, X2, trans]TESTING;
Simplify[Contract[(Transpose[jac].g.jac), jac, 2, 1]/.Trans2rules[X1,X2,trans]]]/;ind=={-1,-1, -1};
ChangeCoords[g_, X1_, X2_, trans_, ind_] :=Module[{jac},
jac = TJacobian[X2, X1, trans];
Simplify[(jac.g.Transpose[jac])/.Trans2rules[X1,X2,trans]]]/;ind=={1,1};
ChangeCoords[g_, X1_, X2_, trans_, ind_] :=Module[{jac},
jac = TJacobian[X1, X2, trans];
Simplify[(Transpose[jac].g.Inverse[jac])/.Trans2rules[X1,X2,trans]]]/;ind=={-1,1};
TJacobian[X1_, X2_, trans_]:= Module[{ru, X1new,res},
ru = Trans2rules[X1, X2, trans];
X1new = X1/.ru;
res = Table[SimpleDeriv[ X1new[[ia]], X2], {ia, Length[X1]}];
Simplify[res]];
Trans2rules[X1_,X2_,trans_]:=Module[{},
Simplify[Solve[trans, X1][[1]]]];

Pullback[\[Omega]_, X_(*big space*), Y_(*hypersurface*), surf_, {-1,-1,-1}]:=Module[{J},
J = TJacobian[X, Y, surf];
Table[Sum[\[Omega][[ia,ib,ic]]J[[ia,iap]]J[[ib,ibp]]J[[ic,icp]], {ia,Length[X]}, {ib, Length[X]}, {ic,Length[X]}], {iap, Length[Y]}, {ibp, Length[Y]}, {icp, Length[Y]}]
];


(* 3d chern-simons, 1-form-valued, d=3*)
ChernSimons[w_, X_]:=Module[{res},
res = w \[Wedge]exterior[w,X] + 2/3 w\[Wedge]w\[Wedge]w ;
Simplify[res]
]/;Length[X]== 3&& Dimensions[w]== {3};

(* 3d chern-simons, matrix-1-form-valued. d=3 *)
ChernSimons[w_, X_]:=Module[{res},
res = TraceFirst[w \[Diamond]MatrixExterior[w,X] + 2/3 w\[Diamond]w\[Diamond]w ];
Simplify[res]
]/;Length[X]== 3&&Dimensions[w] != {3};


FieldStrength[A_, X_] := SimpleDeriv[A,X]-Transpose[SimpleDeriv[A,X]];

(*4d*)
MaxwellStressTensor[Fdd_, Gdd_]:=Module[{Guu, res, Fdu, Fuu, ND},
ND = Length[Gdd];
Guu = IMetric[Gdd];
Fdu = Fdd.Guu//S;
Fuu = Guu.Fdu//S;
res = -1/4 Gdd Sum[Fuu[[ia,ib]]Fdd[[ia,ib]], {ia,ND}, {ib, ND}] + Table[Sum[Fdu[[iu,ig]] Fdd[[iv,ig]], {ig,ND}], {iu,ND}, {iv,ND}];
res];

(*3d *)
CottonTensor[gdd_, X_, sig_:-1]:= (*assumes epsilon_uv = + sqrt(g); 3d*)Module[{Rdd, Rcurve, DdRdd, DdRgdd, paren, res, eddd, guu},
guu = IMetric[gdd];
Rdd = Ricci[gdd,X];
Rcurve = SCurvature[gdd,X];
DdRdd = CoD[Rdd, gdd, X, {-1,-1}];
DdRgdd = CoD[Rcurve*gdd, gdd, X, {-1,-1}];
paren = DdRdd - 1/4 DdRgdd;
eddd =  Sqrt[sig*Det[gdd]]Table[Signature[{ia,ib,ic}], {ia,3}, {ib,3}, {ic,3}];
res = Table[
Sum[eddd[[iu, ibprime, iaprime]]guu[[ia,iaprime]]guu[[ib,ibprime]]paren[[ia,iv,ib]], {ia,3}, {ib,3}, {iaprime,3}, {ibprime,3}], 
{iu, 3}, {iv, 3}];
Simplify[res,TimeConstraint->1]];

Disp[str_] := Module[{},
Print[str, " = ", ToExpression[str]]];


GeodesicEquation[\[CapitalGamma]_, X_, \[Lambda]_]:=Module[{res,v,x, \[CapitalGamma]fix},
x = Table[X[[ia]][\[Lambda]], {ia, Length[X]}];
v = D[x,\[Lambda]];
\[CapitalGamma]fix = \[CapitalGamma]/.Table[X[[ia]]->X[[ia]][\[Lambda]], {ia,Length[X]}];
res = Table[D[x[[ia]],{\[Lambda],2}] + Sum[\[CapitalGamma]fix[[ia,ib,ic]]v[[ib]]v[[ic]], {ib,Length[X]}, {ic,Length[X]}], {ia,Length[X]}]//GSimplify;
res];
InducedMetric[Gdd_, coord_, csign_, X_]:= Module[{nd,\[Gamma]dd},
nd = UnitNormal[Gdd, coord, csign];
\[Gamma]dd =Gdd - Outer[Times, nd,nd];
\[Gamma]dd
];

ExtrinsicCurvature[Gdd_, coord_, csign_,  X_]:=Module[{nd, \[Gamma]dd, delN, delSym, Kdd, \[Gamma]ud},
nd = UnitNormal[Gdd, coord, csign];
\[Gamma]dd =Gdd - Outer[Times, nd,nd];
\[Gamma]ud = IMetric[Gdd].\[Gamma]dd;
delN = CoD[nd, Gdd, X, {-1}];
delN = Simplify[delN, TimeConstraint->.1];
delSym = 1/2(delN + Transpose[delN]);
Kdd = -Transpose[\[Gamma]ud].delSym.\[Gamma]ud;
Kdd];

(* sgn should be +1 for spacelike normal, -1 for timelike normal.  This takes surfaces of constant X[[coord]] *)
UnitNormal[met_, coord_, sgn_] :=Module[{greatertemp, greaterres},
greatertemp =sgn/Sqrt[IMetric[met][[coord,coord]]*sgn];
greaterres = Table[0, {Length[met]}];
greaterres[[coord]] = greatertemp;
GSimplify[greaterres]
];


Unprotect[NonCommutativeMultiply];
NonCommutativeMultiply[tensors__]:=Outer[Times,tensors]//Simplify
Protect[NonCommutativeMultiply];

wedge[forms__]:=Apply[Multinomial,Map[ArrayDepth,List[forms]]]*antisymmetrize[NonCommutativeMultiply[forms]];
exterior[form_, X_]:=(ArrayDepth[form]+1) antisymmetrize[SimpleDeriv[form,X]];
exteriort[form_, X_]:=(ArrayDepth[form]+1) antisymmetrize[SimpleDerivt[form,X]];

(* Assumes first two indices are rectangular matrix *)
MatrixExterior[m_,X_]:= Module[{res, d},
d = Dimensions[m];
res = Table[exterior[m[[ia,ib]],X], {ia, d[[1]]}, {ib, d[[2]]}];
Simplify[res]];
MatrixExteriort[m_,X_]:= Module[{res, d},
d = Dimensions[m];
res = Table[exteriort[m[[ia,ib]],X], {ia, d[[1]]}, {ib, d[[2]]}];
Simplify[res]];

MatrixWedge[w1_, w2_] :=Module[{res,d1, d2},
d1 = Dimensions[w1];
d2 = Dimensions[w2];
res = Table[
Sum[
w1[[ia,ib]]\[Wedge]w2[[ib, ic]],
{ib, d1[[2]]}],
{ia,d1[[1]]}, {ic,d2[[2]]}];
res];

MatrixWedge[a_,b_, c__] := MatrixWedge[a, MatrixWedge[b,c]];

(*Wedge[forms__]:= wedge[forms];
Diamond[forms__]:= MatrixWedge[forms];*)


DiracGamma[rep_, signtt_:-1]:=Module[{res, ia},
res  = Table[{}, {ia,4}];
res[[1]] = DiagonalMatrix[{1,1,-1,-1}];
res[[2]] = ({
 {0, 0, 0, 1},
 {0, 0, 1, 0},
 {0, -1, 0, 0},
 {-1, 0, 0, 0}
});
res[[3]] = ({
 {0, 0, 0, -I},
 {0, 0, I, 0},
 {0, I, 0, 0},
 {-I, 0, 0, 0}
});
res[[4]] = ({
 {0, 0, 1, 0},
 {0, 0, 0, -1},
 {-1, 0, 0, 0},
 {0, 1, 0, 0}
});
If[signtt==-1, res = I res, res = res];
res
]/;rep== "dirac";

DiracGamma[rep_, signtt_:-1]:=Module[{res, ia},
res = DiracGamma["dirac",signtt];
res[[1]] = ({
 {0, 0, 1, 0},
 {0, 0, 0, 1},
 {1, 0, 0, 0},
 {0, 1, 0, 0}
});
If[signtt==-1, res [[1]] = I res[[1]], res[[1]]=res[[1]]];
res
]/;rep== "weyl";
DiracGamma5[rep_, signtt_:-1]:=Module[{gam, res5},
gam = DiracGamma[rep,signtt];
res5 = I *gam[[1]].gam[[2]].gam[[3]].gam[[4]];
res5];
(******  4d SPINORS--this looks general why did I claim 4d only?  **********)
SpinConnection::usage = "SpinConnection[Gdd, \[Eta]bb, Ebu, X]  =  \!\(\*SuperscriptBox[SubscriptBox[\(\[Omega]\), \(\[Mu]\)], \(ab\)]\)";
SpinConnection[Gdd_,\[Eta]bb_,  Ebu_, X_]:=Module[{Guu, \[Eta]aa, Eau, Eua, Ead,DdEad, Ebd, Edb,DdEdb,term1, term2, term3},
Guu = IMetric[Gdd];
\[Eta]aa = IMetric[\[Eta]bb];
Eau = \[Eta]aa.Ebu;
Eua= Transpose[Eau];
Ead = \[Eta]aa.Ebu.Gdd;
DdEad = SimpleDeriv[Ead, X];
Ebd = Ebu.Gdd;
Edb = Transpose[Ebd];
DdEdb = SimpleDeriv[Edb,X];
term1 = 1/2 Table[Sum[ Eua[[iv, ia]](DdEad[[iu,ib,iv]] - DdEad[[iv,ib,iu]]), {iv,Length[X]}], {iu, Length[X]}, {ia, Length[X]}, {ib, Length[X]}]//S;
term2 = Table[term1[[iu, ib,ia]], {iu,Length[X]}, {ia,Length[X]}, {ib,Length[X]}]//S;
term3 = 1/2 Table[Sum[  Eua[[ip,ia]]Eua[[is,ib]](DdEdb[[ip,is,ic]] - DdEdb[[is,ip,ic]])Ead[[ic,iu]], {ip,Length[X]}, {is, Length[X]}, {ic, Length[X]}], {iu, Length[X]}, {ia, Length[X]}, {ib, Length[X]}]//S;
S[term1 - term2 - term3]
];

(*****  END 4d SPINORS   ******)


DAlembertian[f_, met_, X_, sig_:-1]:= 1/Sqrt[sig Det[met]]* Tr[SimpleDeriv[Sqrt[sig Det[met]]IMetric[met].SimpleDeriv[f, X],X]];


dCollect[expr_, f_, proc___]:=Module[{tf,newf},
If[Length[f]==0, tf={f}, tf= f];
newf =Flatten[ Table[{tf[[ia]], tf[[ia]][__], Derivative[__][tf[[ia]]][__]}, {ia,Length[tf]}]];
Collect[expr,newf,proc]];

ChangeDVar[expr_, r1_, r2_, trans12_] := Module[{dr2dr1, rule12, rule21,res, opD},
rule21 = Trans2rules[r2, r1, trans12];
rule12 = Trans2rules[r1, r2, trans12];
dr2dr1 = D[r2/.rule21, r1]/.rule12//S;
opD[f_]:=dr2dr1* D[f,r2];
res = expr/.Derivative[n_][f_][r1]:> Nest[D2, f[r2], n]/.f_[r1]:> f[r2]/.rule12/.D2->opD//S;
res];



Op[\[Zeta]u_, X_]:= Function[\[Zeta]u.SimpleDeriv[#,X]];



TensorMath::invalidIndex = "Invalid index: `1`";
TensorMath::invalidTensor = "Invalid tensor: `1`";
TensorMath::indexMismatch = "Index mismatch error; index `1` has different lengths for different tensors.";
TensorMath::componentOrdering ="Ambiguous or invalid index structure, attempting to add components of type `1` to components of type `2`. Specify ordering using the syntax TensorMath[..., {i1,i2,...}].";
TensorMath::generalError = "Cannot evaluate TensorMath, ended up with the expression: `1`";



(* Get info on upper and lower indices in an expression *)
IndexStructure[T__ComponentTensor] := Module[{alli, upper, lower},
alli = Flatten[Cases[List[T],  ComponentTensor[n_, I_]:>I]];
upper = Cases[alli, _Symbol];
lower = -Cases[alli, -_Symbol];
{upper, lower, alli/.{-s_Symbol->s}}
]

(* ComponentTensor is used internally to add together tensor expressions with labelled indices *)
(*Format[ComponentTensor[components_, I_], StandardForm]:= "<ComponentTensor" <> ToString[I] <> ">";*)
ComponentTensor/: ComponentTensor[C1_, I1_]+ ComponentTensor[C2_, I2_]:=Module[{},
If[!MatchQ[I1, I2],Message[TensorMath::componentOrdering, I1, I2]];
ComponentTensor[ C1+ C2, I1]
];
ComponentTensor/: c_*ComponentTensor[C_, I_]:=ComponentTensor[c*C,I];

Attributes[TensorMath] = {HoldFirst};
TensorMath[expr_, ifree_:Null]:=Module[{res},
res = Hold[expr]/.Times->TTimes;
res = res/.{HoldPattern[Part[a_,b__]]:>ComponentTensor[Hold[a], List[b]/.{TTimes[-1,s_Symbol]->-s}]}//ReleaseHold;
(*Print["res1 = ", res];*)
(* Expand and factor out non-tensors *)
res = res//.t_TTimes:>Distribute[t]//.{
TTimes[a___,b_, c___]:> b*TTimes[a,c]/;FreeQ[b,ComponentTensor],
TTimes[a_]:>a/;FreeQ[a,ComponentTensor],
TTimes[a___,b_*c_,d___]:>c TTimes[a,b,d]/;FreeQ[c,ComponentTensor],
TTimes[]:>1
};
(*Print["res2 = ", res];*)
res =res/.{ComponentTensor[a_,b_]:>TTimes[ComponentTensor[a,b]]}//.{TTimes[a___, TTimes[y___], b___]:>TTimes[a,y,b]}//ReleaseHold;
(*Print["res3 = ", res];*)
res = res/.{TTimes[T__]:>ComponentTensorMultiply[{T}, ifree]};
(*Print["res4 = ", res];*)
If[!MatchQ[Head[res],ComponentTensor], Message[TensorMath::generalError,res];Return[]];

res[[1]]
]

ComponentTensorMultiply[{T__ComponentTensor}, ifreeOrdered_]:= Module[{tensors, istruct, idummy, ifree, i, ilen, badIndices, res, sumIters, tableIters, componentExpression, indices, masks, it, components, matrixIndices, ifreePositioned, tempind},
istruct = IndexStructure[T];
idummy = Intersection[istruct[[1]], istruct[[2]]];
ifree = Select[istruct[[3]], !MemberQ[idummy, #]&]; (*done this way to maintain ordering *)

(* if explicit ordering specified, use it *)
If[MatchQ[Head[ifreeOrdered],List],ifree = Select[ifreeOrdered/.{-s_Symbol->s}, MemberQ[ifree, #]&]]; 
ifreePositioned = If[MemberQ[istruct[[1]], #], #, -#]&/@ifree;

(* Make sure componentTensors are valid *)
badTensors = Cases[{T}, ComponentTensor[t_?(!MatchQ[Head[#],List]&), _]:>t];
If[Length[badTensors]>0, Message[TensorMath::invalidTensor, badTensors[[1]]]; Return[]];
badIndex = Cases[Join[ifree,idummy], Except[_Symbol]];
If[Length[badIndex]>0, Message[TensorMath::invalidIndex, badIndex[[1]]]; Return[]];

(* Dimensions associated to each index *)
ilen = Table[
{i, Union[Dimensions[#[[1]]][[#[[2]]]]&/@
Cases[List[T], ComponentTensor[c_, iall:{ipre___, i|-i, ipost___}]:> {c, Length[List[ipre]]+1}]]},
{i, Join[idummy,ifree]}];
badIndices = Cases[ilen, {i_, {len1_, len2__}}, Infinity, 1];
If[Length[badIndices]>0,Message[TensorMath::indexMismatch, badIndices[[1,1]]]; Return[]];
(* re-format ilen to simplify*)
ilen = ilen/.{{i_, {len_}}:> {i,len}};

(* Wrapping table/sum indexes in tempind prevents accidents if the components depend on a variable of the same name *)
sumIters = Select[ilen, MemberQ[idummy, #[[1]]]&]/.{HoldPattern[{ii_Symbol, ilen_}]:> {tempind[ii], ilen}};
tableIters = Select[ilen, MemberQ[ifree, #[[1]]]&]/.{HoldPattern[{ii_Symbol, ilen_}]:> {tempind[ii], ilen}};

{components, matrixIndices} = Transpose[Cases[List[T], ComponentTensor[c_, i_]:>{c,tempind/@(i/.{-s_Symbol->s})}]];

(* The actual component computation *)
res = Timing[Switch[{Length[idummy]>0,Length[ifree]>0},
{True, True},
Table[Sum[Product[components[[it]][[Sequence@@matrixIndices[[it]]]], {it,Length[components]}], ##]&@@sumIters, ##]&@@tableIters,
{True, False},
Sum[Product[components[[it]][[Sequence@@matrixIndices[[it]]]], {it,Length[components]}], ##]&@@sumIters,
{False, True},
Table[Product[components[[it]][[Sequence@@matrixIndices[[it]]]], {it,Length[components]}], ##]&@@tableIters,
{False,False},
Print["ERROR: No indices in tensormath???"]
]];

(*
Print["idummy = ", idummy];
Print["ifree = ", ifree];
Print["ifreePositioned = ", ifreePositioned];
Print["ilen =", ilen];
Print["sumIters = " , sumIters];
Print["tableIters = " , tableIters];
Print["matrixIndices = ", matrixIndices];
Print["components = ", components];
Print["length = ", Length[components]];
Print[" RES = ", res];
Print["timing = ", res[[1]]];
Print["list ", List[T]];
Print["test = ", Cases[List[T], ComponentTensor[c_, iall:{ipre___, a|-a, ipost___}]\[RuleDelayed] {c, Length[List[ipre]]+1}]];
Print["expr = ", Product[components[[it]][[Sequence@@matrixIndices[[it]]]], {it,Length[components]}]]; 
*)
ComponentTensor[res[[2]], ifreePositioned]
];




EndPackage[]
Print["GREATER has been loaded. Some functions = antisymmetrize, CoD (covariant derivative), Contract, ds2met (covert line element to metric), Epsilon, KillingEquations, SimpleDeriv (partial derivative), LieD (Lie derivative), ChangeCoords, DAlembertian , FieldStrength (F = dA), IndexShift, symmetrize, TJacobian, Trans2rules, SelfContract, CottonTensor, wedge, ExtrinsicCurvature, ShowComponents"];



