(* ::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.                                                         *)
(************************************************************************)



Off[General::"spell",General::"spell1",Remove::"rmnsm",UpSet::"write"];


Unprotect["Global`*"];ClearAll["Global`*"];Remove["Global`*"];Unprotect[$];ClearAll[Wedge,K,$];Protect[$];Unprotect[In,Out];Clear[In,Out];Protect[In,Out];$Line=0;$RecursionLimit=256;$IterationLimit=4096;


Forms[i_]:={};AllScalForms={};matForms[i_]:={};AllMatForms={};AllDifForms={};AllMatrices={};AllSymbols={};FormVars={_Wedge,_d};ZeroVars={_Wedge,_tr};HeadList={};ScalHeadList={};MatFormHeadList={};ZeroHeads={};AllMatHeads={};nodHeads={tr,Bar,Pattern,Blank,Condition,RuleDelayed,SeriesData};$EDCversion=389;


zeroQ[0]=True;zeroQ[x_List]:=And@@(zeroQ/@Union[Flatten[x]]);zeroQ[x_SeriesData]:=If[x[[3]]==={},True,False];zeroQ[x_]:=False;


indepTerms[x_,opt_:{}]:=If[NumericQ[x/._Wedge->1],reWrite[indep$cTerms[Factor[x],opt]],indep$cTerms[x,opt]]


indepCoefs[y_]:=indepTerms[If[FormDegree[y]>0||MatQ[y],Flatten[(#1\[Wedge]$\[Wedge]$&)/@indepTerms[y]//.a_Wedge x_+(z_:Null):>{x,z}]/.Null->Sequence[],y]];


SetAttributes[Bar,{Listable}];Bar[Blank]=Blank;Bar[Bar[x_]]:=x;Bar[Complex[x_,y_]]:=Complex[x,-y];Bar[x_Plus|x_Times|x_Wedge|x_Power|x_Rule|x_Equal]:=Bar/@x;Bar[x_?NumericQ]:=x;Bar[x_RuleDelayed]:=Rule[Bar[x[[1]]],Bar[x[[2]]]]/.Rule->RuleDelayed;Bar[x_SeriesData]:=x/.{First[x]->Bar[First[x]],x[[2]]->Bar[x[[2]]],x[[3]]->Bar/@x[[3]]};Bar[DirectedInfinity[x_]]:=DirectedInfinity[Bar[x]];Bar[HoldPattern[tr[x_]]]:=tr[Bar[x]];Bar[d[x_]]:=d[Bar[x]];Bar[x_Condition]:=x/.x[[1]]->Bar[x[[1]]];Bar[Derivative[x__][y_][z__]]:=If[Union[Bar[{z}]]===Union[{z}],Derivative[Sequence@@Coefficient[{x}.{z},Bar[{z}]]][Bar[y]][z],Derivative[x][Bar[y]][Sequence@@Bar[{z}]]];Bar[h_[y__]/;FreeQ[nodHeads,h]&&FreeQ[h,Derivative]]:=If[Union[Bar[{y}]]===Union[{y}],Bar[h][y],Bar[h][Sequence@@Bar[{y}]]]/;FreeQ[{Integer,Blank,BlankSequence,Pattern,Condition},Head[First[{y}]]];


Qvalue[x_List]:=Length[Select[OrderedQ/@Transpose[{Drop[x,-1],Rest[x]}],#===False&]]


AllCaseTable[x_List]:=With[{rvstr=Reverse[Map[Transpose,x]/.-z_:>z]},If[traceOption===2,Join[Table[{Qvalue[RotateLeft[x,i]],RotateLeft[x,i],i},{i,0,Length[x]-1}],Table[{Qvalue[RotateRight[rvstr,i]],RotateRight[rvstr,i],-i-1},{i,0,Length[x]-1}]],Table[{Qvalue[RotateLeft[x,i]],RotateLeft[x,i],i},{i,0,Length[x]-1}]]]


bestOrd[x_List]:=Block[{tmp$=Sort[AllCaseTable[x]]},tmp$=Select[tmp$,First[#1]===tmp$[[1,1]]&];If[Length[tmp$]>1,First[Sort[Select[tmp$,#1[[2]]===tmp$[[1,2]]&],Abs[Last[#1]]<Abs[Last[#2]]&]],First[tmp$]]];


splitFDprod[n_,FDx_List]:=(Plus@@Take[FDx,n])(Plus@@Drop[FDx,n]);allFDprods[FDx_List]:=Sum[FDx[[i]]*(Plus@@Drop[FDx,i]),{i,Length[FDx]-1}];


traceOption=2;tr[x_]:=FacSimp[Plus@@Transpose[x,{1,1}]]/;MatrixQ[x];tr[x_Plus]:=tr/@x;tr[x_SeriesData]:=(x/.x[[3]]->tr/@x[[3]])+(x/.x[[3]]->{});tr[x__ y_]:=(x tr[y])/;!MatQ[x];tr[d[x_]]:=d[tr[x]];tr[x__\[Wedge]y__]:=(x\[Wedge]tr[Wedge[y]])/;ScalFormQ[x];tr[x_]:=(tr[x]=0)/;AntisymQ[x];tr[Bar[x_]]:=With[{evalHx=tr[x]},Bar[evalHx]/;Head[evalHx]=!=tr];tr[x_Wedge/;Union[Head/@Level[x,1]]==={Bar}]:=With[{evalHx=tr[Bar/@x]},Bar[evalHx]/;Head[evalHx]=!=tr];tr[x_Wedge]:=Block[{inpList=List@@x,resList,FDx},resList=bestOrd[inpList];FDx=FormDegree/@inpList;Which[ Last[resList]>0,(-1)^splitFDprod[ Last[resList],FDx]tr[Wedge@@resList[[2]]], Last[resList]===-1,(-1)^(allFDprods[FDx]+Length[Select[inpList,AntisymQ]])tr[Wedge@@resList[[2]]], Last[resList]<-1,(-1)^(allFDprods[FDx]+splitFDprod[Abs[1+ Last[resList]],FDx]+Length[Select[inpList,AntisymQ]])tr[Wedge@@resList[[2]]]]]/;traceOption>0&&bestOrd[List@@x][[2]]=!=List@@x;tr[x__\[Wedge]y__]:=(tr[x\[Wedge]y]=0)/;(SymQ[Wedge[x]]&&AntisymQ[Wedge[y]]||SymQ[Wedge[y]]&&AntisymQ[Wedge[x]]);tr[x_]:=x matDimension/;!MatQ[x];


toComponents[x_+y_,rul_]:=toComponents[x,rul]+toComponents[y,rul];toComponents[x_.*y_,rul_]:=reWrite[If[MatQ[x*y],x*y,x*y*IdentityMatrix[If[IntegerQ[matDimension],matDimension,Length[Select[(#1[[2]]&)/@rul,MatrixQ][[1]]]]]]/.rul];toComponents[x_\[Wedge]y_,rul_]:=reWrite[If[MatQ[x\[Wedge]y],x\[Wedge]y,x\[Wedge]y*IdentityMatrix[If[IntegerQ[matDimension],matDimension,Length[Select[(#1[[2]]&)/@rul,MatrixQ][[1]]]]]]/.rul];


FormDegree[x_Plus]:=FormDegree[First[x]];FormDegree[x_Times]:=Plus@@FormDegree/@List@@x;FormDegree[x_Wedge]:=Plus@@FormDegree/@List@@x;FormDegree[d[x_]]:=FormDegree[d[x]]=1+FormDegree[x];FormDegree[x_List]:=FormDegree[Last[Union[Flatten[x]]]];FormDegree[tr[x_]]:=FormDegree[x];FormDegree[Bar[x_]]:=FormDegree[Bar[x]]=FormDegree[x];FormDegree[x_SeriesData]:=If[x[[3]]==={},0,FormDegree[x[[3]]]];FormDegree[x_]:=0;


FormCoef[x_List,y_]:=(FormCoef[#,y]&)/@x;FormCoef[x_Plus,y_]:=(FormCoef[#1,y]&)/@x;FormCoef[x_.*y_,y_]:=x;FormCoef[x_,y_Wedge]:=Fold[FormCoef,x,Reverse[List@@y]];FormCoef[x_,y_]:=0/;FreeQ[x/.d[y]->d[$J$],y];FormCoef[x_,y_]:=(((x/.d[y]->d[$J$])-(x/.d[y]->d[$J$]/.y->0))/.y->1/.d[$J$]->d[y])/;EvenQ[FormDegree[y]];FormCoef[x_,y_]:=(((x/.d[y]->d[$J$])-(x/.d[y]->d[$J$]/.y->0))/.d[$J$]->d[y]/.{HoldPattern[Wedge[p___,y,q___]]:>(-1)^Length[Select[{q},OddQ[FormDegree[#]]&]]Wedge[p,q]})/;OddQ[FormDegree[y]];


DeclareMatrixForms[z__]:=Module[{x={{z}},xi,rxi,h,ht,k,min1,tmp=AllMatrices},While[Head[x[[1,1]]]===List,x=First[x]];Unprotect[Transpose];Do[xi=x[[i]];If[Length[xi]<3,Print["Missing symbol for transposed matrix"];Abort[]];If[xi[[3]]===xi[[2]]||xi[[3]]===-xi[[2]],rxi={xi[[2]]},rxi=Rest[xi]];matForms[First[xi]]=Union[matForms[First[xi]],rxi];
AllMatrices=Union[AllMatrices,rxi];h=Head[xi[[2]]];If[First[xi]>0,AllMatForms=Union[AllMatForms,rxi],If[h=!=Symbol,ZeroHeads=Union[{h},ZeroHeads];If[Head[xi[[3]]]=!=Times,ZeroHeads=Union[{Head[xi[[3]]]},ZeroHeads]]]];If[h===Symbol,FormDegree[xi[[2]]]=First[xi];BasicMatQ[xi[[2]]]=True;BasicScalFormQ[xi[[2]]]=False;Transpose[xi[[2]]]=xi[[3]];If[xi[[3]]===-xi[[2]],tr[xi[[2]]]=0,If[xi[[3]]=!=xi[[2]],FormDegree[xi[[3]]]=First[xi];BasicMatQ[xi[[3]]]=True;BasicScalFormQ[xi[[3]]]=False;Transpose[xi[[3]]]=xi[[2]];tr[xi[[3]]]=tr[xi[[2]]]]],FormDegree[_h]=First[xi];BasicMatQ[_h]=True;BasicScalFormQ[_h]=False;ht=Head[xi[[3]]];If[ht===Times,min1=-1;ht=h;tr[h[k__]]=0,min1=1;If[ht=!=h,FormDegree[_ht]=First[xi];BasicMatQ[_ht]=True;BasicScalFormQ[_ht]=False;tr[ht[k__]]=tr[h[k]]]];Transpose[h[k__]]=min1 ht[k];Transpose[ht[k__]]=min1 h[k]],{i,Length[x]}];Protect[Transpose];
AllDifForms=Union[AllScalForms,AllMatForms];AllSymbols=Union[AllDifForms,AllMatrices];MatFormHeadList=Complement[Union[Head/@AllMatForms,MatFormHeadList],{Symbol}];AllMatHeads=Union[ZeroHeads,MatFormHeadList];
HeadList=Union[Head/@AllSymbols];DifFormSymbols=Drop[AllDifForms,-Length[Union[ScalHeadList,MatFormHeadList]]];nodHeads=Union[nodHeads,ZeroHeads];k=Thread[Blank[Union[ScalHeadList,MatFormHeadList]]];FormVars=Flatten[{_Wedge,_d,Union[k,Bar[k]],Union[DifFormSymbols,Bar[DifFormSymbols]],_tr}];ZeroMatSymbols=Drop[matForms[0],-Length[ZeroHeads]];ZeroVars=Flatten[{_Wedge,Union[Thread[Blank[ZeroHeads]],Bar[Thread[Blank[ZeroHeads]]]],Union[ZeroMatSymbols,Bar[ZeroMatSymbols]],_tr}];If[tmp=!={},Print["OK"]]];


DeclareForms[z__]:=(DclFrmNoOK[z];Print["OK"])


DclFrmNoOK[z__]:=Module[{h,x={{z}},xi,rxi,k},While[Head[x[[1,1]]]===List,x=First[x]];Do[xi=x[[i]];rxi=Rest[xi];Forms[First[xi]]=Union[Forms[First[xi]],rxi];
AllScalForms=Union[AllScalForms,rxi];
Do[h=Head[rxi[[j]]];If[h===Symbol,FormDegree[rxi[[j]]]=First[xi];BasicScalFormQ[rxi[[j]]]=True;BasicMatQ[rxi[[j]]]=False;,FormDegree[_h]=First[xi];BasicScalFormQ[_h]=True;BasicMatQ[_h]=False;],{j,Length[rxi]}],{i,Length[x]}];
AllDifForms=Union[AllScalForms,AllMatForms];ScalHeadList=Complement[Union[Head/@AllScalForms,ScalHeadList],{Symbol}];AllSymbols=Union[AllDifForms,AllMatrices];
HeadList=Union[Head/@AllSymbols];DifFormSymbols=Drop[AllDifForms,-Length[Union[ScalHeadList,MatFormHeadList]]];k=Thread[Blank[Union[ScalHeadList,MatFormHeadList]]];FormVars=Flatten[{_Wedge,_d,Union[k,Bar[k]],Union[DifFormSymbols,Bar[DifFormSymbols]],_tr}]];


NoDif[z__]:=(nodHeads=Union[nodHeads,Flatten[{z}]];If[z===interiorProduct||z===int$Prod||z===LieDcartan,Return[],Print["OK"]])


BasicMatQ[Bar[x_]]:=BasicMatQ[x];BasicMatQ[d[x_]]:=BasicMatQ[x];BasicMatQ[x_List]=True;BasicMatQ[x_]:=False;


MatQ[x_Times|x_Wedge|x_Plus]:=Or@@MatQ/@List@@x;MatQ[x_SeriesData]:=Or@@MatQ/@x[[3]];MatQ[x_]:=BasicMatQ[x];


BasicScalFormQ[HoldPattern[tr[x_]]]:=FormDegree[x]>0;BasicScalFormQ[Bar[x_]]:=BasicScalFormQ[x];BasicScalFormQ[d[x_]]:=!BasicMatQ[x];BasicScalFormQ[x_]:=False;


ScalFormQ[x_Times]:=Or@@ScalFormQ/@List@@x;ScalFormQ[x_Wedge|x_Plus]:=And@@ScalFormQ/@List@@x;ScalFormQ[x_]:=BasicScalFormQ[x];


Unprotect[Transpose];Transpose[x_Plus]:=Transpose/@x;Transpose[x_SeriesData]:=x/.x[[3]]->Transpose/@x[[3]];Transpose[x__ y_]:=x Transpose[y]/;!MatQ[x];Transpose[x__\[Wedge]y__]:=(x\[Wedge]Transpose[Wedge[y]])/;ScalFormQ[x];Transpose[x_Wedge]:=(-1)^allFDprods[FormDegree/@(List@@x)]Wedge@@Map[Transpose,Reverse[List@@x]];Transpose[d[x_]]:=d[Transpose[x]];HoldPattern[Transpose[Bar[x_]]]:=Bar[Transpose[x]];Transpose[x_]:=x/;!MatQ[x];Protect[Transpose];


SymQ[Bar[x_]]:=SymQ[x];SymQ[x_]:=Union[Flatten[{reWrite[Transpose[x]-x]}]]==={0};AntisymQ[Bar[x_]]:=AntisymQ[x];AntisymQ[x_]:=Union[Flatten[{reWrite[Transpose[x]+x]}]]==={0};


SetAttributes[d,{Listable}];d[x_Times|x_Wedge]:=d[First[x]]\[Wedge]Rest[x]+(-1)^FormDegree[First[x]]*First[x]\[Wedge]d[Rest[x]];d[x_?NumericQ|x_d]=0;d[matDimension]=0;d[Power[y_,n_]]:=n y^(n-1) d[y]+y^n Log[y]d[n];d[x_Plus]:=d/@x;HoldPattern[d[tr[x_]]]:=With[{evalHx=d[x]},tr[evalHx]/;Head[evalHx]=!=d];HoldPattern[d[Bar[x_]]]:=With[{evalHx=d[x]},Bar[evalHx]/;Head[evalHx]=!=d];d[x_Rule|x_Equal]:=reWrite[d/@x];d[x_SeriesData]:=(x/.x[[3]]->d[x[[3]]])+Wedge[d[First[x]],D[x,First[x]]];d[h_[y__]/;FreeQ[nodHeads,h]]:=Sum[(Derivative[Sequence@@RotateRight[Append[Table[0,{Length[{y}]-1}],1],i]][h][y])d[{y}[[i]]],{i,Length[{y}]}]/;FormDegree[h[y]]===0&&FreeQ[{Integer,Blank,Pattern,Condition},Head[First[{y}]]];


NoDif[interiorProduct];interiorProduct[x_List,0,$z_:Null]=0;interiorProduct[x_List,y_,$z_:Null]:=0/;FormDegree[y]===0;interiorProduct[x_List,y_,$z_:Null]:=reWrite[Which[zeroQ[y/.e[_]->0],If[(y/.d[_]->0)=!=y,er$Wrn[6]];int$Prod[x,y],zeroQ[y/.d[_]->0&&FormDegree[y/.d[k_]:>k]===0],If[$z===Null,er$Wrn[4];Return[],dvarTOe=MapThread[Rule,{d/@$z,e/@Range[Length[$z]]}];eTOdvar=MapThread[Rule,{e/@Range[Length[$z]],d/@$z}];If[(y/.dvarTOe)=!=(y/.dvarTOe/.d[_]->0),er$Wrn[6]];int$Prod[x,y/.dvarTOe]/.eTOdvar],zeroQ[y/.d[_]->0]&&FreeQ[y/.d[k_]:>k,Wedge],Return[],True,tmp=Union[Head/@Select[Union[Flatten[y/.Plus->List/.Times->List/.Wedge->List]],FormDegree[#]===1&]];If[Length[tmp]===1,tmp=tmp[[1]];int$Prod[x,y/.tmp->e]/.e->tmp,If[Length[tmp/.d->Sequence[]]===1,er$Wrn[6];int$Prod[x,y/.tmp->e]/.e->tmp,er$Wrn[5]]]]]/;Head[y]=!=Symbol;


er$Wrn[4]:=stPrint["Must give ordered list of variables (coordinates) as 3rd argument"];


er$Wrn[5]:=stPrint["Error: more than one symbol used for 1-forms: "<>ToString[tmp]];


er$Wrn[6]:=stPrint["To complete the calculation, all exterior derivatives must be expressed in terms of the basis 1-forms!"];


NoDif[int$Prod];int$Prod[x_,0]=0;int$Prod[x_,y_]:=0/;zeroQ[x];int$Prod[x_,y_List]:=int$Prod[x,#]&/@y;int$Prod[x_,y_]:=0/;FormDegree[y]===0;int$Prod[x_,y_Plus]:=int$Prod[x,#]&/@y;int$Prod[x_,u_*y_]:=u*int$Prod[x,y]/;FormDegree[u]===0;int$Prod[x_,e[b_]]:=x[[b]];int$Prod[x_,Wedge[e[k_],p__]]:=x[[k]]*Wedge[p]-Wedge[e[k],int$Prod[x,Wedge[p]]];


NoDif[LieDcartan];LieDcartan[x_,y_,$z_:Null]:=reWrite[d[interiorProduct[x,y,$z]]+interiorProduct[x,d[y],$z]];


newSer$[x_SeriesData,k_]:=SeriesData[First[x],x[[2]],Flatten[Transpose[Prepend[Table[Table[0,{Length[x[[3]]]}],{k-1}],x[[3]]]]],k x[[4]],k x[[5]],k Last[x]]


Wedge[x_]:=x/;Length[{x}]<2&&Head[{x}[[1]]]=!=Pattern


Default[Wedge]:=1;SetAttributes[Wedge,{Flat,OneIdentity}];Wedge[0,y__]=0;Wedge[x__,0]=0;Wedge[x_SeriesData,y_SeriesData]:=Block[{x$,y$,s1$,s2$,re$,x3$,y3$},
  If[Last[x]===Last[y],x$=x;y$=y,x$=newSer$[x,LCM[Last[x],Last[y]]/Last[x]];
    y$=newSer$[y,LCM[Last[x],Last[y]]/Last[y]]];s1$=x$[[-3]]+y$[[-3]];
  s2$=Min[x$[[-2]]+y$[[-3]],x$[[-3]]+y$[[-2]]];
  If[Length[x$[[3]]]<x$[[-2]]-x$[[-3]],
    x3$=Join[x$[[3]],Table[0,{x$[[-2]]-x$[[-3]]-Length[x$[[3]]]}]],
    x3$=x$[[3]]];
  If[Length[y$[[3]]]<y$[[-2]]-y$[[-3]],
    y3$=Join[y$[[3]],Table[0,{y$[[-2]]-y$[[-3]]-Length[y$[[3]]]}]],
    y3$=y$[[3]]];
  Which[s2$===s1$,re$={},s2$-s1$===x$[[-2]]-x$[[-3]],
    re$=Sum[Map[Wedge[x3$[[i]],#]&,Join[Table[0,{i-1}],Drop[y3$,-i+1]]],{i,
          Length[x3$]}],True,
    re$=Sum[Map[Wedge[#,y3$[[i]]]&,Join[Table[0,{i-1}],Drop[x3$,-i+1]]],{i,
          Length[y3$]}]];SeriesData[First[x$],x$[[2]],re$,s1$,s2$,Last[x$]]]/;First[x]===First[y]&&x[[2]]===y[[2]];Wedge[y_,x_SeriesData]:=x/.x[[3]]->Map[Wedge[y,#]&,x[[3]]];Wedge[x_SeriesData,y_]:=x/.x[[3]]->Map[Wedge[#,y]&,x[[3]]];Wedge[x__,y_Plus]:=Plus@@Map[Wedge[x,#]&,List@@y];Wedge[x_Plus,y__]:=Plus@@Map[Wedge[#,y]&,List@@x];Wedge[z__,Times[x_,y_]]:=Times[x,Wedge[z,y]]/;NumericQ[x]||(FormDegree[x]===0&&!MatQ[x]);Wedge[Times[x_,y_],z__]:=Times[x,Wedge[y,z]]/;NumericQ[x]||(FormDegree[x]===0&&!MatQ[x]);Wedge[x_^n_.,y_]:=x^n*y/;FormDegree[x]===0&&!MatQ[x];Wedge[y_,x_^n_.]:=x^n*y/;FormDegree[x]===0&&!MatQ[x];Wedge[x_,y___,x_]:=0/;OddQ[FormDegree[x]]&&ScalFormQ[x];Wedge[x__,y__]:=(-1)^(FormDegree[x]*FormDegree[y])*y\[Wedge]x/;MatQ[x]&&ScalFormQ[y];Wedge[x__]:=Block[{doubL=Transpose[{FormDegree/@{x},{x}}]},Signature[Select[doubL,OddQ[First[#]]&]]Wedge@@Map[Last[#1]&,Sort[doubL]]]/;Union[BasicScalFormQ/@{x}]==={True}&&Map[Last[#1]&,Sort[Transpose[{FormDegree/@{x},{x}}]]]=!={x};Wedge[x_List,y_List]:=If[(FormDegree[x]===0&&!MatQ[Last[Union[Flatten[x]]]])||(FormDegree[y]===0&&!MatQ[Last[Union[Flatten[y]]]]),Dot[x,y],Inner[Wedge,x,y]];Wedge[x_,y_List]:=If[FormDegree[x]===0||(FormDegree[y]===0&&!MatQ[Last[Union[Flatten[y]]]]),x *y,Map[Wedge[x,#]&,y]]/;!MatQ[x];Wedge[x_List,y_]:=If[FormDegree[y]===0||(FormDegree[x]===0&&!MatQ[Last[Union[Flatten[x]]]]),x*y,Map[Wedge[#,y]&,x]]/;!MatQ[y];


simpRules = {x_.*Cos[h_]^2+x_.*Sin[h_]^2:>x};varList={};


SetAttributes[reWrite,{Listable}];reWrite[0]=0;reWrite[x_Equal]:=Equal[reWrite[First[x]-Last[x]],0];reWrite[x_Rule]:=Rule[First[x],reWrite[Last[x]]];reWrite[x_RuleDelayed]:=Rule[First[x],reWrite[x[[2]]]]/.Rule->RuleDelayed;reWrite[x_SeriesData]:=(x/.x[[3]]->reWrite[x[[3]]]);reWrite[x_Condition]:=x/.x[[1]]->reWrite[x[[1]]];reWrite[x_]:=(Collect[x/.simpRules,FormVars,bestFacRul]/.simpRules)/;FormDegree[x]>0;reWrite[x_]:=(Collect[x/.simpRules,ZeroVars,bestFacRul]/.simpRules)/;MatQ[x];reWrite[x_]:=bestFacRul[x/.simpRules]/.simpRules;


FreeALLQ[x_,y_List]:=And@@(FreeQ[x,#]&/@y);


bestFacRul[x_]:=If[varList==={}||FreeALLQ[x,varList],minFacRul[x],varFacRul[x]];varFacRul[x_]:=Collect[Expand[x]/.simpRules,varList,Factor]/.simpRules;minFacRul[x_]:=Factor[Expand[x]/.simpRules]/.simpRules;


Unprotect[SeriesData];If[$VersionNumber<5,Times[x_SeriesData,0]^=0];SeriesData[x1_,x2_,x3_,x4_,x5_,x6_]:= Plus@@If[x2===Infinity,Map[First[#]/x1^((x4+Last[#]-1)/x6)&,Transpose[{x3,Range[Length[x3]]}]]+SeriesData[x1,x2,{},x4,x5,x6],   Map[First[#]*(x1-x2)^((x4+Last[#]-1)/x6)&,Transpose[{x3,Range[Length[x3]]}]]+SeriesData[x1,x2,{},x4,x5,x6]]/;!(FreeQ[x3,x1]||NumericQ[x1]);Protect[SeriesData];


subMatDR[m_,p_,q_]:=Transpose[Drop[Transpose[Drop[m,p]],q]];
subMatTK[m_,p_,q_]:=Transpose[Take[Transpose[Take[m,p]],q]];


FlatOuter[x_?MatrixQ,y_?MatrixQ]:=Partition[Flatten[Transpose[Outer[Times,x,y],{1,3,2,4}]],Dimensions[x][[2]]*Dimensions[y][[2]]];


FlatOuter[x__]:=Fold[FlatOuter,First[{x}],Rest[{x}]]/;Length[{x}]>2;


DirectSum[x_?MatrixQ,y_?MatrixQ]:=Block[{n$1=Dimensions[x],n$2=Dimensions[y],zr$},zr$[__]=0;Join[MapThread[Join,{x,Array[zr$,{n$1[[1]],n$2[[2]]}]}],MapThread[Join,{Array[zr$,{n$2[[1]],n$1[[2]]}],y}]]];


DirectSum[x__]:=Fold[DirectSum,First[{x}],Rest[{x}]]/;Length[{x}]>2;


DeclareMatrixForms[{0,$,$t}];


CHrule=CHrul;CHcoef[x_,0]=1;CHcoef[x_,1]:=tr[x];CHcoef[x_,n_]:=CHcoef[x,n]=reWrite[tr[CHeq[x,n]]/n];


CHeq[x_,n_]:=reWrite[CHcoef[x,n-1]\[Wedge] x+Sum[(-1)^(i-1)CHcoef[x,n-i]\[Wedge]Wedge@@Table[x,{i}],{i,2,n}]];


CHrul[n_,x_:$]:=Module[{cn=CHeq[x,n],cn0,rul,tmp,r1,r2},If[FormDegree[x]>0,Print[x," is not a 0-form matrix"];Return[]]; matDimension=n;rul=First[Solve[cn==cn0,Wedge@@Table[x,{n}]]]/.cn0->tr[cn]/n;rul=reWrite[rul];If[x===$,tmp=MapAt[HoldPattern,rul/.Rule->RuleDelayed,{1,1}];r1=tmp[[1,1]];r2=r1/.$->Pattern[$,Blank[]];tmp={Rule[r2,tmp[[1,2]]/;MatQ[$]&&FormDegree[$]===0&&matDimension===n]}/.Rule->RuleDelayed,tmp=rul]];


detTOtr[n_,A_:$]:=Module[{expr=CHrul[n][[1,2,1]]},(-1)^(1+matDimension) Plus@@Select[List@@expr,MatQ[#1]===False&]/.$->A];


Protect[CHeq,CHrul,CHrule,detTOtr,indepCoefs,FormCoef,subMatTK,subMatDR,FlatOuter,DirectSum];


TrigRules={x_.*Cos[h_]^2+x_.*Sin[h_]^2:>x,x_.*Cot[h_]^2+y_.*Csc[h_]^2:>y/;x+y===0,(Cot[h_]-Csc[h_]) (Cot[h_]+Csc[h_])->-1};HypTrigRules={x_.*Cosh[h_]^2+y_.*Sinh[h_]^2:>x/;x+y===0,x_.*Coth[h_]^2+y_.*Csch[h_]^2:>x/;x+y===0,(Cosh[x_]-Sinh[x_]) (Cosh[x_]+Sinh[x_])->1,(Coth[h_]-Csch[h_]) (Coth[h_]+Csch[h_])->1};imgTrigRules={(Cos[x_]-I Sin[x_]) (Cos[x_]+I Sin[x_])->1,( Cos[x_]+I Sin[x_]) (I Cos[x_]+Sin[x_])->I,(Cos[x_]-I Sin[x_])(-I Cos[x_]+Sin[x_])->-I};RootRules={(-1+x_)(1+x_)/Sqrt[1-x_^2]:>-Sqrt[1-x^2],(-1+x_)(1+x_)/Sqrt[x_^2-1]:>Sqrt[x^2-1],(-x_.+Sqrt[y_])(x_.+Sqrt[y_]):>y-x^2,(x_-Sqrt[y_])(x_+Sqrt[y_]):>-y+x^2};BasicRules=Join[TrigRules,HypTrigRules,RootRules,imgTrigRules];$RGTCversion=389;


SetAttributes[FaSi,{Listable}];FaSi[0]=0;FaSi[x_SeriesData]:=(x/.x[[3]]->FaSi[x[[3]]]);FaSi[x_]:=(#1/.simpRules&)/@If[varList==={}||FreeALLQ[x,varList],Factor[x/.simpRules]/.simpRules,Collect[x/.simpRules,varList,Factor]/.simpRules];


SetAttributes[FacSimp,{Listable}];FacSimp[0]=0;FacSimp[x_Equal]:=Equal[FacSimp[First[x]-x[[2]]],0];FacSimp[x_Rule]:=Rule[First[x],FacSimp[x[[2]]]];FacSimp[x_SeriesData]:=(x/.x[[3]]->FacSimp[x[[3]]]);FacSimp[x_]:=(ReplaceAll[#,simpRules]&/@Collect[x/.simpRules,FormVars,FaSi])/;FormDegree[x]>0;FacSimp[x_]:=(#1/.simpRules&)/@If[varList==={}||FreeALLQ[x,varList],Factor[x/.simpRules]/.simpRules,Collect[x/.simpRules,varList,Factor]/.simpRules];


nonZeroL[x_]:=Select[Union[Flatten[{x}]],!zeroQ[#]&]


nonZeroN[x_]:=Length[Select[Flatten[{x}],!zeroQ[#]&]]


indep$cTerms[x_,opt_:{}]:=redPlus[redTimes[nonZeroL[x],opt]];


redTimes[x_,opt_:{}]:=Block[{tmp$=Select[x,Head[#1]===Times &],rxi},rxi=Complement[x,tmp$]; Union[rxi,Map[numb2[#,opt]&,tmp$]]]


numb1[x_,opt_:{}]:=1/;NumericQ[x];numb1[x_^y_.,opt_:{}]:=1/;MemberQ[opt,x];numb1[x_,opt_:{}]:=x;numb2[x_,opt_:{}]:=Times@@Map[numb1[#,opt]&,Level[x,1]]


redPlus[x_]:=Block[{tmp$=Select[x,Head[#1]===Plus&],cmn$,rst$,i$,negi$},rst$=Complement[x,tmp$];cmn$=Intersection[tmp$,-tmp$];If[cmn$==={},Return[x]];i$=1;While[i$<Length[cmn$],negi$=Select[cmn$,#1===-cmn$[[i$]]&];If[negi$=!={},cmn$=Complement[cmn$,negi$]];i$=i$+1];Union[rst$,Complement[tmp$,-cmn$]]];


RiemSym[h_]:=(h[a_,a_,c_,d_]=0;h[a_,b_,c_,c_]=0;h[a_,b_,c_,d_]:=-h[b,a,c,d]/;a>b;h[a_,b_,c_,d_]:=-h[a,b,d,c]/;c>d;h[a_,b_,c_,d_]:=h[c,d,a,b]/;a>c||(a===c&&b>d));


epsilon[a__]:=Signature[{a}];


zeroTensor[rnk1_]:=Fold[Table,0,List/@Table[Dim,{rnk1}]];


FuncRepRules[z_[y__],g_,ord_:2]:=FuncRepRules[z[y],g,ord]=Block[{tmp$={z[y]->g},indL=indexList[Length[{y}],ord]},Do[tmp$={tmp$,Derivative[Sequence@@indL[[i]]][z][y] ->Expand[FaSi[D[g,Sequence@@Transpose[{{y},indL[[i]]}]]]]},{i,Length[indL]}];Flatten[tmp$]];


indexList[varNo_,ord_]:=Block[{tmp$,indL},tmp$=NestList[RotateRight,Join[{1},Table[0,{varNo-1}]],varNo-1];If[ord>1,Do[indL=tmp$;Do[tmp$=Join[tmp$,Map[Plus[indL[[i]],#]&,indL]],{i,varNo}],{ord-1}]];Union[tmp$]]


M$F[x_]:= Which[ByteCount[x]<40000/Dim,MatrixForm[x],ByteCount[x]<60000,x,True,Map[Short[#,30/Dim]&,x,2]];DclFrmNoOK[1,e[_]];


If[$VersionNumber<6.,st$pr6$[x__]:=StylePrint[x],st$pr6$[x__]:=Print[Style[x]]]


stPrint[x_,y_:Null]:=If[y===Null,st$pr6$[x, FontSize->16,FontSlant->"Italic",FontColor->RGBColor[1, 0, 1]],st$pr6$[x, FontSize->20,FontWeight->"Bold",FontColor->RGBColor[1, 0, 1]]];


timePrint[x_,sT_]:=Print[x<>" computed in ",TimeUsed[]-sT," sec"];


testSER[x_,y_String,z_:""]:=Module[{tmp,k},If[FreeQ[x,SeriesData]||(Dim===2&&TensRnk[x]<4),y<>z,tmp=Sort[Union[Select[Flatten[Which[TensRnk[x]===4,Raise[x,1,2],TensRnk[x]>1,Raise[x,1],True,x]],Head[#]===SeriesData&]],#1[[-2]]/#1[[-1]]<#2[[-2]]/#2[[-1]]&];k=" to order "<>ToString[(tmp[[1,-2]]-1)/tmp[[1,-1]]];y<>If[Length[tmp]>1,k=k<>" or higher",k]]];


Gn$Cf0[x_List,y__]:=Gn$Cf0[#,y]&/@x;Gn$Cf0[x_SeriesData,y_List]:=Gn$Cf0[x,#]&/@y;Gn$Cf0[x_SeriesData,y__]:=x/.(x[[3]]->Map[Coefficient[#,y]&,x[[3]]]);Gn$Cf0[x_,y__]:=Coefficient[x,y];GenCoef[x_,y__]:=FacSimp[Gn$Cf0[x,y]];


fct$rls={(-1+x_Symbol)^b_.*(1+x_Symbol)^b_.:>(-1+x^2)^b,(1-x_Symbol)^b_.*(1+x_Symbol)^b_.:>(1-x^2)^b,(x_Symbol-y_Symbol)^b_.*(x_Symbol+y_Symbol)^b_.:>(x^2-y^2)^b};


er$Wrn[1]:=stPrint["coframe does not span space!"];


RGtensors[gIN_,xIN_,opt___]:=Module[{frameOpt=0,eIN,Ropt=True,Wopt=True,Eopt=True,idMat,eFrame,dxRul,Bmat,Amat,NP$=False,de,k$,StrCon,\[Gamma]Udd,\[Gamma]ddd,g1ddd,g2ddd,Gamddd,rmn,tmp,RUd,Rg,gg,sT,Cflat,EinSp},Clear[GUdd,covD,covDiv,LieD,\[Omega]Ud,RUddd,Rdddd,Rdd,EUd,R,Wdddd,detg,p$D,rtAbsdetg,eta,HStar,Lower,Raise,\[Tau],\[Kappa],\[Rho],\[Sigma],\[Gamma],\[Epsilon],\[Alpha],\[Beta],\[Nu],$\[Pi],\[Lambda],\[Mu],\[CapitalLambda],\[CapitalPsi],\[CapitalPhi],\[CapitalDelta],\[GothicCapitalD],\[Delta],\[Delta]B];DEF$covD;DEF$covDiv;DEF$LieD;DEF$LWR;DEF$RSE;
		
		Which[Head[xIN]=!=List,tmp="Coordinates ",Head[gIN]=!=List,tmp="Metric "];If[Head[tmp]===String,stPrint[tmp<>"not defined"];Return[]];Dim=Length[xIN];If[Dim=!=Length[gIN],stPrint["Metric and coordinate dimensions don't match"];Print["Metric ",Dimensions[gIN],"   ",xIN];Return[]];If[Transpose[gIN]=!=gIN,Print[Map[Short[#,10/Dim]&,gIN,{2}]," is not a symmetric matrix"];Return[]];eFrame=e/@Range[Dim];coordList=xIN;(* frameOpt=0 => coordFrame, frameOpt=1 => explicitFrame, frameOpt=2 => implicitFrame *)
	
			If[Length[{opt}]>0,tmp=Select[Flatten[{opt}],!(FreeQ[#,_e]&&FreeALLQ[#,d[coordList]])&];If[zeroQ[tmp/._e->0]&&Length[tmp]>1,stPrint["Evaluate Clear$dx to clear the coordinate differentials!"];Return[]];
If[zeroQ[tmp/._e->0]&&Length[tmp]===1,If[Head[dxRuleList]===List&&Head[deList]===List,frameOpt=2;eIN=eFrame;If[Length[deList]=!=Dim,stPrint["Length[deList] is not equal to "<>ToString[Dim]<>"!"];Return[]],stPrint["dxRuleList and/or deList have not been defined!"];Return[]]];
		Which[Length[{opt}]===2&&FormDegree[{opt}[[1]]]===1,If[frameOpt=!=2,eIN={opt}[[1]];frameOpt=1;If[Length[eIN]=!=Dim,er$Wrn[1];Return[]]];If[{opt}[[2,1]]===0,Ropt=False];If[{opt}[[2,2]]===0,Wopt=False];If[Length[{opt}[[2]]]===3&&{opt}[[2,3]]===0,Eopt=False],Length[{opt}]===1&&IntegerQ[opt[[1]]],
If[opt[[1]]===0,Ropt=False];If[opt[[2]]===0,Wopt=False];If[Length[{opt}[[1]]]===3&&{opt}[[1,3]]===0,Eopt=False],Length[{opt}]===1&&FormDegree[{opt}[[1]]]===1,If[frameOpt=!=2,eIN=opt;frameOpt=1;If[Length[eIN]=!=Dim,er$Wrn[1];Return[]]],True,stPrint["Wrong optional arguments!"];Return[]]];If[frameOpt===0,eIN=d/@xIN];If[gIN==={{0,1,0,0},{1,0,0,0},{0,0,0,-1},{0,0,-1,0}},If[FreeALLQ[Union[Flatten[{eIN,xIN,{deList,dxRuleList}}]],{\[Tau],\[Kappa],\[Rho],\[Sigma],\[Gamma],\[Epsilon],\[Alpha],\[Beta],\[Nu],$\[Pi],\[Lambda],\[Mu],\[CapitalLambda],\[CapitalPsi],\[CapitalPhi],\[CapitalDelta],\[GothicCapitalD],\[Delta],\[Delta]B}],NP$=True,stPrint["NP symbols used in coordinates or frame: they will not be defined!"]]];
				
				If[!FreeQ[{gIN,eIN,{deList}},SeriesData],tmp=Select[Flatten[{gIN,eIN,{deList}}],Head[#]===SeriesData&][[1,1]];If[MemberQ[xIN,tmp],de[1]=ToExpression["d"<>ToString[tmp]];de[2]=ToExpression["\[ScriptD]"<>ToString[tmp]];DclFrmNoOK[1,\[ScriptD][_],de[1],de[2]];d[tmp]=de[1];\[ScriptD][tmp]=de[2];d[de[1]]=0;If[frameOpt>0&&!FreeQ[{eIN,{deList,dxRuleList}},Derivative[1][d][0]],stPrint["d["<>ToString[tmp]<>"] not set equal to "<>ToString[de[1]]<>".  Evaluate RGtensors[...] again!"];Return[]],d[tmp]=0]];If[frameOpt===1,DclFrmNoOK[1,\[ScriptD][_]];eTO$dx=MapThread[Rule,{eFrame,FacSimp[eIN/.MapThread[Rule,{d/@xIN,\[ScriptD]/@xIN}]]}]];
				
				If[frameOpt>0&&(Union[FormDegree/@eIN]=!={1}||!zeroQ[eIN/.{d[k_]->0,e[_]->0,y_|Bar[y_]:>0/;MemberQ[DifFormSymbols,y]}]),stPrint["coframe has incorrect form"];Return[]];If[Head[varList]=!=List,varList={}];If[Head[simpRules]===Rule||Head[simpRules]===RuleDelayed,simpRules={simpRules}];If[Head[simpRules]=!=List,simpRules={}];idMat=IdentityMatrix[Dim];
				
				gdd=gIN;Print["gdd = ",M$F[gdd]];tmp=eIN.gIN.eIN;k$=If[frameOpt===2,eFrame,d/@xIN];k$=Union[Flatten[Outer[Times,k$,k$]]];LineElement=If[Head[tmp]===SeriesData,tmp/.tmp[[3]]->Collect[tmp[[3]],k$,FaSi],Collect[tmp,k$,FaSi]]//.fct$rls;Print["LineElement = ",Short[LineElement,8Dim]];sT=TimeUsed[];gUU=FixedPoint[FaSi,Inverse[gIN]]//.fct$rls;detg=FixedPoint[FaSi,Det[gIN]]//.fct$rls;Print["gUU = ",M$F[gUU]];timePrint["gUU",sT];
				
				Which[frameOpt===0,Amat=idMat;Bmat=idMat;dxRul={},frameOpt===1&&zeroQ[eIN/.{d[k_]->0,y_|Bar[y_]:>0/;MemberQ[DifFormSymbols,y]}],Bmat=Array[Gn$Cf0[eIN[[#1]],d[xIN[[#2]]]]&,{Dim,Dim}];If[!zeroQ[FaSi[Det[Bmat]]],sT=TimeUsed[];If[ByteCount[eIN]>1500*Dim,stPrint["Computing d[coframe]..."]];Amat=FaSi[Inverse[Bmat]];dxRul=MapThread[Rule,{d[xIN],Amat.eFrame}];tmp=indep$cTerms[Flatten[Outer[Wedge,d[xIN],d[xIN]]]];tmp=MapThread[Rule,{tmp,FacSimp[tmp/.dxRul]}];de=FacSimp[d[eIN]];de=FacSimp[de/.tmp];timePrint["d[coframe]",sT];tmp=d[xIN]/.dxRul;MapThread[Set,{d[xIN],tmp}],er$Wrn[1];Return[]],frameOpt===2,If[!FreeQ[Join[dxRuleList,deList],e[0]],stPrint["The notation e[0] cannot be used; renumber the frame components as e[1]...e[Dim]"];Abort[]];de=deList;dxRul=dxRuleList;tmp=d[xIN]/.dxRul;MapThread[Set,{d[xIN],tmp}];Amat=Array[Gn$Cf0[tmp[[#1]],eIN[[#2]]]&,{Dim,Dim}]];
		
		If[frameOpt===2&&(Length[Union[FreeQ[#,SeriesData]&/@{deList,dxRuleList}]]>1||!FreeQ[gIN,SeriesData]&&(FreeQ[deList,SeriesData]||FreeQ[dxRuleList,SeriesData])),stPrint["RGtensors will run faster if both deList and dxRuleList are expanded in series"]];
		
		rtAbsdetg=PowerExpand[Factor[detg^2]^(1/4)];SetAttributes[HStar,{Listable}];HStar[0]=0;If[NP$,rtAbsdetg=I];HStar[x_SeriesData]:=x/.x[[3]]->HStar/@x[[3]];HStar[x_Plus]:=FacSimp[HStar/@x];HStar[x_*y_]:=x*HStar[y]/;FormDegree[x]===0;HStar[x_]:=x*rtAbsdetg*Wedge@@eFrame/;FormDegree[x]===0;HStar[x_Wedge]:=rtAbsdetg/detg/;Length[x]===Dim;HStar[x_Wedge|x_e]:=HStar[x]=rtAbsdetg*FacSimp[inn[x,Wedge@@eFrame]];eta[]:=eta[]=rtAbsdetg*Array[Signature[{##}]&,Table[Dim,{Dim}]];eta[a__]:=eta[a]=Raise[eta[],a];sT=TimeUsed[];
		
		If[frameOpt>0,d[e[a_]]:=de[[a]];StrCon[c_,a_,a_]=0;StrCon[c_,a_,b_]:=-StrCon[c,b,a]/;a>b;StrCon[c_,a_,b_]:=StrCon[c,a,b]=-Gn$Cf0[de[[c]],e[a]\[Wedge]e[b]];\[Gamma]Udd=Array[StrCon,{Dim,Dim,Dim}];\[Gamma]ddd=gdd.\[Gamma]Udd;If[NP$,d$f$[0]=de[[3]]\[Wedge]e[1]-e[3]\[Wedge]de[[1]];d$f$[2]=de[[2]]\[Wedge]e[4]-e[2]\[Wedge]de[[4]];d$f$[1]=FacSimp[(de[[3]]\[Wedge]e[4]-e[3]\[Wedge]de[[4]]+de[[2]]\[Wedge]e[1]-e[2]\[Wedge]de[[1]])/2];],\[Gamma]Udd=zeroTensor[3];\[Gamma]ddd=\[Gamma]Udd];
	
		(** Tests **)If[frameOpt>0,tmp=(d[Union[Flatten[gdd]]]/.dxRul)/.e[_]->0;If[!zeroQ[tmp],stPrint["Warning: d[] on some symbols not defined!"];Print[indep$cTerms[FacSimp[tmp]]]];If[frameOpt===2,If[(ByteCount[deList]+ByteCount[dxRuleList])/Dim<3000,tmp=FacSimp[FacSimp[d[deList]]];If[zeroQ[tmp],Print["deList closure verified"],stPrint["Warning: deList not closed! d[deList]="];Print[tmp];If[!zeroQ[tmp/.d[_] -> 0],Return[]]],stPrint["Warning: deList closure not tested"]]]];
		
		p$D[x_,y_]:=p$D[x,y]=If[Amat===idMat,D[x,coordList[[y]]],Sum[Amat[[k$,y]]*D[x,coordList[[k$]]],{k$,Dim}]];
		
sT=TimeUsed[];g1ddd[a_,b_,c_]:=g1ddd[b,a,c]/;a>b;g1ddd[a_,b_,c_]:=g1ddd[a,b,c]=p$D[gdd[[a,b]],c];
g2ddd[c_,a_,b_]:=g2ddd[c,b,a]/;a>b;g2ddd[c_,a_,b_]:=g2ddd[c,a,b]=FaSi[(g1ddd[a,c,b]+g1ddd[b,c,a]-g1ddd[a,b,c])/2];If[frameOpt===0,Gamddd=Array[g2ddd,{Dim,Dim,Dim}];
Clear[tmp];k$=gUU.Gamddd;tmp[c_,a_,b_]:=tmp[c,b,a]/;a>b;tmp[c_,a_,b_]:=tmp[c,a,b]=FaSi[k$[[c,a,b]]];GUdd=Array[tmp,{Dim,Dim,Dim}],
Gamddd=FaSi[Array[g2ddd,{Dim,Dim,Dim}]+(\[Gamma]ddd+Transpose[\[Gamma]ddd]-Transpose[\[Gamma]ddd,{3,2,1}])/2];GUdd=If[Union[NumericQ/@Flatten[gUU]]==={True}&&nonZeroN[gUU]===Dim,gUU.Gamddd,FaSi[gUU.Gamddd]]]; timePrint["Gamma",sT];If[frameOpt>0,\[Omega]Ud=Array[Sum[GUdd[[#1,k$,#2]]*e[k$],{k$,Dim}]&,{Dim,Dim}]];If[NP$,NPdefs1];
		
	sT=TimeUsed[];RiemSym[rmn];rmn[a_,b_,c_,d$_]:=rmn[a,b,c,d$]=FaSi[p$D[Gamddd[[a,d$,b]],c]-p$D[Gamddd[[a,c,b]],d$]+Sum[GUdd[[k$,c,b]]*Gamddd[[k$,d$,a]]-GUdd[[k$,d$,b]]*Gamddd[[k$,c,a]]-Gamddd[[a,k$,b]]*\[Gamma]Udd[[k$,c,d$]],{k$,Dim}]];Rdddd=Array[rmn,{Dim,Dim,Dim,Dim}];timePrint["Riemann(dddd)",sT];sT=TimeUsed[];If[zeroQ[Rdddd],stPrint[testSER[Rdddd,"Flat Space","!"],1];Return[]];Rg=gUU . Rdddd;If[Ropt,Clear[rmn];rmn[a_,b_,c_,c_]=0;rmn[a_,b_,c_,d$_]:=-rmn[a,b,d$,c]/;c>d$;rmn[a_,b_,c_,d$_]:=rmn[a,b,c,d$]=FaSi[Rg[[a,b,c,d$]]];RUddd=Array[rmn,{Dim,Dim,Dim,Dim}];timePrint["Riemann(Uddd)",sT];  Rg=RUddd,stPrint["RUddd not computed"]];
				
		sT=TimeUsed[];Rg=Plus@@Transpose[Rg,{1,2,1,3}];Clear[tmp];tmp[a_,b_]:=tmp[b,a]/;a>b;tmp[a_,b_]:=tmp[a,b]=FaSi[Rg[[a,b]]];Rdd=Array[tmp,{Dim,Dim}];RUd=gUU . Rdd;R=FaSi[Plus@@Transpose[RUd,{1,1}]];timePrint["Ricci",sT];sT=TimeUsed[];
				
		If[Wopt||Dim<4,If[Dim>3,If[zeroQ[Rdd],Wdddd=Rdddd,gg=Transpose[Outer[Times,gdd,gdd],{1,3,2,4}];gg=gg-Transpose[gg];Rg=Transpose[Outer[Times,Rdd,gdd],{1,3,2,4}];Rg=Rg-Transpose[Rg];Rg=Rg-Transpose[Rg,{1,2,4,3}];Clear[rmn];RiemSym[rmn];rmn[a_,b_,c_,d$_]:=rmn[a,b,c,d$]=FaSi[Rdddd[[a,b,c,d$]]-Rg[[a,b,c,d$]]/(Dim-2)+R*gg[[a,b,c,d$]]/((Dim-1)*(Dim-2))];Wdddd=Array[rmn,{Dim,Dim,Dim,Dim}]],Wdddd=zeroTensor[4]];timePrint["Weyl",sT];If[zeroQ[Wdddd],Cflat=True;If[Dim=!=3,stPrint[testSER[Wdddd,"Conformally Flat"],1],tmp=CF3d[Wopt];If[zeroQ[tmp],stPrint[testSER[tmp,"Conformally Flat"],1]]]],stPrint["Weyl tensor not computed"]];If[NP$,NPdefs2];
			
					sT=TimeUsed[];If[Eopt,EUd=FaSi[RUd-1/2*R*idMat],stPrint["Einstein tensor not computed"]];If[zeroQ[Rdd],stPrint[testSER[Rdd,"Ricci Flat"],1],If[Eopt,timePrint["Einstein",sT]];If[Eopt&&zeroQ[EUd+(1/2-1/Dim)*R*idMat],stPrint[testSER[Rdd-1/Dim*R*gdd,"Einstein Space"],1];EinSp=True]];If[Dim>2&&Cflat&&EinSp||Dim===2&&zeroQ[FaSi[covD[R]]],stPrint[testSER[Rdd-1/Dim*R*gdd,"Space of Constant Curvature","!"],1]];"All tasks completed in "]//Timing//Apply[#2<>ToString[#1]<>If[$VersionNumber>5.5," seconds",""] &, 
  If[#[[2]]===Return[]||#[[2]]===Null,#/.#[[2]] -> 
      "Aborted after ",#]/. Second ->"seconds"]&//Print


CF3d[x_]:=If[x,stPrint["Testing for 3-dim conformal flatness..."];Block[{tmp$=covD[Rdd]-Outer[Times,gdd,covD[R]/2/(Dim-1)]},FaSi[tmp$-Transpose[tmp$,{1,3,2}]]],stPrint["3-dim conformal flatness not tested"];Return[1]];


NPdefs1:=($\[Pi]=GenCoef[d$f$[1],e[1]\[Wedge]e[2]\[Wedge]e[3]];\[Mu]=GenCoef[d$f$[1],e[1]\[Wedge]e[3]\[Wedge]e[4]];\[Tau]=-GenCoef[d$f$[1],e[1]\[Wedge]e[2]\[Wedge]e[4]];\[Rho]=-GenCoef[d$f$[1],e[2]\[Wedge]e[3]\[Wedge]e[4]];tmp$=FacSimp[d$f$[0]+\[Tau] e[1]\[Wedge]e[3]\[Wedge]e[4]+\[Rho] e[1]\[Wedge]e[2]\[Wedge]e[3]];\[Epsilon]=GenCoef[tmp$,e[1]\[Wedge]e[2]\[Wedge]e[3]]/2;\[Beta]=GenCoef[tmp$,e[1]\[Wedge]e[3]\[Wedge]e[4]]/2;\[Sigma]=-GenCoef[tmp$,e[1]\[Wedge]e[2]\[Wedge]e[4]];\[Kappa]=-GenCoef[tmp$,e[2]\[Wedge]e[3]\[Wedge]e[4]];tmp$=FacSimp[d$f$[2]-$\[Pi] e[2]\[Wedge]e[3]\[Wedge]e[4]-\[Mu] e[1]\[Wedge]e[2]\[Wedge]e[4]];\[Gamma]=-GenCoef[tmp$,e[1]\[Wedge]e[2]\[Wedge]e[4]]/2;\[Alpha]=-GenCoef[tmp$,e[2]\[Wedge]e[3]\[Wedge]e[4]]/2;\[Lambda]=GenCoef[tmp$,e[1]\[Wedge]e[2]\[Wedge]e[3]];\[Nu]=GenCoef[tmp$,e[1]\[Wedge]e[3]\[Wedge]e[4]];Print["Spin coefficients \[Tau], \[Kappa], \[Rho], \[Sigma], \[Gamma], \[Epsilon], \[Alpha], \[Beta], \[Nu], $\[Pi], \[Lambda], \[Mu] defined"];Bar[e[1]]=e[1];Bar[e[2]]=e[2];Bar[e[4]]=e[3];Bar[e[3]]=e[4];Bar[x_]:=x;\[CapitalDelta][x_]:=FacSimp[p$D[x,1]];\[GothicCapitalD][x_]:=FacSimp[p$D[x,2]];\[Delta][x_]:=FacSimp[-p$D[x,4]];\[Delta]B[x_]:=FacSimp[-p$D[x,3]];Print["Differential Operators \[GothicCapitalD][x_], \[CapitalDelta][x_], \[Delta][x_], \[Delta]B[x_] defined"];)


NPdefs2:=(\[CapitalPhi][0, 0]=Rdd[[2,2]]/2;\[CapitalPhi][0, 1]=-Rdd[[2,4]]/2;\[CapitalPhi][1, 0]=-Rdd[[2,3]]/2;\[CapitalPhi][2, 2]=Rdd[[1,1]]/2; \[CapitalPhi][2, 1]=-Rdd[[3,1]]/2;\[CapitalPhi][1, 2]=-Rdd[[4,1]]/2;\[CapitalPhi][0, 2]=Rdd[[4,4]]/2;\[CapitalPhi][2, 0]=Rdd[[3,3]]/2; \[CapitalPhi][1, 1]=FaSi[(Rdd[[1,2]]+Rdd[[3,4]])/4];\[CapitalLambda]=-R/24;tmp$="NP variables \[CapitalPhi][i,j], \[CapitalLambda] ";If[Head[Wdddd]=!=Symbol,\[CapitalPsi][0]=-Wdddd[[2,4,2,4]];\[CapitalPsi][1]=Wdddd[[2,4,2,1]];\[CapitalPsi][2]=-Wdddd[[2,4,3,1]]; \[CapitalPsi][3]=Wdddd[[2,1,3,1]];\[CapitalPsi][4]=-Wdddd[[3,1,3,1]];tmp$=tmp$<>"and \[CapitalPsi][i] "]; Print[tmp$<>"defined"];);


inn[y_Wedge,x_]:=inn[Rest[y],inn[First[y],x]];inn[e[a_],x_+y_]:=inn[e[a],x]+inn[e[a],y];inn[e[a_],x_*y_]:=x*inn[e[a],y]/;FormDegree[x]===0;inn[e[a_],x_]:=0/;FormDegree[x]===0;inn[e[a_],x_Wedge]:=inn[e[a],First[x]]*Rest[x]+(-1)^FormDegree[First[x]]*First[x]\[Wedge]inn[e[a],Rest[x]];inn[e[a_],e[b_]]:=gUU[[a,b]];


TensRnk[x_]:=If[Head[x]===List,Length[Dimensions[x]],0]


MoveIndex[x_,i_,g_]:=Block[{rnk1=TensRnk[x],indL},Which[i===1,g.x,i===rnk1,x.g,True,indL=Range[rnk1]/.{i->1,1->i};Transpose[g.Transpose[x,indL],indL]]];


DEF$LWR:=(Lower[x_,i__]:=Lower[x,i]=FacSimp[FoldList[MoveIndex[#1,#2,gdd]&,x,Flatten[{i}]][[-1]]];)


DEF$RSE:=(Raise[x_,i__]:=Raise[x,i]=FacSimp[FoldList[MoveIndex[#1,#2,gUU]&,x,Flatten[{i}]][[-1]]];)


DEF$LWR;DEF$RSE;


multiDot[x_,y_,indPr__]:=Module[{rnk1=TensRnk[x],rnk2=TensRnk[y],indL=Sort/@Transpose[{indPr}],ordPr=Sort[{indPr}]},If[!And@@MapThread[SameQ[Dimensions[x][[#1]],Dimensions[y][[#2]]]&,Transpose[ordPr]],er$Wrn[3];Return[]];If[Union[First[indL]]=!=First[indL]||indL[[1,-1]]>rnk1||Union[indL[[2]]]=!=indL[[2]]||Last[indL[[2]]]>rnk2,er$Wrn[2];Return[]];Which[Length[ordPr]===1,FacSimp[dotOne[x,y,indPr]],rnk1===rnk2&&rnk1===Length[ordPr],FacSimp[FacSimp[dotAll[x,y,ordPr]]],True,OKindCon[dotOne[x,y,First[ordPr]],newPrList[ReplaceAll[#,{a_,b_}:>{a,b+rnk1}]&/@ordPr]]]];


dotOne[x_,y_,{i_,j_}]:=Block[{rnk1=TensRnk[x],rnk2=TensRnk[y],tmp$},If[(FormDegree[x]>0||MatQ[Last[Union[Flatten[x]]]])&&(FormDegree[y]>0||MatQ[Last[Union[Flatten[y]]]]),tmp$=Wedge,tmp$=Dot];tmp$[Transpose[x,Join[Range[i-1],RotateRight[Range[i,rnk1]]]],Transpose[y,Join[RotateLeft[Range[j]],Range[j+1,rnk2]]]]];


dotAll[x_,y_,indPr_]:=Module[{s$,allInd,sumInd,tmp$},allInd=s$/@Range[Length[indPr]];allInd={allInd,allInd/.Map[Rule[#[[1]],#[[2]]]&,RotateLeft/@indPr]};sumInd=Transpose[{First[allInd],Dimensions[x]}];If[(FormDegree[x]>0||MatQ[Last[Union[Flatten[x]]]])&&(FormDegree[y]>0||MatQ[Last[Union[Flatten[y]]]]),tmp$=Wedge,tmp$=Times];Sum@@{Unevaluated[tmp$[x[[Sequence@@First[allInd]]],y[[Sequence@@allInd[[2]]]]]],Sequence@@sumInd}];


Contract[x_,indPr__]:=Module[{rnk1=TensRnk[x],indL=Flatten[{indPr}]},If[!And@@Map[SameQ[Dimensions[x][[#[[1]]]],Dimensions[x][[#[[2]]]]]&,{indPr}],er$Wrn[3];Return[]];If[Length[Union[indL]]<Length[indL]||Length[indL]>rnk1||Sort[indL][[-1]]>rnk1,er$Wrn[2];Return[],OKindCon[x,Sort[Sort/@{indPr}]]]];


newPrList[x_]:=Rest[x]-1/.k_:>k-1/;k>x[[1,2]]-1;


OKindCon[x_,indPr_]:=FacSimp[ConOne[x,First[indPr]]]/;Length[indPr]===1;OKindCon[x_,indPr_]:=OKindCon[ConOne[x,First[indPr]],newPrList[indPr]];


ConOne[x_,{i_,j_}]:=Block[{rnk1=TensRnk[x],tr1List,tr2List},tr1List=(Range[rnk1]/.{j->i})/.k_:>k-1/;k>j;tr2List=Join[Range[i-1],RotateRight[Range[i,rnk1-1]]];Apply[Plus,Transpose[Transpose[x,tr1List],tr2List],{rnk1-2}]];


er$Wrn[2]:=stPrint["Wrong tensor rank or number(s) of indices!"];er$Wrn[3]:=stPrint["Summed indices have unequal dimensions!"];


DEF$covD:=(covD[x_List,y_List]:=covD[x,y]=c$D[x,y];covD[x_,y_List]:=If[Length[y]>0,er$Wrn[2],covD[x]];covD[x_]:=covD[x]=c$D[x];)


c$D[x_,UpList_:{}]:=Module[{k=TensRnk[x]},If[Length[UpList]>0&&(Length[Flatten[UpList]]>k||Sort[Flatten[UpList]][[-1]]>k),er$Wrn[2];Return[]];If[Head[x]===List,Array[cDcomp[x,k,UpList],Table[Dim,{k+1}]],FaSi[Array[p$D[x,#]&,{Dim}]]]];


cDcomp[x_,rnk1_,UpList_:{}][y__]:=Block[{i$,j$,y$=Drop[{y},-1],z$=Last[{y}]},FaSi[p$D[x[[Sequence@@y$]],z$]+Sum[If[MemberQ[UpList,i$],GUdd[[y$[[i$]],z$,j$]]x[[Sequence@@ReplacePart[y$,j$,i$]]],-GUdd[[j$,z$,y$[[i$]]]]x[[Sequence@@ReplacePart[y$,j$,i$]]]],{j$,Dim},{i$,rnk1}]]];


DEF$covDiv:=(covDiv[x_,UpList_]:=covDiv[x,UpList]=Module[{k=TensRnk[x],divInd,j$,z$,FlUpL=Flatten[UpList]},If[Length[FlUpL]>k||Sort[FlUpL][[-1]]>k||Head[x]=!=List,er$Wrn[2];Return[]]; If[Length[UpList]>1,divInd=Flatten[Select[UpList,Head[#]===List&]][[1]],divInd=FlUpL[[1]]];
If[k===1,FaSi[Sum[p$D[x[[z$]],z$],{z$,Dim}]+Sum[GUdd[[z$,z$,j$]]x[[j$]],{j$,Dim},{z$,Dim}]],
Array[cDIVcomp[x,k,divInd,FlUpL],Table[Dim,{k-1}]]]];)


cDIVcomp[x_,rnk1_,divInd_,UpList_:{}][y__]:=Block[{z$,y$,i$,j$},y$=Insert[{y},z$,divInd];FaSi[Sum[p$D[x[[Sequence@@y$]],z$],{z$,Dim}]+Sum[If[MemberQ[UpList,i$],GUdd[[y$[[i$]],z$,j$]]x[[Sequence@@ReplacePart[y$,j$,i$]]],-GUdd[[j$,z$,y$[[i$]]]]x[[Sequence@@ReplacePart[y$,j$,i$]]]],{j$,Dim},{i$,rnk1},{z$,Dim}]]];


DEF$LieD:=(LieD[x_List,y_List,z_List]:=LieD[x,y,z]=LieDerv[x,y,z];LieD[x_List,y_]:=LieD[x,y]=LieDerv[x,y];)


LieDerv[h_List,x_,UpList_:{}]:=Module[{k=TensRnk[x],hjdU,hUjd,i$,ASgam=GUdd-Transpose[GUdd,{1,3,2}]},
If[Length[Dimensions[h]]=!=1,stPrint["First argument not a vector!"];Return[]];If[Length[UpList]>0&&(Length[Flatten[UpList]]>k||Sort[Flatten[UpList]][[-1]]>k),er$Wrn[2];Return[]];hjdU=Array[p$D[h,#1]&,{Dim}];hUjd=Transpose[hjdU]+ASgam.h;hjdU=Transpose[hUjd];FaSi[If[FormDegree[x]>0,LieDcartan[h,x,coordList],h.Array[p$D[x,#1]&,{Dim}]]+If[Head[x]===List,Sum[Which[i$===1,If[MemberQ[UpList,i$],-hUjd.x,hjdU.x],
i$===k,If[MemberQ[UpList,i$],-x.hjdU,x.hUjd],
True,Transpose[Inner[Times,x,If[MemberQ[UpList,i$],-hjdU,hUjd],Plus,i$],OrdN[Join[Range[i$-1],RotateRight[Range[i$,k]]]]]],{i$,k}],0]]];


OrdN[x_List]:=Flatten[Position[x,#]&/@Sort[x]]


Laplacian[Xx^2+Yy^2,{Xx,Yy}];Unprotect[Laplacian];Laplacian[x_,UpList_:{}]:=covDiv[covD[x,UpList].gUU,Join[UpList,{{TensRnk[x]+1}}]]/;UpList==={}||Union[Head/@Flatten[UpList]]==={Integer};s1$=DownValues[Laplacian];s2$=Position[s1$,UpList];If[s2$[[1,1]]=!=1,DownValues[Laplacian]=Join[{s1$[[s2$[[1,1]]]]},Drop[s1$,{s2$[[1,1]]}]]];Clear[s1$,s2$];Off[Laplacian::argtu];Protect[Laplacian];


Grad2Norm[x_,y_:Null]:=FaSi[If[y===Null,gUU.covD[x].covD[x],gUU.covD[x].covD[y]]]/;Head[x]=!=List&&Head[y]=!=List;


Bianchi[0]:=Block[{tmp$=covD[Rdddd]},indep$cTerms[FaSi[indep$cTerms[tmp$+Transpose[tmp$,{1,2,4,5,3}]+Transpose[tmp$,{1,2,5,3,4}]]]]];


Bianchi[1]:=Block[{tmp$=covD[Rdd]},If[Head[RUddd]===Symbol,RUddd=Raise[Rdddd,1]];indep$cTerms[FaSi[indep$cTerms[covDiv[RUddd,{1}]+tmp$-Transpose[tmp$,{1,3,2}]]]]];


Bianchi[2]:=indep$cTerms[covDiv[EUd,{1}]];


Clear$dx:=Which[Head[coordList]===Symbol,Print["OK"];Return[],Union[Head /@ d /@ coordList,{d,Symbol}]==={d,Symbol},If[dxRuleList==={},Clear[dxRuleList,deList];If[Length[varList]>0,stPrint["varList cleared"]];stPrint["deList and dxRuleList also cleared"],Print["OK"]];Return[],True,Block[{var$,z$,i$=1},While[i$<=Dim,var$=coordList[[i$]];z$="d"<>ToString[var$];If[Head[d[var$]]=!=d,(d[var$])=.];If[FormDegree[ToExpression[z$]]===1,ToExpression[z$<>"=."];d[var$]=ToExpression[z$]];i$=i$+1]];If[Length[varList]===0&&Head[dxRuleList]===Symbol&&Head[deList]===Symbol,
Print["OK"]];If[Length[varList]>0,stPrint["varList cleared"]];If[Head[dxRuleList]===List&&Head[deList]===List,stPrint["deList and dxRuleList also cleared"]];Clear[eTO$dx,dxRuleList,deList];varList={};If[Union[Head/@d/@e/@Range[Dim]]=!={d},(d[e[a_]])=.]];


metric[LinEl_,x_List]:=FacSimp[Which[FormDegree[x]===0&&Length[indepTerms[d[x]/.e[_]->0]]===Length[x],Array[D[LinEl/2, d[x[[#1]]], d[x[[#2]]]]&,{Length[x],Length[x]}],x===e/@Range[Length[x]]||(indepTerms[d[x]]==={}&&Length[indepTerms[x/.e[_]->0]]===Length[x]),Array[D[LinEl/2,x[[#1]],x[[#2]]]&,{Length[x],Length[x]}],FormDegree[x]===0&&Length[indepTerms[d[x]/.e[_]->0]]<Length[x],stPrint["The metric components are given with respect to the coframe {e[1],...,e[n}}"];Array[D[LinEl/2,e[#1],e[#2]]&,{Length[x],Length[x]}],True,stPrint["The second argument must be either the list of coordinates, or the coframe list {e[1],...,e[n}}"]]];


dot$2[x_,y_]:=multiDot[x,y,{3,1},{4,2}];tr$2[x_]:=Contract[x,{1,3},{2,4}];


Plebanski[x_,y_]:=Module[{R2d4=Transpose[Outer[Times,x,y]+Outer[Times,y,x],{1,3,2,4}],Rg=Transpose[Outer[Times,tr[y.gUU]x+tr[x.gUU]y,gdd],{1,3,2,4}],gg=Transpose[Outer[Times,gdd,gdd],{1,3,2,4}],Rsq=FaSi[x.gUU.y+y.gUU.x],trRsq,R2g,Rscal},If[!zeroQ[FaSi[x-Transpose[x]]]||!zeroQ[FaSi[y-Transpose[y]]],stPrint["Not a symmetric tensor!"];Return[]];R2d4=R2d4-Transpose[R2d4];Rg=Rg-Transpose[Rg];Rg=Rg-Transpose[Rg,{1,2,4,3}];gg=gg-Transpose[gg];trRsq=Contract[Rsq.gUU,{1,2}];R2g=Transpose[Outer[Times,Rsq,gdd],{1,3,2,4}];R2g=R2g-Transpose[R2g];R2g=R2g-Transpose[R2g,{1,2,4,3}];Rscal=Contract[gUU.Rg.gUU/(2(Dim-1)),{1,3},{2,4}];FaSi[R2d4- (Rg-R2g)/(Dim-2)+(Rscal-trRsq)/((Dim-1)(Dim-2)) gg]];


Plebanski[x_]:=Module[{R2d4=Transpose[Outer[Times,x,x],{1,3,2,4}],Rg=Transpose[Outer[Times,x,gdd],{1,3,2,4}],gg=Transpose[Outer[Times,gdd,gdd],{1,3,2,4}],Rsq=FaSi[x.gUU.x],trRsq,R2g,Rscal=Contract[x.gUU,{1,2}]},If[!zeroQ[FaSi[x-Transpose[x]]],stPrint["Not a symmetric tensor!"];Return[]];R2d4=R2d4-Transpose[R2d4];Rg=Rg-Transpose[Rg];Rg=Rg-Transpose[Rg,{1,2,4,3}];gg=gg-Transpose[gg];trRsq=Contract[Rsq.gUU,{1,2}];R2g=Transpose[Outer[Times,Rsq,gdd],{1,3,2,4}];R2g=R2g-Transpose[R2g];R2g=R2g-Transpose[R2g,{1,2,4,3}];FaSi[R2d4- (Rscal Rg-R2g)/(Dim-2)+(Rscal^2-trRsq)/((Dim-1)(Dim-2)) gg]];


ClassifyMetric[x_]:=Which[Union[Flatten[Dimensions[x]]]=!={4},stPrint["Classification defined for 4-d tensors only!"];Return[],Head[x]=!=List||Dimensions[x]=!={4,4}&&Dimensions[x]=!={4,4,4,4},stPrint["Undefined or incorrect-rank tensor!"];Return[],Dimensions[x]==={4,4,4,4},If[zeroQ[FaSi[x+Transpose[x,{1,2,4,3}]]]&&zeroQ[FaSi[x+Transpose[x,{2,1,3,4}]]]&&zeroQ[FaSi[x-Transpose[x,{3,4,1,2}]]]&&zeroQ[Contract[x.gUU,{2,4}]],ClassifyWeyl[x],stPrint["Not a Weyl-type tensor!"];Return[]],Dimensions[x]==={4,4},If[zeroQ[FaSi[x-Transpose[x]]],tmp$=FaSi[Plus@@Transpose[x.gUU/2,{1,1}]];If[zeroQ[tmp$],ClassifyRic09[x],ClassifyRic10[x,tmp$]],stPrint["Not a symmetric tensor!"];Return[]]];


indx$Ruls[1]={{#1->1,#2->2},{#1->1,#2->3},{#1->1,#2->4},{#1->2,#2->3},{#1->2,#2->4},{#1->3,#2->4}};indx$Ruls[2]=indx$Ruls[1]/.{#1->#3,#2->#4};indx$Ruls[12]=Array[Join[indx$Ruls[1][[#1]],indx$Ruls[2][[#2]]]&,{6,6}];


ClassifyWeyl[x_]:=If[zeroQ[x],Print["Type O"];Return[{0}],Module[{gg=Transpose[Outer[Times,gdd,gdd],{1,3,2,4}],j$,ID$4},ID$4= (Raise[gg-Transpose[gg],3,4]-j$ eta[3,4])/4;gg=indep$cTerms[Collect[PowerExpand[Union[Flatten[dot$2[ID$4,ID$4]-ID$4]]],j$,FaSi]];Which[zeroQ[FaSi[gg/.j$->I]],j$=I,zeroQ[FaSi[gg/.j$->1]],j$=1,True,stPrint["Cannot perform classification; try different coordinates or frame!"];Return[]];Off[General::pspec];ID$4=ID$4[[#1,#2,#3,#4]]/.indx$Ruls[12];CW$[1]=FaSi[-ID$4.(Raise[x,3,4][[#1,#2,#3,#4]]/.indx$Ruls[12])];On[General::pspec];CW$[2]=FaSi[2CW$[1].CW$[1]];Is$2=tr[CW$[2]];Js$3=FaSi[Plus@@FaSi[Transpose[2/3CW$[1].CW$[2],{1,1}]]];Which[zeroQ[CW$[1]],Print["Type O"];Return[{0}],zeroQ[CW$[2]],Print["Type N"];Return[{4}],zeroQ[{Is$2,Js$3}],Print["Type III"];Return[{31}],zeroQ[FaSi[Is$2^3-27Js$3^2]],CW$[20]=FaSi[CW$[2]-2/3 Is$2 ID$4];If[zeroQ[FaSi[Is$2*CW$[20]-3 Js$3*CW$[1]]],Print["Type D"];Return[{22}],Print["Type II"];Return[{211}]],True,Which[zeroQ[Is$2],Print["Type I, invariant I = 0"],zeroQ[Js$3],Print["Type I, invariant J = 0"],True,Print["Type I - G"]];Return[{1111}]]]];


testSolve[x_,y_]:=Block[{$ol=Union[Solve[x==0,y]]},If[$ol==={},Return[False]];If[zeroQ[FaSi[x/.$ol]],True,False]];


ClassifyRic09[x_]:=Module[{ID$4=IdentityMatrix[4],trR,al$p,tmp$,b4$,c6$},SM$[1]=FaSi[x.gUU];SM$[2]=FaSi[SM$[1].SM$[1]];SM$[3]=FaSi[SM$[2].SM$[1]];trR[2]=tr[SM$[2]/2];trR[3]=tr[SM$[3]/2];trR[4]=FaSi[Plus@@FaSi[Transpose[1/2 SM$[2].SM$[2],{1,1}]]];Print["Traceless"];al$p=zeroQ/@{FaSi[trR[2]^2-trR[4]],trR[3],trR[2]};Which[al$p==={True,True,True},tmp$="Quadruple 0",al$p==={True,True,False},tmp$="Double 0",al$p==={True,False,False},tmp$="Single 0",True,tmp$=""];Clear[al$p];If[StringLength[tmp$]>0,Print[tmp$<>" eigenvalue"]];(* Type IV **)If[zeroQ[trR/@{2,3,4}],Which[zeroQ[SM$[1]],Print["Segre [(111,1)] (Vacuum)"]; Return["IV-M1(0)"],zeroQ[SM$[2]],Print["Segre [(11,2)] (Null EM field)"]; Return["IV-M2(00)"],zeroQ[SM$[3]],Print["Segre [(1,3)]"]; Return["IV-M3(000)"],True,Print["Segre [(4)]??"]; Return["IV-G(0000)"]]];(* Type IID **)If[zeroQ[{trR[3],FaSi[-trR[2]^2+2trR[4]]}],Which[testSolve[Union[Flatten[SM$[2]-al$p[2]ID$4]],al$p[2]],Print["Segre [(11)(1,1)] (Non-Null EM field)"]; Return["IID-M2(01)"],testSolve[Union[Flatten[SM$[3]-al$p[1]SM$[2]-al$p[1]^2SM$[1]+al$p[1]^3ID$4]],al$p[1]],Print["Segre [(11),2]"]; Return["IID-M3(111)"],True,Print["Segre [2,2]??"]; Return["IID-G(0101)"]]];
b4$=FaSi[-7trR[2]^2+6trR[4]];c6$=FaSi[-2trR[2]^3+3trR[3]^2];
(* Type III **)If[zeroQ[{b4$,c6$}],Which[testSolve[Union[Flatten[SM$[2]+2al$p[1]SM$[1]-3al$p[1]^2ID$4]],al$p[1]],Print["Segre [(111),1] or [1(11,1)] (radiation-type perfect fluid)"]; Return["III-M2(11)"],testSolve[Union[Flatten[SM$[3]+al$p[1]SM$[2]-5al$p[1]^2SM$[1]+3al$p[1]^3ID$4]],al$p[1]],Print["Segre [1(1,2)]"]; Return["III-M3(111)"],True,Print["Segre [1,3]"]; Return["III-G(0111)"]]];
(* Type II **)If[zeroQ[FaSi[b4$^3+(2c6$-3trR[2]b4$)^2]],If[testSolve[Union[Flatten[SM$[3]+al$p[1] SM$[2]-(2al$p[1]^2+2al$p[1]al$p[2]+al$p[2]^2)SM$[1]+al$p[1]al$p[2](2al$p[1]+al$p[2])ID$4]],al$p/@{1,2}],Print["Segre [11(1,1)] or [(11)1,1] or [(11),Z\!\(\*OverscriptBox[\"Z\", \"_\"]\)]"]; Which[testSolve[Union[Flatten[SM$[3]-al$p[2]SM$[1]]],al$p[2]],Return["II-M3(010)"],testSolve[Union[Flatten[SM$[3]+al$p[1]SM$[2]-2al$p[1]^2SM$[1]]],al$p[1]],Return["II-M3(110)"],True,Return["II-M3(111)"]],Print["Segre [11,2]"]; Which[zeroQ[{trR[3],FaSi[-trR[2]^2+trR[4]]}],Return["II-G(0100)"],zeroQ[{trR[2],trR[4]}],Return["II-G(0010)"],zeroQ[{trR[2],trR[3]}],Return["II-G(0001)"],zeroQ[trR[2]],Return["II-G(0011)"],zeroQ[trR[3]],Return["II-G(0101)"],zeroQ[FaSi[-trR[2]^2+trR[4]]],Return["II-G(0110)"],True,Return["II-G(0111)"]]]];(* Type I **)Print["Segre [111,1] or [11,Z\!\(\*OverscriptBox[\"Z\", \"_\"]\)]"];Which[zeroQ[{trR[2],trR[3]}],Return["I-G(0001)"],zeroQ[{trR[2],trR[4]}],Return["I-G(0010)"],zeroQ[trR[2]],Return["I-G(0011)"],zeroQ[trR[3]],Return["I-G(0101)"],zeroQ[FaSi[-trR[2]^2+trR[4]]],Return["I-G(0110)"],True,Return["I-G(0111)"]]];


ClassifyRic10[x_,y_]:=Module[{ID$4=IdentityMatrix[4],trR,al$p,tmp$,b4$,c6$},SM$[1]=FaSi[x.gUU];SM$[2]=FaSi[SM$[1].SM$[1]];SM$[3]=FaSi[SM$[2].SM$[1]];trR[1]=y;trR[2]=tr[SM$[2]/2];trR[3]=tr[SM$[3]/2];trR[4]=FaSi[Plus@@FaSi[Transpose[1/2 SM$[2].SM$[2],{1,1}]]];al$p=zeroQ/@FaSi[{4trR[1]^4-12trR[1]^2trR[2]+3trR[2]^2+8trR[1]trR[3]-3trR[4],2trR[1]^3-3trR[1]trR[2]+trR[3],2trR[1]^2-trR[2],trR[1]}];Which[al$p==={True,True,True,True},tmp$="Quadruple 0",al$p==={True,True,True,False},tmp$="Triple 0",al$p==={True,True,False,False},tmp$="Double 0",al$p==={True,False,False,False},tmp$="Single 0",True,tmp$=""];Clear[al$p];If[StringLength[tmp$]>0,Print[tmp$<>" eigenvalue"]];(* Type IV **)If[zeroQ[FaSi[{-trR[1]^2+2trR[2],-trR[1]^3+4trR[3],-trR[1]^4+8trR[4]}]],Which[zeroQ[FaSi[SM$[1]-trR[1]/2 ID$4]],Print["Segre [(111,1)] (Lamda Term)"]; Return["IV-M1(1)"],zeroQ[FaSi[SM$[2]-trR[1] SM$[1]+trR[1]^2/4 ID$4]],Print["Segre [(11,2)]"]; Return["IV-M2(11)"],zeroQ[FaSi[SM$[3]-3/2trR[1]SM$[2]+3/4trR[1]^2SM$[1]-1/8trR[1]^3ID$4]],Print["Segre [(1,3)]"]; Return["IV-M3(111)"],True,Print["Segre [(4)]??"]; Return["IV-G(1111)"]]];(* Type IID **)If[zeroQ[FaSi[{trR[1]^3-3trR[1]trR[2]+2trR[3],trR[1]^4-2trR[1]^2trR[2]-trR[2]^2+2trR[4]}]],Which[testSolve[Union[Flatten[SM$[2]-(al$p[1]+al$p[2])SM$[1]+al$p[1]al$p[2]ID$4]],al$p/@{1,2}],Print["Segre [(11)(1,1)]"]; If[testSolve[Union[Flatten[SM$[2]-al$p[2]SM$[1]]],al$p[2]],Return["IID-M2(10)"],Return["IID-M2(11)"]],testSolve[Union[Flatten[SM$[3]-(2al$p[1]+al$p[2])SM$[2]+(2al$p[1]al$p[2]+al$p[1]^2)SM$[1]-al$p[1]^2al$p[2]ID$4]],al$p/@{1,2}],Print["Segre [(11),2]"]; Which[testSolve[Union[Flatten[SM$[3]-al$p[2]SM$[2]]],al$p[2]],Return["IID-M3(100)"],testSolve[Union[Flatten[SM$[3]-2al$p[1]SM$[2]+al$p[1]^2SM$[1]]],al$p[1]],Return["IID-M3(110)"],testSolve[Union[Flatten[SM$[3]-3al$p[1]^2SM$[1]+2al$p[1]^3ID$4]],al$p[1]],Return["IID-M3(011)"],testSolve[Union[Flatten[SM$[3]+3al$p[2]SM$[2]-4al$p[2]^3ID$4]],al$p[2]],Return["IID-M3(101)"],True,Return["IID-M3(111)"]],True,Print["Segre [2,2]??"]; Return["IID-G(1111)"]]];
b4$=FaSi[-4trR[1]^4+16trR[1]^2trR[2]-7trR[2]^2-12trR[1]trR[3]+6trR[4]];c6$=FaSi[4trR[1]^6-24trR[1]^4trR[2]+39trR[1]^2trR[2]^2-8trR[2]^3+12trR[1]^3trR[3]-36trR[1]trR[2]trR[3]+12trR[3]^2];
(* Type III **)
If[zeroQ[{b4$,c6$}],Which[testSolve[Union[Flatten[SM$[2]-(al$p[1]+al$p[2])SM$[1]+al$p[1]al$p[2]ID$4]],al$p/@{1,2}],Print["Segre [(111),1] or [1(11,1)] (Perfect fluid)"]; Which[testSolve[Union[Flatten[SM$[2]-al$p[2]ID$4]],al$p[2]],Return["III-M2(01)"],testSolve[Union[Flatten[SM$[2]-al$p[2]SM$[1]]],al$p[2]],Return["III-M2(10)"],True,Return["III-M2(11)"]],testSolve[Union[Flatten[SM$[3]-(2al$p[1]+al$p[2])SM$[2]+(2al$p[1]al$p[2]+al$p[1]^2)SM$[1]-al$p[1]^2al$p[2]ID$4]],al$p/@{1,2}],Print["Segre [1(1,2)]"]; Which[testSolve[Union[Flatten[SM$[3]-al$p[2]SM$[2]]],al$p[2]],Return["III-M3(100)"],testSolve[Union[Flatten[SM$[3]-2al$p[1]SM$[2]+al$p[1]^2SM$[1]]],al$p[1]],Return["III-M3(110)"],testSolve[Union[Flatten[SM$[3]-3al$p[1]^2SM$[1]+2al$p[1]^3ID$4]],al$p[1]],Return["III-M3(011)"],testSolve[Union[Flatten[SM$[3]+3al$p[2]SM$[2]-4al$p[2]^3ID$4]],al$p[2]],Return["III-M3(101)"],True,Return["III-M3(111)"]],True,Print["Segre [1,3]"]; Return["III-G(1111)"]]];(* Type II **)If[zeroQ[FaSi[4b4$^3+(c6$-3(-trR[1]^2+2trR[2])b4$)^2]],If[testSolve[Union[Flatten[SM$[3]-(al$p[1]+al$p[2]+al$p[3])SM$[2]+(al$p[1]al$p[2]+al$p[2]al$p[3]+al$p[3]al$p[1])SM$[1]-al$p[1]al$p[2]al$p[3]ID$4]],al$p/@{1,2,3}],Print["Segre [11(1,1)] or [(11)1,1] or [(11),Z\!\(\*OverscriptBox[\"Z\", \"_\"]\)]"]; Which[testSolve[Union[Flatten[SM$[3]-al$p[2]ID$4]],al$p[2]],Return["II-M3(001)"],(*Impossible*)testSolve[Union[Flatten[SM$[3]-al$p[2]SM$[1]]],al$p[2]],Return["II-M3(010)"],testSolve[Union[Flatten[SM$[3]-al$p[2]SM$[2]]],al$p[2]],Return["II-M3(100)"],(*Impossible*)testSolve[Union[Flatten[SM$[3]-(al$p[1]^2+al$p[1]al$p[2]+al$p[2]^2)SM$[1]+al$p[1]al$p[2](al$p[1]+al$p[2])ID$4]],al$p/@{1,2}],Return["II-M3(011)"],testSolve[Union[Flatten[SM$[3]-(al$p[1]-al$p[2])SM$[2]+al$p[1]al$p[2]^2ID$4]],al$p/@{1,2}],Return["II-M3(101)"],testSolve[Union[Flatten[SM$[3]-2al$p[2]SM$[2]+al$p[2]^2SM$[1]]],al$p[2]],Return["II-M3(110)"],True,Return["II-M3(111)"]],Print["Segre [11,2]"]; Which[zeroQ[FaSi[{-4trR[1]^3+trR[3],-2trR[1]^2+trR[2]}]],Return["II-G(1001)"],zeroQ[FaSi[{-2trR[1]^2+trR[2],8trR[1]^4-8trR[1]trR[3]+3trR[4]}]],Return["II-G(1010)"],zeroQ[FaSi[{4trR[1]^4-4trR[1]^2trR[2]-trR[2]^2+trR[4],2trR[1]^3-3trR[1]trR[2]+trR[3]}]],Return["II-G(1100)"],zeroQ[FaSi[2trR[1]^2-trR[2]]],Return["II-G(1011)"],zeroQ[FaSi[2trR[1]^3-3trR[1]trR[2]+trR[3]]],Return["II-G(1101)"],zeroQ[FaSi[4trR[1]^4-12trR[1]^2trR[2]+3trR[2]^2+8trR[1]trR[3]-3trR[4]]],Return["II-G(1110)"],True,Return["II-G(1111)"]]]];(* Type I **)Print["Segre [111,1] or [11,Z\!\(\*OverscriptBox[\"Z\", \"_\"]\)]"];Which[zeroQ[FaSi[{-4trR[1]^3+trR[3],-2trR[1]^2+trR[2]}]],Return["I-G(1001)"],zeroQ[FaSi[{-2trR[1]^2+trR[2],8trR[1]^4-8trR[1]trR[3]+3trR[4]}]],Return["I-G(1010)"],zeroQ[FaS[2trR[1]^2-trR[2]]],Return["I-G(1011)"],zeroQ[FaSi[2trR[1]^3-3trR[1]trR[2]+trR[3]]],Return["I-G(1101)"],zeroQ[FaSi[4trR[1]^4-12trR[1]^2trR[2]+3trR[2]^2+8trR[1]trR[3]-3trR[4]]],Return["I-G(1110)"],True,Return["I-G(1111)"]]];


$PrePrint=If[ByteCount[#]>200000,Short[#,100],#]&;


Protect[cDcomp,CF3d,ConOne,c$D,dot$2,LieDerv,dotAll,dotOne,epsilon,FacSimp,FaSi,GenCoef,indepTerms,indep$cTerms,indexList,indx$Ruls,inn,metric,M$F,MoveIndex,newPrList,numb1,numb2,OKindCon,redPlus,redTimes,RiemSym,zeroTensor,testSER,testSolve,tr$2,OrdN];


On[General::"spell",General::"spell1",UpSet::"write"];



