(* ::Package:: *)

(* ::Subsubtitle:: *)
(* ccgrg.m*)


(* ::Text:: *)
(*Andrzej Woszczyna*)
(*Copernicus Center for Interdisciplinary Studies, Krak\[OAcute]w*)


(* ::Text:: *)
(*in collaboration with: *)
(*Wojciech Czaja,  Krzysztof Glod, Zdzislaw Golda, Rados\[LSlash]aw Kycia, Andrzej Odrzywo\[LSlash]ek, Piotr Plaszczyk, Lech Sokolowski,  Sebastian Szybka*)


(* ::Section::Plain::Closed:: *)
(*Start*)


(* ::Subsection::Plain::Closed:: *)
(*begin*)


BeginPackage["ccgrg`"];


(* ::Subsection::Plain::Closed:: *)
(*names*)


(* ::Subsubsection::Plain::Closed:: *)
(*auxiliary*)


swap;MF;TF;


(* ::Subsubsection::Plain::Closed:: *)
(*general*)


(* ::Text:: *)
(*package version:*)


ccgrg;version;ccver;


(* ::Text:: *)
(*series:*)


trunc;pertOrder;ord\[DoubleDagger];pertVar;pert\[DoubleDagger];


(* ::Text:: *)
(*simplification:*)


simplification;smp;simp;


denominatorlist; collectdenominators;


(* ::Subsubsection::Plain::Closed:: *)
(*memory*)


(* memoization controll *)
tasama;erasable;memRegistry;clr;
cacheview;retreat;mviewall;fun\[DoubleDagger];associated;


(* ::Subsubsection::Plain::Closed:: *)
(*tensor and geometry*)


(* tensor identifier *)
id\[DoubleDagger];


(* coordinates *)
x\[DoubleDagger];coordinateExt;


(* index *)
klamra;
dim;all;all2;all3;var;varL;
idx;idX;indeX;idxx;idXX;


(* conditions *)
sub(*podstawienia*);
auxi(*pomocnicze*);
restrict(*warunki*);
unspecified;RealFunctionsAll;verify;addassumption;


(* metrics *)
g\[DoubleDagger];metricExt;metricDet;gAux;auxiliarymetric;chess;toMatrix;
sqrtdet;metricsign;


(* vector *)
vectorsquared;unitvector;vctr;


(* tensor *)
functionName;
tensorExt;tExt;tensorAux;tensorQ;
setrank;rank;rankT;countUp;countUp2;zerosUp;zerosUp2;
span;tabular;tensor;


(* connection *)
\[CapitalGamma];\[CapitalGamma]Aux;auxiliaryconnection;
pion;


(* antisymmetric *)
\[Eta]\[DoubleDagger];etaExt;


(* hypersurfaces *)
h\[DoubleDagger];projectionExt;\[Chi]\[DoubleDagger];IImetricExt;


(* curvature *)
Ric;tRicciR;tEinsteinG;G\[DoubleDagger];Riem;tRiemannR;tWeylC;tWeylCdual;sEinsteinG;


(* derivatives  *)
orderD;
covariantDAll;covariantderivative;partialD;partialDS;covariantD;\[EmptyDownTriangle];


(* invariants *)
S\[DoubleDagger];tPlebanskiS;CMinvR1;CMinvR2;CMinvR3;CMinvM1;CMinvM2;CMinvM3;CMinvM4;CMinvM5;CMinvW1;CMinvW2;


(* geometrical procedures *)
open;antisymmetric;hypersurfaces;curvature;CMinvariants;


(* geodesic *)
geo;geodesic;


(* ::Subsubsection::Plain::Closed:: *)
(*general relativity*)


Schwarzschild;Friedman;LTB; 


velocityU;fourvelocity;timelikeEigenvector;timelikeV;TimelikeEigenvector1;TimelikeEigenvector2;TimelikeEigenvectorG;


energydensity;isotropicpressure;


positiveEM;


(* ::Subsubsection::Plain::Closed:: *)
(*vector field*)


valuesG;
gradU;exp\[Theta];scalar\[Theta];rot\[Omega];shear\[Sigma];\[Omega]\[Omega];\[Sigma]\[Sigma];
hydrodynamics;vectorfield;flow;


tProjectiveShear;t\[Sigma];
sProjectiveShear;s\[Sigma];
tProjectiveExpansion;t\[Theta];
sProjectiveExpansion;s\[Theta];
tProjectiveVorticity;t\[Omega];
sProjectivelVorticity;s\[Omega];

t\[Theta]cov;t\[Omega]cov;


(* ::Subsection::Plain::Closed:: *)
(*usage*)


(* ::Subsubsection::Plain::Closed:: *)
(*opening sesion*)


open::usage = "open[x, g]
open[x, g, simplification->Simplify]"


(* ::Subsubsection::Plain::Closed:: *)
(*tensor definition*)


tensorExt::usage ="
T[i_,j_]:=\!\(\*StyleBox[\"tensorExt\",\nFontSlant->\"Italic\"]\)[\!\(\*SubscriptBox[\(T\), \(cov\)]\)][i,j]

T[i_,j_]:=\!\(\*StyleBox[\"tensorExt\",\nFontSlant->\"Italic\"]\)[\!\(\*SubscriptBox[\(T\), \(cov\)]\),valence->{\"d\",\"d\"}][i,j]

T[i_,j_]:=\!\(\*StyleBox[\"tensorExt\",\nFontSlant->\"Italic\"]\)[\!\(\*SubscriptBox[\(T\), \(mix\)]\),valence->{\"u\",\"d\"}][i,j]"


rank::usage ="\!\(\*StyleBox[\"rank\",\nFontSlant->\"Italic\"]\)[T] gives the tensor rank of T"


tensorQ::usage ="\!\(\*StyleBox[\"tensorQ\",\nFontSlant->\"Italic\"]\)[T] gives True if T is a tensor"


id\[DoubleDagger]::usage ="tensor identifier: T[i,0,m...] =\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"id\[DoubleDagger]\",\nFontSlant->\"Italic\"]\)"


idx::usage ="Tcov[i_/;\!\(\*StyleBox[\"idx\",\nFontSlant->\"Italic\"]\)[i]]:=expr[x][i] 

covariant index condition:  \!\(\*StyleBox[\"idx\",\nFontSlant->\"Italic\"]\)[i] gives True for Integer i: 0<i<=dim";


idX::usage ="T[i_/;\!\(\*
StyleBox[\"idX\",\nFontSlant->\"Italic\"]\)[i]]:=\!\(\*
StyleBox[\"tensorExt\",\nFontSlant->\"Italic\"]\)[Tcov][i] 
 
index range condition: \!\(\*StyleBox[\"idX\",\nFontSlant->\"Italic\"]\)[i] gives True if for Integer i:  0<Abs[i]<=dim";


indeX::usage ="T[i_,j_,k_]/;\!\(\*
StyleBox[\"indeX\",\nFontSlant->\"Italic\"]\)[i,j,k]:=\!\(\*
StyleBox[\"tensorExt\",\nFontSlant->\"Italic\"]\)[Tcov][i,j,k]

multiple index range condition: \!\(\*StyleBox[\"indeX\",\nFontSlant->\"Italic\"]\)[i,j,k] gives True if for each Integer i,j,k:  0<=Abs[i_]<=dim.";


tabular::usage ="R=\!\(\*
StyleBox[\"tabular\",\nFontSlant->\"Italic\"]\)[tRicciR],  then R[\!\(\*
StyleBox[\"all\",\nFontSlant->\"Italic\"]\),\!\(\*
StyleBox[\"all\",\nFontSlant->\"Italic\"]\)] and R[-\!\(\*
StyleBox[\"all\",\nFontSlant->\"Italic\"]\),\!\(\*
StyleBox[\"all\",\nFontSlant->\"Italic\"]\)] are matrixes";


(* ::Subsubsection::Plain::Closed:: *)
(*metric*)


x\[DoubleDagger]::usage = "coordinates:
x\[DoubleDagger][1],x\[DoubleDagger][2]...(covariant), x\[DoubleDagger][-1],x\[DoubleDagger][-2]...(contravariant)"


g\[DoubleDagger]::usage = "Extended metric tensor: g\[DoubleDagger][m,n] =\!\(\*SubscriptBox[\(g\), \(mn\)]\),  g\[DoubleDagger][-m,-n] =\!\(\*SuperscriptBox[\(g\), \(mn\)]\),  g\[DoubleDagger][m,-n] =\!\(\*SuperscriptBox[SubscriptBox[\(\[Delta]\), \(m\)], \(n\)]\)
g\[DoubleDagger] is short name for metricExt. (\[DoubleDagger] is used only with the single-letter tensor names.)";


metricExt::usage = "Extended metric tensor: 
	\!\(\*
StyleBox[\"metricExt\",\nFontSlant->\"Italic\"]\)[m,n] = \!\(\*SubscriptBox[\(g\), \(mn\)]\) 
	\!\(\*
StyleBox[\"metricExt\",\nFontSlant->\"Italic\"]\)[-m,-n] = \!\(\*SuperscriptBox[\(g\), \(mn\)]\)
	\!\(\*
StyleBox[\"metricExt\",\nFontSlant->\"Italic\"]\)[m,-n] = \!\(\*SuperscriptBox[SubscriptBox[\(\[Delta]\), \(m\)], \(n\)]\)
\!\(\*
StyleBox[\"metricExt\",\nFontSlant->\"Italic\"]\) is equivalent to g\[DoubleDagger] ";


metricDet::usage ="metric tensor determinant";


metricsign::usage ="Sign[\!\(\*
StyleBox[\"metricDet\",\nFontSlant->\"Italic\"]\)]";


sqrtdet::usage ="\!\(\*
StyleBox[\"sqrtdet\",\nFontSlant->\"Italic\"]\)=\!\(\*SqrtBox[\(\*
StyleBox[\"metricDet\",\nFontSlant->\"Italic\"]\\\ Sign[\*
StyleBox[\"metricDet\",\nFontSlant->\"Italic\"]]\)]\)"


metricSet::usage ="Schwarzschild, Friedman, LTB";


(* ::Subsubsection::Plain::Closed:: *)
(*connection*)


\[CapitalGamma]::usage = "Christoffel symbol: \[CapitalGamma][s,m,n] = \!\(\*SubscriptBox[\(\[CapitalGamma]\), \(smn\)]\),  \[CapitalGamma][-s,m,n] = \!\(\*SubscriptBox[SuperscriptBox[\(\[CapitalGamma]\), \(s\)], \(mn\)]\)";


(* ::Subsubsection::Plain::Closed:: *)
(*curvature*)


tRicciR::usage = "Ricci tensor:  \!\(\*
StyleBox[\"tRicciR\",\nFontSlant->\"Italic\"]\)[m, n] = \!\(\*SubscriptBox[\(R\), \(mn\)]\), \!\(\*
StyleBox[\"tRicciR\",\nFontSlant->\"Italic\"]\)[-m,-n] = \!\(\*SuperscriptBox[\(R\), \(mn\)]\), \!\(\*
StyleBox[\"tRicciR\",\nFontSlant->\"Italic\"]\)[m,-n] = \!\(\*SuperscriptBox[SubscriptBox[\(R\), \(m\)], \(n\)]\)";


sRicciR::usage = "Ricci scalar";


tEinsteinG::usage = "Einstein tensor:  \!\(\*
StyleBox[\"tEinsteinG\",\nFontSlant->\"Italic\"]\)[m, n] = \!\(\*SubscriptBox[\(G\), \(mn\)]\), \!\(\*
StyleBox[\"tEinsteinG\",\nFontSlant->\"Italic\"]\)[-m,-n]= \!\(\*SuperscriptBox[\(G\), \(mn\)]\), \!\(\*
StyleBox[\"tEinsteinG\",\nFontSlant->\"Italic\"]\)[m,-n] = \!\(\*SuperscriptBox[SubscriptBox[\(G\), \(m\)], \(n\)]\)";
G\[DoubleDagger]::usage = "G\[DoubleDagger] is short name for tEinsteinG. (\[DoubleDagger] is used only with the single-letter tensor names.)";


tRiemannR::usage ="Riemann tensor: \!\(\*
StyleBox[\"tRiemannR\",\nFontSlant->\"Italic\"]\)[i,j,m,n] = \!\(\*SubscriptBox[\(R\), \(ijmn\)]\), 
	\!\(\*
StyleBox[\"tRiemannR\",\nFontSlant->\"Italic\"]\)[-i, j,m,n] = \!\(\*SubscriptBox[SuperscriptBox[\(R\), \(i\)], \(jmn\)]\),   \!\(\*
StyleBox[\"tRiemannR\",\nFontSlant->\"Italic\"]\)[i,-j, m,n] = \!\(\*SubscriptBox[SuperscriptBox[SubscriptBox[\(R\), \(i\)], \(j\)], \(mn\)]\), ...
	\!\(\*
StyleBox[\"tRiemannR\",\nFontSlant->\"Italic\"]\)[-i,-j,m,n] = \!\(\*SubscriptBox[SuperscriptBox[\(R\), \(ij\)], \(mn\)]\),  \!\(\*
StyleBox[\"tRiemannR\",\nFontSlant->\"Italic\"]\)[i,-j,-m,n] = \!\(\*SubscriptBox[SuperscriptBox[SubscriptBox[\(R\), \(i\)], \(jm\)], \(n\)]\), ...
	(...)
	\!\(\*
StyleBox[\"tRiemannR\",\nFontSlant->\"Italic\"]\)[-i,-j,-m,-n] = \!\(\*SuperscriptBox[\(R\), \(ijmn\)]\) ";


tWeylC::usage ="Weyl tensor: \!\(\*
StyleBox[\"tWeylC\",\nFontSlant->\"Italic\"]\)[i,j,m,n] =\!\(\*SubscriptBox[\(C\), \(ijmn\)]\), 
	\!\(\*
StyleBox[\"tWeylC\",\nFontSlant->\"Italic\"]\)[-i, j,m,n] = \!\(\*SubscriptBox[SuperscriptBox[\(C\), \(i\)], \(jmn\)]\),   \!\(\*
StyleBox[\"tWeylC\",\nFontSlant->\"Italic\"]\)[i,-j, m,n] = \!\(\*SubscriptBox[SuperscriptBox[SubscriptBox[\(C\), \(i\)], \(j\)], \(mn\)]\), ...
	\!\(\*
StyleBox[\"tWeylC\",\nFontSlant->\"Italic\"]\)[-i,-j,m,n] = \!\(\*SubscriptBox[SuperscriptBox[\(C\), \(ij\)], \(mn\)]\),   \!\(\*
StyleBox[\"tWeylC\",\nFontSlant->\"Italic\"]\)[i,-j,-m,n] = \!\(\*SubscriptBox[SuperscriptBox[SubscriptBox[\(C\), \(i\)], \(jm\)], \(n\)]\), ...
	(...)
	\!\(\*
StyleBox[\"tWeylC\",\nFontSlant->\"Italic\"]\)[-i,-j,-m,-n] = \!\(\*SuperscriptBox[\(C\), \(ijmn\)]\) ";


tWeylCdual::usage ="dual Weyl tensor: \!\(\*
StyleBox[\"tWeylCdual\",\nFontSlant->\"Italic\"]\)[i,j,m,n] =\!\(\*SubscriptBox[\(\!\(\*SuperscriptBox[\(C\), \(*\)]\)\), \(ijmn\)]\)";


CMinvR1::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(R\), \(1\)]\)";
CMinvR2::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(R\), \(2\)]\)";
CMinvR3::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(R\), \(3\)]\)";

CMinvM1::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(M\), \(1\)]\)";
CMinvM2::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(M\), \(2\)]\)";
CMinvM3::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(M\), \(3\)]\)";
CMinvM4::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(M\), \(4\)]\)";
CMinvM5::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(M\), \(5\)]\)";

CMinvW1::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(W\), \(1\)]\)";
CMinvW2::usage ="Carminati-McLenaghan invariant \!\(\*SubscriptBox[\(W\), \(2\)]\)";


(* ::Subsubsection::Plain::Closed:: *)
(*derivativeses*)


covariantD::usage = "\!\(\*
StyleBox[\"covariantD\",\nFontSlant->\"Italic\"]\)[T][{i,j},{k,l,m}] = \!\(\*
StyleBox[\"covariantD\",\nFontSlant->\"Italic\"]\)[T][i,j,k,l,m] =\!\(\*SubscriptBox[\(T\), \(ij; klm\)]\)";
\[EmptyDownTriangle]::usage ="\[EmptyDownTriangle] is short name for covariantD";


partialD::usage = "\!\(\*
StyleBox[\"partialD\",\nFontSlant->\"Italic\"]\)[T][{i,j},{k,l,m}] = \!\(\*
StyleBox[\"partialD\",\nFontSlant->\"Italic\"]\)[T][i,j,k,l,m] =\!\(\*SubscriptBox[\(T\), \(ij, klm\)]\)";


LieD::usage="\!\(\*
StyleBox[\"LieD\",\nFontSlant->\"Italic\"]\)[u][T][i,j] = \!\(\*SubscriptBox[\(\[ScriptCapitalL]\), \(u\)]\)\!\(\*SubscriptBox[\(T\), \(ij\)]\)";


LaplaceOperator::usage ="LaplaceOperator@\[Phi][x]";


(* ::Subsubsection::Plain::Closed:: *)
(*vectors and hypersurfaces*)


vectorsquared::usage ="\!\(\*
StyleBox[\"vectorsquared\",\nFontSlant->\"Italic\"]\)[v] = \!\(\*SubscriptBox[\(v\), \(j\)]\)\!\(\*SuperscriptBox[\(v\), \(j\)]\)";


unitvector::usage = "\!\(\*
StyleBox[\"unitvector\",\nFontSlant->\"Italic\"]\)[v][k] = \!\(\*FractionBox[SuperscriptBox[\(v\), \(k\)], SqrtBox[\(Abs[\*SubscriptBox[\(v\), \(j\)]\\\ \*SuperscriptBox[\(v\), \(j\)]]\)]]\)";


h\[DoubleDagger]::usage = "Induced metrics and the projection tensor on hypersurfaces orthogonal to v;

h\[DoubleDagger][v][i,j] = \!\(\*SubscriptBox[\(h\), \(ij\)]\) = \!\(\*SubscriptBox[\(g\), \(ij\)]\)- \!\(\*FractionBox[\(\*SubscriptBox[\(v\), \(i\)] \*SubscriptBox[\(v\), \(j\)]\), \(\*SubscriptBox[\(v\), \(k\)] \*SuperscriptBox[\(v\), \(k\)]\)]\),      h\[DoubleDagger][v][i,-j] = \!\(\*SuperscriptBox[SubscriptBox[\(h\), \(i\)], \(j\)]\)= \!\(\*SuperscriptBox[SubscriptBox[\(\[Delta]\), \(i\)], \(j\)]\)- \!\(\*FractionBox[\(\*SubscriptBox[\(v\), \(i\)] \*SuperscriptBox[\(v\), \(j\)]\), \(\*SubscriptBox[\(v\), \(k\)] \*SuperscriptBox[\(v\), \(k\)]\)]\),      h\[DoubleDagger][v][-i,-j]= \!\(\*SuperscriptBox[\(h\), \(ij\)]\)= \!\(\*SuperscriptBox[\(g\), \(ij\)]\)- \!\(\*FractionBox[\(\*SuperscriptBox[\(v\), \(j\)] \*SuperscriptBox[\(v\), \(j\)]\), \(\*SubscriptBox[\(v\), \(k\)] \*SuperscriptBox[\(v\), \(k\)]\)]\);

h\[DoubleDagger] is short name for projectionExt. (\[DoubleDagger] is used only with the single-letter tensor names.)";

projectionExt::usage ="\!\(\*
StyleBox[\"projectionExt\",\nFontSlant->\"Italic\"]\) is equivalent to h\[DoubleDagger]";


\[Chi]\[DoubleDagger]::usage = "Second fundamental form on the hypersurfaces orthogonal to v;

\[Chi]\[DoubleDagger][v][i,j] = h\[DoubleDagger][v][i,-m] h\[DoubleDagger][v][j,-n] \[EmptyDownTriangle][v][m,n];

\!\(\*SubscriptBox[\(\[Chi]\), \(ij\)]\) = \!\(\*SuperscriptBox[SubscriptBox[\(h\), \(i\)], \(m\)]\)\!\(\*SuperscriptBox[SubscriptBox[\(h\), \(j\)], \(n\)]\) \!\(\*SubscriptBox[\(v\), \(\(:\)\(mn\)\)]\) ";


t\[Theta]::usage ="t\[Theta][u_vector]
reletivistic expansion tensor:  t\[Theta][u][m,n]=\!\(\*SubscriptBox[\(\[Theta]\), \(mn\)]\)=\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[SubscriptBox[\(h\), \(m\)], \(p\)]\)\!\(\*SuperscriptBox[SubscriptBox[\(h\), \(n\)], \(q\)]\)(\!\(\*SubscriptBox[\(u\), \(p; q\)]\)+\!\(\*SubscriptBox[\(u\), \(q; p\)]\))
equivalent names: tProjectiveExpansion=t\[Theta]" ;


s\[Theta]::usage ="s\[Theta][u_vector]
reletivistic expansion rate:  s\[Theta][u]=\[Theta]=\!\(\*SuperscriptBox[SubscriptBox[\(\[Theta]\), \(m\)], \(m\)]\)
equivalent names: sProjectiveExpansion=s\[Theta]" ;


t\[Omega]::usage ="t\[Omega][u_vector]
reletivistic vorticity tensor:  t\[Omega][u][m,n]=\!\(\*SubscriptBox[\(\[Omega]\), \(mn\)]\)=\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[SubscriptBox[\(h\), \(m\)], \(p\)]\)\!\(\*SuperscriptBox[SubscriptBox[\(h\), \(n\)], \(q\)]\)(\!\(\*SubscriptBox[\(u\), \(p; q\)]\)-\!\(\*SubscriptBox[\(u\), \(q; p\)]\))
equivalent names: tProjectiveVorticity=t\[Omega];" ;


t\[Sigma]::usage ="t\[Sigma][u_vector]
reletivistic shear tensor:  t\[Sigma][u][m,n]=\!\(\*SubscriptBox[\(\[Sigma]\), \(mn\)]\)=\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[SubscriptBox[\(h\), \(m\)], \(p\)]\)\!\(\*SuperscriptBox[SubscriptBox[\(h\), \(n\)], \(q\)]\)(\!\(\*SubscriptBox[\(u\), \(p; q\)]\)+\!\(\*SubscriptBox[\(u\), \(q; p\)]\)-\!\(\*FractionBox[\(2\), \(3\)]\)\[Theta] \!\(\*SubscriptBox[\(h\), \(pq\)]\))
equivalent names: tProjectiveShear=t\[Sigma];" ;


s\[Sigma]::usage ="s\[Sigma][u_vector]
reletivistic shear scalar:  s\[Sigma][u]=\!\(\*SqrtBox[\(\*FractionBox[\(1\), \(2\)] \*SubscriptBox[\(\[Sigma]\), \(mn\)] \*SuperscriptBox[\(\[Sigma]\), \(mn\)]\)]\)
equivalent names: sProjectiveShear=s\[Sigma];" ;


s\[Omega]::usage ="s\[Omega][u_vector]
reletivistic vorticity scalar:  s\[Omega][u]=\!\(\*SqrtBox[\(\*FractionBox[\(1\), \(2\)] \*SubscriptBox[\(\[Omega]\), \(mn\)] \*SuperscriptBox[\(\[Omega]\), \(mn\)]\)]\)
equivalent names: sProjectiveVorticiyy=s\[Omega];" ;


(* ::Subsubsection::Plain::Closed:: *)
(*geodesic*)


geodesic::usage ="\!\(\*
StyleBox[\"geodesic\",\nFontSlant->\"Italic\"]\)[s]/@\!\(\*
StyleBox[\"all\",\nFontSlant->\"Italic\"]\)";


(* ::Subsubsection::Plain::Closed:: *)
(*assumptions*)


addassumption::usage ="\!\(\*
StyleBox[\"addassumption\",\nFontSlant->\"Italic\"]\)[condition] adds new condition to \!\(\*
StyleBox[\"$Assumptions\",\nFontSlant->\"Italic\"]\)";


positiveEM::usage ="finds conditions for Einstein tensor G to be positively defined";


timelikeEigenvector::usage ="\!\(\*
StyleBox[\"timelikeEigenvector\",\nFontSlant->\"Italic\"]\)[T]  attempts to find timelike eigenvector of the tensor T."


(* ::Subsubsection::Plain::Closed:: *)
(*security*)


erasable::usage ="\!\(\*
StyleBox[\"erasable\",\nFontSlant->\"Italic\"]\)[tensor] adds new tensor to \!\(\*
StyleBox[\"memRegistry\",\nFontSlant->\"Italic\"]\) - see section \!\(\*
StyleBox[\"Memory content\",\nFontSlant->\"Italic\"]\)";


cacheview::usage ="\!\(\*
StyleBox[\"cacheview\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"[\",\nFontSlant->\"Italic\"]\)symbol\!\(\*
StyleBox[\"]\",\nFontSlant->\"Italic\"]\) views memorized values directly assigned to symbol";


associated::usage ="\!\(\*
StyleBox[\"associated\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"[\",\nFontSlant->\"Italic\"]\)symbol\!\(\*
StyleBox[\"]\",\nFontSlant->\"Italic\"]\) views memorized values associated with symbol\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)in whatever way";


retreat::usage ="\!\(\*
StyleBox[\"retreat\",\nFontSlant->\"Italic\"]\)[symbol] clears out values listed by \!\(\*
StyleBox[\"cacheview\",\nFontSlant->\"Italic\"]\)[symbol]

\!\(\*
StyleBox[\"retreat\",\nFontSlant->\"Italic\"]\)[symbol,\!\(\*
StyleBox[\"associated\",\nFontSlant->\"Italic\"]\)] clears out all values shown by \!\(\*
StyleBox[\"associated\",\nFontSlant->\"Italic\"]\)[symbol]";


(* ::Subsubsection::Plain::Closed:: *)
(*auxiliary*)


verify::usage ="\!\(\*
StyleBox[\"verify\",\nFontSlant->\"Italic\"]\)[identity_,indexes__]:=Timing[DeleteDuplicates[Flatten[Simplify[Outer[identity,indexes]]]]]

\!\(\*
StyleBox[\"verify\",\nFontSlant->\"Italic\"]\) attempts to checks whether identity is true or false for specified list of indexes. ";


denominatorlist::usage ="\!\(\*
StyleBox[\"denominatorlist\",\nFontSlant->\"Italic\"]\)=1/Select[DeleteDuplicates[Denominator/@Level[#,Infinity]],!NumericQ[#]&]&;";


(* ::Subsubsection::Plain::Closed:: *)
(*perturbations*)


trunc::usage="\!\(\*
StyleBox[\"trunc\",\nFontSlant->\"Italic\"]\)[] = \!\(\*
StyleBox[\"trunc\",\nFontSlant->\"Italic\"]\)[ord->ord\[DoubleDagger]]"


simp::usage ="\!\(\*
StyleBox[\"simp\",\nFontSlant->\"Italic\"]\)[] = \!\(\*
StyleBox[\"simp\",\nFontSlant->\"Italic\"]\)[mtd->simplification,ord->ord\[DoubleDagger]]
\!\(\*
StyleBox[\"simp\",\nFontSlant->\"Italic\"]\)[mtd->Simplify]
\!\(\*
StyleBox[\"simp\",\nFontSlant->\"Italic\"]\)[ord->0]";


pert\[DoubleDagger]::usage ="The perturbation parameter declared by user in the metric tensor definition.";
ord\[DoubleDagger]::usage ="The perturbation order declared by user in the metric tensor definition.";


(* ::Subsection::Plain::Closed:: *)
(*messages*)


Weyl::lowdim="Space dimension less than 3."


covariantD::argsqnce="The list of indexes shorter than the tensor rank.";
covariantD::indrng="Index out of range.";
covariantD::tnsr="Tensor expected.";
covariantD::zero="0-th order derivative";


LieD::vect="The underlying field `1` is not a vector.";
LieD::argx="The length of the sequence of indexes is different from the rank of the differentiated tensor `1`.";


partialD::argsqnce="The list of indexes shorter than the tensor rank.";
partialD::cntr="No contravariant indices allowed in partial differentiation.";


tensorAux::tnsrank="Tensor rank not determined.";


tensorExt::indrng="Index out of range.";
tensorExt::invval="Invalid tensor valence.";


indeX::indrng="Index out of range.";


span::tnsrank="Tensor rank not determined.";


(* ::Subsubsection::Plain:: *)
(*help*)


(* ccgrg := "{g\[DoubleDagger], \[CapitalGamma], tRicciR, sRicciR, tEinsteinG, tRiemannR, tWeylC, tWeylCdual, covariantD}" *)


(* ::Section::Plain::Closed:: *)
(*Context*)


Begin["`Private`"];


(* ::Section::Plain::Closed:: *)
(*Auxiliary commands*)


(* ::Subsubsection::Plain::Closed:: *)
(*version*)


(* displays the Package version *)
ccver="ver-1.04";
version:=Module[{tam},
	tam=FileNameJoin[{$UserBaseDirectory,"Applications","ccgrg.m"}];
	{ccver,DateString[FileDate[tam]],FileByteCount[tam]}
];


(* ::Subsubsection::Plain::Closed:: *)
(*constants*)


fura=16;


If[!ValueQ[sub],sub={}];
If[!ValueQ[auxi],auxi={}];
If[!ValueQ[restrict],restrict={}];
If[!ValueQ[memRegistry],memRegistry={}];


(* ::Subsubsection::Plain::Closed:: *)
(*general shortcuts*)


MF=MatrixForm;
TF=TraditionalForm;


(* ::Section::Plain::Closed:: *)
(*Functions and arguments*)


(* ::Subsection::Plain::Closed:: *)
(*Functions*)


(* ::Subsubsection::Plain::Closed:: *)
(*simplification*)


Options[simp]={mtd->smp,ord->ord\[DoubleDagger]};
simp[OptionsPattern[]]:=Module[{method,orger,sqq1,sqq2},
	method=OptionValue[mtd];
	order=OptionValue[ord];
	sqq1:=method[#//.sub]&; 
	sqq2:=Normal[method[Normal[#]//.sub//trunc[order]]]&; 
	If[ord\[DoubleDagger]==0,sqq1,(addassumption[0<pert\[DoubleDagger]];sqq2)]
];


denominatorlist=1/Select[DeleteDuplicates[Denominator/@Level[#,Infinity]],!NumericQ[#]&]&;
collectdenominators=Collect[#,denominatorlist[#]]&;


(* ::Subsubsection::Plain::Closed:: *)
(*function name*)


functionName[ff_]:=Last[Select[NestList[Head,ff,fura],#1=!=Symbol&]];


(* ::Subsubsection::Plain::Closed:: *)
(*unspecified functions*)


(**)
unspecified[variables_,expr_]:=Module[{},
Select[DeleteDuplicates[Flatten[Level[expr,{0,\[Infinity]}]]],!NumericQ[#]&&(Head[#]==Symbol||!MemberQ[Attributes[Evaluate[Head[#]]],Protected])&]//Sort//Quiet
];


(* ::Subsubsection::Plain::Closed:: *)
(*argument operations*)


(* ::Text:: *)
(*f[x,y] -> f[{x,y}]*)


klamra:=Head[#]@Level[#,1]&


swap[a_,b_]:=If[a=!=b,Permute[#,Cycles[{{a,b}}]],#]&


(* ::Subsubsection::Plain::Closed:: *)
(*verify*)


(* ::Text:: *)
(*Attempts to checks whether identity is true or false for specified list of indexes*)


verify[identity_,indexes__]:=Timing[DeleteDuplicates[Flatten[Simplify[Outer[identity,indexes]]]]]


(* ::Subsection::Plain::Closed:: *)
(*Perturbation calculus*)


(* ::Subsubsection::Plain::Closed:: *)
(*series expansion*)


trunc[n_:ord\[DoubleDagger]]:=Series[#,{pert\[DoubleDagger],0,n}]&;          


(* ::Subsubsection::Plain::Closed:: *)
(*perturbation order*)


(* ::Text:: *)
(*Detects perturbation variable in the metric form*)


pertVar[g2_]:=Module[{parpert,LL1,LL2},
	parpert[f1_]:=Module[{},
		LL1=Apply[List,#]&@ f1;
		If[!ArrayQ[LL1]&&!NumberQ[LL1]&&Head[LL1]=!=Symbol,First[LL1]]
		];
	LL2=Select[DeleteDuplicates[Flatten[Map[parpert,g2,{2}]]],#=!=Null&];
	If[LL2=!={},#&@@LL2]
];


(* ::Text:: *)
(*Detects perturbation order*)


pertOrder[g2_]:=Module[{elementOrder,order},
	elementOrder[f1_]:=Module[{L1,L2},
		L1=Apply[List,#]&@ f1;
		L2=If[!ArrayQ[L1]&&!NumberQ[L1]&&Head[L1]=!=Symbol,Take[L1,{-2,-1}],0];
		Divide@@L2
		];
	order=Max[Max[Flatten[Map[elementOrder,g2,{2}]]]-1,0]
];


(* ::Section::Plain::Closed:: *)
(*Memory content, and related*)


(* ::Subsection::Plain::Closed:: *)
(*assumptions*)


(* ::Subsubsection::Plain::Closed:: *)
(*addassumption*)


(* ::Text:: *)
(*Adds new assumption to actual calue of $Assumptions*)


(* adds new assumption to $Assumptions *)
addassumption[condition_/;condition=!=False]:=
Module[{restrict},
	restrict=$Assumptions&&condition;
	$Assumptions=If[Head[restrict]===And,And@@DeleteDuplicates[List@@restrict],restrict]
];


(* ::Subsection::Plain::Closed:: *)
(*memory*)


(* ::Subsubsection::Plain::Closed:: *)
(*erasable*)


(* ::Text:: *)
(*erasable adds function name to the tensor registry*)


(* adds function name to the tensor registry *)
erasable[symbol_]:=Module[{},memRegistry=DeleteDuplicates[Append[memRegistry,symbol]];]


(* ::Subsubsection::Plain::Closed:: *)
(*cacheview*)


(* ::Text:: *)
(*Displays the list of memorized components:*)


(* list of memorized components *)
mviewall[function_]:=Module[{begin,y1,y2,y3,y4,y5,y6,k,zobacz,w},
	y1=StringToStream[ToString[function//Definition]];
	zobacz:=Find[y1,ToString[function]];
	Module[{},
	k=1;
	Label[begin];w[k]=zobacz;If[w[k]=!=EndOfFile, (k++;Goto[begin])]];
		y2=Most[Table[w[i],{i,k}]];
	If[Select[y2,Head[#]==String&]=={},y6={},
		y3=Select[y2,Head[#]==String&]//Most//Quiet;
		y4=First[StringSplit[#,"="]]&/@ y3//Quiet;
		y5=ToExpression[StringReplace[#,ToString[function]->"fun\[DoubleDagger]"]]&/@ y4//Quiet;
		y6=If[Length[Select[y5,#=!=$Failed&]]>0,
		{#1[[0]]/.fun\[DoubleDagger]->function,List@@#1}&/@ Select[y5,#=!=$Failed&],{}]
]];


(* ::Text:: *)
(*Search for the memorized functions values:*)


cacheview[function_]:=Module[{hd,mb},
(* search for symbols in the list of memorized components  *)
	hd=Append[Level[#,Infinity,Heads->True],#]&;
	mb[f1_]:=MemberQ[#,f1]& @hd@First[#]&;
	Select[mviewall[functionName[function]],mb[function][#]&&DeleteDuplicates[#[[2]]]=!={0}&]
];


(* ::Text:: *)
(*Search for symbols in the list of memorized values:*)


associated[tensor_]:=Module[{X1,X2,X3},
	X1=Flatten[cacheview/@memRegistry,1]//DeleteDuplicates;
	X2=Position[X1,tensor]//DeleteDuplicates;
	X3=X1[[#]]&@@#&/@X2
];


(* ::Subsubsection::Plain::Closed:: *)
(*retreat*)


clr[u_/;vctr[u]]:=Module[{},Clear[u];
(* clears out the projection tensor and second fundamental form *)
	h\[DoubleDagger][u]=.;\[Chi]\[DoubleDagger][u]=.;
	h\[DoubleDagger][u][0,0]={-1};
	\[Chi]\[DoubleDagger][u][0,0]={-1};
]//Quiet;


(* ::Text:: *)
(*Clears out specitied cahed components:*)


(* clears out indicated memorized components *)
retreat[function_,type_:cacheview]:=Module[{},clr[function];
	DeleteDuplicates[Apply[#1[Sequence@@#2]=.&,#]&/@type[function]];
];


(* ::Section::Plain::Closed:: *)
(*Tensor procedures*)


(* ::Subsubsection::Plain::Closed:: *)
(*metric form*)


(* finds the metric tensor for a quadratic form *)
toMatrix[forma_,variables_]:=
	Module[{dimension,independent,g,gg,cf1,cf2,linsys,Dx}, 
	Dx=Dt/@variables;dimension=Length[variables];
	gg=Array[g,{dimension,dimension}];g[i_,k_]:=g[k,i]/;k<i; 
	independent=DeleteDuplicates[Flatten[gg]];
	cf1=Simplify[forma]-Simplify[Dx.gg.Dx];
	cf2=DeleteCases[Flatten[CoefficientList[cf1,Dx]],0];
	linsys=#==0&/@ cf2;
	gg=gg/.Flatten[Solve[linsys,independent]];
	Simplify[gg]
];


(* ::Subsubsection::Plain::Closed:: *)
(*tensor rank and index ranges*)


(* tensor identifier *)
id\[DoubleDagger]={};


rank[f_]:=Module[{X},
	X=Flatten[Position[Apply[f[##1]&,zerosUp,1],id\[DoubleDagger]]]/.{a_}:>a/.{}:>0//Quiet;
	If[NumberQ[X],X,Infinity]];


zerosUp=ConstantArray[0,#]& /@ Range[fura];
countUp[f_]:=Module[{X},
	X=Flatten[Position[Apply[ValueQ[f[##1]]&,zerosUp,1],True]]/.{a_}:>a/.{}:>0//Quiet;
	If[NumberQ[X],X,Infinity]];


zerosUp2=Prepend[ConstantArray[0,#]& /@ Range[fura],{}];
countUp2[f2_]:=Position[Apply[ValueQ[f2[##1]]&,zerosUp2,1],True]-1//.{a_}:>a


vctr:=TrueQ[countUp[#]==1]&;


setrank[tensorname_,rnk_Integer]:=Module[{},
	Set[Evaluate[tensorname@@ConstantArray[0,rnk]],id\[DoubleDagger]]//Quiet;
	erasable[functionName[tensorname]];
];


(* checks if function is recognized as tensor *)
tensorQ=Module[{zl,wsad},
	zl=countUp[Evaluate[#]];
	wsad=If[zl=!=0,#@@ConstantArray[0,zl]];
	If[zl==0||wsad=={},True,False]
]&;


pion[T1_]:=Module[{sa,perm,active},
	sa[q1_,q2__]:=Boole[ValueQ[T1[q1,q2]]]; 
	SetAttributes[sa,Listable];
	perm=Reverse[Sort[Permutations[ReplacePart[ConstantArray[1,countUp[T1]],1->all2]]]];
	active=Apply[sa,perm,{1}] ;
	Cases[#,Except[0]]&/@(Times[all2,#]&/@active)
	];


(* ::Subsubsection::Plain::Closed:: *)
(*tensor indexes*)


var[nazwa_,rz\:0105dtensora_]:=ToExpression[ToString[nazwa]<>ToString[#]]&/@ Range[rz\:0105dtensora];
varL[nazwa_,rz\:0105dtensora_]:=ToExpression[ToString[nazwa]<>ToString[#]<>"_"]&/@ Range[rz\:0105dtensora];


(* strong index condition *)
idxx[i_Integer]:=0<i<=dim;
idXX[i_Integer]:=0<Abs[i]<=dim;


(* strong index condition *)
idx:=If[And@@(#1\[Element]Integers&&0<#1<=dim&)/@Flatten[{##1}],True,Message[indeX::indrng]]&;
idX:=If[And@@(#1\[Element]Integers&&0<Abs[#1]<=dim&)/@Flatten[{##1}],True,Message[indeX::indrng]]&;


(* weak index condition*)
indeX:=If[And@@(#1\[Element]Integers&&0<=Abs[#1]<=dim&)/@Flatten[{##1}],True,Message[indeX::indrng]]&;


(* ::Subsubsection::Plain::Closed:: *)
(*metric tensor extension*)


(* Extended metrics 8x8 *)
chess:=Module[{m1,m2,dim},dim=Length[#];
m1={Join[#,Reverse[IdentityMatrix[dim]]],
    Join[Reverse[IdentityMatrix[dim]],Reverse[Reverse/@Inverse[#]]]};
m2=Flatten/@Thread[m1]//simp[]//Simplify]&;


(* ::Subsubsection::Plain::Closed:: *)
(*tensor definitions*)


Options[tensorExt]={valence->1,metrics->metricExt};
tensorExt[Tensor6_,OptionsPattern[]][index6__]/;Times@index6==0:=id\[DoubleDagger];
tensorExt[Tensor6_,OptionsPattern[]][index6__]/;If[indeX[{index6}],True,Message[tensorExt::indrng];False]:=
	Module[{mtensor,numarg6,tval,tExt},
	Clear[mtensor,numarg6,tval,tExt];
	numarg6=Length[{index6}];
	tval=OptionValue[valence]/.{d->1,u->-1}/.{"d"->1,"u"->-1};
	mtensor=OptionValue[metrics];
	tExt[functionname_/;functionname=!=\[CapitalGamma]&&!NumberQ[functionname],
		index_List/;DeleteDuplicates[Flatten[index]]=={0}||Complement[index,all2]=={} ]:=
		If[DeleteDuplicates[Flatten[List[index]]]=={0},id\[DoubleDagger],
			Module[{p,f,metric,numarg,mobileIdx,fixedIdx,whereMobile,whereFixed,dummyIndex,dummyIdxBlank,skipIdx,term,skladnik,zakres,fkt},
				Clear[p,f,metric,numarg,mobileIdx,fixedIdx,whereMobile,whereFixed,dummyIndex,dummyIdxBlank,skipIdx,term,skladnik,zakres,fkt];
				numarg=Length[index];
				whereMobile=Position[Sign/@index tval,-1]//Flatten;
				whereFixed=Position[Sign/@index tval,1]//Flatten;
				mobileIdx=Part[index ,whereMobile];
				fixedIdx=Part[index ,whereFixed];
				dummyIndex=var[p,numarg][[whereMobile]];
				dummyIdxBlank=varL[p,numarg][[whereMobile]];
				skipIdx=var[p,numarg][[whereFixed]];
				If[Length[mobileIdx]==0,functionname@@(index tval),
					(fkt=If[tval===1,1,tval[[whereMobile]]];
					term=Times@@(metric@@#&/@Thread[{mobileIdx,-dummyIndex}])f@@(tval var[p,numarg]);(**)
					Clear[skladnik];Evaluate[skladnik@@dummyIdxBlank]:=Evaluate[term/.Thread[skipIdx->fixedIdx]];
					zakres=fkt ConstantArray[all,Length[dummyIndex]];Plus@@Flatten[Outer[skladnik,Sequence@@zakres]]/.metric->mtensor/.f->functionname//simp[])
				](* If *)
			](* Module *)
		];(* If *)
	If[(Head[tval]===List)&&(Abs[Times@@Flatten[tval]]===1)&&(numarg6==Length[tval])|| tval===1,tExt[Tensor6,Flatten[{index6}]],(Message[tensorExt::invval])]
];


span[T_,Tcov_,metrictensor_: metricExt]:=Module[{rozmiar},
	T[q___/;Length[{q}]==countUp[Tcov]]:=T[q]=tensorExt[Tcov,metrics->metrictensor][q];(* erasable[T] *)
];


(* ::Subsubsection::Plain::Closed:: *)
(*auxiliary tensor*)


tensorAux[T_,TAux_]:=
Module[{X,funk,INdex,zakresy,array1,array2,rk,ins,values,numericComponents,symbolicComponents,original,allcomponents,multiple,subst,c},
	Clear[TAux]; erasable[TAux];
	If[!NumberQ[countUp[T]],(Message[tensorAux::tnsrank];Abort[];)];
	INdex=pion[T];//Quiet;
    array1=tabular[T][Sequence@@INdex];
	array2=tabular[TAux][Sequence@@INdex];
	rk=Length[INdex];
	ins[liczba_]:=
		Module[{L1,L2},
			L1=Position[array1,liczba,rk];
			L2=array2[[Sequence@@#]]&/@L1;
			Set[#,liczba]&/@L2//Quiet;
	];
	allcomponents=DeleteDuplicates[Flatten[array1]];
	numericComponents=Select[allcomponents,NumericQ];
	ins/@numericComponents;
	symbolicComponents=Complement[allcomponents,numericComponents];
	multiple=Select[Position[array1,#,rk]&/@ symbolicComponents,Length[#]>1&];
	values=array2[[Sequence@@#1]]&;
	original=Map[values,multiple,{2}]//Sort;
	subst=ConstantArray[Last[#],Length[#]]&/@original;
	Set@@#&/@ Thread[{Flatten[original],Flatten[subst]}];
];


(* ::Subsubsection::Plain::Closed:: *)
(*tabular forms*)


 tabular[tensor_]:=simp[][Outer[tensor,##]]&; (*simp[]*)


(* ::Subsubsection::Plain::Closed:: *)
(*antisymetric tensor*)


metricsign=Sign[metricDet]//Simplify;
sqrtdet=Sqrt[metricsign metricDet]//Simplify;


antisymmetric:=Module[{etacov},
	etacov[a_,b__]:=metricsign sqrtdet Signature[{a,b}];   (* Soko\[LSlash]owski  3.25 *)
	etaExt[i_,j__]/;indeX[i,j]:=tensorExt[etacov][i,j];\[Eta]\[DoubleDagger]=etaExt;erasable[etaExt];
];


(* ::Subsubsection::Plain::Closed:: *)
(*vector tools*)


vectorsquared[vector_/;vctr[vector]]:=Sum[metricExt[-p,-q]vector[p]vector[q],{p,dim},{q,dim}];


unitvector[vector_/;vctr[vector]][m_Integer/;Abs[m]<=dim]:=
Module[{denominator,original,final},
	denominator=Sqrt[Sign[vectorsquared[vector]]vectorsquared[vector]]//Simplify; (* Simplify *) 
	If[denominator=!=0,
		original[k_/;0<k<=dim]:=vector[k]/denominator;
		original[0]=id\[DoubleDagger];
	];
	span[final,original];
	final[m]
];


(* ::Section::Plain::Closed:: *)
(*Differential geometry*)


(* ::Subsection::Plain::Closed:: *)
(*Curvature*)


(* ::Subsubsection::Plain::Closed:: *)
(*Riemann curvature*)


(* ::Text:: *)
(*Connection, Ricci, Riemann, Weyl:*)


curvature:=Module[{p1,p2},
	(* connection *)
	\[CapitalGamma][k_?idX,i_?idx,j_?idx]:= \[CapitalGamma][k,i,j]=simp[][1/2 Sum[g\[DoubleDagger][k,-s](-D[g\[DoubleDagger][i,j],x\[DoubleDagger][-s]]+D[g\[DoubleDagger][i,s],x\[DoubleDagger][-j]]+D[g\[DoubleDagger][j,s],x\[DoubleDagger][-i]]),{s,1,dim}]];
	\[CapitalGamma][k_?idX,j_?idx,i_?idx]/;i<j:=\[CapitalGamma][k,i,j];\[CapitalGamma][0,0,0]={-1};erasable[\[CapitalGamma]];
	tensorAux[\[CapitalGamma],\[CapitalGamma]Aux];\[CapitalGamma]Aux[0,0,0]={-1};
	(* Ricci tensor *)
	erasable[Ric]; 
	Ric[0,0]={-1};
	Ric[j_,i_]/;i<j:=Ric[i,j];
	Ric[i_,j_]:=Ric[i,j]=simp[][Sum[D[\[CapitalGamma][-k,i,j],x\[DoubleDagger][-k]],{k,1,dim}]-D[Sum[\[CapitalGamma][-q,i,q],{q,1,dim}],x\[DoubleDagger][-j]]+Sum[\[CapitalGamma][-p,i,j]Sum[\[CapitalGamma][-q,p,q],{q,1,dim}],{p,1,dim}]-Sum[\[CapitalGamma][-p,i,q]\[CapitalGamma][-q,j,p],{p,1,dim},{q,1,dim}]];
	erasable[tRicciR];
	tRicciR[i_,j_]/;Times@@{i,j}==0=id\[DoubleDagger];
	tRicciR[i_,j_]/;indeX[i,j]:=tRicciR[i,j]=simp[][tensorExt[Ric][i,j]];
	sRicciR[]:=simp[][Sum[tRicciR[q,-q],{q,1,dim}]];
	(* Einstein tensor *)
	erasable[tEinsteinG];
	tEinsteinG[i_,j_]/;Times@@{i,j}==0=id\[DoubleDagger];
	tEinsteinG[i_,j_]/;indeX[i,j]:=tEinsteinG[i,j]=simp[][tRicciR[i,j]-1/2 g\[DoubleDagger][i,j]sRicciR[]];G\[DoubleDagger]=tEinsteinG;
	sEinsteinG[]:=simp[][Sum[tEinsteinG[q,-q],{q,1,dim}]];
	(* Riemann tensor *)
	erasable[Riem];
	Riem[0,0,0,0]={-1};
	Riem[i_,j_,m_,n_]/;i==j:=Riem[i,j,m,n]=0;
	Riem[i_,j_,m_,n_]/;m==n:=Riem[i,j,m,n]=0;
	Riem[j_,i_,m_,n_]/;i<j:=-Riem[i,j,m,n];
	Riem[i_,j_,n_,m_]/;m<n:=-Riem[i,j,m,n];
	Riem[m_,n_,i_,j_]/;m+n<i+j:=Riem[i,j,m,n];
	Riem[i_,j_,m_,n_]:=Riem[i,j,m,n]=simp[][D[\[CapitalGamma][j,i,m],x\[DoubleDagger][-n]]-D[\[CapitalGamma][j,i,n],x\[DoubleDagger][-m]]+Sum[\[CapitalGamma][-p,i,n]\[CapitalGamma][p,j,m],{p,1,dim}]-Sum[\[CapitalGamma][-p,i,m]\[CapitalGamma][p,j,n],{p,1,dim}]];
	erasable[tRiemannR];
	tRiemannR[i_,j_,m_,n_]/;Times@@{i,j,m,n}==0=id\[DoubleDagger];
	tRiemannR[i_,j_,m_,n_]/;indeX[i,j,m,n]:=tRiemannR[i,j,m,n]=simp[][tensorExt[Riem][i,j,m,n]];
	(* Weyl tensor *)
	erasable[tWeylC];
	tWeylC[i_,j_,m_,n_]/;Times@@{i,j,m,n}==0=id\[DoubleDagger];
	tWeylC[i_,k_,l_,m_]/;If[indeX[i,k,l,m]&&dim>2,True,Message[Weyl::lowdim];False]:=tWeylC[i,k,l,m]=simp[][tRiemannR[i,k,l,m]+1/(dim-2) (tRicciR[i,m]g\[DoubleDagger][k,l]+tRicciR[k,l]g\[DoubleDagger][i,m]-tRicciR[k,m]g\[DoubleDagger][i,l]-tRicciR[i,l]g\[DoubleDagger][k,m])+1/((dim-1)(dim-2)) sRicciR[](g\[DoubleDagger][i,l] g\[DoubleDagger][k,m]-g\[DoubleDagger][i,m] g\[DoubleDagger][k,l])];
	erasable[tWeylCdual];
	tWeylCdual[i_,j_,m_,n_]/;Times@@{i,j,m,n}==0=id\[DoubleDagger];
	tWeylCdual[i_,j_,m_,n_]/;If[indeX[i,j,m,n]&&dim>2,True,Message[Weyl::lowdim];False]:=tWeylCdual[i,j,m,n]=simp[][1/2 Sum[etaExt[i,j,p1,p2]tWeylC[-p1,-p2,m,n],{p1,dim},{p2,dim}]];
];


(* ::Subsubsection::Plain::Closed:: *)
(*Carminati - McLenaghan invariants*)


(* Carminati-McLenaghan invariants *)
CMinvariants:=Module[{},
	tPlebanskiS[a_Integer/;Abs[a]<=dim,b_Integer/;Abs[b]<=dim]:=tPlebanskiS[a,b]=simp[][(tRicciR[a,b]-1/dim sRicciR[] g\[DoubleDagger][a,b])];erasable[tPlebanskiS];
	S\[DoubleDagger]=tPlebanskiS;
	CMinvR1:=simp[][Sum[1/4 S\[DoubleDagger][-q1,q2]S\[DoubleDagger][-q2,q1],{q1,dim},{q2,dim}]];
	CMinvR2:=simp[][Sum[-1/8 S\[DoubleDagger][-q1,q2]S\[DoubleDagger][-q2,q3]S\[DoubleDagger][-q3,q1],{q1,dim},{q2,dim},{q3,dim}]];
	CMinvR3:=simp[][Sum[1/16 S\[DoubleDagger][-q1,q2]S\[DoubleDagger][-q2,q3]S\[DoubleDagger][-q3,q4]S\[DoubleDagger][-q4,q1],{q1,dim},{q2,dim},{q3,dim},{q4,dim}]];
	CMinvM3/;If[dim>2,True,Message[Weyl::lowdim];False]:=Module[{F3},
		F3[f_]:=simp[][Sum[1/16 S\[DoubleDagger][-q2,-q3]S\[DoubleDagger][q5,q6]f[q1,q2,q3,q4]f[-q1,-q5,-q6,-q4],{q1,dim},{q2,dim},{q3,dim},{q4,dim},{q5,dim},{q6,dim}]];
		F3[tWeylC]+F3[tWeylCdual]];
	CMinvM4/;If[dim>2,True,Message[Weyl::lowdim];False]:=Module[{F4},
		F4[f_]:=simp[][Sum[-1/32 S\[DoubleDagger][-q1,-q7] S\[DoubleDagger][-q5,-q6]S\[DoubleDagger][-q3,q4]f[q1,q3,-q4,-q2]f[q2,q5,q6,q7],{q1,dim},{q2,dim},{q3,dim},{q4,dim},{q5,dim},{q6,dim},{q7,all}]];
		F4[tWeylC]+F4[tWeylCdual]];
	CMinvW1/;If[dim>2,True,Message[Weyl::lowdim];False]:=simp[][Sum[1/8 (tWeylC[a1,a2,a3,a4]+I tWeylCdual[a1,a2,a3,a4])tWeylC[-a1,-a2,-a3,-a4],{a1,dim},{a2,dim},{a3,dim},{a4,dim}]];
	CMinvW2/;If[dim>2,True,Message[Weyl::lowdim];False]:=simp[][Sum[-1/16 (tWeylC[q1,q2,-q3,-q4]+I tWeylCdual[q1,q2,-q3,-q4])tWeylC[q3,q4,-q5,-q6]tWeylC[q5,q6,-q1,-q2],{q1,all},{q2,all},{q3,all},{q4,all},{q5,all},{q6,all}]];
	CMinvM1/;If[dim>2,True,Message[Weyl::lowdim];False]:=simp[][Sum[1/8 S\[DoubleDagger][-q1,-q4]S\[DoubleDagger][-q2,-q3](tWeylC[q1,q2,q3,q4]+I tWeylCdual[q1,q2,q3,q4]),{q1,all},{q2,all},{q3,all},{q4,all}]];
	CMinvM2/;If[dim>2,True,Message[Weyl::lowdim];False]:=simp[][Sum[1/16 S\[DoubleDagger][-q2,-q3]S\[DoubleDagger][q5,q6](tWeylC[q1,q2,q3,q4]tWeylC[-q1,-q5,-q6,-q4]-tWeylCdual[q1,q2,q3,q4]tWeylCdual[-q1,-q5,-q6,-q4])+I/8 S\[DoubleDagger][-q2,-q3]S\[DoubleDagger][q5,q6]tWeylCdual[q1,q2,q3,q4]tWeylC[-q1,-q5,-q6,-q4],{q1,all},{q2,all},{q3,all},{q4,all},{q5,all},{q6,all}]];
	CMinvM5/;If[dim>2,True,Message[Weyl::lowdim];False]:=simp[][Sum[1/32 S\[DoubleDagger][-q2,-q3]S\[DoubleDagger][-q5,-q6](tWeylC[-q1,-q7,-q8,-q4]+I tWeylCdual[-q1,-q7,-q8,-q4])(tWeylC[q1,q2,q3,q4]tWeylC[q7,q5,q6,q8]+tWeylCdual[q1,q2,q3,q4]tWeylCdual[q7,q5,q6,q8]),{q1,all},{q2,all},{q3,all},{q4,all},{q5,all},{q6,all},{q7,all},{q8,all}]];
];


(* ::Subsection::Plain::Closed:: *)
(*Derivatives*)


(* ::Subsubsection::Plain::Closed:: *)
(*covariant derivative: covariantD*)


(* ::Text:: *)
(*n^th - order covariant derivative*)


(* n-order covariant derivative *)
covariantderivative:=Module[{covariantDTensor,P,sg,sd3,P1,lt},

covariantDTensor[Tensor_,Index_]:=
Module[{coord,lowerAll,firstOrder,aff,affAll,diffAll,pDim,p,q},
	P[0,_]=0;sg=Sign/@Index;
	rankT=rank[Tensor];
	orderD=Length[Index]-rankT;
	lowerAll[tensorT_,Lq_]:=lowerAll[tensorT,Lq]=
	Module[{tt},
		 firstOrder[tensorT1_,Lq1_]:=firstOrder[tensorT1,Lq1]=
		 Module[{qDim},
			qDim=var[q,Length[Lq1]];
			aff[tensorT2_,wsk_]:=(Sum[\[CapitalGamma]Aux[-s,wsk,Last[qDim]]tensorT2 @@(Most[qDim]/.wsk->s),{s,dim}]);
			affAll[tensorT2_,Lq2_]:= Total@(aff[tensorT2,#]&/@Most[qDim]);
			diffAll[tensorT2_,Lq2_]:=P[tensorT2@@Most[qDim],coord[-Last[qDim]]];
			diffAll[tensorT1,Lq1]-If[rankT>0,affAll[tensorT1,Lq1],0]//.Thread[qDim->Lq1]
		]; 
		tt[0]=tensorT; 
		pDim=var[p,Length[Lq]];
		Last[Table[tt[k][Sequence@@varL[p,rankT+k]]=firstOrder[tt[k-1],var[p,rankT+k]],{k,orderD}]]//.Thread[pDim->Lq]
	];
	Clear[sd3];Evaluate[sd3@@varL[k,rankT+orderD]]:=Evaluate[lowerAll[Tensor,var[k,rankT+orderD]]];
	tensorExt[sd3][Index]//.{coord->x\[DoubleDagger],gAux->g\[DoubleDagger],\[CapitalGamma]Aux->\[CapitalGamma]}/.P->D//simp[]
];

covariantDAll[tensor_,index_List/;Complement[index,all3]=={}]:=
Module[{P1lt,rk},
	P1[k_]:=tensorExt[partialD[tensor]][k];
	rk=rank[tensor];lt=Length[index];
	Which[
		lt<rk, (Message[covariantD::argsqnce];Abort[];),
		lt==rk,(Message[covariantD::zero];tensor@@index),
		lt>rk&&rk>0, covariantDTensor[tensor,index],
		rk==0&&lt==1,tensorExt[Sum[metricExt[#,-k] P1[k],{k,all}]&][index[[1]]],
		rk==0&&lt>1, covariantDTensor[P1,index]
	]
];

covariantD[Tensor6_][index6__]/;tensorQ[Tensor6]&&Times@Flatten[{index6}]==0:=id\[DoubleDagger];
covariantD[Tensor6_][index6__]/;
	If[tensorQ[Tensor6],
		If[TrueQ[Complement[Select[Flatten[List[index6]],NumericQ],all3]=={}],True,Message[covariantD::indrng];False],
		Message[covariantD::tnsr];False
	]:=covariantD[Tensor6][index6]= covariantDAll[Tensor6,Flatten[List[index6]]
];

\[EmptyDownTriangle]=covariantD;erasable[covariantD];
];


(* ::Subsubsection::Plain::Closed:: *)
(*partial derivative: partialD*)


(* ::Text:: *)
(*Partial derivative with syntax identical to covariant derivative*)


(* partial derivative with syntax identical to covariant derivative *)
partialDS[f_,pdindex_List/;Complement[Flatten[pdindex],all2]=={}]:=
Module[{pdindex1,pdindex2,fpdx,pdrankf,pdinddiv,fAux},
	Clear[pdrankf,pdindex1,pdindex2,pdinddiv];
	fAux=Evaluate[f];pdrankf=rank[fAux];
	fpdx=Flatten[pdindex];
	If[Length[fpdx]<pdrankf,(Message[partialD::argsqnce];Abort[];),
	(pdinddiv={Take[fpdx,pdrankf],Drop[fpdx,pdrankf]};
	If[Complement[Flatten[pdinddiv[[2]]],all]=!={},(Message[partialD::cntr];(* Abort[]; *)),
	D[If[countUp2[fAux]==={},fAux,fAux@@pdinddiv[[1]]],Sequence@@coordinateExt/@Flatten[-pdinddiv[[2]]]]]
	)]
];
partialD[f6_][pindex6__]:=partialDS[f6,Flatten[List[pindex6]]];


(* ::Subsubsection::Plain::Closed:: *)
(*Lie derivative: LieD*)


(* Lie derivative *)
LieD[U_][T_][i___Integer]/;
If[rank[U]==1,True,Message[LieD::vect,U]]&&If[Length[{i}]==rank[T],True,Message[LieD::argx,T]]:=
	Module[{},
		LieD[U][T][Sequence@@ConstantArray[0,rank[T]]]=id\[DoubleDagger];
		Sum[U[-j]D[T[i],x\[DoubleDagger][-j]],{j,dim}]+Total[MapIndexed[Sign[#1]Sum[T@@ReplacePart[{i},First[#2]->Sign[#1]j](D[U[#1],x\[DoubleDagger][-#2]]&@@Sort[{#1,-Sign[#1]j},#1<#2/.j->1&]),{j,dim}] (* Sum *)&,{i}](* MapIndexed *)] (* Total *)
];(* Module *) 
(* Krzysztof G\[LSlash]\[OAcute]d *)


(* ::Subsection::Plain::Closed:: *)
(*Hypersurfaces*)


(* ::Subsubsection::Plain::Closed:: *)
(*hypersurfaces *)


hypersurfaces:=Module[{hcov,\[Chi]cov},
	hcov[v_][i_,j_]:=g\[DoubleDagger][i,j]-v[i]v[j]/vectorsquared[v]//simp[];
	hcov[v_][j_,i_]/;i<j:=hcov[v][j,i]=hcov[v][i,j];
	h\[DoubleDagger][v_][i_,j_]/;indeX[i,j]:=h\[DoubleDagger][v][i,j]=tensorExt[hcov[v]][i,j]//simp[];
	erasable/@{hcov,h\[DoubleDagger]};

	\[Chi]cov[v_][m_, n_]:=Sum[covariantD[v][i,j]h\[DoubleDagger][v][-i,m]h\[DoubleDagger][v][-j,n],{i,all},{j,all}]//simp[];
	\[Chi]\[DoubleDagger][v_][i_,j_]/;indeX[i,j]:=\[Chi]\[DoubleDagger][v][i,j]=tensorExt[\[Chi]cov[v]][i,j]//simp[];
	erasable/@{\[Chi]cov,\[Chi]\[DoubleDagger]};

	tProjectiveExpansion=t\[Theta];erasable[t\[Theta]];
	sProjectiveExpansion=s\[Theta];erasable[s\[Theta]];
	t\[Theta]cov[u_][i_,j_]:=t\[Theta]cov[u][i,j]=1/2 Sum[h\[DoubleDagger][u][i,-q]h\[DoubleDagger][u][j,-p](\[Chi]\[DoubleDagger][u][p,q]+\[Chi]\[DoubleDagger][u][q,p]),{p,all},{q,all}]//simp[];
	t\[Theta][u_][i_,j_]/;indeX[i,j]:=t\[Theta][u][i,j]=tensorExt[t\[Theta]cov[u]][i,j]//simp[];
	erasable[t\[Theta]cov];

	tProjectiveVorticity=t\[Omega];erasable[t\[Omega]];
	sProjectivelVorticity=s\[Omega];erasable[s\[Omega]];
	\[Omega]\[Omega][u_]:=Outer[t\[Omega][u],all,-all];
	t\[Omega]cov[u_][i_,j_]:=t\[Omega]cov[u][i,j]=1/2 Sum[h\[DoubleDagger][u][i,-q]h\[DoubleDagger][u][j,-p](\[Chi]\[DoubleDagger][u][p,q]-\[Chi]\[DoubleDagger][u][q,p]),{p,all},{q,all}]//simp[mtd->Simplify,ord->n][ord\[DoubleDagger]+1];
	t\[Omega][u_][i_,j_]/;indeX[i,j]:=t\[Omega][u][i,j]=tensorExt[t\[Omega]cov[u]][i,j]//simp[mtd->Simplify,ord->n][ord\[DoubleDagger]+1];
	s\[Omega][u_]:=Sqrt[Tr[\[Omega]\[Omega][u].\[Omega]\[Omega][u]]/2]//simp[mtd->Simplify,ord->n][ord\[DoubleDagger]+1];
	erasable[t\[Omega]cov];

	tProjectiveShear=t\[Sigma];erasable[t\[Sigma]];
	sProjectiveShear=s\[Sigma];erasable[s\[Sigma]];
	\[Sigma]\[Sigma][u_]:=Outer[t\[Sigma][u],all,-all];
	s\[Theta][u_]:=Outer[t\[Theta][u],all,-all]//Tr//simp[]; erasable[s\[Theta]];
	t\[Sigma][u_][i_,j_]/;indeX[i,j]:=t\[Theta][u][i,j]-1/3 s\[Theta][u]h\[DoubleDagger][u][i,j]//simp[mtd->Simplify,ord->n][ord\[DoubleDagger]+1];
	s\[Sigma][u_]:=Sqrt[Tr[\[Sigma]\[Sigma][u].\[Sigma]\[Sigma][u]]/2]//simp[mtd->Simplify,ord->n][ord\[DoubleDagger]+1];
];


(* ::Subsection::Plain::Closed:: *)
(*Geodesic*)


(* ::Subsubsection::Plain::Closed:: *)
(*geodesic*)


(* geodesic equation: explicit form *)
geo[s_][i_]:=D[x\[DoubleDagger][-i],{s,2}]+Sum[\[CapitalGamma][-i,j,k]D[x\[DoubleDagger][-j],s]D[x\[DoubleDagger][-k],s],{j,dim},{k,dim}]==0; 
geodesic[s_][i_]:=Equal@@@Flatten[Simplify[Solve[geo[s][i],D[x\[DoubleDagger][-i],{s,2}]]]];


(* ::Section::Plain::Closed:: *)
(*General relativity*)


(* ::Subsection::Plain::Closed:: *)
(*GRG-sample metric tensors*)


(* ::Subsubsection::Plain::Closed:: *)
(*Schwarzschild*)


Schwarzschild:=
Module[{},
	Clear[Global`\[ScriptT], Global`\[ScriptR], Global`\[CurlyTheta], Global`\[CurlyPhi], Global`M, Global`\[ScriptR]];
	Global`x={Global`\[ScriptT], Global`\[ScriptR], Global`\[CurlyTheta], Global`\[CurlyPhi]}; 
	Global`g=DiagonalMatrix[{-(1-(2Global`M)/Global`\[ScriptR]),(1-(2Global`M)/Global`\[ScriptR])^-1,Global`\[ScriptR]^2,Global`\[ScriptR]^2 Sin[Global`\[CurlyTheta]]^2}]; 
	$Assumptions=Global`\[ScriptR]>0;
	open[Global`x,Global`g]
	(* Global`g//MatrixForm *)
];


(* ::Subsubsection::Plain::Closed:: *)
(*Friedman*)


Friedman:=
Module[{},
	Clear[Global`\[Eta], Global`\[ScriptX], Global`\[ScriptY], Global`\[ScriptZ], Global`a, K, Global`x];
	Global`x={Global`\[Eta], Global`\[ScriptX], Global`\[ScriptY],Global`\[ScriptZ]};
	Global`g=Global`a[Global`\[Eta]]^2 DiagonalMatrix[{-1,(1+K/4 (Global`x.Global`x-Global`\[Eta]^2))^-2,(1+K/4 (Global`x.Global`x-Global`\[Eta]^2))^-2,(1+K/4 (Global`x.Global`x-Global`\[Eta]^2))^-2}]; (* Roberson-Walker *)
	$Assumptions=Global`a[Global`\[Eta]]>0&&1+K/4 (Global`x.Global`x-Global`\[Eta]^2)>0; 
	open[Global`x,Global`g]
	(* Global`g//MatrixForm *)
];


(* ::Subsubsection::Plain::Closed:: *)
(*Lemaitre-Tolman-Bondi*)


LTB:=
Module[{},
	Clear[Global`\[ScriptT], Global`\[ScriptR], Global`\[CurlyTheta], Global`\[CurlyPhi]];
	Global`x={Global`\[ScriptT], Global`\[ScriptR], Global`\[CurlyTheta], Global`\[CurlyPhi]};
	Global`g=DiagonalMatrix[{-1,((D[Global`a[Global`\[ScriptT], Global`\[ScriptR]], Global`\[ScriptR]])^2/(1+2 Global`\[ScriptF][Global`\[ScriptR]])),Global`a[Global`\[ScriptT],Global`\[ScriptR]]^2,Global`a[Global`\[ScriptT],Global`\[ScriptR]]^2 Sin[Global`\[CurlyTheta]]^2}]; 
	$Assumptions=Global`r>0&&1+2 Global`\[ScriptF][Global`\[ScriptR]]>0&&0<Global`\[CurlyPhi]<2\[Pi]&&0<Global`\[CurlyTheta]<\[Pi];
	open[Global`x,Global`g,simplification->Simplify]
	(* Global`g//MatrixForm *)
];


(* ::Subsection::Plain::Closed:: *)
(*GRG dedicated*)


(* ::Subsubsection::Plain::Closed:: *)
(*energy conditions*)


(* ::Text:: *)
(*Positively defined Einstein tensor G:*)


(* positively defined Einstein tensor G *)
positiveEM:=
	Module[{M1,M2},
	M1=Outer[tEinsteinG,all,all]//simp[];
	M2=M1/.pert\[DoubleDagger]->0//simp[];
	valuesG=Eigenvalues[M2];
	Simplify[And@@Simplify[(0<=#1&)/@ valuesG]]
];


(* ::Subsubsection::Plain::Closed:: *)
(*timelike eigenvector *)


(* ::Text:: *)
(*Attempts to find timelike eigenvector of the Einstein tensor G*)


(* timelike eigenvector of Einstein tensor G *)
timelikeEigenvector[Tensor_,Vector_,podstawieniaWW_:{}]:=
Module[{M,gg,eigvG,prenorm,Prenorm,Prenorm0,sign,timelike,W,single},
	gg=Outer[g\[DoubleDagger],all,all];
	M=Outer[Tensor,all,-all]/.podstawieniaWW//simp[mtd->Simplify];
	eigvG=Eigenvectors[M]//simp[mtd->Simplify];
	prenorm[i_Integer]:=eigvG[[i]].eigvG[[i]]/eigvG[[i]].gg.eigvG[[i]]//simp[mtd->Simplify];
	Prenorm=Array[prenorm,dim]/.podstawieniaWW;
	Prenorm0=Limit[Prenorm,pert\[DoubleDagger]->0];
	sign=2Boole/@ simp[mtd->Simplify][Positive/@Prenorm0]-1;
	single=Select[Tally[sign],Last[#]==1&]//Flatten//First;
	timelike=#&@@Flatten[Position[sign,single]];
	W[n_,i_]:=PowerExpand[simp[mtd->Simplify][eigvG[[n]]/Sqrt[sign[[n]] eigvG[[n]].gg.eigvG[[n]]]]][[-i]]//simp[mtd->Simplify];
	Vector[i_]:=Vector[i]=\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(dim\)]\(g\[DoubleDagger][i, j] W[timelike, \(-j\)]\)\)//simp[mtd->Simplify]; 
	erasable[Vector];Vector[0]=id\[DoubleDagger];
];


(* ::Section::Plain::Closed:: *)
(*Open*)


(* ::Subsection::Plain::Closed:: *)
(*open >>*)


Options[open]={simplification->Together};


open[x1_,g1_,OptionsPattern[]]:=Module[{gg},

tasama=(g1===g2&&OptionValue[simplification]===smp);
If[!ValueQ[memRegistry],memRegistry={}];
If[!tasama,(Clear[dim,all,all2]; retreat/@memRegistry;)];

(* conditions *)
If[!ValueQ[sub],sub={}];
If[!ValueQ[auxi],auxi={}];
If[!ValueQ[restrict],restrict={}];
RealFunctionsAll=unspecified[x1,g1]\[Element] Reals;
restrict=x1\[Element]Reals&&RealFunctionsAll; 
addassumption[restrict];

(* dimension *)
dim=Length[x1];all=Range[dim];
all2=RotateLeft[Complement[Range[-dim,dim],{0}],dim];
all3=Range[-dim,dim];

(* expansion *)
pert\[DoubleDagger] = pertVar[g1];
ord\[DoubleDagger]=Max[0,pertOrder[g1]];

(* simplification *)
smp=OptionValue[simplification];

(* coordinates *)
Table[x\[DoubleDagger][-i]=x1[[i]],{i,dim}];coordinateExt=x\[DoubleDagger];SetAttributes[x\[DoubleDagger],Listable];erasable[x\[DoubleDagger]];

(* extended metric tensor 8x8 *)
Table[g\[DoubleDagger][i,j]=chess[simp[mtd->Simplify][g1]][[i,j]],{i,all2},{j,all2}]; erasable[g\[DoubleDagger]];g\[DoubleDagger][0,0]=id\[DoubleDagger];
metricExt=g\[DoubleDagger];SetAttributes[g\[DoubleDagger],Listable];
metricDet=Det[Outer[g\[DoubleDagger],all,all]]//Simplify;
tensorAux[g\[DoubleDagger],gAux];gAux[0,0]={-1};


(* antisymmetric symbol *)antisymmetric;
(* Riemann curvature *)curvature;
(* nth-odrer covariant derivative *)covariantderivative;
(* hypersurfaces *)hypersurfaces;
(* Carminati-McLenaghan invariants *)CMinvariants;
	
g2=g1;
Column[{If[tasama, "continuation:  ","new metric:  "],"  ", MF[g1]," ","simplification method: " OptionValue[simplification]}]

];


(* ::Section::Plain::Closed:: *)
(*End*)


End[];
Print[Style["ccgrg - Copernicus Center General Relativity Package for Mathematica 8/9", Italic, FontFamily->"Helvetica"]]
EndPackage[];