(* ::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]]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]]]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$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]]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]]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"];