xPand.m

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



xAct`xPand`$Version={"0.4.3",{2018,04,28}};


xAct`xPand`$xTensorVersionExpected={"1.1.3",{2018,2,28}};
xAct`xPand`$xPertVersionExpected={"1.0.6",{2018,2,28}};


(* xPand: Cosmological perturbations about homogeneous space-times *)

(* Copyright (C) 2012-2018 Cyril Pitrou, Xavier Roy and Obinna Umeh *)

(* This program is free software; you can redistribute it and/or
 modify it under the terms of the GNU General Public License as
 published by the Free Software Foundation; either version 2 of
 the License,or (at your option) any later version.

This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 General Public License for more details.

You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place-Suite 330, Boston, MA 02111-1307,
  USA. 
*)


(* :Title: xPand *)

(* :Author: Cyril Pitrou, Xavier Roy and Obinna Umeh *)

(* :Summary: Computer algebra for cosmological perturbations around homogenous space-times *)

(* :Brief Discussion:
   - Background 3+1 splitting of space-time;
   - Spatial sections are homogeneous (Friedmann-Lemaitre or Bianchi of all types) *)
  
(* :Context: xAct`xPand` *)

(* :Package Version: 0.4.3 *)

(* :Copyright: Cyril Pitrou, Xavier Roy and Obinna Umeh (2012-2018) *)

(* :History: see xPand.History *)

(* :Keywords: *)

(* :Source: xPand.nb *)

(* :Warning: *)

(* :Mathematica Version: 7.0 and later *)

(* :Limitations: *)


If[Unevaluated[xAct`xCore`Private`$LastPackage]===xAct`xCore`Private`$LastPackage,xAct`xCore`Private`$LastPackage="xAct`xPand`"];


Off[General::nostdvar]
Off[General::nostdopt]
BeginPackage["xAct`xPand`",{"xAct`xPert`","xAct`xTensor`","xAct`xPerm`","xAct`xCore`","xAct`ExpressionManipulation`"}]


If[Not@OrderedQ@Map[Last,{$xTensorVersionExpected,xAct`xTensor`$Version}],Throw@Message[General::versions,"xTensor",xAct`xTensor`$Version,$xTensorVersionExpected]]

If[Not@OrderedQ@Map[Last,{$xPertVersionExpected,xAct`xPert`$Version}],Throw@Message[General::versions,"xPert",xAct`xPert`$Version,$xPertVersionExpected]]


Print[xAct`xCore`Private`bars];
Print["Package xAct`xPand`  version ",$Version[[1]],", ",$Version[[2]]];
Print["CopyRight (C) 2012-2018, Cyril Pitrou, Xavier Roy and Obinna Umeh under the General Public License."];


Off[General::shdw]
xAct`xPand`Disclaimer[]:=Print["These are points 11 and 12 of the General Public License:\n\nBECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM `AS IS\.b4 WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.\n\nIN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES."]
On[General::shdw]


If[xAct`xCore`Private`$LastPackage==="xAct`xPand`",
Unset[xAct`xCore`Private`$LastPackage];
Print[xAct`xCore`Private`bars];
Print["These packages come with ABSOLUTELY NO WARRANTY; for details type Disclaimer[]. This is free software, and you are welcome to redistribute it under certain conditions. See the General Public License for details."];
Print[xAct`xCore`Private`bars]];


$CovDFormat="Prefix";
Message[General::nostdvar,"$CovDFormat","Prefix"];


(*** VERSIONS ***)

$Version::usage="$Version is a global variable giving the version of the package xPand in use.";

$xTensorVersionExpected::usage="$xTensorVersionExpected is a global variable giving the oldest possible version of the package xTensor which is required by the version of the package xPand in use.";

$xPertVersionExpected::usage="$xPertVersionExpected is a global variable giving the oldest possible version of the package xPert which is required by the version of the package xPand in use.";


(*** BOOLEANS ***)

BackgroundFieldMethod::usage = "BackgroundFieldMethod is a boolean option, set by default to 'False'. If set to 'True', then only the first order of the metric perturbation is kept, and it therefore stands for the total metric perturbation. 

This option should be placed at the begining of the notebook, right after loading the xPand package.";
(* Recently implemented. Does not work correctly. CP needs to fix that. To be checked... *)

$CacheRules::usage = "If set to 'True', the rules are remembered. This is potentially dangerous, so default is 'False'.";
(* This is obsolete. This will have to be erased. *)

$ConformalTime::usage = "If set to 'True' (default value), then the prime ' stands for the Lie derivative along the vector 'n' normal to the spatial slices. If the time coordinate is taken to be the proper time of the Eulerian observers (identified by the normal vector), then the Lie derivative along 'n' corresponds to the partial time derivative with respect to the conformal time, namely: \!\(\*SubscriptBox[\(\[ScriptCapitalL]\), \(n\)]\) .. = \[PartialD]/\[PartialD]\[Eta] .. . 

The second Label-Index (LI) of projected tensors represents the number of Lie derivatives along 'n' that are applied to them. It also corresponds to the number of partial derivatives with respect to conformal time. 

In the cases for which we prefer to use the cosmic time, we can choose to employ the second LI not for the number of \!\(\*SubscriptBox[\(\[ScriptCapitalL]\), \(n\)]\), but rather for the number of 1/a \!\(\*SubscriptBox[\(\[ScriptCapitalL]\), \(n\)]\), where 'a' is the scale factor. This is done by setting '$ConformalTime' to 'False'. The label index then indicates a derivative with respect to cosmic time.

The output format of the second LI is a prime for a conformal time derivative, and a dot for a cosmic time derivative.";

$FirstOrderVectorPerturbations::usage = "Boolean option to decide if first order vector modes should be included. Default is 'True'.";

$FirstOrderTensorPerturbations::usage = "Boolean option to decide if first order tensor modes should be included. Default is 'True'.";

$OpenConstantsOfStructure::usage = "Boolean value to decide if the constants of structure should be parameterized as in Ellis & MacCallum (1969). If 'False', then the constants of structure are kept general. If 'True', then the parameterization is: \!\(\*SubscriptBox[SuperscriptBox[\(C\), \(i\)], \(\(jk\)\(\\\ \)\)]\)= \!\(\*SubscriptBox[\(\[Epsilon]\), \(jkl\)]\)\!\(\*SuperscriptBox[\(N\), \(li\)]\) + \!\(\*SubscriptBox[\(A\), \([j\)]\)\!\(\*SubscriptBox[SuperscriptBox[\(\[Delta]\), \(i\)], \(\(k\)\(]\)\)]\).";

$SortCovDAutomatic::usage = "Boolean value to decide whether or not xTensor should commute automatically the induced derivatives. Default is 'True' in order to optimize the canonicalization of expression. Note that some adhoc commutation relations defined in '$CommutecdRules' are also applied in order to enforce the transverse properties of vector and tensor fields, on the one hand, and the appearance of Laplacians, on the other hand.";


(*** ERROR MESSAGES ***)

(** SetSlicing **)

SetSlicing::dimension = "The dimension of the manifold is `1`; it is too small to perform a slicing.";

SetSlicing::invalidnormalvector = "The normal vector `1` used to perform the background slicing should have only one index, belonging to the VBundle of the manifold. To avoid issues, it is possible to not predefine the normal vector and let the function SetSlicing define it automatically.";

(********** The following error message is not yet used. **********)

SetSlicing::ninds = "The number of indices is insufficient to define the (N-1)+1 splitting. Use at least `1` indices.";

SetSlicing::noambientmetric = "The ambient metric `1` was not previously defined as a metric. Define an ambient metric first, and then call SetSlicing.";

(********** The following error message does not seem correct. It should be modified. **********)

SetSlicing::signature = "The signature of the metric `1` should be -1 in order to perform a (N-1)+1 background splitting along a time-like vector. This comes from the fact that in the current version, the norm of the time-like vectors is fixed to -1.";

(** DefProjectedTensor **)

(********** The following error message is not used anymore. **********)

DefProjectedTensor::invalidindices = "No label-index needs to be given for the definition of projected tensors. They are automatically added by subroutines.";

DefProjectedTensor::notdownindices = "In order to be defined, the projected tensor `1` has to be written with down indices.";

(** DefProjectedTensorProperties **)

DefProjectedTensorProperties::symmetrictensors="The current version of DefProjectedTensorProperties can only add properties to tensors that are, at least, fully symmetric.";

(** ExtractComponents **)

ExtractComponents::invalidprojector = "One of the projectors in `1` is not valid. The possible projectors are: 'Time' or '*NameOfNormalVector*', and 'Space' or '*NameOfInducedMetric*'.";

(** RulesSplitCovDsOfTensor **)

RulesSplitCovDsOfTensor::maxcovdnumber = "The number of CovDs was larger than `1`. The last (optional) argument of RulesSplitCovDsOfTensor needs to be set to a larger value in order to split correctly all covariant derivatives.";

(** ToxPand **)

(********** The following two error messages are not yet edited. The first one is not referred to, and the utility of the second one has to be discussed. **********)

ToxPand::invalidindices = "Incorrect indices for tensor `1`";

ToxPand::invalidconffactor = "A conformal factor should have two label indices, one for the perturbation order which is always zero, and one for the number of Lie derivatives along u. And in the argument of SplitGaugeChange, you need only to put the name of the conformal factor (this is not a[LI[0],LI[0]] but just a).";

xPand::makeboxes = "For anisotropic manifolds, the second label-index can be only interpreted as a Lie derivative for indices down.";


(*** FUNCTIONS ***)

SetSlicing::usage = "SetSlicing[g,n,nNorm,h,cd,{cdPost,cdPre},SpaceTimeType]. 

Given a background manifold of dimension 'N' and an ambient metric 'g' on that manifold, this function performs a (N-1)+1 slicing. 

More precisely, it builds the induced metric 'h' on the hypersurfaces, the associated covariant derivative 'cd', and the vector 'n' normal to these hypersurfaces. It also sets the type 'SpaceTimeType' for the hypersurfaces. The latter argument can take a value among {Anisotropic, BianchiB, BianchiA, BianchiI, FLCurved, FLFlat, Minkowski}. The norm 'nNorm' of the vector 'n' is set by default to -1 if omitted.";

$CommutecdRules::usage = "Set of remembered rules for the commutation of the induced derivatives of the hypersurfaces. This is used to handle the successive operations of two, three or four induced derivatices.";

Conformal::usage = "Conformal[metricfinal][metric1,metric2][expr] performs a conformal transformation of 'expr' from 'metric1' to 'metric2' and expresses the result in terms of 'metricfinal' (along with its associated covariant derivative and curvature tensors). If needed, the function first reformulates 'expr' in terms of 'metric1'. The two metrics 'metric1' and 'metric2' have to be conformally related by DefConformalMetric. 

Conformal[metric1,metric2][expr] provides the final result with respect to the ambient metric, which is the first metric defined on the manifold.";

ConformalWeight::usage = "ConformalWeight[tensor[inds]] specifies the power of the conformal factor applied in the conformal transformation of the tensor 'tensor'. It corresponds to the number of down indices, minus the number of up indices, plus the quantity ConformalWeight[tensor] (without the indices). The latter is set by default to zero. 

For instance, in the setting ConformalWeight[tensor] = 0, we have: ConformalWeight[tensor[\[Alpha],\[Beta]]] = -2, and ConformalWeight[tensor[-\[Alpha],-\[Beta]]] = 2.";

DefConformalMetric::usage = "DefConformalMetric[g,S] defines a metric conformally related to 'g', with the conformal factor 'S'. The new metric is named 'gS2', and its inverse is given by 'Inv[gS2]'. If other metrics are conformally related to 'g', then 'gS2' is related to them by transitivity of the conformal transformation. 

When calling the function SetSlicing, the conformal factor is chosen to be the scale factor 'a[h]' (or, equivalently, 'ah'). The constructed metric 'gah2' then corresponds to the 'real' background metric (i.e. the one that would be obtained from the physical metric of the manifold).";

DefMetricFields::usage = "DefMetricFields[g,gpert,h,parameter] defines the perturbed metric 'gpert' from its background value 'g'. The argument 'h' is the induced metric on the spatial hypersurfaces, and the optional argument 'parameter' (set by default to \[Epsilon]) is the symbol used in the perturbative expansions. Several tensors that are used to decompose 'gpert' are defined when calling this function.";

DefMatterFields::usage = "DefMatterFields[uf,ufpert,h,parameter] defines the perturbed vector field 'ufpert' from its background value 'uf'. The argument 'h' is the induced metric on the spatial hypersurfaces, and the optional argument 'parameter' (set by default to \[Epsilon]) is the symbol used in the perturbative expansions. Several tensors that are used to decompose 'ufpert' are defined when calling this function.";

DefProjectedTensor::usage = "DefProjectedTensor[Name[inds],h,options] defines a projected tensor on the spatial hypersurfaces. 

'Name' refers to the name of the tensor, and 'inds' (which has to be with down indices) defines its rank. The argument 'h' is the induced metric on the spatial hypersurfaces. 

There exist 3 types of options. 

The first type is: 'TensorProperties', by default set to 'Transverse', 'Traceless' and 'SymmetricTensor'. These can be overwritten with 'TensorProperties->ListOfProperties', where 'ListOfProperties' is a list whose elements can be: 'Transverse', 'Traceless' and 'SymmetricTensor'. 

The second option is: 'SpaceTimesOfDefinition', by default set to 'Background' and 'Perturbed'. These can be also overwritten with 'SpaceTimesOfDefinition->SpaceTimesList', where 'SpaceTimesList' is a list whose elements can be: 'Background' and 'Perturbed'. 

The last option is: 'PrintAs', by default set to 'Identity'. This can be overwritten with 'PrintAs->String', in order to specify the output form of the tensor." ;

ExtractComponents::usage = "ExtractComponents[expr,h,ListOfProjectors,ListOfFreeIndices] projects 'expr' along its time component(s) (i.e. along the vector normal to the hypersurfaces) and along its space component(s) (i.e. onto the hypersurfaces), being defined by the background slicing associated with the induced metric 'h'. 

'ListOfProjectors' is a list whose elements are among: 'Time' (or, equivalently, '*NameOfNormalVector*') and 'Space' (or, equivalently, '*NameOfInducedMetric*'). The length of the list has to equal the rank of the tensors of 'expr'. The free indices (in canonical order) of 'expr' are then projected according to the 'ListOfProjectors'. 

'ListOfFreeIndices' is the list of free indices on which the ListOfProjectors is to be applied. If this is omitted, then it is the list of free indices in alphabetic order. 

When the argument 'ListOfProjectors' is omitted, all the projections for rank-0 to rank-2 tensors are displayed in a table.";

(********** The following command, ExtractOrder, will have to be re-edited (a check of the command is in order). **********)

ExtractOrder::usage = "ExtractOrder[expr,order] displays the expression at order 'order' of 'expr'.";

IndicesUp::usage = "IndicesUp[expr] raises the indices of all the tensors of 'expr' using the ambient metric.";
IndicesDown::usage = "IndicesDown[expr] lowers the indices of all the tensors of 'expr' using the ambient inverse metric.";



(********** The following command, SplitMatter, will have to be re-edited. **********)

SplitMatter::usage = "SplitMatter[uf,ufpert,normuf,h,gauge,order,tiltspecification] builds the rules to split the perturbed four-velocity 'ufpert' (with index up) in a given 'gauge', and according to the background slicing associated with the induced metric 'h'. The rules are computed to the order 'order'. 

'normuf' is the squared norm of the four-velocity 'uf' (usually taken as -1), and the optional argument 'tiltspecification' (set by default to 'NotTilted') allows for a tilt between the fluid four-velocity and the vector normal to the hypersurfaces. 

The perturbations of all the fields that depend on 'uf' are also computed.";

(********** The following command, SplitMetric, will have to be re-edited. **********)

SplitMetric::usage = "SplitMetric[g,gpert,h,gauge] builds the rules to split the perturbed metric 'gpert' in a given gauge and according to the background slicing associated with the induced metric 'h'.";


(********* This command was added in version 0.4.2. More documentation is needed and syntax can still evolve ********)
SplitNormalVector::usage = "SplitNormalVector[h,order], where h is an induced metric, and order is the order up to which the perturbations of the normal vector should be computed in function of the perturbed metric. This outputs a list of rules, one for each order of perturbation of the normal vector. All rules are expressed in function of the perturbed metric";


$MyRules::usage = "List of global rules that are remembered by xPand.";

(*ConstantsDecompositionRules::usage = "ConstantsDecompositionRules[h] provides the rules to transform the constants of structure according to the Bianchi type of the homogeneous hypersurfaces associated with the induced metric 'h'. This option is only valid for three-dimensional hypersurfaces. 

The parametrization follows Ellis & MacCallum (1969), that is, the constants of structure are decomposed as: Subscript[C^i, jk ]= Subscript[\[Epsilon], jkl]N^li + Subscript[A, [j]Subscript[\[Delta]^i, k]].";*)

ToBianchiType::usage = "ToBianchiType[h][expr] applies to expr the necessary rules to express the constants of structure according to the Bianchi type of the homogeneous hypersurfaces associated with the induced metric 'h'. This option is only valid for three-dimensional hypersurfaces. 

The parametrization follows Ellis & MacCallum (1969), that is, the constants of structure are decomposed as: \!\(\*SubscriptBox[SuperscriptBox[\(C\), \(i\)], \(\(jk\)\(\\\ \)\)]\)= \!\(\*SubscriptBox[\(\[Epsilon]\), \(jkl\)]\)\!\(\*SuperscriptBox[\(N\), \(li\)]\) + \!\(\*SubscriptBox[\(A\), \([j\)]\)\!\(\*SubscriptBox[SuperscriptBox[\(\[Delta]\), \(i\)], \(\(k\)\(]\)\)]\).";

(*ConstantsOfStructureRules::usage = "ConstantsOfStructureRules[h]. 

Given an induced metric h corresponding to a non trivial Bianchi type, ConstantsOfStructureRules[h] gives all the rules to express the Riemann, the Ricci and the Ricci scalar in function of the constants of structure of the spatial sections.";*)

ToConstantsOfStructure::usage = "ToConstantsOfStructure[h][expr]. 

Given an induced metric h corresponding to a non trivial Bianchi type, ToConstantsOfStructure[h][expr] applies to expr all the rules to express the Riemann, the Ricci and the Ricci scalar in function of the constants of structure of the spatial sections.";

(*(********** The following command, RulesCovDsOfTensor, is not yet defined. **********)

RulesCovDsOfTensor::usage = "";*)


(* ************** *)

ToCosmicTime::usage = "Transforms an expression where the second label index of projected tensor means \!\(\*SubscriptBox[\(\[ScriptCapitalL]\), SuperscriptBox[\(n\), \(\[Alpha]\)]]\), that is a derivative with respect to conformal time, into an expression where the second index means '1/a *\!\(\*SubscriptBox[\(\[ScriptCapitalL]\), SuperscriptBox[\(n\), \(\[Alpha]\)]]\)', that is a derivative with respect to cosmic time. The global variable $ConformalTime is set to False in the process.";

ToConformalTime::usage = "Transforms an expression where the second label index of projected tensor means '1/a * \!\(\*SubscriptBox[\(\[ScriptCapitalL]\), SuperscriptBox[\(n\), \(\[Alpha]\)]]\)', that is a derivative with respect to cosmic time, into an expression where the second index means '\!\(\*SubscriptBox[\(\[ScriptCapitalL]\), SuperscriptBox[\(n\), \(\[Alpha]\)]]\)', that is a derivative with respect to conformal time. The global variable $ConformalTime is set to True in the process.";

(********** The following command, RulesVelocitySpatial, is not yet edited. **********)

RulesVelocitySpatial::usage = "RulesVelocitySpatial[h,uf,duf,NormVectorSquare,gauge,order,tiltedvalue] gives the standard rules for the spatial velocity in a given gauge. h is the induced metric, uf and duf are the fluid velocity and its perturbation. NormVectorSquare is the norm squared of velocity, and gauge is the gauge choice. Order is the order up to which the perturbations are computed. The last argument, which can be '\"Tilted\"' or '\"NotTiled\"' is by default set to \"NotTilted\" and is optional.";

SplitGaugeChange::usage = "SplitGaugeChange[expr,ListPairs,\[Xi],h,order] performs a gauge transformation of order 'order' on 'expr', with a conformal transformation of conformal factor 'a[h]', and splits the result according to the background slicing associated with the induced metric 'h'. 

'\[Xi]' is the general vector field used for the gauge transformation. It must have a Label-Index (LI) for the order of perturbation (i.e. it has to be of the form \[Xi][LI[order],index]). 'ListPairs' is the list of rules used for the background slicing. 

The theory and formulas for the power expansion of gauge transformation was developed by Bruni et al. in: Class. Quantum Grav. 14, 2585 (1997).";

(********** The following command, SplitFieldsAndGaugeChange, is not yet edited. **********)

SplitFieldsAndGaugeChange::usage = "SplitFieldsAndGaugeChange[expr,g,dg,uf,duf,h,order,tiltspecification] performs a gauge transformation on expr where expr must involve metric perturbations (or curvature tensor associated with this metric) or matter field perturbation such as the fluid velocity, and splits the results according to the most general gauge choice. \nh is the induced metric, and refers to the background n+1 splitting. g is the metric and dg is the metric perturbation. uf is the fluid velocity and duf is the fluid velocity perturbation (when index is up). \n\[Xi] is the general vector field used for the gauge transformation which must have a label index for the order of perturbation (it is of the form \[Xi][LI[order],Index]), and order is the order of the gauge transformation. 'tiltspecification' is optional and by default is \"NotTilted\". \nThe theory and formulas for the power expansion of changes of gauge was developed by Bruni et al. in reference Class. Quantum Grav. 14, 2585 (1997).";

(*SplitPerturbationsOld::usage = "This is the old version of SplitPerturbations. It is obsolete, and it will have to be removed.";*)

SplitPerturbations::usage="SplitPerturbations[expr,ListPairs,h] splits the expression 'expr', using the list of rules 'ListPairs', according to the background slicing associated with the induced metric 'h'.";

(********** The following command, ToMetric, will have to be re-edited after its modification. **********)

ToMetric::usage = "ToMetric[expr,metric] formulates 'expr' in terms of 'metric', and its associated covariant derivative and curvature tensors.";

ToxPand::usage = "ToxPand[expr,gpert,uf,ufpert,h,gauge,order,tiltspecification]. 

This function first performs on 'expr' a conformal transformation with the scale factor 'a[h]' (except for a Minkowski 'SpaceType'); then it perturbs the result up to the order 'order' by means of xPert tools; and finally, it splits the result using the rules associated with the 'gauge', and according to the background slicing associated with the induced metric 'h'. 'tiltspecification' is optional and by default is \"NotTilted\" 

All the perturbed fields are derived from the perturbed metric 'gpert' and the perturbed fluid four-velocity 'ufpert'.

The syntax for multifluid is xPand[expr,gpert,{{uf1,uf1pert},{uf2,uf2pert},...},h,gauge,order,tiltspecification]. 
";

ToxPandFromRules::usage = "ToxPandFromRules[expr,ListPairs,h,order]. 

This function first performs on 'expr' a conformal transformation with the scale factor 'a[h]' (except for a Minkowski 'SpaceType'); then it perturbs the result up to the order 'order' by means of xPert tools; and finally, it splits the result using the list of rules 'ListPairs' and according to the background slicing associated with the induced metric 'h'. 

Compared to the function ToxPand, the present function requires to give the rules for the splitting of the perturbations.";

VisualizeTensor::usage = "VisualizeTensor[expr,h] displays in a table the projections of 'expr' along its time components (i.e. along the vector normal to the hypersurfaces) and along its space components (i.e. onto the hypersurfaces), being defined by the background slicing associated with the induced metric 'h'. 

'expr' has to be an expression involving tensors of rank 2 only.";
(* VisualizeTensor will be put as private, since the function ExtractComponents already handle this. *)


(*** LISTS OF FIELDS ***)

$ListFieldsBackgroundOnly::usage="$ListFieldsBackgroundOnly[h], with 'h' being the induced metric on the hypersurfaces, is the list of fields that live on the background only. It is automatically defined when calling the function SetSlicing.";

$ListFieldsPerturbedOnly::usage="$ListFieldsPerturbedOnly[h], with 'h' being the induced metric on the hypersurfaces, is the list of fields that live on the perturbed manifold only. It contains the metric perturbations. It is automatically defined when calling the function DefMetricFields.";

$ListFieldsBackgroundAndPerturbed::usage="$ListFieldsBackgroundAndPerturbed[h,uf], with 'h' and 'uf' being respectively the induced metric on the hypersurfaces and the fluid four-velocity, is the list of fields that live on the background manifold and that can be perturbed. It is automatically defined when calling the function DefMatterFields.";


(*** PREDEFINED FIELDS: Geometrical quantities ***)

a::usage = "a[h] is the scale factor of the spatial hypersurfaces whose induced metric is given by 'h'. It is automatically defined when calling the function SetSlicing. 

The command a[h] also builds the symbol 'ah', so one can use a[h] or, equivalently, ah. 

The default printed form of a[h][] is: 'a'.";

H::usage = "H[h] is the Hubble factor of the spatial hypersurfaces whose induced metric is given by 'h'. It is automatically defined when calling the function SetSlicing. 

The command H[h] also builds the symbol 'Hh', so one can use H[h] or, equivalently, Hh. 

The default printed form of H[h][] is: '\[ScriptCapitalH]'.";

Connection::usage = "Connection[h][index,index,index] are the coefficients of the connection associated with the induced metric 'h'. They are automatically defined when calling the function SetSlicing. 

The command Connection[h] builds the symbol 'Connectionh', so one can use Connection[h][index,index,index] or, equivalently, Connectionh[index,index,index].";

CS::usage = "CS[h][index,index,index] are the constants of structure on the hypersurfaces associated with the induced metric 'h'. They are automatically defined when calling the function SetSlicing. 

The command CS[h] builds the symbol 'CSh', so one can use CS[h][index,index,index] or, equivalently, CSh[index,index,index]. 

The default printed form of CS[h][-index,-index,-index] is: '\!\(\*SubscriptBox[\(C\), \(index\\\ index\\\ index\)]\)'.";

nt::usage = "nt[h][index,index] is the symmetric tensor part involved in the decomposition of the constants of structure: \!\(\*SubscriptBox[SuperscriptBox[\(C\), \(i\)], \(\(jk\)\(\\\ \)\)]\)= \!\(\*SubscriptBox[\(\[Epsilon]\), \(jkl\)]\)\!\(\*SuperscriptBox[\(N\), \(li\)]\) + \!\(\*SubscriptBox[\(A\), \([j\)]\)\!\(\*SubscriptBox[SuperscriptBox[\(\[Delta]\), \(i\)], \(\(k\)\(]\)\)]\), in Ellis and MacCallum (1969). It is automatically defined when calling the function SetSlicing for Bianchi space-times of type other then I. The presence as an argument of the induced metric 'h' is reminiscent of the choice of the background slicing. 

The command nt[h] builds the symbol 'nth', so one can use nt[h][index,index] or, equivalently, nth[index,index]. 

The default printed form of nt[h][-index,-index] is: '\!\(\*SubscriptBox[\(N\), \(index\\\ index\)]\)'.";

av::usage = "av[h][index] is the 'vector' part involved in the decomposition of the constants of structure: \!\(\*SubscriptBox[SuperscriptBox[\(C\), \(i\)], \(\(jk\)\(\\\ \)\)]\)= \!\(\*SubscriptBox[\(\[Epsilon]\), \(jkl\)]\)\!\(\*SuperscriptBox[\(N\), \(li\)]\) + \!\(\*SubscriptBox[\(A\), \([j\)]\)\!\(\*SubscriptBox[SuperscriptBox[\(\[Delta]\), \(i\)], \(\(k\)\(]\)\)]\), in Ellis & MacCallum (1969). It is automatically defined when calling the function SetSlicing for Bianchi space-times of type other than I. The presence as an argument of the induced metric 'h' is reminiscent of the choice of the background slicing.

The command av[h] builds the symbol 'avh', so one can use av[h][index] or, equivalently, avh[index]. 

The default printed form of av[h][-index] is: '\!\(\*SubscriptBox[\(A\), \(index\)]\)'.";

\[ScriptK]::usage = "\[ScriptK][h] is the (constant) curvature parameter of the isotropic and homogeneous hypersurfaces associated with the induced metric 'h'. It is automatically defined when calling the function SetSlicing for Friedmann-Lemaitre cosmologies. 

The command \[ScriptK][h] also builds the symbol '\[ScriptK]h', so one can use \[ScriptK][h] or, equivalently, \[ScriptK]h. 

The default printed form of \[ScriptK][h][] is: \[ScriptK].";


(*** PREDEFINED FIELDS: Scalar field and fluid quantities ***)

\[CurlyPhi]::usage = "\[CurlyPhi][LI[p],LI[q]] is a general scalar field, automatically defined when calling the function DefMatterFields. 

The first Label-Index (LI) gives the order 'p' of '\[CurlyPhi]' (hence, for p = 0, '\[CurlyPhi]' corresponds to the background value of the scalar field). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of '\[CurlyPhi]'. 

The commands \[CurlyPhi][] and \[CurlyPhi][LI[p]] are defined to be equal to \[CurlyPhi][LI[0],LI[0]] and \[CurlyPhi][LI[p],LI[0]], respectively.";

\[Rho]::usage = "\[Rho][uf][LI[p],LI[q]] is the energy density of the fluid, automatically defined when calling the function DefMatterFields. It depends on the fluid four-velocity 'uf'. 

The first Label-Index (LI) gives the order 'p' of '\[Rho]' (hence, for p = 0, '\[Rho]' corresponds to the background value of the energy density). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of '\[Rho]'. 

The commands \[Rho][uf][] and \[Rho][uf][LI[p]] are defined to be equal to \[Rho][uf][LI[0],LI[0]] and \[Rho][uf][LI[p],LI[0]], respectively. 

The command \[Rho][uf] builds the symbol '\[Rho]uf', so one can use \[Rho][uf][LI[p],LI[q]] or, equivalently, \[Rho]uf[LI[p],LI[q]].";
	
P::usage = "P[uf][LI[p],LI[q]] is the pressure of the fluid, automatically defined when calling the function DefMatterFields. It depends on the fluid four-velocity 'uf'. 

The first Label-Index (LI) gives the order 'p' of 'P' (hence, for p = 0, 'P' corresponds to the background value of the pressure). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'P'. 

The commands P[uf][] and P[uf][LI[p]] are defined to be equal to P[uf][LI[0],LI[0]] and P[uf][LI[p],LI[0]], respectively. 

The command P[uf] builds the symbol 'Puf', so one can use P[uf][LI[p],LI[q]] or, equivalently, Puf[LI[p],LI[q]].";

(* \[CapitalPi]::usage = "\[CapitalPi][uf][LI[p],LI[q],index,index] is the anisotropic stress tensor of the fluid, automatically defined when calling the function DefMatterFields. It depends on the fluid four-velocity 'uf'. 

The first Label-Index (LI) gives the order 'p' of '\[CapitalPi]' (hence, for p = 0, '\[CapitalPi]' corresponds to the (vanishing) background value of the anisotropic stress). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of '\[CapitalPi]'. 

The commands \[CapitalPi][uf][index,index] and \[CapitalPi][uf][LI[p],index,index] are defined to be equal to \[CapitalPi][uf][LI[0],LI[0],index,index] and \[CapitalPi][uf][LI[p],LI[0],index,index], respectively. 

The command \[CapitalPi][uf] builds the symbol '\[CapitalPi]uf', so one can use \[CapitalPi][uf][LI[p],LI[q],index,index] or, equivalently, \[CapitalPi]uf[LI[p],LI[q],index,index]."; *)


(*** PREDEFINED FIELDS: Perturbed fluid four-velocity ***)

V0::usage = "V0[h,uf][LI[p],LI[q]] is the time part of the perturbed fluid four-velocity 'ufpert' (with index up): ufpert = V0 \[Times] '*NameOfNormalVector*' + Vspat. It is automatically defined when calling the function DefMatterFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of 'V0' (hence, if the fluid is not tilted in the background, 'V0' equals 1 for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'V0'. 

The commands V0[h,uf][] and V0[h,uf][LI[p]] are defined to be equal to V0[h,uf][LI[0],LI[0]] and V0[h,uf][LI[p],LI[0]], respectively. 

The command V0[h,uf] builds the symbol 'V0huf', so one can use V0[h,uf][LI[p],LI[q]] or, equivalently, V0huf[LI[p],LI[q]]. 

The default printed form of V0[h,uf][] is: 'V0'.";

Vspat::usage = "Vspat[h,uf][LI[p],LI[q],index] is the spatial part of the perturbed fluid four-velocity 'ufpert' (with index up): ufpert = V0 \[Times] '*NameOfNormalVector*' + Vspat. It is automatically defined when calling the function DefMatterFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of 'Vspat' (hence, if the fluid is not tilted in the background, 'Vspat' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Vspat'. 

The commands Vspat[h,uf][index] and Vspat[h,uf][LI[p],index] are defined to be equal to Vspat[h,uf][LI[0],LI[0],index] and Vspat[h,uf][LI[p],LI[0],index], respectively. 

The command Vspat[h,uf] builds the symbol 'Vspathuf', so one can use Vspat[h,uf][LI[p],LI[q],index] or, equivalently, Vspathuf[LI[p],LI[q],index]. 

The default printed form of Vspat[h,uf][-index] is: '\!\(\*SubscriptBox[\(\[ScriptCapitalV]\), \(index\)]\)'.";

Vs::usage = "Vs[h,uf][LI[p],LI[q]] is the scalar part of the scalar-vector decomposition of 'Vspat', where 'Vspat' is the spatial part of the perturbed fluid four-velocity 'ufpert' (with index up). It is automatically defined when calling the function DefMatterFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of 'Vs' (hence, if the fluid is not tilted in the background, 'Vs' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Vs'. 

The commands Vs[h,uf][] and Vs[h,uf][LI[p]] are defined to be equal to Vs[h,uf][LI[0],LI[0]] and Vs[h,uf][LI[p],LI[0]], respectively. 

The command Vs[h,uf] builds the symbol 'Vshuf', so one can use Vs[h,uf][LI[p],LI[q]] or, equivalently, Vshuf[LI[p],LI[q]]. 

The default printed form of Vs[h,uf][] is: 'V'.";

Vv::usage = "Vv[h,uf][LI[p],LI[q],index] is the vector part of the scalar-vector decomposition of 'Vspat', where 'Vspat' is the spatial part of the perturbed fluid four-velocity 'ufpert' (with index up). It is automatically defined when calling the function DefMatterFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of 'Vv' (hence, if the fluid is not tilted in the background, 'Vv' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Vv'. 

The commands Vv[h,uf][index] and Vv[h,uf][LI[p],index] are defined to be equal to Vv[h,uf][LI[0],LI[0],index] and Vv[h,uf][LI[p],LI[0],index], respectively. 

The command Vv[h,uf] builds the symbol 'Vvhuf', so one can use Vv[h,uf][LI[p],LI[q],index] or, equivalently, Vvhuf[LI[p],LI[q],index]. 

The default printed form of Vv[h,uf][-index] is: '\!\(\*SubscriptBox[\(V\), \(index\)]\)'.";


(*** PREDEFINED FIELDS: Perturbation of the metric ***)

(** Time component of the metric perturbation **)

\[Phi]::usage = "\[Phi][h][LI[p],LI[q]] is the first Bardeen potential, corresponding to the 00 component of the metric perturbation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of the perturbation '\[Phi]' (hence '\[Phi]' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of '\[Phi]'. 

The commands \[Phi][h][] and \[Phi][h][LI[p]] are defined to be equal to \[Phi][h][LI[0],LI[0]] and \[Phi][h][LI[p],LI[0]], respectively. 

The command \[Phi][h] builds the symbol '\[Phi]h', so one can use \[Phi][h][LI[p],LI[q]] or, equivalently, \[Phi]h[LI[p],LI[q]].";

(** 0i components of the metric perturbation **)

Bs::usage = "Bs[h][LI[p],LI[q]] is the scalar part of the scalar-vector decomposition of the 0i components of the metric perturbation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of the perturbation 'Bs' (hence 'Bs' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Bs'. 

The commands Bs[h][] and Bs[h][LI[p]] are defined to be equal to Bs[h][LI[0],LI[0]] and Bs[h][LI[p],LI[0]], respectively. 

The command Bs[h] builds the symbol 'Bsh', so one can use Bs[h][LI[p],LI[q]] or, equivalently, Bsh[LI[p],LI[q]]. 

The default printed form of Bs[h][] is: 'B'.";

Bv::usage = "Bv[h][LI[p],LI[q],index] is the vector part of the scalar-vector decomposition of the 0i components of the metric perturbation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of the perturbation 'Bv' (hence 'Bv' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Bv'. 

The commands Bv[h][index] and Bv[h][LI[p],index] are defined to be equal to Bv[h][LI[0],LI[0],index] and Bv[h][LI[p],LI[0],index], respectively. 

The command Bv[h] builds the symbol 'Bvh', so one can use Bv[h][LI[p],LI[q],index] or, equivalently, Bvh[LI[p],LI[q],index]. 

The default printed form of Bv[h][-index] is: '\!\(\*SubscriptBox[\(B\), \(index\)]\)'.";

(** Spatial components of the metric perturbation **)

\[Psi]::usage = "\[Psi][h][LI[p],LI[q]] is the second Bardeen potential, present in the spatial components of the metric perturbation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of the perturbation '\[Psi]' (hence '\[Psi]' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of '\[Psi]'. 

The commands \[Psi][h][] and \[Psi][h][LI[p]] are defined to be equal to \[Psi][h][LI[0],LI[0]] and \[Psi][h][LI[p],LI[0]], respectively. 

The command \[Psi][h] builds the symbol '\[Psi]h', so one can use \[Psi][h][LI[p],LI[q]] or, equivalently, \[Psi]h[LI[p],LI[q]].";

Es::usage = "Es[h][LI[p],LI[q]] is the scalar part of the scalar-vector-tensor decomposition of the spatial components of the metric perturbation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of the perturbation 'Es' (hence 'Es' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Es'. 

The commands Es[h][] and Es[h][LI[p]] are defined to be equal to Es[h][LI[0],LI[0]] and Es[h][LI[p],LI[0]], respectively. 

The command Es[h] builds the symbol 'Esh', so one can use Es[h][LI[p],LI[q]] or, equivalently, Esh[LI[p],LI[q]]. 

The default printed form of Es[h][] is: 'E'.";

Ev::usage = "Ev[h][LI[p],LI[q],index] is the vector part of the scalar-vector-tensor decomposition of the spatial components of the metric perturbation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of the perturbation 'Ev' (hence 'Ev' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Ev'. 

The commands Ev[h][index] and Ev[h][LI[p],index] are defined to be equal to Ev[h][LI[0],LI[0],index] and Ev[h][LI[p],LI[0],index], respectively. 

The command Ev[h] builds the symbol 'Evh', so one can use Ev[h][LI[p],LI[q],index] or, equivalently, Evh[LI[p],LI[q],index]. 

The default printed form of Ev[h][-index] is: '\!\(\*SubscriptBox[\(E\), \(index\)]\)'.";

Et::usage = "Et[h][LI[p],LI[q],index,index] is the tensor part of the scalar-vector-tensor decomposition of the spatial components of the metric perturbation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of the perturbation 'Et' (hence 'Et' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Et'. 

The commands Et[h][index,index] and Et[h][LI[p],index,index] are defined to be equal to Et[h][LI[0],LI[0],index,index] and Et[h][LI[p],LI[0],index,index], respectively. 

The command Et[h] builds the symbol 'Eth', so one can use Et[h][LI[p],LI[q],index,index] or, equivalently, Eth[LI[p],LI[q],index,index]. 

The default printed form of Et[h][-index,-index] is: '\!\(\*SubscriptBox[\(E\), \(index\\\ index\)]\)'.";


(*** PREDEFINED FIELDS: Gauge vector, its decomposition and perturbation parameter ***)

\[Xi]::usage = "\[Xi][h][LI[order],index] is the general vector field used in a gauge transformation. It is automatically defined when calling the function SplitFieldsAndGaugeChange. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The Label-Index (LI) gives the order 'order' of the vector field '\[Xi]' (hence '\[Xi]' is null for p = 0). 

The command \[Xi][h] builds the symbol '\[Xi]h', so one can use \[Xi][h][LI[order],index] or, equivalently, \[Xi]h[LI[order],index].";

T::usage = "T[h][LI[p],LI[q]] is the factor of the time component of the vector field '\[Xi]' used in a gauge transformation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of 'T' (hence 'T' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'T'. 

The commands T[h][] and T[h][LI[p]] are defined to be equal to T[h][LI[0],LI[0]] and T[h][LI[p],LI[0]], respectively. 

The command T[h] builds the symbol 'Th', so one can use T[h][LI[p],LI[q]] or, equivalently, Th[LI[p],LI[q]]. 

The default printed form of T[h][] is: 'T'.";

Ls::usage = "Ls[h][LI[p],LI[q]] is the scalar part of the scalar-vector decomposition of the spatial component of the vector field '\[Xi]' used in a gauge transformation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of 'Ls' (hence 'Ls' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Ls'. 

The commands Ls[h][] and Ls[h][LI[p]] are defined to be equal to Ls[h][LI[0],LI[0]] and Ls[h][LI[p],LI[0]], respectively. 

The command Ls[h] builds the symbol 'Lsh', so one can use Ls[h][LI[p],LI[q]] or, equivalently, Lsh[LI[p],LI[q]]. 

The default printed form of Ls[h][] is: 'L'.";

Lv::usage = "Lv[h][LI[p],LI[q],index] is the vector part of the scalar-vector decomposition of the spatial component of the vector field '\[Xi]' used in a gauge transformation. It is automatically defined when calling the function DefMetricFields. It depends on the choice of the background slicing, and hence on the induced metric 'h'. 

The first Label-Index (LI) gives the order 'p' of 'Lv' (hence 'Lv' is null for p = 0). The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of 'Lv'. 

The commands Lv[h][index] and Lv[h][LI[p],index] are defined to be equal to Lv[h][LI[0],LI[0],index] and Lv[h][LI[p],LI[0],index], respectively. 

The command Lv[h] builds the symbol 'Lvh', so one can use Lv[h][LI[p],LI[q],index] or, equivalently, Lvh[LI[p],LI[q],index]. 

The default printed form of Lv[h][-index] is: '\!\(\*SubscriptBox[\(L\), \(index\)]\)'.";

\[Epsilon]::usage = "\[Epsilon] is the default perturbation parameter. It is automatically defined when calling either the function DefMetricFields or the function DefMatterFields.";



\[ScriptCapitalN]0::usage = "\[ScriptCapitalN]0[h][LI[p],LI[q]] is the time part of the perturbation of the normal vector to constant time hypersurfaces (whith index down). It is automatically defined when calling SetSlicing.
In the ADM formalism, it is equal to the lapse function.

The first Label-Index (LI) gives the order 'p' of the perturbation '\[ScriptCapitalN]0'. The second LI represents the number 'q' of Lie derivatives along the vector normal to the hypersurfaces acting on the 'p'-order of '\[ScriptCapitalN]'. Hence \[ScriptCapitalN]0[h][LI[0],LI[0]] is the background lapse and is equal to 1.

The commands \[ScriptCapitalN]0[h][index] and \[ScriptCapitalN]0[h][LI[p],index] are defined to be equal to \[ScriptCapitalN]0[h][LI[0],LI[0],index] and \[ScriptCapitalN]0[h][LI[p],LI[0],index], respectively. 

The command \[ScriptCapitalN]0[h] builds the symbol '\[ScriptCapitalN]0h', so one can use \[ScriptCapitalN]0[h][LI[p],LI[q],index] or, equivalently, \[ScriptCapitalN]0h[LI[p],LI[q],index]. 

The default printed form of \[ScriptCapitalN]0[h][-index] is: '\!\(\*SubscriptBox[\(\[ScriptCapitalN]0\), \(index\)]\)'.";


\[ScriptCapitalN]::usage = "\[ScriptCapitalN][h][index] is the normal vector to constant time hypersurfaces. It is automatically defined when calling SetSlicing. Though on the background spacetime it is equal to \!\(\*SuperscriptBox[\(n\), \(\[Mu]\)]\), where \!\(\*SuperscriptBox[\(n\), \(\[Mu]\)]\) is the vector normal to h on the background space-time, this definition allows manipulation of the normal vector, like perturbation, conformal transformations, which would be impossible directly with \!\(\*SuperscriptBox[\(n\), \(\[Mu]\)]\) since many definitions are already automatized for it (e.g. normalization to -1).";


d\[ScriptCapitalN]::usage = "d\[ScriptCapitalN][h][LI[p],index] is the perturbation of the normal vector to constant time hypersurfaces. 'p' is the order of perturbation.

The command d\[ScriptCapitalN][h] builds the symbol 'd\[ScriptCapitalN]h', so one can use d\[ScriptCapitalN][h][LI[p],index] or, equivalently, d\[ScriptCapitalN]h[LI[p]index].";

\[ScriptCapitalN]i::usage = "\[ScriptCapitalN]i[h][index] is the shift vector (in the ADM formalism) of the normal vector to constant time hypersurfaces. It is defined as \[ScriptCapitalN]i[h][upindex] = -Projector[h][ \[ScriptCapitalN]0[h] \[ScriptCapitalN][h][upindex] ].
It is only meaningful when used with and upper index. Its main purpose is the computation of perturbations for the shift vector in the ADM formalism.

The command \[ScriptCapitalN]i[h] builds the symbol '\[ScriptCapitalN]ih', so one can use \[ScriptCapitalN]i[h][upindex] or, equivalently, \[ScriptCapitalN]ih[upindex]. 
";


(*** RESERVED WORDS AND PROTECTED NAMES ***)

SpaceTimesOfDefinition::usage = "";
TensorProperties::usage = "";


(*** TYPES OF GAUGE ***)

$ListOfGauges::usage = "Possible gauges are: AnyGauge, FluidComovingGauge, ScalarFieldComovingGauge FlatGauge, IsoDensityGauge, NewtonGauge and SynchronousGauge."


(*** TYPES OF MANIFOLD ***)

(*Anisotropic::usage="General type of spatial hypersurfaces for an anisotropic manifold. The boolean function '$OpenConstantsOfStructure' is automatically set to 'False'. For a four-dimensional manifold, the user can choose instead the particular types: BianchiB, BianchiA or BianchiI.";

BianchiB::usage = "Type of spatial hypersurfaces corresponding to a general Bianchi space-time. 

If '$OpenConstantsOfStructure' is set to 'True', then the constants of structure Subscript[C^i, jk] are expressed in terms of N^ij and A^i, according to: Subscript[C^i, jk ]= Subscript[\[Epsilon], jkl]N^li + Subscript[A, [j]Subscript[\[Delta]^i, k]].";

BianchiA::usage = "Type of spatial hypersurfaces corresponding to a Bianchi type A space-time. 

If '$OpenConstantsOfStructure' is set to 'True', then the constants of structure Subscript[C^i, jk] are expressed in terms of N^ij, according to: Subscript[C^i, jk ]= Subscript[\[Epsilon], jkl]N^li.";

BianchiI::usage = "Type of spatial hypersurfaces corresponding to a Bianchi type I space-time.";
FLCurved::usage = "Type of spatial hypersurfaces corresponding to a curved Friedmann-Lemaitre manifold (that is, with spatial curvature).";
FLFlat::usage = "Type of spatial hypersurfaces corresponding to a flat Friedmann-Lemaitre manifold (that is, without spatial curvature).";
Minkowski::usage = "Type of spatial hypersurfaces corresponding to a Minkowski manifold (that is, a flat FL space-time without expansion).";
*)

SpaceType::usage = "SpaceType[h], with 'h' being the induced metric on the hypersurfaces, displays the type of hypersurfaces chosen for the SetSlicing. It can be: Anisotropic, BianchiB, BianchiA, BianchiI, FLCurved, FLFlat or Minkowski.";

$ListOfSpaceTypes::usage = "Possible space types are: Anisotropic, BianchiB, BianchiA, BianchiI, FLCurved, FLFlat and Minkowski."


Begin["xAct`xPand`Private`"]


(** Why a scalar head inside a scalar function? **)
(** We choose that the Scalar head should be removed**)
Unprotect[xAct`xTensor`NoScalar];
xAct`xTensor`NoScalar[f_?ScalarFunctionQ[expr_]]:=f[NoScalar@expr]
Protect[xAct`xTensor`NoScalar];

(** To avoid complaints from MakeRule when the rule has already been defined and the lhs is zero. **)
(*xAct`xTensor`MakeRule[{lhs_?(Evaluate[#]===0&),rhs_,conditions___},options:OptionsPattern[]]:={};*)

(* We create our own function based on MakeRule, but separating the special problematic case. This avoids overloading a definition from xTensor.
It does not change the way xPand works, but that way, the code is clearly separated from xTensor. 
This is a much more polite way to modify the behaviour of xTensor. This method was advised by Leo Stein**)

SetAttributes[BuildRule,HoldFirst]
BuildRule[{lhs_?(Evaluate[#]===0&),rhs_,conditions___},options:OptionsPattern[]]:={};
BuildRule[{lhs_,rhs_,conditions___}]:=MakeRule[{lhs,rhs,conditions}];
BuildRule[{lhs_,rhs_,conditions___},options:OptionsPattern[]]:=MakeRule[{lhs,rhs,conditions},options];


(*** DEFAULT OPTIONS AND PROTECTED NAMES ***)

BackgroundFieldMethod = False;
$DebugInfoQ=False;
$CommutecdRules = {};
$ConformalTime = True;
$FirstOrderVectorPerturbations=True;
$FirstOrderTensorPerturbations=True;
$OpenConstantsOfStructure=True;
$PrePrint=ScreenDollarIndices;
$SortCovDAutomatic=True;

Off[RuleDelayed::rhs];


Block[{Print},DefInertHead[ProtectMyScalar,LinearQ->True]];

FixProtectScalar:=ProtectMyScalar[expr_]:>ProtectMyScalar[xAct`xTensor`Private`MathInputExpand[expr]]


(*** COLORATION OF THE PERTURBATION ORDERS ***)

DerCharacter:=If[$ConformalTime,"\[Prime]","."];
(* The output format of the Lie derivative with respect to the vector normal to the hypersurfaces is a prime when one considers the conformal time, and a dot when one considers the cosmic time. *)

PerturbationOrderColor[1]:=RGBColor[0.85,0,0];
PerturbationOrderColor[2]:=RGBColor[0,0.6,0];
PerturbationOrderColor[3]:=RGBColor[0,0,0.85];
PerturbationOrderColor[4]:=RGBColor[0.6,0.4,0.2];
PerturbationOrderColor[5]:=RGBColor[0.7,0,0.7];
PerturbationOrderColor[_]:=RGBColor[1,0.6,0];



xTensorFormStop[Tensor]

MakeBoxes[Tens_?xTensorQ[LI[p_?((IntegerQ[#]&&#>=0)&)],LI[q_?((IntegerQ[#]&&#>=0)&)],inds___],StandardForm]:=
xAct`xTensor`Private`interpretbox[Tens[LI[p],LI[q],inds],
If[(AnisotropyBool[SpaceType@InducedMetricOf[Tens]]||BianchiBool[SpaceType@InducedMetricOf[Tens]])&&Cases[{inds},_?UpIndexQ]=!={}&&q>=1,

Message[xPand::makeboxes];

RowBox[{OverscriptBox["\[Null]",RowBox[{"(",ToString[p],")"}]],"\[NegativeThinSpace]",OverscriptBox[MakeBoxes[Tens[inds],StandardForm],ToString[q]]}],

Switch[p,
0,RowBox[{OverscriptBox[MakeBoxes[Tens[inds],StandardForm],StringJoin@Table[DerCharacter,{i,1,q}]]}],
_,RowBox[{OverscriptBox["\[Null]",RowBox[{StyleBox["("<>ToString[p]<>")",FontColor->PerturbationOrderColor[p]]}]],"\[NegativeThinSpace]",OverscriptBox[MakeBoxes[Tens[inds],StandardForm],StringJoin@Table[DerCharacter,{i,1,q}]]}]
]
]
];

xTensorFormStart[Tensor]


(*** CONVENIENT FUNCTIONS TO BUILD RULES ***)

PatternLeft[(left_:> right_),ElementsToBePatterned_List]:=(left/.((#->Pattern[#,_])&/@ ElementsToBePatterned)):>right
PatternLeft[(left_-> right_),ElementsToBePatterned_List]:=(left/.((#->Pattern[#,_])&/@ ElementsToBePatterned))->right


(*** HANDLING EXPRESSIONS ***)

collect[expr_]:=Collect[expr,$PerturbationParameter,Identity]

(** fix by Jolyon Bloomfield. No Idea if this works or not. This should correct the bug found by Adam Solomon in August 2014 **)
fixScalar:=Scalar[expr_]:>Scalar[xAct`xTensor`Private`MathInputExpand[expr]]
org[expr_]:=Collect[ContractMetric[expr],$PerturbationParameter,ToCanonical[#/.fixScalar]&]
(* Old org before the fix*)
(* org[expr_]:=Collect[ContractMetric[expr],$PerturbationParameter,ToCanonical[#/.fixScalar]&]*)


TCnoCM[expr_]:=ToCanonical[expr,UseMetricOnVBundle->None]


(*** PRIVATE BOOLEAN FUNCTIONS ***)

(** Miscellaneous **)

$DefInfoQ=True;
(* If set to 'False', the information messages are not printed. *)

(** Testing expressions **)

AnyIndicesListQ[inds_List]:=Cases[inds,AnyIndices[_]]=!={}
(* If the list of indices contains the head 'AnyIndices', then AnyIndicesListQ[inds] returns 'True'; otherwise it returns 'False'. *)

DefTensorQ[symb_]:=Cases[$Tensors,symb]==={symb}
(* If 'symb' has been defined as a tensor, then DefTensorQ[symb] returns 'True'; otherwise it returns 'False'. *)

GaugeQ[strg_]:=Cases[{"AnyGauge","FluidComovingGauge","ScalarFieldComovingGauge","FlatGauge","IsoDensityGauge","NewtonGauge","SynchronousGauge"},strg]==={strg}
(* If 'strg' matches one of the element in the list, then GaugeQ[strg] returns 'True'; otherwise it returns 'False'. *)

InducedMetricQ[symb_]:=If[MetricQ[symb],InducedFrom[symb]=!=Null,False]
(* If 'symb' is not defined as a metric, then InducedMetricQ[symb] returns 'False'. Otherwise, it checks whether 'symb' is an induced metric. 
If this is the case, then InducedMetricQ[symb] gives 'True'; otherwise it gives 'False'. *)

SpaceTimeQ[strg_]:=Cases[{"Anisotropic", "BianchiB","BianchiA","BianchiI","FLCurved","FLFlat","Minkowski"},strg]==={strg}
(* If 'strg' matches one of the element of the list, then SpaceTimeQ[strg] returns 'True'; otherwise it returns 'False'. *)

TensorNullQ[tens_]:=If[Cases[SlotsOfTensor[tens],Labels]=!={},Print["** Warning: the function TensorNullQ is only suited to test tensors defined without label-indices. ",tens],With[{inds=DummyIn/@SlotsOfTensor[tens]},tens@@inds===0]]
(* TensorNullQ[tens] returns 'True' if 'tens' is a null tensor and 'False' otherwise. This function is not suited to test tensors defined with label-indices. *)

(** Type of manifolds **)

BianchiBool[Spacetype_]:=(Spacetype==="BianchiB"||Spacetype==="BianchiA"||Spacetype==="BianchiI")
FlatSpaceBool[Spacetype_]:=(Spacetype==="BianchiI"||Spacetype==="FLFlat"||Spacetype==="Minkowski")
NoAnisotropyBool[Spacetype_]:=(Spacetype==="FLCurved"||Spacetype==="FLFlat"||Spacetype==="Minkowski")
AnisotropyBool[Spacetype_]:=(Spacetype==="Anisotropic")
(* The name 'AnisotropyBool' should be changed. *)


(*** BUILDING SYMBOLS: Geometrical quantities ***)

a[symb_]:=SymbolJoin[a,symb];
H[symb_]:=SymbolJoin[H,symb];

Connection[symb_]:=SymbolJoin[Connection,symb];
CS[symb_]:=SymbolJoin[CS,symb];
nt[symb_]:=SymbolJoin[nt,symb];
av[symb_]:=SymbolJoin[av,symb];

K[symb_]:=SymbolJoin[K,symb];
\[ScriptK][symb_]:=SymbolJoin[\[ScriptK],symb];


(*** BUILDING SYMBOLS: Fluid quantities ***)

\[Rho][symb_]:=SymbolJoin[\[Rho],symb];
P[symb_]:=SymbolJoin[P,symb];
(* \[CapitalPi][symb_]:=SymbolJoin[\[CapitalPi],symb]; *)


(*** BUILDING SYMBOLS: Perturbed fluid four-velocity ***)

V0[h_?InducedMetricQ,velocity_]:=SymbolJoin[V0,h,velocity];
Vspat[h_?InducedMetricQ,velocity_]:=SymbolJoin[Vspat,h,velocity];
Vs[h_?InducedMetricQ,velocity_]:=SymbolJoin[Vs,h,velocity];
Vv[h_?InducedMetricQ,velocity_]:=SymbolJoin[Vv,h,velocity];


(*** BUILDING SYMBOLS: Perturbation of the metric ***)

\[Phi][h_?InducedMetricQ]:=SymbolJoin[\[Phi],h];

Bs[h_?InducedMetricQ]:=SymbolJoin[Bs,h];
Bv[h_?InducedMetricQ]:=SymbolJoin[Bv,h];

\[Psi][h_?InducedMetricQ]:=SymbolJoin[\[Psi],h];
Es[h_?InducedMetricQ]:=SymbolJoin[Es,h];
Ev[h_?InducedMetricQ]:=SymbolJoin[Ev,h];
Et[h_?InducedMetricQ]:=SymbolJoin[Et,h];


(*** BUILDING SYMBOLS: Gauge vector and decomposition ***)

\[Xi][h_?InducedMetricQ]:=SymbolJoin[\[Xi],h];

T[h_?InducedMetricQ]:=SymbolJoin[T,h];
Ls[h_?InducedMetricQ]:=SymbolJoin[Ls,h];
Lv[h_?InducedMetricQ]:=SymbolJoin[Lv,h];


$ListFieldsBackgroundOnly[h_?InducedMetricQ]:={{a[h][],"a"}, {H[h][],"\[ScriptCapitalH]"}};

$ListFieldsPerturbedOnly[h_?InducedMetricQ]:=With[{M=ManifoldOfCovD@CovDOfMetric@First@InducedFrom@h},
Block[{\[Mu]},
{\[Mu]}=GetIndicesOfVBundle[Tangent@M,1];
{{\[Phi][h][],"\[Phi]"},{Bs[h][],"B"},{Bv[h][-\[Mu]],"B"},{\[Psi][h][],"\[Psi]"},{Es[h][],"E"},{Ev[h][-\[Mu]],"E"},{Et[h][-\[Mu],-\[Nu]],"E"},{T[h][],"T"},{Ls[h][],"L"},{Lv[h][-\[Mu]],"L"}}
]
];

(* The former implementation of $ListFieldsBackgroundAndPerturbed failed to treat separately different fluids. *)
(* Obinna proposed the replace by the version below. *)
(* We thank Enrico Pajer for reporting this bug. *)

$ListFieldsBackgroundAndPerturbed[h_?InducedMetricQ,velocity_]:=With[{M=ManifoldOfCovD@CovDOfMetric@First@InducedFrom@h},
Block[{\[Mu],\[Nu]},
{\[Mu],\[Nu]}=GetIndicesOfVBundle[Tangent@M,2];
{{\[CurlyPhi][],"\[CurlyPhi]"},{\[Rho][velocity][],ToString@StringJoin[{"\[Rho]"},ToString[velocity]]},{P[velocity][],ToString@StringJoin[{"P"},ToString[velocity]]},{Vspat[h,velocity][-\[Mu]],ToString@StringJoin[{"\[ScriptCapitalV]",ToString[velocity]}]},{V0[h,velocity][],ToString@StringJoin[{"V0",ToString[velocity]}]},{Vs[h,velocity][],ToString@StringJoin[{"V",ToString[velocity]}]},{Vv[h,velocity][-\[Mu]],ToString@StringJoin[{"V",ToString[velocity]}]}}
]
];


(*** DEFAULT VALUES: DefProjectedTensor & DefProjectedTensorProperties ***)

Options[DefProjectedTensor]={PrintAs->Identity,TensorProperties->{"SymmetricTensor","Traceless","Transverse"},SpaceTimesOfDefinition->{"Background","Perturbed"}};

DefProjectedTensorQ[Name_,h___]:=False;
PropertiesList[Name_]:={};
InducedMetricOf[Name_]:={};


(*** MODULE: DefProjectedTensor ***)

DefProjectedTensor[Name_[inds___],h_?InducedMetricQ,options___?OptionQ]:=Catch@Module[{IndsNoLI,p,q,PrAs,SpaTimeDef,TensProp},

With[{M=ManifoldOfCovD@CovDOfMetric[First@InducedFrom@h]},

If[DefProjectedTensorQ[Name,h],

If[$DefInfoQ,
Throw@Print["** DefProjectedTensor: The projection properties on the hypersurfaces associated with the induced metric ", h," have already been defined for the tensor ", Name,"."],
Throw[Null]
];

];

If[DefProjectedTensorQ[Name],

If[$DefInfoQ,
Print["** DefProjectedTensor: Projection properties for the tensor ", Name," have been defined for another slicing. New projection properties on the hypersurfaces associated with the induced metric ", h," are now added."]
];

];

(* If projection properties for the tensor 'Name' have been defined for another slicing, the command 'Name[inds]', with 'inds' being abstract index, is substituted by 'Name[LI[0],LI[0],inds]' (refer to the function DefProjectedTensorProperties for this rule). The following command then serves to remove the possible label indices. *)
IndsNoLI=Select[{inds},Not@LIndexQ[#]&];

If[(Length[IndsNoLI]>= 1)&&(Not@AnyIndicesListQ[IndsNoLI])&&(Select[IndsNoLI,UpIndexQ] =!={}),
Throw[Message[DefProjectedTensor::notdownindices,Name]]];

(********** The following Message is not yet defined. **********)

If[(Length[IndsNoLI]>= 2)&&(AnyIndicesListQ[IndsNoLI]),
Throw[Message[DefProjectedTensor::invalidanyindices,Name]]];

PrAs=PrintAs/.CheckOptions[options]/.Options[DefProjectedTensor];
SpaTimeDef=SpaceTimesOfDefinition/.CheckOptions[options]/.Options[DefProjectedTensor];
TensProp=TensorProperties/.CheckOptions[options]/.Options[DefProjectedTensor];


If[DefProjectedTensorQ[Name]||DefTensorQ[Name],

If[$DefInfoQ,
Print["** DefProjectedTensor: The tensor ", Name," already exists. The projection properties on the hypersurfaces associated with the induced metric ", h," are now defined."]],
DefTensor[Name[LI[p],LI[q],IndsNoLI/.List->Sequence],M,Symmetric[IndsNoLI],PrintAs->PrAs]

];


(* Definition of the projection properties for the tensor 'Name'. *)
If[Not@AnyIndicesListQ[IndsNoLI],
DefProjectedTensorProperties[Name,IndsNoLI/.List->Sequence,h,TensProp,SpaTimeDef],
DefProjectedTensorPropertiesAnyIndices[Name,h,TensProp,SpaTimeDef]];

]
]

SetNumberOfArguments[DefProjectedTensor,{2,Infinity}]
Protect[DefProjectedTensor];


(*** MODULE: DefProjectedTensorProperties ***)

DefProjectedTensorProperties[Name_,inds___?DownIndexQ,h_?InducedMetricQ,Properties_List,Spacetimes_List]:=Catch@Module[{prot,Lengthindices},

With[{

g=First@InducedFrom@h,
n=Last@InducedFrom@h},

With[{

cd=CovDOfMetric[h],
M=ManifoldOfCovD[CovDOfMetric[g]],

SymmetricBool=(Cases[Properties,"SymmetricTensor"]==={"SymmetricTensor"}),
TracelessBool=(Cases[Properties,"Traceless"]==={"Traceless"}),
TransverseBool=(Cases[Properties,"Transverse"]==={"Transverse"}),

BackgroundBool=(Cases[Spacetimes,"Background"]==={"Background"}),
PerturbedBool=(Cases[Spacetimes,"Perturbed"]==={"Perturbed"}),

ToCan=ToCanonical[#,UseMetricOnVBundle->None]&
},

If[Not[SymmetricBool]&&(Length[{inds}]>=2),
Throw@Message[DefProjectedTensorProperties::symmetrictensors]];

DefProjectedTensorQ[Name,h]^=True;

InducedMetricOf[Name]^=h;

If[DefProjectedTensorQ[Name]===False,

PropertiesList[Name]^=Join[Properties,Which[
Length[{inds}]===0,{"Scalar"},
Length[{inds}]===1,{"Vector"},
Length[{inds}]>=2,{"Tensor"}]];


Name[indices___?AIndexQ]:=Name[LI[0],LI[0],indices]
/;(Length[{indices}]===Length[{inds}]);

Name[LI[p_?((IntegerQ[#]&&#>=0)&)],indices___?AIndexQ]:=Name[LI[p],LI[0],indices]
/;(Length[{indices}]===Length[{inds}]);


(* For tensors of rank larger than or equal to 2, the following rules provide the traceless property. *)
If[(Length[{inds}]>=2)&&TracelessBool,

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,Dum_,indices2___,-Dum_,indices3___]:=0 /;( Length[Join[{indices1},{indices2},{indices3}]]+2===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,-Dum_,indices2___,Dum_,indices3___]:=0 /;( Length[Join[{indices1},{indices2},{indices3}]]+2===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,Dum_,indices2___,-Dum_,indices3___]:=
Module[{Dummy1,Dummy2},

Dummy1=DummyIn[Tangent[M]];
Dummy2=DummyIn[Tangent[M]];

If[$ConformalTime,1,1/a[h][]]*(LieD[n[Dummy1]][ToCan[Name[LI[p],LI[q-1],indices1,Dum,indices2,-Dum,indices3]]]

+2*Name[LI[p],LI[q-1],indices1,-Dummy1,indices2,-Dummy2,indices3]ExtrinsicK[h][Dummy1,Dummy2])
 
]/;( Length[Join[{indices1},{indices2},{indices3}]]+2===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,-Dum_,indices2___,Dum_,indices3___]:=
Module[{Dummy1,Dummy2},

Dummy1=DummyIn[Tangent[M]];
Dummy2=DummyIn[Tangent[M]];

If[$ConformalTime,1,1/a[h][]]*(LieD[n[Dummy1]][ToCan[Name[LI[p],LI[q-1],indices1,Dum,indices2,-Dum,indices3]]]

+2 * Name[LI[p],LI[q-1],indices1,-Dummy1,indices2,-Dummy2,indices3]ExtrinsicK[h][Dummy1,Dummy2])
 
]/;( Length[Join[{indices1},{indices2},{indices3}]]+2===Length[{inds}]);

];


PerturbationOrder[Name[LI[0],LI[q_],indices___]]^:=Catch@Module[{},

If[Length[{indices}]=!=Length[{inds}],
Throw[Print["** Warning: The number of indices for the tensor ", Name," is incorrect."]]];

If[Cases[AIndexQ/@{indices},False]=!={},
Throw[Print["** Warning: The indices of the tensor ", Name," have to be abstract indices."]]];

If[(IntegerQ[q]&&q>=0),
0,
(*Throw[Print["** Warning: The second label-index has to be a positive integer."]]]*)
If[NumericQ[q],Throw[Print["** Warning: The second label-index has to be a positive integer."]]]]

];


(* TO DO: Here, we have to comment on the commutation between perturbing and Lie deriving. This commutation stems from the fact that the Lie derivative is a background Lie derivative. *)
If[PerturbedBool,

(* For any perturbed tensors, we have: *)
PerturbationOrder[Name[LI[p_],LI[q_],indices___]]^:=Catch@Module[{},

If[Length[{indices}]=!=Length[{inds}],
Throw[Print["** Warning: The number of indices for the tensor ", Name," is incorrect."]]];

If[Cases[AIndexQ/@{indices},False]=!={},
Throw[Print["** Warning: The indices of the tensor ", Name," have to be abstract indices."]]];

If[(IntegerQ[p]&&p>=1)&&(IntegerQ[q]&&q>=0),
p,
(*Throw[Print["** Warning: The label-indices have to be positive integers."]]]*)
If[NumericQ[p]&&NumericQ[q],Throw[Print["** Warning: The label-indices have to be positive integers."]]]]

];
(* (the rule for 'p = 0' is defined above) *)

(* and for perturbed scalar quantities, we have: *)
Perturbation[Name[LI[p_],LI[q_]]]^:=Catch@If[(IntegerQ[p]&&p>=0)&&(IntegerQ[q]&&q>=0),
Name[LI[p+1],LI[q]],
Throw[Print["** Warning: The label-indices have to be positive integers."]]];

Perturbation[Name[LI[p_],LI[q_]],PertOrder_]^:=Catch@Module[{},

If[Not[IntegerQ[PertOrder]],
Throw[Print["** Warning: The order of the perturbation has to be an integer."]]];

If[(IntegerQ[p]&&p>=0)&&(IntegerQ[q]&&q>=0),
Name[LI[p+PertOrder],LI[q]],
Throw[Print["** Warning: The label-indices have to be positive integers."]]]

];

];


(* For tensors living on the background only, we have: *)
If[BackgroundBool&&Not[PerturbedBool],

Perturbation[Name[LI[0],LI[q_],indices___],PertOrder_]^:=Catch@Module[{},

If[Length[{indices}]=!=Length[{inds}],
Throw[Print["** Warning: The number of indices for the tensor ", Name," is incorrect."]]];

If[Cases[AIndexQ/@{indices},False]=!={},
Throw[Print["** Warning: The indices of the tensor ", Name," have to be abstract indices."]]];

If[Not[IntegerQ[q]&&q>=0],
Throw[Print["** Warning: The second label-index has to be a positive integer."]]];

If[PertOrder===0,
Name[LI[0],LI[q],indices],
If[(IntegerQ[PertOrder]&&PertOrder>=1),
0,
Throw[Print["** Warning: The order of the perturbation has to be a positive integer."]]]]

];

Perturbation[Name[LI[0],LI[q_],indices___]]^:=Catch@Module[{},

If[Length[{indices}]=!=Length[{inds}],
Throw[Print["** Warning: The number of indices for the tensor ", Name," is incorrect."]]];

If[Cases[AIndexQ/@{indices},False]=!={},
Throw[Print["** Warning: The indices of the tensor ", Name," have to be abstract indices."]]];

If[IntegerQ[q]&&q>=0,
0,
Throw[Print["** Warning: The second label-index has to be a positive integer."]]]

];

];


(* For pure perturbations (i.e. for tensors without background values), we have: *)
If[Not[BackgroundBool]&&PerturbedBool,

(* This is the 0.4.0 version way to do it. We define a delayed 0 value for the background*)
(* However this is problematic when we want to perturb. Because then when we perturbe this quantity*)
(* we write Perturbed[quantity] but Mathematica will read Perturbed[0], and so this leads to 0.*)
(* Instead we should append this to a list of rules, that should be applied in SplitPerturbations *)
(* The easiest way is to define a set of global rules and to append a new rule each time there is a tensor vanishing on the background *)
	
(* Old implementation *)
(*Name[LI[0],LI[q_?((IntegerQ[#]&&#\[GreaterEqual]0)&)],indices___?AIndexQ]:=0
/;(Length[{indices}]===Length[{inds}]);*)

(* New implementation *)
Lengthindices=Length[{inds}];
$RulesVanishingBackgroundFields[h]=Append[$RulesVanishingBackgroundFields[h],Name[LI[0],LI[q_?((IntegerQ[#]&&#>=0)&)],indices___?AIndexQ]:>0
/;(Length[{indices}]===Lengthindices) ];

];

DefProjectedTensorQ[Name]^=True;

];


PropertiesList[Name]^=
Join[PropertiesList[Name],Which[
Length[{inds}]===1&&TransverseBool,{"SVT-Vector associated with the induced metric " h ""},Length[{inds}]>=2&&SymmetricBool&&TracelessBool&&TransverseBool,{"SVT-Tensor associated with the induced metric " h ""},Length[{inds}]>=0,{}]];
(* 'SVT-Vector' is a vector satisfying the SVT-decomposition properties (hence it is transverse). 
'SVT-Tensor' is a tensor satisfying the SVT-decomposition properties (hence it is symmetric, traceless and transverse).*)


(* The Lie derivative for tensors with indices down is represented by the second label-index. *)
Name/:LieD[n[Dum_]][Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices___?DownIndexQ]]:= If[$ConformalTime,1,a[h][]]*Name[LI[p],LI[q+1],indices]
/;(Length[{indices}]===Length[{inds}]);


(* For tensors of rank larger than or equal to 1, *)
If[Length[{inds}]>=1,

(* we bypass OrthogonalToVectorQ: *)
Name/:OrthogonalToVectorQ[n][Name]=True;

(* Projection conditions *)
(* The contraction with the vector normal to the hypersurfaces gives 0 (since the tensor is projected): *)

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,Dum_,indices2___]n[-Dum_]:=0 /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,-Dum_,indices2___]n[Dum_]:=0 /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

(* and the action of the projector onto the hypersurfaces is an invariant operation: *)
(* However this is left to post processing that's why we load this rule into $Rulecdh *)

$Rulecdh[h]=Append[$Rulecdh[h],
Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,Dum1_,indices2___]h[-Dum1_,-Dum2_]:>Name[LI[p],LI[0],indices1,-Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}])];

$Rulecdh[h]=Append[$Rulecdh[h],
Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,-Dum1_,indices2___]h[Dum1_?UpIndexQ,Dum2_?UpIndexQ]:>Name[LI[p],LI[0],indices1,Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}])];


Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,Dum1_,indices2___]h[-Dum1_,Dum2_?UpIndexQ]:=Name[LI[p],LI[0],indices1,Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,Dum1_,indices2___]h[Dum2_?UpIndexQ,-Dum1_]:=Name[LI[p],LI[0],indices1,Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,-Dum1_,indices2___]h[Dum1_,-Dum2_]:=Name[LI[p],LI[0],indices1,-Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,-Dum1_,indices2___]h[-Dum2_,Dum1_]:=Name[LI[p],LI[0],indices1,-Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

(* When there are Lie derivatives, h contracted automatically only when the free index is down. Otherwise it requires ContractMetric to perform contraction *)

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,Dum1_,indices2___]h[-Dum1_,-Dum2_]:=Name[LI[p],LI[q],indices1,-Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,Dum1_,indices2___]h[-Dum2_,-Dum1_]:=Name[LI[p],LI[q],indices1,-Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,-Dum1_,indices2___]h[Dum1_,-Dum2_]:=Name[LI[p],LI[q],indices1,-Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,-Dum1_,indices2___]h[-Dum2_,Dum1_]:=Name[LI[p],LI[q],indices1,-Dum2,indices2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);


(* g converted to h when contracted with a projected tensor. It is made automatic *)

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,Dum1_,indices2___]g[-Dum1_,Dum2_]:=Name[LI[p],LI[q],indices1,Dum1,indices2]h[-Dum1,Dum2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,Dum1_,indices2___]g[Dum2_,-Dum1_]:=Name[LI[p],LI[q],indices1,Dum1,indices2] h[Dum2,-Dum1] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,-Dum1_,indices2___]g[Dum1_,Dum2_]:=Name[LI[p],LI[q],indices1,-Dum1,indices2] h[Dum1,Dum2] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,-Dum1_,indices2___]g[Dum2_,Dum1_]:=Name[LI[p],LI[q],indices1,-Dum1,indices2] h[Dum2,Dum1] /;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

(* If the tensor belongs (at least) to the perturbed manifold, and if it is tranverse, *)
If[PerturbedBool && TransverseBool,

(* (transverse property) *)
Name/:cd[Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[0],indices1___,-Dum_,indices2___]]:=
0/;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

Name/:cd[-Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[0],indices1___,Dum_,indices2___]]:=
0/;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

(* then the first rule: D^aSubscript[\[ScriptCapitalL], n]Subscript[T, a... ]= Subscript[\[ScriptCapitalL], n](D^aSubscript[T, a...]) - Subscript[\[ScriptCapitalL], n](g^ab)Subscript[D, b]Subscript[Subscript[T, a], ...] + (g^ab)[Subscript[D, b],Subscript[\[ScriptCapitalL], n]]Subscript[T, a...] *)
Name/:cd[Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___?DownIndexQ,-Dum_,indices2___?DownIndexQ]]:=Module[{Dummy1,Dummy2},

Dummy1=DummyIn[Tangent[M]];
Dummy2=DummyIn[Tangent[M]];

ToCan@ContractMetric[

If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][cd[Dum][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]]

- ToCan[If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][g[Dum,Dummy2]]cd[-Dummy2][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]//MetricToProjector[#,h] &]

+g[Dum,Dummy2]ToCan[cd[-Dummy2][Name[LI[p],LI[q],indices1,-Dum,indices2]]-If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][cd[-Dummy2][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]]]

]
]/;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

(* and the second rule: Subscript[D, a]g^abSubscript[\[ScriptCapitalL], n]Subscript[Subscript[T, b], ... ]= Subscript[\[ScriptCapitalL], n](D^aSubscript[T, a...]) - Subscript[\[ScriptCapitalL], n](g^ab)Subscript[D, b]Subscript[Subscript[T, a], ...] + (g^ab)[Subscript[D, b],Subscript[\[ScriptCapitalL], n]]Subscript[T, a...] *)
(* XR remark 09/2015: careful that the index up means that it is raised after the Lie derivative is taken. Thus the rhs is the same as the previous rule.*)
Name/:cd[-Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___?DownIndexQ,Dum_,indices2___?DownIndexQ]]:=Module[{Dummy1,Dummy2},

Dummy1=DummyIn[Tangent[M]];
Dummy2=DummyIn[Tangent[M]];

ToCan@ContractMetric[

If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][cd[Dum][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]]

- ToCan[If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][g[Dum,Dummy2]]cd[-Dummy2][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]//MetricToProjector[#,h] &]

+ g[Dum,Dummy2]ToCan[cd[-Dummy2][Name[LI[p],LI[q],indices1,-Dum,indices2]]-If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][cd[-Dummy2][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]]]

]

]/;( Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

(* If the indices are not down, we separate them. In practice this happens almost never... *)
(* This is for the case where we have at least one derivative, for which the meaning is only when the index is down.*)
Name/:cd[Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,indup_?UpIndexQ,indices3___,-Dum_,indices2___]]:=Module[{Dummy},Dummy=DummyIn[Tangent[M]];

g[indup,Dummy]cd[Dum][Name[LI[p],LI[q],indices1,-Dummy,indices3,-Dum,indices2]]

]/;( Length[Join[{indices1},{indices2},{indices3}]]+2===Length[{inds}]);

Name/:cd[Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,-Dum_,indices3___,indup_?UpIndexQ,indices2___]]:=Module[{Dummy},Dummy=DummyIn[Tangent[M]];

g[indup,Dummy]cd[Dum][Name[LI[p],LI[q],indices1,-Dum,indices3,-Dummy,indices2]]

]/;( Length[Join[{indices1},{indices2},{indices3}]]+2===Length[{inds}]);

Name/:cd[-Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,indup_?UpIndexQ,indices3___,Dum_,indices2___]]:=Module[{Dummy},Dummy=DummyIn[Tangent[M]];

g[indup,Dummy]cd[-Dum][Name[LI[p],LI[q],indices1,-Dummy,indices3,Dum,indices2]]

]/;( Length[Join[{indices1},{indices2},{indices3}]]+2===Length[{inds}]);

Name/:cd[-Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,Dum_,indices3___,indup_?UpIndexQ,indices2___]]:=Module[{Dummy},Dummy=DummyIn[Tangent[M]];

g[indup,Dummy]cd[-Dum][Name[LI[p],LI[q],indices1,Dum,indices3,-Dummy,indices2]]

]/;( Length[Join[{indices1},{indices2},{indices3}]]+2===Length[{inds}]);



];
];


(* The following rule serves to express the three-covariant derivative of (homogeneous) background quantities in terms of the connection components. For instance, 'Subscript[D, i]Subscript[\[Omega], j]' is replaced by: '- Subscript[\[CapitalGamma]^a, ij]Subscript[\[Omega], a]'. *)
If[BackgroundBool,

Name/:cd[Dummy2_][Name[LI[0],LI[q_?((IntegerQ[#]&&#>=0)&)],indices___]]:=Module[{Dummy1},

Dummy1=DummyIn[Tangent[M]];

ToCan[Plus@@((-Connection[h][Dummy1,Dummy2,#]ReplaceIndex[Evaluate[Name[LI[0],LI[q],indices]],#->-Dummy1])&/@ {indices})]

]/;(Length[{indices}]===Length[{inds}]);

];
(* Confer the function SetSlicing for the definition of 'Connection' and its properties. *)
(* Is the above relation correct for space-time indices? To check! *)
(* What is the need for 'Evaluate'? I forgot. *)


(* In principle, the following rule needs not to be used. *)
(* However, it makes sure that, if something went wrong in our algorithm, and the Lie derivatives acts on a tensor with up indices, then this converted to down indices.*)
(* This ensures that the algorithm for splitting the covariant derivatives into induced derivatives and Lie derivative always works. *)
Name/:LieD[n[Dummy1_]][Name[LI[p_?((IntegerQ[#]&&#>=0)&)],LI[q_?((IntegerQ[#]&&#>=0)&)],indices1___?AIndexQ,IndexUp_?UpIndexQ,indices2___?AIndexQ]]:=Module[{Dum},

Dum=DummyIn[Tangent[M]];

g[IndexUp,Dum]LieD[n[Dummy1]][Name[LI[p],LI[q],indices1,-Dum,indices2]]+ LieD[n[Dummy1]][g[IndexUp,Dum]] Name[LI[p],LI[q],indices1,-Dum,indices2]

]/;(Length[Join[{indices1},{indices2}]]+1===Length[{inds}]);

]
]
]


(*** MODULE: DefProjectedTensorPropertiesAnyIndices ***)

DefProjectedTensorPropertiesAnyIndices[Name_,h_?InducedMetricQ,Properties_List,Spacetimes_List]:=Catch@Module[{},

With[{

g=First@InducedFrom@h,
n=Last@InducedFrom@h,

cd=CovDOfMetric[h]},

With[{

M=ManifoldOfCovD[CovDOfMetric[g]],

SymmetricBool=(Cases[Properties,"SymmetricTensor"]==={"SymmetricTensor"}),
TracelessBool=(Cases[Properties,"Traceless"]==={"Traceless"}),
TransverseBool=(Cases[Properties,"Transverse"]==={"Transverse"}),

BackgroundBool=(Cases[Spacetimes,"Background"]==={"Background"}),
PerturbedBool=(Cases[Spacetimes,"Perturbed"]==={"Perturbed"}),

ToCan=ToCanonical[#,UseMetricOnVBundle->None]&
},

DefProjectedTensorQ[Name,h]^=True;

InducedMetricOf[Name]^=h;

If[DefProjectedTensorQ[Name]===False,

PropertiesList[Name]^=
Join[Properties,{"Scalar","Vector","Tensor"}];

Name[indices___?AIndexQ]:=Name[LI[0],LI[0],indices];

Name[LI[p_?((IntegerQ[#]&&#>=0)&)],indices___?AIndexQ]:=Name[LI[p],LI[0],indices];


(* The following rule provides the symmetry property. *)
If[SymmetricBool,
SymmetryGroupOfTensor[Name[LI[p_],LI[q_],indices__]]^:=Symmetric[2+Range@Length@{indices},Cycles];
];


(* The following rules provide the traceless property. *)
If[TracelessBool,

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,Dum_,indices2___,-Dum_,indices3___]:=0 ;

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,-Dum_,indices2___,Dum_,indices3___]:=0 ;

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,Dum_,indices2___,-Dum_,indices3___]:=
Module[{Dummy1,Dummy2},

Dummy1=DummyIn[Tangent[M]];
Dummy2=DummyIn[Tangent[M]];

If[$ConformalTime,1,1/a[h][]]*(LieD[n[Dummy1]][ToCan[Name[LI[p],LI[q-1],indices1,Dum,indices2,-Dum,indices3]]]

+2 * Name[LI[p],LI[q-1],indices1,-Dummy1,indices2,-Dummy2,indices3]ExtrinsicK[h][Dummy1,Dummy2])
 
];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,-Dum_,indices2___,Dum_,indices3___]:=
Module[{Dummy1,Dummy2},

Dummy1=DummyIn[Tangent[M]];
Dummy2=DummyIn[Tangent[M]];

If[$ConformalTime,1,1/a[h][]]*(LieD[n[Dummy1]][ToCan[Name[LI[p],LI[q-1],indices1,Dum,indices2,-Dum,indices3]]]

+2 * Name[LI[p],LI[q-1],indices1,-Dummy1,indices2,-Dummy2,indices3]ExtrinsicK[h][Dummy1,Dummy2])
 
];

];


PerturbationOrder[Name[LI[0],LI[q_],indices___]]^:=Catch@Module[{},

If[Cases[AIndexQ/@{indices},False]=!={},
Throw[Print["** Warning: The indices of the tensor ", Name," have to be abstract indices."]]];

If[(IntegerQ[q]&&q>=0),
0,
(*Throw[Print["** Warning: The second label-index has to be a positive integer."]]]*)
If[NumericQ[q],Throw[Print["** Warning: The second label-index has to be a positive integer."]]]]

];


(* TO DO: We here have to comment on the commutation between perturbing and Lie deriving. This commutation stems from the fact that the Lie derivative is a background Lie derivative. *)
If[PerturbedBool,

(* For any perturbed tensors, we have: *)
PerturbationOrder[Name[LI[p_],LI[q_],indices___]]^:=Catch@Module[{},

If[Cases[AIndexQ/@{indices},False]=!={},
Throw[Print["** Warning: The indices of the tensor ", Name," have to be abstract indices."]]];

If[(IntegerQ[p]&&p>=1)&&(IntegerQ[q]&&q>=0),
p,
(*Throw[Print["** Warning: The label-indices have to be positive integers."]]]*)
If[NumericQ[p]&&NumericQ[q],Throw[Print["** Warning: The label-indices have to be positive integers."]]]]

];
(* The rule for 'p = 0' is defined at the beginning of the module. *)

(* and for perturbed scalar quantities, we have: *)
Perturbation[Name[LI[p_],LI[q_]]]^:=Catch@If[(IntegerQ[p]&&p>=0)&&(IntegerQ[q]&&q>=0),
Name[LI[p+1],LI[q]],
Throw[Print["** Warning: The label-indices have to be positive integers."]]];

Perturbation[Name[LI[p_],LI[q_]],PertOrder_]^:=Catch@Module[{},

If[Not[IntegerQ[PertOrder]],
Throw[Print["** Warning: The order of the perturbation has to be an integer."]]];

If[(IntegerQ[p]&&p>=0)&&(IntegerQ[q]&&q>=0),
Name[LI[p+PertOrder],LI[q]],
Throw[Print["** Warning: The label-indices have to be positive integers."]]]

];

];


(* For tensors living on the background only, we have: *)
If[BackgroundBool&&Not[PerturbedBool],

Perturbation[Name[LI[0],LI[q_],indices___],PertOrder_]^:=Catch@Module[{},

If[Cases[AIndexQ/@{indices},False]=!={},
Throw[Print["** Warning: The indices of the tensor ", Name," have to be abstract indices."]]];

If[Not[IntegerQ[q]&&q>=0],
Throw[Print["** Warning: The second label-index has to be a positive integer."]]];

If[PertOrder===0,
Name[LI[0],LI[q],indices],
If[(IntegerQ[PertOrder]&&PertOrder>=1),
0,
Throw[Print["** Warning: The order of the perturbation has to be a positive integer."]]]]

];

Perturbation[Name[LI[0],LI[q_],indices___]]^:=Catch@Module[{},

If[Cases[AIndexQ/@{indices},False]=!={},
Throw[Print["** Warning: The indices of the tensor ", Name," have to be abstract indices."]]];

If[IntegerQ[q]&&q>=0,
0,
Throw[Print["** Warning: The second label-index has to be a positive integer."]]]

];

Name/:Name[LI[p_?((IntegerQ[#]&&#>=1)&)],LI[q_],indices___]:=0

];


(* For perturbations (i.e. for tensors without background values), we have: *)
If[Not[BackgroundBool]&&PerturbedBool,

(* Again we do not Set it to 0 but rather add a rule which is used in SplitPerturbations*)
(* So when we perturbed a quantity this avoids to do Perturbed[0]. *)
(*Name[LI[0],LI[q_?((IntegerQ[#]&&#\[GreaterEqual]0)&)],indices___?AIndexQ]:=0*)
$RulesVanishingBackgroundFields[h]=Append[$RulesVanishingBackgroundFields[h],Name[LI[0],LI[q_?((IntegerQ[#]&&#>=0)&)],indices___?AIndexQ]:>0];
];

DefProjectedTensorQ[Name]^=True;

];


PropertiesList[Name]^=
Join[PropertiesList[Name],If[SymmetricBool&&TracelessBool&&TransverseBool,{"SVT-Vector associated with the induced metric "h "","SVT-Tensor associated with the induced metric "h ""},{}]];


(* The Lie derivative for tensors with indices down is represented by the second label-index. *)
Name/:LieD[n[Dum_]][Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices___?DownIndexQ]]:=If[$ConformalTime,1,a[h][]]*Name[LI[p],LI[q+1],indices];


(* We bypass OrthogonalToVectorQ *)
Name/:OrthogonalToVectorQ[n][Name]=True;

(* The contraction with the vector normal to the hypersurfaces gives 0 (since the tensor is projected): *)

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,Dum_,indices2___]n[-Dum_]:=0;

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,-Dum_,indices2___]n[Dum_]:=0;

(* and the action of the projector onto the hypersurfaces is an invariant operation: *)

$Rulecdh[h]=Append[$Rulecdh[h],
Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,Dum1_,indices2___]h[-Dum1_,-Dum2_]:>Name[LI[p],LI[0],indices1,-Dum2,indices2] ];

$Rulecdh[h]=Append[$Rulecdh[h],
Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,-Dum1_,indices2___]h[Dum1_?UpIndexQ,Dum2_?UpIndexQ]:>Name[LI[p],LI[0],indices1,Dum2,indices2]];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,Dum1_,indices2___]h[-Dum1_,Dum2_?UpIndexQ]:=Name[LI[p],LI[0],indices1,Dum2,indices2] ;

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,Dum1_,indices2___]h[Dum2_?UpIndexQ,-Dum1_]:=Name[LI[p],LI[0],indices1,Dum2,indices2];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,-Dum1_,indices2___]h[Dum1_,-Dum2_]:=Name[LI[p],LI[0],indices1,-Dum2,indices2];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[0],indices1___,-Dum1_,indices2___]h[-Dum2_,Dum1_]:=Name[LI[p],LI[0],indices1,-Dum2,indices2];

(* When there are Lie derivatives, h contracted automatically only when the free index is down. Otherwise it requires ContractMetric to perform contraction *)

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,Dum1_,indices2___]h[-Dum1_,-Dum2_]:=Name[LI[p],LI[q],indices1,-Dum2,indices2];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,Dum1_,indices2___]h[-Dum2_,-Dum1_]:=Name[LI[p],LI[q],indices1,-Dum2,indices2];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,-Dum1_,indices2___]h[Dum1_,-Dum2_]:=Name[LI[p],LI[q],indices1,-Dum2,indices2];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,-Dum1_,indices2___]h[-Dum2_,Dum1_]:=Name[LI[p],LI[q],indices1,-Dum2,indices2];


(*g to h  is automatic *)

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,Dum1_,indices2___]g[-Dum1_,Dum2_]:=Name[LI[p],LI[q],indices1,Dum1,indices2]h[-Dum1,Dum2];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,Dum1_,indices2___]g[Dum2_,-Dum1_]:=Name[LI[p],LI[q],indices1,Dum1,indices2] h[Dum2,-Dum1];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,-Dum1_,indices2___]g[Dum1_,Dum2_]:=Name[LI[p],LI[q],indices1,-Dum1,indices2] h[Dum1,Dum2];

Name/:Name[LI[p_?((IntegerQ[#] && #>=0) &)],LI[q_?((IntegerQ[#] && #>=0) &)],indices1___,-Dum1_,indices2___]g[Dum2_,Dum1_]:=Name[LI[p],LI[q],indices1,-Dum1,indices2] h[Dum2,Dum1];





(* If the tensor belongs (at least) to the perturbed manifold, and if it is transverse, *)
If[PerturbedBool&& TransverseBool,

(* (transverse property) *)
Name/:cd[Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[0],indices1___,-Dum_,indices2___]]:=0;

Name/:cd[-Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[0],indices1___,Dum_,indices2___]]:=0;

(* then the first rule: D^aSubscript[\[ScriptCapitalL], n]Subscript[T, a... ]= Subscript[\[ScriptCapitalL], n](D^aSubscript[T, a...]) - Subscript[\[ScriptCapitalL], n](g^ab)Subscript[D, b]Subscript[Subscript[T, a], ...] + (g^ab)[Subscript[D, b],Subscript[\[ScriptCapitalL], n]]Subscript[T, a...] *)
Name/:cd[Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___?DownIndexQ,-Dum_,indices2___?DownIndexQ]]:=Module[{Dummy1,Dummy2},

Dummy1=DummyIn[Tangent[M]];
Dummy2=DummyIn[Tangent[M]];

ToCan@ContractMetric[

If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][cd[Dum][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]]

- ToCan[If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][g[Dum,Dummy2]]cd[-Dummy2][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]//MetricToProjector[#,h] &]

+g[Dum,Dummy2]ToCan[cd[-Dummy2][Name[LI[p],LI[q],indices1,-Dum,indices2]]-If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][cd[-Dummy2][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]]]

]
];

(* and the second rule: Subscript[D, a]Subscript[\[ScriptCapitalL], n]Subscript[T^a, ... ]= Subscript[\[ScriptCapitalL], n](D^aSubscript[T, a...]) - Subscript[\[ScriptCapitalL], n](g^ab)Subscript[D, b]Subscript[Subscript[T, a], ...] + (g^ab)[Subscript[D, b],Subscript[\[ScriptCapitalL], n]]Subscript[T, a...] *)
Name/:cd[-Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___?DownIndexQ,Dum_,indices2___?DownIndexQ]]:=Module[{Dummy1,Dummy2},

Dummy1=DummyIn[Tangent[M]];
Dummy2=DummyIn[Tangent[M]];

ToCan@ContractMetric[

If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][cd[Dum][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]]

- ToCan[If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][g[Dum,Dummy2]]cd[-Dummy2][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]//MetricToProjector[#,h] &]

+ g[Dum,Dummy2]ToCan[cd[-Dummy2][Name[LI[p],LI[q],indices1,-Dum,indices2]]-If[$ConformalTime,1,1/a[h][]]*LieD[n[Dummy1]][cd[-Dummy2][Name[LI[p],LI[q-1],indices1,-Dum,indices2]]]]

]

];


(* If the indices are not down, we separate them. In practice this happens almost never... *)
Name/:cd[Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,indup_?UpIndexQ,indices3___,-Dum_,indices2___]]:=Module[{Dummy},Dummy=DummyIn[Tangent[M]];

g[indup,Dummy]cd[Dum][Name[LI[p],LI[q],indices1,-Dummy,indices3,-Dum,indices2]]

];

Name/:cd[Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,-Dum_,indices3___,indup_?UpIndexQ,indices2___]]:=Module[{Dummy},Dummy=DummyIn[Tangent[M]];

g[indup,Dummy]cd[Dum][Name[LI[p],LI[q],indices1,-Dum,indices3,-Dummy,indices2]]

];

Name/:cd[-Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,indup_?UpIndexQ,indices3___,Dum_,indices2___]]:=Module[{Dummy},Dummy=DummyIn[Tangent[M]];

g[indup,Dummy]cd[-Dum][Name[LI[p],LI[q],indices1,-Dummy,indices3,Dum,indices2]]

];

Name/:cd[-Dum_][Name[LI[p_?((IntegerQ[#] && #>=1) &)],LI[q_?((IntegerQ[#] && #>=1) &)],indices1___,Dum_,indices3___,indup_?UpIndexQ,indices2___]]:=Module[{Dummy},Dummy=DummyIn[Tangent[M]];

g[indup,Dummy]cd[-Dum][Name[LI[p],LI[q],indices1,Dum,indices3,-Dummy,indices2]]

];

];


(* The following rule serves to express the three-covariant derivative of (homogeneous) background quantities in terms of the connection components. For instance, 'Subscript[D, i]Subscript[\[Omega], j]' is replaced by: '- Subscript[\[CapitalGamma]^a, ij]Subscript[\[Omega], a]'. *)
If[BackgroundBool,

Name/:cd[Dummy2_][Name[LI[0],LI[q_?((IntegerQ[#]&&#>=0)&)],indices___]]:=Module[{Dummy1},

Dummy1=DummyIn[Tangent[M]];

ToCan[Plus@@((-Connection[h][Dummy1,Dummy2,#]ReplaceIndex[Evaluate[Name[LI[0],LI[q],indices]],#->-Dummy1])&/@ {indices})]

];

];
(* Confer the function SetSlicing for the definition of 'Connection' and its properties. *)
(* Is the above relation correct for space-time indices? To check! *)
(* What is the need for 'Evaluate'? *)


(* In principle, the following rule needs not to be used. *)
(* Indeed, the algorithm of SPlitPerturbation is designed so that there is never any Lie derivative acting on tensors with up indices.*)
(* If that was the case because of a strange use of xPand, this rule ensure that the Lie derivative of a tensor with up indices is converted into a Lie derivative with down indices. This is to make sure that the algo always works.*)
Name/:LieD[n[Dummy1_]][Name[LI[p_?((IntegerQ[#]&&#>=0)&)],LI[q_?((IntegerQ[#]&&#>=0)&)],Indices1___?AIndexQ,IndexUp_?UpIndexQ,Indices2___?AIndexQ]]:=Module[{Dum},

Dum=DummyIn[Tangent[M]];

g[IndexUp,Dum]LieD[n[Dummy1]][Name[LI[p],LI[q],Indices1,-Dum,Indices2]]+ LieD[n[Dummy1]][g[IndexUp,Dum]] Name[LI[p],LI[q],Indices1,-Dum,Indices2]

];

]
]
]


TensorHasAnyIndicesQ[tens_]:=AnyIndicesListQ[SlotsOfTensor[tens]](*(Length@Cases[SlotsOfTensor[tens],AnyIndices[_]]>0)*)


(* Note that only down indices can be passed to this function *)
IsElementInList[el_,list_]:=Length@Cases[list,el]>=1;

(*TranverseTensorQ[tens_]:=xTensorQ[tens]&&IsElementInList[Transverse,PropertiesList[tens]];*)
ScalarTensorQ[tens_]:=xTensorQ[tens]&&IsElementInList["Scalar",PropertiesList[tens]];
VectorTensorQ[tens_]:=xTensorQ[tens]&&IsElementInList["Vector",PropertiesList[tens]];
TensorTensorQ[tens_]:=xTensorQ[tens]&&IsElementInList["Tensor",PropertiesList[tens]];


(* Rules used to build Automaticrules. We first use MakeRule, and then we use these rule to put patterns on the Lable indices of the left hand side.*)

PatternTensorLeftScalar[(lhs_:>rhs_),TensDummy_?ScalarTensorQ[inds___]]:=Block[{TScal},
(Evaluate[lhs/.(TensDummy->PatternTest[Pattern[TensDummy,_],ScalarTensorQ])]:>(rhs))/.(TensDummy->TScal)
]


PatternTensorLeftVector[(lhs_:>rhs_),TensDummy_?VectorTensorQ[inds___]]:=Block[{TVect},
(Evaluate[lhs/.(TensDummy->PatternTest[Pattern[TensDummy,_],VectorTensorQ])]:>(rhs))/.(TensDummy->TVect)
]


PatternTensorLeftTensor[(lhs_:>rhs_),TensDummy_?TensorTensorQ[inds___]]:=Block[{TTens},
(Evaluate[lhs/.(TensDummy->PatternTest[Pattern[TensDummy,_],TensorTensorQ])]:>(rhs))/.(TensDummy->TTens)
]


SetSlicing[g_?MetricQ,u_,normu_:-1,h_,cd_,{cdpost_String,cdpre_String},SpaceTimeType_?SpaceTimeQ]:=
Module[{n,p,q,DummyS,DummyV,DummyT,ui,indsdimminustwo,indsdim,dim,prot},Catch[

If[Not@MetricQ[g],Throw[Message[SetSlicing::noambientmetric,g]]];

If[SignDetOfMetric@g=!=-1,Throw[Message[SetSlicing::signature,g]]];

With[{Manifold=ManifoldOfCovD@CovDOfMetric[g],CD=CovDOfMetric[g]},

dim=DimOfManifold[Manifold];

If[NumericQ[dim]&&dim<2,Throw[Message[SetSlicing::dimension,dim]]];

If[BianchiBool[SpaceTimeType],
If[Not@NumericQ[dim]||dim=!=4,Print["** Warning: You define a Bianchi spacetime for a dimension which is not 4. The Ellis&MacCallum-1969 parameterization of the constants of structure as \!\(\*SubscriptBox[SuperscriptBox[\(C\), \(i\)], \(jk\)]\)=\!\(\*SubscriptBox[\(\[Epsilon]\), \(jkl\)]\)\!\(\*SuperscriptBox[\(n\), \(li\)]\) + \!\(\*SubscriptBox[SuperscriptBox[\(\[Delta]\), \(i\)], \([j\)]\) \!\(\*SubscriptBox[\(a\), \(\(k\)\(]\)\)]\) is not possible. The option $OpenConstantsOfStructure is set to False. **"];
$OpenConstantsOfStructure=False;
]];

With[{ind1=DummyIn[Tangent[Manifold]],
ind2=DummyIn[Tangent[Manifold]],
ind3=DummyIn[Tangent[Manifold]],
ind4=DummyIn[Tangent[Manifold]],
ind5=DummyIn[Tangent[Manifold]],
ind6=DummyIn[Tangent[Manifold]],
ind7=DummyIn[Tangent[Manifold]],
i1=DummyIn[Tangent[Manifold]],
i2=DummyIn[Tangent[Manifold]],
i3=DummyIn[Tangent[Manifold]],
i4=DummyIn[Tangent[Manifold]],
i5=DummyIn[Tangent[Manifold]],
dummy=DummyIn[Tangent[Manifold]],
ah=a[h],Hh=H[h],
avh=av[h],nth=nt[h],
CSh=CS[h],Connectionh=Connection[h]},

If[Not@DefTensorQ[u],DefTensor[u[ind1],{Manifold},PrintAs->"\!\("<>ToString[u]<>"\&-\)"];,
If[Length[SlotsOfTensor[u]]=!=1||VBundleOfIndex@DummyIn[First@SlotsOfTensor@u]=!=Tangent[Manifold],Throw[Message[SetSlicing::invalidnormalvector,u]]];

];

Off[DefMetric::old];(* Annoying message turned off*)
(* contribution of Obinna. But removed by Cyril *)
(*DefMetric[1,h[-ind1,-ind2],cd,{cdpost,cdpre},InducedFrom\[Rule]{g,u},If[$ConformalTime,If[FlatSpaceBool[SpaceTimeType],PrintAs\[Rule]"\[Delta]",PrintAs\[Rule]"\!\("<>ToString[h]<>"\&-\)"],PrintAs\[Rule]"\!\("<>ToString[h]<>"\&-\)"]];*)

DefMetric[1,h[-ind1,-ind2],cd,{cdpost,cdpre},InducedFrom->{g,u},PrintAs->"\!\("<>ToString[h]<>"\&-\)"];

(* Maybe not so useful...*)
AutomaticRules[h,BuildRule[{g[ind1,ind2]h[-ind2,-ind3],h[ind1,-ind3]}]];
AutomaticRules[h,BuildRule[{g[-ind1,-ind2]h[ind2,ind3],h[-ind1,ind3]}]];

On[DefMetric::old];

(* We remove automatic Leibniz rule when there is a Scalar Head. This is to ensure that the induced derivative does not spoil an 'InducedDecomposition'*)
prot=Unprotect[cd];
cd[a_][Scalar[expr_]]=.;
Protect[prot];

$Rulecdh[h1_]:={h1[-a_,b_]cd[a_][expr1_]:>cd[b][expr1],
h1[a_,b_]cd[-a_][expr1_]:>cd[b][expr1],
h1[b_,-a_]cd[a_][expr1_]:>cd[b][expr1],
h1[b_,a_]cd[-a_][expr1_]:>cd[b][expr1],
h1[-a_,b_]cd[c_]@cd[a_][expr1_]:>cd[c]@cd[b][expr1],
h1[a_,b_]cd[c_]@cd[-a_][expr1_]:>cd[c]@cd[b][expr1],
h1[b_,-a_]cd[c_]@cd[a_][expr1_]:>cd[c]@cd[b][expr1],
h1[b_,a_]cd[c_]@cd[-a_][expr1_]:>cd[c]@cd[b][expr1]};

$RulesVanishingBackgroundFields[h1_]:={};

(** All of this is now included in xAct 1.0.5 so we comment these lines **)
(** I keep it in comments to remember it**)
(*cd[i1_][x_ y_]=.;
cd[_?GIndexQ][expr_]=.;*)
(* and we redefine them *)
(*cd[i1_][x_  g[ind1_,ind2_]]:=h[ind1,ind2] cd[i1][x](*/;OrthogonalToVectorQ[u][x]*);
cd[i1_][x_  delta[ind1_,ind2_]]:=h[ind1,ind2] cd[i1][x](*/;OrthogonalToVectorQ[u][x]*);
cd[i1_][g[ind1_,ind2_]]:=0;*)

(*cd[i1_][x_  y_]:=cd[i1][x]y+cd[i1][y]x/;And[OrthogonalToVectorQ[u][x],OrthogonalToVectorQ[u][y]];

cd[_?GIndexQ][expr_]:=Throw@Message[Validate::error,"Induced derivative acting on non-projected expression."]/;Not@OrthogonalToVectorQ[u][expr];

cd[i1_][x_  y_]:=cd[i1][Expand@GradNormalToExtrinsicK@Expand[InducedDecomposition[x,{h,u}]InducedDecomposition[y,{h,u}]]]/;Head@x=!=g&&Head@y=!=g&&OrthogonalToVectorQ[u][x y];*)
(********)
(* ******* End of Patch *)


SpaceType[h]^=SpaceTimeType;

(* The default positionof indices for the extrinsic curvature and the acceleration is down*)
(SlotsOfTensor[#]^:={-Tangent[Manifold],-Tangent[Manifold]})&/@{ExtrinsicK[h]};
(SlotsOfTensor[#]^:={-Tangent[Manifold]})&/@{Acceleration[u]};

(* The acceleration should vanish for homogeneous spacetimes.*)
Acceleration[u][ind1_] = 0;
AutomaticRules[u,BuildRule[{u[ind1]u[-ind1],normu}]];
AutomaticRules[u,BuildRule[{u[-ind1]g[ind1,ind2],u[ind2]}]];
AutomaticRules[u,BuildRule[{g[ind1,ind2]u[-ind2]u[-ind1],normu}]];
If[IntegerQ@dim&&dim>=2,
indsdim=GetIndicesOfVBundle[Tangent@Manifold,dim,{ind5}];
AutomaticRules[epsilon[g],BuildRule[Evaluate[{epsilon[g]@@indsdim u[-indsdim[[1]]]h[-indsdim[[2]],ind5],ReplaceIndex[Evaluate[epsilon[g]@@indsdim],indsdim[[2]]-> ind5] u[-indsdim[[1]]]}]]];
];


(* We define the normal vector to constant time hypersurfaces (it should be equal to the normal vector on the background spacetime.). This defines also the perturbation of the normal vector, and the fields related to it (lapse scalar and shift vector in ADM formalism). *)
DefNormalFields[h];

(* We define a[h] the scale factor associated to h, and the Hubble factor related.*)
(DefProjectedTensor[#[[1]],h,SpaceTimesOfDefinition->{"Background"},TensorProperties->{"Traceless","Transverse","SymmetricTensor"},PrintAs->#[[2]]])&/@$ListFieldsBackgroundOnly[h];

ah::usage=a::usage;
Hh::usage=H::usage;

(* Patch to have a nice output for the Hubble factor *)
Hh/:PrintAs[Hh]=.;
PrintAs[Hh]^:=If[$ConformalTime,"\[ScriptCapitalH]","H"];

Unprotect[NoScalar];
(* We know it is safe to remove the scalar heads on scale factors and Hubble factors, because we shall never replace a Hubble factor nor a scale factor by something else which is not a pure scalar field...*)
NoScalar[Power[Scalar[Hh[LI[0],LI[0]]],n_Integer]]:=Power[Hh[LI[0],LI[0]],n];
NoScalar[Power[Scalar[ah[LI[0],LI[0]]],n_Integer]]:=Power[ah[LI[0],LI[0]],n];
Protect[NoScalar];

If[SpaceTimeType==="Minkowski",ah[LI[0],LI[0]]=1;Hh[LI[0],LI[0]]=0;,
ah[LI[0],LI[1]]:=Hh[LI[0],LI[0]]ah[LI[0],LI[0]];
ah[LI[0],LI[q_?((IntegerQ[#]&&#>=2)&)]]:=NoScalar@org@Nest[LieD[u[ind1]][#]&,ah[LI[0],LI[0]],q];
DefConformalMetric[g,ah];
];

(* Obvious. SHould be automatic in xPert*)
cd[_][$PerturbationParameter]=0;

(* The vector used for the background slicing is indeed background so should not be perturbed *)
(* However we want to be able to perturb the normal vector to constan time hypersurfaces.*)
(* In order to do so, we will use the vector \[ScriptCapitalN][h][\[Mu]] which is equal to the background normal vector n^\[Mu].
Howeverm it is ano rule attached to it, so it can be perturbed. Its perturbation is then d\[ScriptCapitalN][LI[p],\[Mu]] *)
(* See the section devoted to the perturbation of the normal vector for more explanations*)
Perturbation[u[ind1_],n_]^:=0/;n>=1;


(* Lie Derivatives of metric with up indices should be automatic*)
AutomaticRules[g,BuildRule[Evaluate[{LieD[u[dummy]][g[ind1,ind2]],LieD[u[dummy]][g[ind1,ind2]]//MetricToProjector[#,h]&}],MetricOn->None]];
AutomaticRules[g,BuildRule[Evaluate[{LieD[u[dummy]][g[-ind1,-ind2]],LieD[u[dummy]][g[-ind1,-ind2]]//MetricToProjector[#,h]&}],MetricOn->None]];


DefTensor[Connectionh[-ind1,-ind2,-ind3],{Manifold},Antisymmetric[{-ind1,-ind3}],OrthogonalTo->{u[ind1],u[ind2],u[ind3]},ProjectedWith->{h[ind1,-ind4],h[ind2,-ind5],h[ind3,-ind6]},PrintAs->StringJoin["\[CapitalGamma]",ToString[h]]] ;

DefTensor[CSh[-ind1,-ind2,-ind3],{Manifold},Antisymmetric[{-ind2,-ind3}],OrthogonalTo->{u[ind1],u[ind2],u[ind3]},ProjectedWith->{h[ind1,-ind4],h[ind2,-ind5],h[ind3,-ind6]},PrintAs->StringJoin["C"(*,ToString[h]*)]] ;

If[Not[FlatSpaceBool[SpaceTimeType]&&dim===4&&BianchiBool[SpaceTimeType]],

DefTensor[nth[-ind1,-ind2],{Manifold},Symmetric[{-ind1,-ind2}],OrthogonalTo->{u[ind1],u[ind2]},ProjectedWith->{h[ind1,-ind4],h[ind2,-ind5]},PrintAs->StringJoin["N"(*,ToString[h]*)]] ;

DefTensor[avh[-ind1],{Manifold},OrthogonalTo->{u[ind1]},ProjectedWith->{h[ind1,-ind4]},PrintAs->StringJoin["A"(*,ToString[h]*)]] ;

If[SpaceTimeType==="BianchiA",
avh[ind1_]=0;,
AutomaticRules[avh,BuildRule[Evaluate[{avh[ind1]nth[-ind1,-ind2],0}]]];
];
];

(* Properties of constancy of these constants...*)
(* So far the only way we found is to be exhaustive *)

(* We first distinguish the case of Flat Friedmann and Curved Friedmann where the connections have a simple expression*)
Which[
FlatSpaceBool[SpaceTimeType],
CSh[ind1_,ind2_,ind3_]=0;
Riemann[cd][i1_,i2_,i3_,i4_]=0;
Ricci[cd][i1_,i2_]=0;
RicciScalar[cd][]=0;,

SpaceTimeType==="FLCurved",
DefTensor[\[ScriptK][h][],{Manifold},PrintAs->StringJoin["\[ScriptK]"(*,ToString[h]*)]];

CD[ind1_][\[ScriptK][h][]]:=0(*-1*(-2/3)*u[ind1]\[ScriptK][h][]CD[-ind2][u[ind2]]*); (* This is because Subscript[\[ScriptCapitalL], u](\[ScriptK]^(1/2)Subscript[\[Epsilon]^\[Mu], \[Nu]\[Sigma]]) must be zero since Subscript[\[ScriptCapitalL], u](Subscript[C^\[Mu], \[Nu]\[Sigma]])=0*)
(* Given that we first work on the conformal metric where the trace of the extrinsic curvature is 0, then this point does not matter at all for us.*)


(*In case the user uses a Lie derivative on spatila curvature instead of a covariant derivative, we also specify that it should vanish. This was suggested by Adam Solomon.*)
Evaluate[\[ScriptK][h]]/:LieD[u[ind1_]][\[ScriptK][h][]]:=0;(*Module[{dum},dum=DummyIn[Tangent[Manifold]]u[-dum]CD[dum][\[ScriptK][h][]]]*)
cd[ind1_][\[ScriptK][h][]]=0;
If[dim===4,CSh[ind1_,ind2_,ind3_]:=2Sqrt[\[ScriptK][h][]]epsilon[h][ind1,ind2,ind3]];
,


(BianchiBool@SpaceTimeType)&&(Not@FlatSpaceBool@SpaceTimeType),
(* Lie derivatives of the constants of structure in general Bianchi cases It is 0 when indices are in position Subscript[(C^i), jk] .*)
(* We are exhaustive. Maybe there is  simpler way to write it.*)
AutomaticRules[CSh,BuildRule[Evaluate[{LieD[u[ind1]][CSh[ind2,-ind3,-ind4]],0}],MetricOn->None]];

AutomaticRules[CSh,BuildRule[Evaluate[{LieD[u[ind1]][CSh[-ind2,-ind3,-ind4]],ToCanonical@ContractMetric[LieD[u[ind1]][g[-ind2,-ind5]]CSh[ind5,-ind3,-ind4]//MetricToProjector[#,h]&]}],MetricOn->None]];

AutomaticRules[CSh,BuildRule[Evaluate[{LieD[u[ind1]][CSh[ind2,-ind3,ind4]],ToCanonical@ContractMetric[LieD[u[ind1]][g[ind4,ind6]]CSh[ind2,-ind3,-ind6]//MetricToProjector[#,h]&]}],MetricOn->None]];

AutomaticRules[CSh,BuildRule[Evaluate[{LieD[u[ind1]][CSh[ind2,ind3,-ind4]],ToCanonical@ContractMetric[LieD[u[ind1]][g[ind3,ind5]]CSh[ind2,-ind5,-ind4]//MetricToProjector[#,h]&]}],MetricOn->None]];

AutomaticRules[CSh,BuildRule[Evaluate[{LieD[u[ind1]][CSh[-ind2,ind3,-ind4]],ToCanonical@ContractMetric[LieD[u[ind1]][g[ind3,ind5]]g[-ind6,-ind2]CSh[ind6,-ind5,-ind4]+g[ind3,ind5]LieD[u[ind1]][g[-ind6,-ind2]]CSh[ind6,-ind5,-ind4]//MetricToProjector[#,h]&]}],MetricOn->None]];

AutomaticRules[CSh,BuildRule[Evaluate[{LieD[u[ind1]][CSh[-ind2,-ind4,ind3]],ToCanonical@ContractMetric[LieD[u[ind1]][g[ind3,ind5]]g[-ind6,-ind2]CSh[ind6,-ind4,-ind5]+g[ind3,ind5]LieD[u[ind1]][g[-ind6,-ind2]]CSh[ind6,-ind4,-ind5]//MetricToProjector[#,h]&]}],MetricOn->None]];

AutomaticRules[CSh,BuildRule[Evaluate[{LieD[u[ind1]][CSh[-ind2,ind3,ind4]],
ToCanonical@ContractMetric[LieD[u[ind1]][g[ind3,ind5]]g[-ind6,-ind2]g[ind4,ind7]CSh[ind6,-ind5,-ind7]
+g[ind3,ind5]LieD[u[ind1]][g[-ind6,-ind2]]g[ind4,ind7]CSh[ind6,-ind5,-ind7]
+g[ind3,ind5]g[-ind6,-ind2]LieD[u[ind1]][g[ind4,ind7]]CSh[ind6,-ind5,-ind7]//MetricToProjector[#,h]&]}],MetricOn->None]];

AutomaticRules[CSh,BuildRule[Evaluate[{LieD[u[ind1]][CSh[ind2,ind3,ind4]],
ToCanonical@ContractMetric[LieD[u[ind1]][g[ind3,ind5]]g[ind4,ind6]CSh[ind2,-ind5,-ind6]
+g[ind3,ind5]LieD[u[ind1]][g[ind4,ind6]]CSh[ind2,-ind5,-ind6]//MetricToProjector[#,h]&]}],MetricOn->None]];
];


(* -------- *)
(* We now define the connection and the Riemann tensor of the induced metric. See draft *)
IndexSet[Connectionh[i1_,i2_,i3_],1/2(CSh[i2,i1,i3]-CSh[i1,i3,i2]+CSh[i3,i1,i2])];

If[SpaceTimeType==="FLCurved",
IndexSet[Riemann[cd][i1_,i2_,i3_,i4_],\[ScriptK][h][](h[i1,i3]h[i2,i4]-h[i1,i4]h[i2,i3])];
IndexSet[Ricci[cd][i1_,i2_],\[ScriptK][h][]h[i1,i2](h[i3,-i3]-1)];
IndexSet[RicciScalar[cd][],\[ScriptK][h][]h[i1,-i1](h[i3,-i3]-1)];,

(* Rieman, Ricci and RicciScalar tensors. We build rules to express them in terms of constants of structure. *)

ConstantsOfStructureRules[h]=Flatten@Join[BuildRule[Evaluate[{Riemann[cd][i1,i2,i3,i4], ToCanonical@ContractMetric[(-Connectionh[i5,i3,i2]Connectionh[i1,i4,-i5]+Connectionh[i5,i4,i2]Connectionh[i1,i3,-i5]-CSh[i5,i3,i4]Connectionh[i1,-i5,i2])(*/.ConstantsDecompositionRules[h]*)]}]],
BuildRule[Evaluate[{Ricci[cd][i1,i3], ToCanonical@ContractMetric[(h[-i2,-i4](-Connectionh[i5,i3,i2]Connectionh[i1,i4,-i5]+Connectionh[i5,i4,i2]Connectionh[i1,i3,-i5]-CSh[i5,i3,i4]Connectionh[i1,-i5,i2]))(*/.ConstantsDecompositionRules[h]*)]}]],
BuildRule[Evaluate[{RicciScalar[cd][], ToCanonical@ContractMetric[(h[-i2,-i4]h[-i1,-i3](-Connectionh[i5,i3,i2]Connectionh[i1,i4,-i5]+Connectionh[i5,i4,i2]Connectionh[i1,i3,-i5]-CSh[i5,i3,i4]Connectionh[i1,-i5,i2]))(*/.ConstantsDecompositionRules[h]*)]}]]
];

];


(* We replace the extrinsic curvature with an adhoc tensor which can have derivatives label indices*)
DefProjectedTensor[K[h][-ind1,-ind2],h,TensorProperties->{"SymmetricTensor","Traceless"},SpaceTimesOfDefinition->{"Background"},PrintAs->"K"];
ExtrinsicK[h][i1_,i2_]:=If[$ConformalTime,1,a[h][]]*K[h][LI[0],LI[0],i1,i2];
(*ExtrinsicK[h][LI[0],LI[q_],i1_,i2_]:=K[h][LI[0],LI[q],i1,i2];*)
(*PrintAs[Kh]^=printextrinsicK;*)


If[NoAnisotropyBool[SpaceTimeType],
ExtrinsicK[h][i1_,i2_]=0;Evaluate[K[h]][LI[0],LI[0],i1_,i2_]=0;];


(* Commutation of Lie and Partial derivatives. Very important here. See draft for the formula which is implemented...*)
Unprotect[LieD];
(* A revoir... This is to commute Lie and induced derivative when the indices are up. Basically it just ensure that indices should be put down.*)

(* The code below was reported not to work when there was a contraction of indices*)
(* Thank to Enrico Pajer for finding this.*)

(* This was the first patch but it did not work*)
(*LieD[u[ind1_]][cd[ind2_][expr1_]]:=LieD[u[ind1]][IndicesDown[cd[ind2][expr1] ] ]/;Length[IndicesOf[AIndex,Up][cd[ind2][expr1]]]=!=0&&OrthogonalToVectorQ[u][expr1]&&Abs[u[ind1]u[-ind1]]===1;*)

(* Instead we define a new rule for this very special case of a Laplacian...*)
(* This sound like a stupid pathc rather than something deeply thought*)
(* We basically introduce the metric by hand in order to lower the up index in the Laplacian. Then it goes through the Lie Derivative*)
(*  So we treat this special case by hand.*)

LieD[u[ind1_]][cd[ind2_][cd[-ind2_][expr1_]]]:=Module[{dum},dum=DummyIn[Tangent[Manifold]];
ContractMetric[LieD[u[ind1]][g[dum,ind2]]cd[-dum][cd[-ind2][expr1]]
+g[dum,ind2]LieD[u[ind1]][cd[-dum][cd[-ind2][expr1]]]]
];

LieD[u[ind1_]][cd[-ind2_][cd[ind2_][expr1_]]]:=Module[{dum},dum=DummyIn[Tangent[Manifold]];
ContractMetric[LieD[u[ind1]][g[dum,ind2]]cd[-dum][cd[-ind2][expr1]]
+g[dum,ind2]LieD[u[ind1]][cd[-dum][cd[-ind2][expr1]]]]
];


LieD[u[ind1_]][cd[ind2_][expr1_]]:=LieD[u[ind1]][IndicesDown[cd[ind2][expr1] ] ]/;Length[IndicesOf[Free,Up][cd[ind2][expr1]]]=!=0&&OrthogonalToVectorQ[u][expr1]&&Abs[u[ind1]u[-ind1]]===1;

LieD[u[ind1_]][cd[ind2_?DownIndexQ][expr1_]]:=Module[{dum},dum=DummyIn[Tangent[Manifold]];
With[{frees=FindFreeIndices[expr1]},
ToCanonical[
(cd[ind2][LieD[u[ind1]][expr1]]
+$ExtrinsicKSign *Plus@@(
(-cd[#][ExtrinsicK[h][ind2,dum]]ReplaceIndex[expr1,#->-dum]
+cd[dum][ExtrinsicK[h][#,ind2]]ReplaceIndex[expr1,#->-dum]
-cd[ind2][ExtrinsicK[h][dum,#]]ReplaceIndex[expr1,#->-dum])&/@frees)
)
,UseMetricOnVBundle->None]]]/;Length[IndicesOf[Free,Up][expr1]]===0&&OrthogonalToVectorQ[u][expr1]&&Abs[u[ind1]u[-ind1]]===1;
Protect[LieD];


If[(BianchiBool@SpaceTimeType)&&(Not@FlatSpaceBool@SpaceTimeType),
(* Automatic rules for the gradients of the constants of structure. I should define one function only to define these rules since it is also called in the DefSTFTensor... Here I keep the way it was implemented in the first version*)
AutomaticRules[CSh,
BuildRule[Evaluate[{cd[ind1][CSh[ind2,ind3,ind4]],ToCanonical[
(* Minus the connection for down indices *)
-Plus@@((Connectionh[dummy,ind1,#] ReplaceIndex[CSh[ind2,ind3,ind4],#->-dummy])&/@{ind2,ind3,ind4})]
}]]];

(* From this, the n+1 splitting of Covariant derivatives is obvious. *)
AutomaticRules[CSh,
BuildRule[Evaluate[{CD[ind1][CSh[ind2,ind3,ind4]],
ToCanonical[ContractMetric[ToInducedDerivative[IndicesDown[CD[ind1][CSh[ind2,ind3,ind4]]],CD,cd]//GradNormalToExtrinsicK]]}]]];

(* TODO Finish correctly this part *)
(* It should work as it is, but it might not be optimal in terms of computatiomal time *)
(* When doing the 1+3 splitting of the covariant derivative acting on the 3D riemann tensor, there are Lie derivatives. *)
(* These are handled because the 3D Riemann tensor is eventually replaced by constants of structure and for the constants of structure we have defined what the Lie derivatives should give. *)
(* So all rules are there, but this might not be optimal in terms of compyutational time because the optimal way would be to rpecompute the Lie derivative on 3D riemann, without having to rely on constants of structure.*)

AutomaticRules[Evaluate[Riemann[cd]],
BuildRule[Evaluate[{cd[ind1][Riemann[cd][ind2,ind3,ind4,ind5]],ToCanonical[
(* Minus the connection for down indices *)
-Plus@@((Connectionh[dummy,ind1,#] ReplaceIndex[Evaluate[Riemann[cd][ind2,ind3,ind4,ind5]],#->-dummy])&/@{ind2,ind3,ind4,ind5})]
}]]];

AutomaticRules[Evaluate[Ricci[cd]],
BuildRule[Evaluate[{cd[ind1][Ricci[cd][ind2,ind3]],ToCanonical[
(* Minus the connection for down indices *)
-Plus@@((Connectionh[dummy,ind1,#] ReplaceIndex[Evaluate[Ricci[cd][ind2,ind3]],#->-dummy])&/@{ind2,ind3})]
}]]];

AutomaticRules[Evaluate[RicciScalar[cd]],
BuildRule[Evaluate[{cd[ind1][RicciScalar[cd][]],0}]]];


(* 3+1 splitting of CD on Riemann and Ricci. Should never be used in general unless the user wants to split CD[Riemann]. Very unlikely. *)
AutomaticRules[Evaluate[Riemann[cd]],
BuildRule[Evaluate[{CD[ind1][Riemann[cd][ind2,ind3,ind4,ind5]],
ToCanonical[ContractMetric[ToInducedDerivative[IndicesDown[CD[ind1][Riemann[cd][ind2,ind3,ind4,ind5]]],CD,cd]//GradNormalToExtrinsicK]]}]]];

AutomaticRules[Evaluate[Ricci[cd]],
BuildRule[Evaluate[{CD[ind1][Ricci[cd][ind2,ind3]],
ToCanonical[ContractMetric[ToInducedDerivative[IndicesDown[CD[ind1][Ricci[cd][ind2,ind3]]],CD,cd]//GradNormalToExtrinsicK]]}]]];

AutomaticRules[Evaluate[RicciScalar[cd]],
BuildRule[Evaluate[{CD[ind1][RicciScalar[cd][]],
ToCanonical[ContractMetric[ToInducedDerivative[IndicesDown[CD[ind1][RicciScalar[cd][]]],CD,cd]//GradNormalToExtrinsicK]]}]]];

(* We also need to treat the Lie Derivative which comes from thsi 3+1 splitting *)
(* Lie D with down indices is treated first*)
(* This is done by expressing the Riemann tensor in terms of constanst of structure. *)
AutomaticRules[Evaluate[Riemann[cd]],BuildRule[Evaluate[{LieD[u[ind1]][Riemann[cd][-ind2,-ind3,-ind4,-ind5]],ToCanonical@ContractMetric[LieD[u[ind1]][Riemann[cd][-ind2,-ind3,-ind4,-ind5]/.ConstantsOfStructureRules[h]]] }],MetricOn->None]];

AutomaticRules[Evaluate[Ricci[cd]],BuildRule[Evaluate[{LieD[u[ind1]][Ricci[cd][-ind2,-ind3]],ToCanonical@ContractMetric[LieD[u[ind1]][Ricci[cd][-ind2,-ind3]/.ConstantsOfStructureRules[h]]] }],MetricOn->None]];

AutomaticRules[Evaluate[RicciScalar[cd]],BuildRule[Evaluate[{LieD[u[ind1]][RicciScalar[cd][]],ToCanonical@ContractMetric[LieD[u[ind1]][RicciScalar[cd][]/.ConstantsOfStructureRules[h]]] }],MetricOn->None]];


(* Extension to any position of indices inside a Lie Derivative*)
Evaluate[Riemann[cd]]/:LieD[u[dummy_]][Evaluate[Riemann[cd]][indices1___,indup_?UpIndexQ,indices2___]]:=
Module[{dum},dum=DummyIn[Tangent[Manifold]];(*Print["** Warning: a Lie derivative with up indices has been converted (safely) to a Lie derivative with down indices **"];*)g[indup,dum]LieD[u[dummy]][Evaluate[Riemann[cd]][indices1,-dum,indices2]]+Identity[LieD[u[dummy]][g[indup,dum]]Evaluate[Riemann[cd]][indices1,-dum,indices2]]];

Evaluate[Ricci[cd]]/:LieD[u[dummy_]][Evaluate[Ricci[cd]][indices1___,indup_?UpIndexQ,indices2___]]:=
Module[{dum},dum=DummyIn[Tangent[Manifold]];(*Print["** Warning: a Lie derivative with up indices has been converted (safely) to a Lie derivative with down indices **"];*)g[indup,dum]LieD[u[dummy]][Evaluate[Ricci[cd]][indices1,-dum,indices2]]+Identity[LieD[u[dummy]][g[indup,dum]]Evaluate[Ricci[cd]][indices1,-dum,indices2]]];



];

(* TODO MAJOR PROBLEM HERE. IF the tensores are not Symmetric or not transverse or not traceless etc... the properties are also applied so this will give wrong results. TODO TODO. SO we should remove the option that these tensors are tranverse. I will patch that. In general users use only transverse vectors though.*)

(* We gather many many rules for the commutation of induced derivatives. These are in general made to make sure that the transverse conditions of vectors and tensors are used. It is also made to gather Laplacians.
*)
DummyS[Sym_]:=SymbolJoin[DS,Sym];
DummyV[Sym_]:=SymbolJoin[DV,Sym];
DummyT[Sym_]:=SymbolJoin[DT,Sym];

Block[{Print},
DefProjectedTensor[Evaluate[DummyS[h]][],h,SpaceTimesOfDefinition->{"Perturbed"}];
DefProjectedTensor[Evaluate[DummyV[h]][-ind1],h,SpaceTimesOfDefinition->{"Perturbed"},TensorProperties->{(*"Transverse"*)}];
DefProjectedTensor[Evaluate[DummyT[h]][-ind1,-ind2],h,TensorProperties->{"SymmetricTensor"(*,(*"Transverse",*)"Traceless"*)},SpaceTimesOfDefinition->{"Perturbed"}];
];
(* Here we consider a general tensor which is not necessarily traceless or transverse. But it is symmetric.*)

$CommutecdRules=Join[$CommutecdRules,Flatten@Join[
PatternLeft[#,{n,q}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[h]][LI[n],LI[q]]]&/@BuildRule[Evaluate[{cd[-ind1][cd[ind2][cd[ind3][cd[ind1][Evaluate[DummyS[h]][LI[n],LI[q]]]]]],org[CommuteCovDs[CommuteCovDs[cd[-ind1][cd[ind2][cd[ind3][cd[ind1][Evaluate[DummyS[h]][LI[n],LI[q]]]]]],cd,{ind2,-ind1}],cd,{ind3,-ind1}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[h]][LI[n],LI[q]]]&/@BuildRule[Evaluate[{cd[ind2][cd[-ind1][cd[ind3][cd[ind1][Evaluate[DummyS[h]][LI[n],LI[q]]]]]],org[CommuteCovDs[cd[ind2][cd[-ind1][cd[ind3][cd[ind1][Evaluate[DummyS[h]][LI[n],LI[q]]]]]],cd,{ind3,-ind1}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[h]][LI[n],LI[q]]]&/@BuildRule[Evaluate[{cd[-ind1][cd[ind1][cd[ind3][cd[ind2][Evaluate[DummyS[h]][LI[n],LI[q]]]]]],org[CommuteCovDs[CommuteCovDs[CommuteCovDs[CommuteCovDs[cd[-ind1][cd[ind1][cd[ind3][cd[ind2][Evaluate[DummyS[h]][LI[n],LI[q]]]]]],cd,{ind3,ind1}],cd,{ind3,-ind1}],cd,{ind2,ind1}],cd,{ind2,-ind1}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[h]][LI[n],LI[q]]]&/@BuildRule[Evaluate[{cd[-ind2][cd[-ind1][cd[ind2][Evaluate[DummyS[h]][LI[n],LI[q]]]]],org[CommuteCovDs[cd[-ind2][cd[-ind1][cd[ind2][Evaluate[DummyS[h]][LI[n],LI[q]]]]],cd,{-ind1,-ind2}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[h]][LI[n],LI[q]]]&/@BuildRule[Evaluate[{cd[-ind2][cd[ind2][cd[-ind1][Evaluate[DummyS[h]][LI[n],LI[q]]]]],org[CommuteCovDs[CommuteCovDs[cd[-ind2][cd[ind2][cd[-ind1][Evaluate[DummyS[h]][LI[n],LI[q]]]]],cd,{-ind1,ind2}],cd,{-ind1,-ind2}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[h]][LI[n],LI[q],ind2]]&/@BuildRule[Evaluate[{cd[-ind2][cd[-ind1][cd[ind1][Evaluate[DummyV[h]][LI[n],LI[q],ind2]]]],org[CommuteCovDs[CommuteCovDs[cd[-ind2][cd[-ind1][cd[ind1][Evaluate[DummyV[h]][LI[n],LI[q],ind2]]]],cd,{-ind1,-ind2}],cd,{ind1,-ind2}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[h]][LI[n],LI[q],ind2]]&/@BuildRule[Evaluate[{cd[-ind2][cd[-ind1][cd[ind1][Evaluate[DummyV[h]][LI[n],LI[q],ind2]]]],org[CommuteCovDs[CommuteCovDs[cd[-ind2][cd[-ind1][cd[ind1][Evaluate[DummyV[h]][LI[n],LI[q],ind2]]]],cd,{-ind1,-ind2}],cd,{ind1,-ind2}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[h]][LI[n],LI[q],ind3]]&/@BuildRule[Evaluate[{cd[-ind3][cd[ind1][cd[ind2][Evaluate[DummyV[h]][LI[n],LI[q],ind3]]]],org[CommuteCovDs[CommuteCovDs[cd[-ind3][cd[ind1][cd[ind2][Evaluate[DummyV[h]][LI[n],LI[q],ind3]]]],cd,{ind1,-ind3}],cd,{ind2,-ind3}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[h]][LI[n],LI[q],ind3]]&/@BuildRule[Evaluate[{cd[ind1][cd[-ind3][cd[ind2][Evaluate[DummyV[h]][LI[n],LI[q],ind3]]]],org[CommuteCovDs[cd[ind1][cd[-ind3][cd[ind2][Evaluate[DummyV[h]][LI[n],LI[q],ind3]]]],cd,{ind2,-ind3}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[h]][LI[n],LI[q],ind2]]&/@BuildRule[Evaluate[{cd[-ind3][cd[-ind1][cd[ind3][Evaluate[DummyV[h]][LI[n],LI[q],ind2]]]],org[CommuteCovDs[cd[-ind3][cd[-ind1][cd[ind3][Evaluate[DummyV[h]][LI[n],LI[q],ind2]]]],cd,{ind3,-ind1}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[h]][LI[n],LI[q],ind1]]&/@BuildRule[Evaluate[{cd[-ind1][cd[-ind2][Evaluate[DummyV[h]][LI[n],LI[q],ind1]]],org[CommuteCovDs[cd[-ind1][cd[-ind2][Evaluate[DummyV[h]][LI[n],LI[q],ind1]]],cd,{-ind2,-ind1}]]}]]),

PatternLeft[#,{n,q}]&/@(PatternTensorLeftTensor[#,Evaluate[DummyT[h]][LI[n],LI[q],ind2,ind3]]&/@BuildRule[Evaluate[{cd[-ind2][cd[-ind1][Evaluate[DummyT[h]][LI[n],LI[q],ind2,ind3]]],org[CommuteCovDs[cd[-ind2][cd[-ind1][Evaluate[DummyT[h]][LI[n],LI[q],ind2,ind3]]],cd,{-ind1,-ind2}]]}]])

]];

Block[{Print},
UndefTensor[Evaluate[DummyS[h]]];
UndefTensor[Evaluate[DummyV[h]]];
UndefTensor[Evaluate[DummyT[h]]];
];


(* If there is a up index on the Lie derivative of extrinsic curvature, it is converted automatically.
That way, the n+1 splitting of Cov Ds on extrinsic curvatrue is completely automatic. It is the only tensor for which it is completely automatic*)
If[Not@NoAnisotropyBool[SpaceTimeType],
Evaluate[K[h]]/:LieD[u[dummy_]][Evaluate[K[h]][LI[n_],LI[q_],indices1___,indup_?UpIndexQ,indices2___]]:=
Module[{dum},dum=DummyIn[Tangent[Manifold]];(*Print["** Warning: a Lie derivative with up indices has been converted (safely) to a Lie derivative with down indices **"];*)g[indup,dum]LieD[u[dummy]][Evaluate[K[h]][LI[n],LI[q],indices1,-dum,indices2]]+Identity[LieD[u[dummy]][g[indup,dum]]Evaluate[K[h]][LI[n],LI[q],indices1,-dum,indices2]]];
];

]
]
]
]

SetNumberOfArguments[SetSlicing,6]
Protect[SetSlicing];

(* Default value. It is overwritten if there is curvature *)
ConstantsOfStructureRules[_]:={};

ToConstantsOfStructure[h_?InducedMetricQ]:=Function[expr,expr/.ConstantsOfStructureRules[h]];



DefinedPerturbationParameter[x_]:=False;

DefMetricFields[g_?MetricQ,dg_,h_?InducedMetricQ,PerturbParameter_:\[Epsilon]]:=Module[{n},
With[{M=ManifoldOfCovD@CovDOfMetric@g,cd=CovDOfMetric@h},
Block[{\[Mu],\[Nu]},
{\[Mu],\[Nu]}=GetIndicesOfVBundle[Tangent@M,2];
If[Not[DefTensorQ[dg]],
DefMetricPerturbation[g, dg, Evaluate[If[DefinedPerturbationParameter[$PerturbationParameter],$PerturbationParameter,PerturbParameter]]];
DefinedPerturbationParameter[$PerturbationParameter]=True;
PrintAs[dg] ^="\[Delta]"<>ToString[g];
dg[LI[n_],\[Mu]_,\[Nu]_]:=0/;n>1&&BackgroundFieldMethod;
dg[\[Mu]_?AIndexQ,\[Nu]_?AIndexQ]:=dg[LI[0],\[Mu],\[Nu]];,
If[$DefInfoQ,Print["** Warning: Metric perturbation already defined. This cannot be redefined without undefining it. **"];];
];

(* Defining some tensors we shall need anyway *)
(DefProjectedTensor[#[[1]],h,TensorProperties->{"Traceless","Transverse","SymmetricTensor"},SpaceTimesOfDefinition->{"Perturbed"},PrintAs->#[[2]]])&/@$ListFieldsPerturbedOnly[h];

Evaluate[\[Phi][h]]::usage=\[Phi]::usage;

Evaluate[Bs[h]]::usage=Bs::usage;
Evaluate[Bv[h]]::usage=Bv::usage;

Evaluate[\[Psi][h]]::usage=\[Psi]::usage;
Evaluate[Es[h]]::usage=Es::usage;
Evaluate[Ev[h]]::usage=Ev::usage;
Evaluate[Et[h]]::usage=Et::usage;

Evaluate[T[h]]::usage=T::usage;
Evaluate[Ls[h]]::usage=Ls::usage;
Evaluate[Lv[h]]::usage=Lv::usage;

(* If we want a nice output for the perturbation parameter,,*)
MakeBoxes[PerturbParameter,StandardForm]:=
StyleBox[ToString[$PerturbationParameter],FontColor->RGBColor[0.3, 0.8, 0.8]]

]
]
]

SetNumberOfArguments[DefMetricFields,{3,4}];
Protect[DefMetricFields];


DefMatterFields[uf_,duf_,h_?InducedMetricQ, PerturbParameter_:\[Epsilon]]:=Module[{n,ord},
With[{g=First@InducedFrom@h,u=Last@InducedFrom@h},
With[{M=ManifoldOfCovD@CovDOfMetric@g,cd=CovDOfMetric@h},
Block[{\[Mu],\[Nu]},
{\[Mu],\[Nu]}=GetIndicesOfVBundle[Tangent@M,2];


(DefProjectedTensor[#[[1]],h,TensorProperties->{"Traceless","Transverse","SymmetricTensor"},SpaceTimesOfDefinition->{"Background","Perturbed"},PrintAs->#[[2]]])&/@$ListFieldsBackgroundAndPerturbed[h,uf];

Evaluate[\[Rho][uf]]::usage=\[Rho]::usage;
Evaluate[P[uf]]::usage=P::usage;

Evaluate[Vspat[h,uf]]::usage=Vspat::usage;
Evaluate[V0[h,uf]]::usage=V0::usage;
Evaluate[Vs[h,uf]]::usage=Vs::usage;
Evaluate[Vv[h,uf]]::usage=Vv::usage;

If[Not[DefTensorQ[uf]],
DefTensor[uf[\[Mu]],M];,
If[$DefInfoQ,Print["** Warning: Tensor ",uf , "is already defined. It cannot be redefined without undefining it. **"];];
];

If[Not[DefTensorQ[duf]],
DefTensorPerturbation[duf[LI[ord],\[Mu]],uf[\[Mu]],M,PrintAs->"\[Delta]"<>PrintAs[uf]];,
If[$DefInfoQ,Print["** Warning: Tensor ",duf , "is already defined. It cannot be redefined without undefining it. **"];];
];

If[Not@DefinedPerturbationParameter[$PerturbationParameter],$PerturbationParameter=PerturbParameter;DefinedPerturbationParameter[$PerturbationParameter]=True];

duf[\[Mu]_?AIndexQ]:=duf[LI[0],\[Mu]];

]
]
]
]

SetNumberOfArguments[DefMatterFields,{3,4}]
Protect[DefMatterFields];


BV1:=Boole@$FirstOrderVectorPerturbations;
BT1:=Boole@$FirstOrderTensorPerturbations;

SplitMetric[g_?MetricQ,dg_,h_?InducedMetricQ,gauge_?GaugeQ]:=Module[{n,ind1,ind2},
With[{u=Last@InducedFrom@h,M=ManifoldOfCovD@CovDOfMetric@g,cd=CovDOfMetric@h,ShBool=Evaluate@BianchiBool[SpaceType[h]]},

{ind1,ind2}=GetIndicesOfVBundle[Tangent@M,2];
With[{i1=ind1,i2=ind2},

If[Not[DefTensorQ[Et[h]]]||Not[DefTensorQ[dg]],Print["** Warning: The perturbed metric, or the fields required to parameterize its splitting were not previously defined **"];
Print["** DefMetricFields is called to build the perturbation of the metric and the fields needed for future splitting **"];
DefMetricFields[g, dg,h]];

(* First the rules at first order. We separate them from higher order rules because the user can choose not to use vectors and tensor at first order.*)


(** JMM explained that I should use a combination of Module and With, and also build first the lhs and rhs that I gather in {lhs,rhs}, and apply the Rule head on it with Rule@@ **)
(** This avoids the use  of PatternLeft and ensure that the indices inside the rule belong to the manifold and thus this avoids conflict with SCreenDollarIndices**)

Join[(*PatternLeft[#,{i1,i2}]&/@*){Rule@@Switch[gauge ,
"NewtonGauge",
{dg[LI[1],i1_,i2_],-u[i1]u[i2]2\[Phi][h][LI[1],LI[0]]
-2\[Psi][h][LI[1],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+BT1 2Et[h][LI[1],LI[0],i1,i2]
-BV1(u[i1]Bv[h][LI[1],LI[0],i2]+u[i2]Bv[h][LI[1],LI[0],i1])},
"FluidComovingGauge",
{dg[LI[1],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[1],LI[0]]
-2\[Psi][h][LI[1],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+BT1 2Et[h][LI[1],LI[0],i1,i2]
-BV1(u[i1]Bv[h][LI[1],LI[0],i2]+u[i2]Bv[h][LI[1],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[1],LI[0]]+u[i2]cd[i1]@Bs[h][LI[1],LI[0]])]},
"ScalarFieldComovingGauge",
{dg[LI[1],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[1],LI[0]]
-2\[Psi][h][LI[1],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+BT1 2Et[h][LI[1],LI[0],i1,i2]
-BV1(u[i1]Bv[h][LI[1],LI[0],i2]+u[i2]Bv[h][LI[1],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[1],LI[0]]+u[i2]cd[i1]@Bs[h][LI[1],LI[0]])]},
"SynchronousGauge",
{dg[LI[1],i1_,i2_],Identity[
-2\[Psi][h][LI[1],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+BT1 2Et[h][LI[1],LI[0],i1,i2]
+2cd[i1]@cd[i2]@Es[h][LI[1],LI[0]]
+BV1(cd[i1]@Ev[h][LI[1],LI[0],i2]+cd[i2]@Ev[h][LI[1],LI[0],i1])]},
"FlatGauge",
{dg[LI[1],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[1],LI[0]]
+BT1 2Et[h][LI[1],LI[0],i1,i2]
-BV1(u[i1]Bv[h][LI[1],LI[0],i2]+u[i2]Bv[h][LI[1],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[1],LI[0]]+u[i2]cd[i1]@Bs[h][LI[1],LI[0]])]},
"AnyGauge",
{dg[LI[1],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[1],LI[0]]
-2\[Psi][h][LI[1],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+BT1 2Et[h][LI[1],LI[0],i1,i2]
+2cd[i1]@cd[i2]@Es[h][LI[1],LI[0]]
+BV1(cd[i1]@Ev[h][LI[1],LI[0],i2]+cd[i2]@Ev[h][LI[1],LI[0],i1])
-BV1(u[i1]Bv[h][LI[1],LI[0],i2]+u[i2]Bv[h][LI[1],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[1],LI[0]]+u[i2]cd[i1]@Bs[h][LI[1],LI[0]])]},
"IsoDensityGauge",
{dg[LI[1],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[1],LI[0]]
-2\[Psi][h][LI[1],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+BT1 2Et[h][LI[1],LI[0],i1,i2]
-BV1(u[i1]Bv[h][LI[1],LI[0],i2]+u[i2]Bv[h][LI[1],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[1],LI[0]]+u[i2]cd[i1]@Bs[h][LI[1],LI[0]])]}
]},

(* And then the rules at order larger that 1*)
(*PatternLeft[#,{i1,i2}]&/@*){Rule@@Switch[gauge ,
"NewtonGauge",
{dg[LI[n_?(#>=2&)],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[n],LI[0]]
-2\[Psi][h][LI[n],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+2Et[h][LI[n],LI[0],i1,i2]
-(u[i1]Bv[h][LI[n],LI[0],i2]+u[i2]Bv[h][LI[n],LI[0],i1])]},
"FluidComovingGauge",
{dg[LI[n_?(#>=2&)],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[n],LI[0]]
-2\[Psi][h][LI[n],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+2Et[h][LI[n],LI[0],i1,i2]
-(u[i1]Bv[h][LI[n],LI[0],i2]+u[i2]Bv[h][LI[n],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[n],LI[0]]+u[i2]cd[i1]@Bs[h][LI[n],LI[0]])]},
"ScalarFieldComovingGauge",
{dg[LI[n_?(#>=2&)],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[n],LI[0]]
-2\[Psi][h][LI[n],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+2Et[h][LI[n],LI[0],i1,i2]
-(u[i1]Bv[h][LI[n],LI[0],i2]+u[i2]Bv[h][LI[n],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[n],LI[0]]+u[i2]cd[i1]@Bs[h][LI[n],LI[0]])]},
"SynchronousGauge",
{dg[LI[n_?(#>=2&)],i1_,i2_],Identity[
-2\[Psi][h][LI[n],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+2Et[h][LI[n],LI[0],i1,i2]
+2cd[i1]@cd[i2]@Es[h][LI[n],LI[0]]
+cd[i1]@Ev[h][LI[n],LI[0],i2]+cd[i2]@Ev[h][LI[n],LI[0],i1]]},
"FlatGauge",
{dg[LI[n_?(#>=2&)],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[n],LI[0]]
+2Et[h][LI[n],LI[0],i1,i2]
-(u[i1]Bv[h][LI[n],LI[0],i2]+u[i2]Bv[h][LI[n],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[n],LI[0]]+u[i2]cd[i1]@Bs[h][LI[n],LI[0]])]},
"AnyGauge",
{dg[LI[n_?(#>=2&)],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[n],LI[0]]
-2\[Psi][h][LI[n],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+2Et[h][LI[n],LI[0],i1,i2]
+2cd[i1]@cd[i2]@Es[h][LI[n],LI[0]]
+cd[i1]@Ev[h][LI[n],LI[0],i2]+cd[i2]@Ev[h][LI[n],LI[0],i1]
-(u[i1]Bv[h][LI[n],LI[0],i2]+u[i2]Bv[h][LI[n],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[n],LI[0]]+u[i2]cd[i1]@Bs[h][LI[n],LI[0]])]},
"IsoDensityGauge",
{dg[LI[n_?(#>=2&)],i1_,i2_],Identity[-u[i1]u[i2]2\[Phi][h][LI[n],LI[0]]
-2\[Psi][h][LI[n],LI[0]](h[i1,i2]+If[ShBool,K[h][LI[0],LI[0],i1,i2]/H[h][LI[0],LI[0]],0])
+2Et[h][LI[n],LI[0],i1,i2]
-(u[i1]Bv[h][LI[n],LI[0],i2]+u[i2]Bv[h][LI[n],LI[0],i1])
-(u[i1]cd[i2]@Bs[h][LI[n],LI[0]]+u[i2]cd[i1]@Bs[h][LI[n],LI[0]])]}]
}]
]
]
];
SetNumberOfArguments[SplitMetric,4]
Protect[SplitMetric];


IndicesDown[expr_]:= Fold[SeparateMetric[First@$Metrics][#1,#2]&,expr,Select[IndicesOf[Up][expr],Not@LIndexQ[#]&]]
IndicesUp[expr_]:= Fold[SeparateMetric[First@$Metrics][#1,#2]&,expr,Select[IndicesOf[Down][expr],Not@LIndexQ[#]&]]


IndicesDown[0]:=0 ;(* This is to avoid bugs...*)
IndicesUp[0]:=0 ;


ConformalMetricName[g_?MetricQ,1]:=g;
ConformalMetricName[g_?MetricQ,conffactor_]:=SymbolJoin[g,conffactor,2];


DefConformalMetric[g_?MetricQ,conffactor_]:=Module[{n,q},Catch@With[{M=ManifoldOfCovD@CovDOfMetric@g,CD=CovDOfMetric@g},
With[{i1=DummyIn[Tangent[M]],i2=DummyIn[Tangent[M]],sy1=SymbolOfCovD[CD][[1]],sy2=SymbolOfCovD[CD][[2]],metlist=Select[$Metrics,InducedFrom[#]===Null&]},


If[Not@DefTensorQ[conffactor],DefTensor[conffactor[LI[n],LI[q]],{M},PrintAs->ToString[conffactor]]];


Off[DefMetric::old];(* Annoying message turned off*)
If[Not[DefTensorQ[ConformalMetricName[g,conffactor]]],


DefMetric[-1,ConformalMetricName[g,conffactor][-i1,-i2],SymbolJoin[CD,conffactor,2],{":",StringJoin[sy2,ToString[conffactor],ToString[2]]},PrintAs->StringJoin["[",PrintAs[g],"\!\("<>PrintAs[conffactor]<>"\^2\)","]"(*ToString[conffactor],ToString[2]*)],ConformalTo->{g[-i1,-i2],conffactor[LI[0],LI[0]]^2}];

];
On[DefMetric::old];

(* We have to ensure it is safe*)
Perturbation[conffactor[LI[0],LI[q_],indices___],n_]^:=0/;n>=1;


Off[ConformalRules::unknown];
(* We use the error sent by ConformalRules to chekc whether or not the metric in the list metlist is conformallyr elated to the metric g.
If it is the case we enforce transitivity of the conformal relations.*)
If[Catch@ConformalRules[g,#]=!=Null,
SetConformalTo[SymbolJoin[g,conffactor,2][-i1,-i2], {#[-i1, -i2],ConformalFactor[g,#]* conffactor[LI[0],LI[0]]^2}]]&/@metlist;
On[ConformalRules::unknown];

]
]
]

SetNumberOfArguments[DefConformalMetric,2]
Protect[DefConformalMetric];


ConfHead[_,_][delta[\[Mu]_,\[Nu]_]]:=delta[\[Mu],\[Nu]](* Because I know that when there is delta function in an expression, it is always with one index up and one down...so this should be fine.*)

 
ConfHead[metric1_?MetricQ,metric2_?MetricQ][ConfHead[metric2_?MetricQ,metric1_?MetricQ][expr_]]:=expr (* Not necessary *)
ConfHead[metric1_?MetricQ,metric1_?MetricQ][expr_]:=expr
ConfHead[metric2_?MetricQ,metric3_?MetricQ][ConfHead[metric1_?MetricQ,metric2_?MetricQ][expr_]]:=ConfHead[metric1,metric3][expr]


(* Thanks to Jolyon Bloomfield and Leo Stein, the definition below should be much more general. *)
(* The main reason is that the delta tensor is greedy and wants to contract through expressions like
ConfHead[...][f[Scalar[phi[]]]].*)
ConfHead/:IsIndexOf[ConfHead[_,_][_],_,delta]:=False;


$BoolBasicConformalWeight=True;

WeightOfIndicesList[indices_List]:=With[{aindex=Select[indices,Not[LIndexQ[#]]&]},Length@Select[aindex,DownIndexQ]-Length@Select[aindex,UpIndexQ]]

(* Conformal weight of a tensor *)
ConformalWeight[tens_?xTensorQ]:=0;
ConformalWeight[tens_?xTensorQ[indices___]]:=ConformalWeight[tens]+WeightOfIndicesList[{indices}]

ConformalWeight[f_?ScalarFunctionQ]:=0;

MyChangeChristoffel[expr_,cd_,cd_]:=expr

MyChangeChristoffel[expr_,cd2list_List,cd1_]:=Fold[MyChangeChristoffel[#1,#2,cd1]&,expr,cd2list]

MyChangeChristoffel[expr_,cd2_,cd1_]:=With[{vb=Tangent[ManifoldOfCovD[cd1]]},With[{chr1=Head[(Christoffel[cd1])[DummyIn[vb],-DummyIn[vb],-DummyIn[vb]]],chr2=Head[(Christoffel[cd2])[DummyIn[vb],-DummyIn[vb],-DummyIn[vb]]],chr21=Head[(Christoffel@@Sort[{cd2,cd1}])[DummyIn[vb],-DummyIn[vb],-DummyIn[vb]]],sign=Order[cd2,cd1]},
expr/.chr2[i1_,i2_,i3_]:>chr1[i1,i2,i3]+sign*chr21[i1,i2,i3]
]
]


ExistInertHead[head_]:=Length@Cases[$InertHeads,head]>0

RulesConf[metric1_?MetricQ,metric2_?MetricQ]:=(
Module[{cd1,cd2,confa2,confa,M,res,inds},cd1=CovDOfMetric[metric1];cd2=CovDOfMetric[metric2];


confa2=ConformalFactor[metric2,metric1];
confa=Sqrt[ConformalFactor[metric2,metric1]]/.Sqrt[x_^n_?EvenQ]:>x^(n/2);
(*In principle this should give a or 1/a depending if we go from metric1 to metric2 or metric2 to metric1*)

M=ManifoldOfCovD[cd1];
inds=DummyIn/@Table[Tangent[M],{Range[4]}];
With[{i1=inds[[1]],i2=inds[[2]],i3=inds[[3]],i4=inds[[4]]},

(* Once confheads are put on expression (as a result of a formal conformal transformation) then we remove them by expressing what they mean in function of the original tensors and the scale factor *)
res=
{RuleDelayed@@Hold[ConfHead[metric1,metric2][(Riemann@cd1)[i1_,i2_,i3_,i4_]],confa^(WeightOfIndicesList[{i1,i2,i3,i4}]-2)(Riemann@cd2)[i1,i2,i3,i4]],
RuleDelayed@@Hold[ConfHead[metric1,metric2][(Ricci@cd1)[i1_,i2_]],confa^(WeightOfIndicesList[{i1,i2}]-2)(Ricci@cd2)[i1,i2]],
RuleDelayed@@Hold[ConfHead[metric1,metric2][(RicciScalar@cd1)[]],(RicciScalar@cd2)[]],
RuleDelayed@@Hold[ConfHead[metric1,metric2][(Christoffel@cd1)[i1_,i2_,i3_]],confa^(WeightOfIndicesList[{i1,i2,i3}]-1)*(Christoffel@cd2)[i1,i2,i3]],
RuleDelayed@@Hold[ConfHead[metric1,metric2][(Determinant[metric1,AIndex])[]],(* This is removed because now xTensor is patched confa2^DimOfManifold[M]. Thanks to Leo Stein.*)(Determinant[metric2,AIndex])[]],

(* This line below is not working well.  The problem should be considered later when xTensor knows how to handle the epsilon of a frozen metric. So this really works only when metric1 is the ambient metric... *)
RuleDelayed@@Hold[ConfHead[metric1,metric2][(epsilon@metric1)[inds__?(Length[{#}]===DimOfManifold[M]&)]],(*confa2^(DimOfManifold[M]/2)*)confa^(WeightOfIndicesList[{inds}])(epsilon@metric1)[inds]],

(* Not really satisfactory but minimalist for scalar functions *)
(* Following Leo Stein suggestion, we allow the scalar function to have several arguments *)
ConfHead[metric1,metric2][f_?ScalarFunctionQ[arg___]]:>Simplify[confa2^((ConformalWeight[f])/2),Assumptions->confa>0]f[arg],

ConfHead[metric1,metric2][tens_?xTensorQ[indss___]]:>Simplify[confa2^(ConformalWeight[tens[indss]]/2),Assumptions->confa>0]tens[indss]
};

res
]
]
)


RemoveInducedDerivative[expr_,cd_]:=Module[{res},With[{h=MetricOfCovD@cd},
If[InducedFrom@h===Null,expr,
With[{g=First@InducedFrom@h},
With[{CD=CovDOfMetric@g},
res=(expr//.cd[ind_][Expr___]:>Projector[h][CD[ind][Expr]])/.Projector[h]->ProjectWith[h];
If[res=!=expr,Print["** Warning: you are using ToMetric or Conformal with induced covariant derivatives. \nThese induced derivatives are first expressed in function of the covariant derivative form which they are induced since we do not know very well how to handle that. **"];];
res
]
]
]
]
];

RemoveAllInducedDerivatives[expr_]:=With[{InducedMetrics=Select[$Metrics,InducedFrom[#]=!=Null&]},
Fold[RemoveInducedDerivative[#1,CovDOfMetric[#2]]&,expr,InducedMetrics]
];


ToMetric[expr_,metric1_?MetricQ]:=If[InducedFrom@metric1=!=Null,expr,
Module[{res,preexpression},
Off[ConformalRules::unknown];
With[{cd1=CovDOfMetric[metric1],$CovDsNotInduced=Select[Rest@$CovDs,InducedFrom[MetricOfCovD[#]]===Null&]},
With[{$CovDsNotInducedRelatedTocd1=Select[$CovDsNotInduced,(Catch@ConformalRules[metric1,MetricOfCovD[#]]=!=Null)&]},

preexpression=(RemoveAllInducedDerivatives[expr]//ProjectorToMetric//EinsteinToRicci//WeylToRiemann//ContractMetric//ToCanonical);

res=ChangeCovD[#,$CovDsNotInduced,cd1]&@
ChristoffelToGradConformal[#,$CovDsNotInducedRelatedTocd1,cd1]&@
MyChangeChristoffel[#,$CovDsNotInducedRelatedTocd1,cd1]&@
ChangeCovD[#,$CovDsNotInduced,cd1]&@
ChangeCurvature[#,$CovDsNotInduced,cd1]&@preexpression;
Off[ConformalRules::unknown];

Fold[(#1/.ConformalRules[MetricOfCovD[#2],metric1])&,res,$CovDsNotInducedRelatedTocd1]
]
]
]
];

ToMetric[expr_]:=ToMetric[expr,First@$Metrics];

SetNumberOfArguments[ToMetric,{1,2}]
Protect[ToMetric];


InverseMetricQ[x_?xTensorQ]:=With[{tid=TensorID@x},(Length@tid>0)&&(tid[[1]]===xAct`xTensor`Private`InvMetric)]
InverseMetricQ[_]:=False


SeparateIndicesDownOfInverseMetric[invmetric_?InverseMetricQ][expr_]:=Fold[SeparateMetric[First@$Metrics][#1,#2]&,expr,IndicesOf[Down,invmetric][expr]];
SeparateIndicesDownOfInverseMetric[_][expr_]:=expr


Conformal[metricbase_?MetricQ][metric1_?MetricQ,metric2_?MetricQ][expr_]:=Module[{cdb,cd1,cd2,res,res2,(*oldpre,*)resbis,exprnoproj,M,i1,i2,beforeputtingconfheads,IDInvMetric,res3,firstmetric},
(* The conflict with ScreenDollarIndicea has now been solved. So there is no need to redefine temporarily $PrePrint*)
(*oldpre=$PrePrint;$PrePrint=Identity;*)

(* we define the Covds associated with the metric. The starting metric is metric1, the conformally transformed metric is metric2, and metricbase i the base metric for raising and lowering indiced. It might be one of the other two, but it might not be...*)
cdb=CovDOfMetric[metricbase];
cd1=CovDOfMetric[metric1];Off[ConformalFactor::"unknown"];
cd2=CovDOfMetric[metric2];


exprnoproj=expr//ProjectorToMetric;

M=ManifoldOfCovD[cd1];
i1=DummyIn[Tangent[M]];
i2=DummyIn[Tangent[M]];

IDInvMetric=SeparateIndicesDownOfInverseMetric[Inv[metric1]];

beforeputtingconfheads=IDInvMetric[(IndicesDown@ToMetric[exprnoproj,metric1])
/.Scalar[ex_]:>Scalar[IndicesDown[ex]]/.sf_?ScalarFunctionQ[args___]:>sf@@IndicesDown/@{args}]
/.Scalar[ex_]:>Scalar[IDInvMetric[ex]]/.sf_?ScalarFunctionQ[args___]:>sf@@IDInvMetric/@{args};
(* Above we make sure that IndicesDown and SeparateIndicesDownOfInverseMetric goes inside the scalar Head*)

(*We use ToMetric to have only references to the metric1 and its associated CovD and curvature tensors *)
(*Print["beforeputtingconfheads ",beforeputtingconfheads];*)

(* Then we place the ConfHead on every expression to perform formally the conformal transformation. *)
res=(beforeputtingconfheads
(* Dirty case of scalar functions *)
/.f_?ScalarFunctionQ[ex___]:>ConfHead[metric1,metric2][f[ex]]
(* tensors *)
/.tens_?xTensorQ[inds___]:>ConfHead[metric1,metric2][tens[inds]]
(* Covariant derivatives *)
/.cd1[i1_?DownIndexQ]:>cd2[i1]
(* Now that we have ConfHead everywhere we need to remove the head by specifying the change*)
(* First Obvious rules for metric and inverse metric*)
/.ConfHead[metric1,metric2][metric1[i1_?DownIndexQ,i2_?DownIndexQ]]:>metric2[i1,i2]
/.ConfHead[metric1,metric2][Inv[metric1][i1_?UpIndexQ,i2_?UpIndexQ]]:>Inv[metric2][i1,i2]);

(* And then all other rules to remove the ConfHead*)


(* A bug with Scalar expressions. I have added a NoSCalar . Check if this works fine for the rest. It seems so.*)
(* It seems to be fine so I should post the modification*)
(* This avoids the appearance of sums inside SymmetryOf which were resulting from a removed Scalar head.*)
resbis=NoScalar[res//.RulesConf[metric1,metric2]];


(* So here we have conformally transformed the expression, but now we want to express it in function of the original metric and orginal covD etc...
Indeed at that point, we still have the second metric, and the Riemann of the second metric for instance. SO we use again ToMetric*)

(*$PrePrint=oldpre;*)

On[ConformalFactor::"unknown"];
Off[ToCanonical::"cmods"];

res2=ToCanonical@ContractMetric@NoScalar[ToMetric[resbis,metricbase]];



(* Cyril Pitrou : March 2018. I have added an attempt to patch the problems I have with multiples frozen metrics coming from multiple conformal transformations. Does not affect xPand as usuall in xPand we never use the conformally related metrics and always use the base metric *)
(* I think it works... The idea is that since for frozen metrics some simplifications to Dirac do not occur when they should, I switch to active metric, let xTensor put the Dirac, and switch back to the desired metric. *)
(* It happened that some definitions in xTensor were missing for frozen metrics. As of 28/03/2018 version of xAct this should be fixed and the patch I am adding below should not be necessary...*)
firstmetric=First@$Metrics;
If[metricbase=!=firstmetric,
res3=ContractMetric[((ToCanonical@NoScalar@ContractMetric[res2/.ConformalRules[metricbase,firstmetric]])/.ConformalRules[firstmetric,metricbase])];,
res3=res2;];

On[ToCanonical::"cmods"];
res3
]

(* In case the base metric is unspecified, it is the base metric of course...*)
Conformal[metric1_?MetricQ,metric2_?MetricQ][expr_]:=Conformal[First@$Metrics][metric1,metric2][expr]




ToCosmicTime[expr_,h_?InducedMetricQ]:=
If[$ConformalTime===False,expr,

$ConformalTime=False;
Module[{cd1=CovDOfMetric[h],i1,
n=Last@InducedFrom@h,exprtemp,res},
i1=DummyIn[Tangent[ManifoldOfCovD[cd1]]];

exprtemp=expr(*/.Tens_?((DefProjectedTensorQ[#,h])&)[LI[p_],LI[q_],inds___]\[RuleDelayed]ConformalToCosmic[Tens[LI[p],LI[q],inds]]*);

res=PostProcess[h]@NoScalar@SameDummies@ToCanonical@ContractMetric[
exprtemp
/.H[h][LI[0],LI[q_?IntegerQ]]:>ToCanonical[Nest[LieD[n[i1]][#]&,a[h][]H[h][LI[0],LI[0]],q],UseMetricOnVBundle->None]

/.K[h][LI[0],LI[q_?IntegerQ],i1_,i2_]:>ToCanonical[Nest[LieD[n[i1]][#]&,a[h][]K[h][LI[0],LI[0],i1,i2],q],UseMetricOnVBundle->None]

/.tens_?((DefProjectedTensorQ[#,h]&&#=!=K[h]&&#=!=a[h]&&#=!=H[h])&)[LI[p_],LI[q_?IntegerQ],inds___]:>(ToCanonical[Nest[LieD[n[i1]][#]&,tens[LI[p],LI[0],inds],q],UseMetricOnVBundle->None])
];
(*$ConformalTime=True;*)
res
]
]

ToConformalTime[expr_,h_?InducedMetricQ]:=
If[$ConformalTime===True,expr,

$ConformalTime=True;
Module[{cd1=CovDOfMetric[h],i1,
n=Last@InducedFrom@h,exprtemp,res},
i1=DummyIn[Tangent[ManifoldOfCovD[cd1]]];

exprtemp=expr(*/.Tens_?((DefProjectedTensorQ[#,h])&)[LI[p_],LI[q_],inds___]\[RuleDelayed]CosmicToConformal[Tens[LI[p],LI[q],inds]]*);

res=PostProcess[h]@NoScalar@SameDummies@ToCanonical@ContractMetric[
exprtemp
/.H[h][LI[0],LI[q_?IntegerQ]]:>ToCanonical[Nest[(*It is this factor of 1/a which was forgotten *)1/a[h][]LieD[n[i1]][#]&,1/a[h][]*H[h][LI[0],LI[0]],q],UseMetricOnVBundle->None]

/.K[h][LI[0],LI[q_?IntegerQ],i1_,i2_]:>ToCanonical[Nest[(*It is this factor of 1/a which was forgotten and lead to the bug *)1/a[h][]LieD[n[i1]][#]&,1/a[h][]K[h][LI[0],LI[0],i1,i2],q],UseMetricOnVBundle->None]

/.tens_?((DefProjectedTensorQ[#,h]&&#=!=K[h]&&#=!=a[h]&&#=!=H[h])&)[LI[p_],LI[q_?IntegerQ],inds___]:>(ToCanonical[Nest[1/a[h][]LieD[n[i1]][#]&,tens[LI[p],LI[0],inds],q],UseMetricOnVBundle->None])

];
(*$ConformalTime=False;*)
res
]
]



(* This function creates a list of Covds of a tensor. That is tens,CD[tens],CD[CD[tens]] etc... *)
TablePatternsCovDs[tens_[inds___],CD_,n_:10]:=
With[{M=ManifoldOfCovD@CD},
Join[{tens[inds]},Table[Fold[(CD[#2]@#1)&,tens[inds],#]&[Take[Pattern[#,_]&/@DummyIn/@Table[Tangent[M],{Range[n]}],i]],{i,1,n}]]
];


PutFreeIndicesDown[expr_]:=Fold[ReplaceIndex[#1,#2->ChangeIndex[#2]]&,expr,IndicesOf[Free,Up][expr]];
PutFreeIndicesUp[expr_]:=(Fold[ReplaceIndex[#1,#2->ChangeIndex[#2]]&,expr,IndicesOf[Free,Down][expr]])


PatternSymbolQ[expr_Symbol]:=(expr===Pattern);
PatternTestSymbolQ[expr_Symbol]:=(expr===PatternTest);

PutPatternIndicesDown[expr_]:=Fold[#1/.pat_?PatternSymbolQ[#2,st_]:>-pat[#2,st]/.a_?PatternTestSymbolQ[-pat_?PatternSymbolQ[#2,st_],fun_]:>-a[pat[#2,st],fun]&,expr,IndicesOf[Free,Up][RemovePatterns[expr]]];

PutPatternIndicesUp[expr_]:=Fold[#1/.pat_?PatternSymbolQ[ChangeIndex[#2],st_]:>-pat[Evaluate@ChangeIndex[#2],st]/.a_?PatternTestSymbolQ[-pat_?PatternSymbolQ[Evaluate@ChangeIndex[#2],st_],fun_]:>-a[pat[Evaluate@ChangeIndex[#2],st],fun]&,expr,IndicesOf[Free,Down][RemovePatterns[expr]]];


RemovePatterns[expr_]:=(expr/.hp_HoldPattern:>First@hp)/.a_PatternTest:>First@a/.a_Pattern:>First@a
RemovePatternTests[expr_]:=(expr/.hp_HoldPattern:>First@hp)/.a_PatternTest:>First@a


ConstantsDecompositionRules[h_?InducedMetricQ]:=
With[{g=First@InducedFrom@h},
With[{M=ManifoldOfCovD@CovDOfMetric[g]},
With[{i=DummyIn[Tangent[M]],dim=DimOfManifold@M,j=DummyIn[Tangent[M]],k=DummyIn[Tangent[M]],\[Mu]=DummyIn[Tangent[M]]},
If[$OpenConstantsOfStructure&&Not[TensorNullQ[CS[h]]]&&(SpaceType[h]==="BianchiB"||SpaceType[h]==="BianchiA")&&dim===4,

BuildRule[Evaluate[{CS[h][i,j,k],epsilon[h][j,k,-\[Mu]](nt[h][\[Mu],i])+av[h][j]h[i,k]-av[h][k]h[i,j]}],MetricOn->All],{}]
]
]
]

SetNumberOfArguments[ConstantsDecompositionRules,1]
Protect[ConstantsDecompositionRules];

ToBianchiType[h_?InducedMetricQ]:=Function[expr,expr/.ConstantsDecompositionRules[h]]


CommuteCDSafe[expr_,cd_]:=Module[{restemp,counter,res,restemp0},
(*Print["CommuteCDSafe",$SortCovDAutomatic,$DebugInfoQ];*)
If[$SortCovDAutomatic,
If[$DebugInfoQ,Print["Preliminar Canonicalisation. "];];
restemp0=SameDummies@ToCanonical@expr;
SortCovDsStart[cd];If[$DebugInfoQ,Print["First Canonicalisation with automatic sorting of Cov Ds.  ",Length[restemp0]," terms."];];
counter=0;
restemp=Map[(If[$DebugInfoQ,counter=counter+1;If[Mod[counter,10]===0(*||counter\[GreaterEqual]4180*),Print["We canonicalize term ",counter," ",#];];];SameDummies@ToCanonical@ContractMetric[#])&,restemp0];
Block[{Print},SortCovDsStop[cd];];,
Block[{Print},Off[Unset::norep];SortCovDsStop[cd];On[Unset::norep];];
restemp=expr;
];
(*Entering the Cov Ds sorting defined by BackgroundSLicing. Made to gather the Laplacian near the perturbations, and to use the transversailty conditions.*)
If[$DebugInfoQ,Print["Entering the commutation of Cov Ds..."];];
counter=0;
res=FixedPoint[(counter=counter+1;If[$DebugInfoQ,Print["We commute the Cov Ds for the ",counter," time."];];SameDummies@ToCanonical@ContractMetric[#/.$CommutecdRules])&,restemp,10];
If[$DebugInfoQ,Print["Induced Derivatives were commuted ",counter," times"];];
If[counter>=9,Print["** Warning, the induced derivatives are commuting endlessly. Stopped after 10 iterations to avoid an infinite loop. This is anormal behaviour. **"];];
res
]


(*SplitCovDs[expr_,CD_,cd_]:=FixedPoint[FixedPoint[GradNormalToExtrinsicK[#]&,ToInducedDerivative[#,CD,cd]&],expr];*)


CheckSTFTensors[expr_,h_,exceptlist_List]:=Module[{tens},With[{listpb=Cases[Expand@expr,tens_?((And@@(Function[t,t=!=#]/@exceptlist))&&xTensorQ[#]&&Not@DefProjectedTensorQ[#,h]&)[inds___]->tens,Infinity]},Print["** Warning: the tensor ",#," was not defined with DefProjectedTensor. The rules necessary for its splitting were thus not defined. **"]&/@(DeleteDuplicates@listpb)];];

(* Warning message to update \[Rule] the slicing should be taken into account. *)

CheckSTFTensors[expr_,h_]:=CheckSTFTensors[expr,h,{}];


FullToInducedDerivative[expr_,CD_,cd_]:=FixedPoint[(ToInducedDerivative[#,CD,cd]//ProjectorToMetric//GradNormalToExtrinsicK)&,expr];

PutDownCD[expr_,supercd_,cd_]:=
With[{indmetric=MetricOfCovD[cd]},With[{vector=Last@InducedFrom[indmetric]},expr/.supercd[ind_][expr1:(_?xTensorQ[___]|cd[_][_]|LieD[_][_])]:>IndicesDown[supercd[ind][expr1]]/;OrthogonalToVectorQ[vector][expr1]
]
];

FullToInducedDerivativeAndCDDown[expr_,CD_,cd_]:=FixedPoint[(ToInducedDerivative[PutDownCD[#,CD,cd],CD,cd]//ProjectorToMetric//GradNormalToExtrinsicK)&,expr];


ContractMetricOutSideProjector[expr_Plus,cd_]:=ContractMetricOutSideProjector[#,cd]&/@expr;
ContractMetricOutSideProjector[rest_. ih_?InertHeadQ[inside_],cd_]:=ContractMetric[rest]ih[inside];
ContractMetricOutSideProjector[rest_.  LieD[u_][ex_],cd_]:=ContractMetric[rest]LieD[u][ex];
ContractMetricOutSideProjector[rest_.  cd_[a_][ex_],cd_]:=ContractMetric[rest]cd[a][ex];
ContractMetricOutSideProjector[rest_,cd_]:=ContractMetric[rest];


(* Determine if expr has (at least) n nested CDs acting on inner, which will probably be a pattern expression (it can even be _) *)
ContainsDerOrderQ[expr_,CD_?CovDQ,inner_,n_Integer?NonNegative]:=!FreeQ[expr,Nest[CD[_],inner,n]];
(* Find the maximum order of CD that expr has acting on inner (to simply count the maximum order of CD, use _ for inner) *)
MaxDerOrder[expr_,CD_?CovDQ,inner_]:=-1+NestWhile[#+1&,0,ContainsDerOrderQ[expr,CD,inner,#]&];
(* Written by Leo Stein *)


(* If the background field method is used, this rule ensures that tensor with a label index for order strictly larger than 1 vanish.*)
$BackgroundFieldRule:=If[BackgroundFieldMethod,
{tens_?xTensorQ[inds___]:>0/;PerturbationOrder[tens[inds]]>1},{}]


RulesCovDsOfTensor[expr_,replacerule_,rulesprojected_,h_?InducedMetricQ]:=Module[{dummiesup,tableleft,tableright,testpatternlistnmax,freedownleft,dummiesdown,leftupindices,leftupindicesnopattern,tableleftnopattern,Listlhsrhsnoreplace,temp,tempdown,rhs,rulepert,leftupindicesnopatterntests,tobecan,oncecan,RulesCovDs},Catch@With[{g=First@InducedFrom@h,u=Last@InducedFrom@h},
With[{CD=CovDOfMetric[g],cd=CovDOfMetric[h],M=ManifoldOfCovD@CovDOfMetric[g]},
With[{CDmax=MaxDerOrder[expr,CD,First@replacerule]},

(* In this function we will precompute the rules for the CovDs of tensor for which the rule is given in the argument replacerule.
For instance if we compute the perturbation of the Ricci tensor, there will be several terms with two covariant derivatives of the perturbed metric. By precomputing the rule once instead of several times, we shall save some computing time. *)

dummiesup=DummyIn/@Table[Tangent[M],{Range[CDmax]}];
dummiesdown=ChangeIndex/@dummiesup;

Off[Pattern::patvar];

freedownleft=IndicesOf[Free,Down][RemovePatterns@First@replacerule];

leftupindices=PutPatternIndicesUp[First@replacerule//ProjectorToMetric];
leftupindicesnopattern=RemovePatterns[leftupindices];
leftupindicesnopatterntests=RemovePatternTests[leftupindices];

(* We build the left hand side *)
tableleft=Join[{leftupindices},Table[Fold[(CD[Pattern[#2,Blank[]]]@#1)&,leftupindices,#]&[Take[dummiesup,i]],{i,1,CDmax}]];

tableleftnopattern=RemovePatterns/@tableleft;

Listlhsrhsnoreplace=Reverse@Transpose[{tableleft,tableleftnopattern}];

RulesCovDs=(
tempdown=IndicesDown[Last@#];
rulepert=leftupindicesnopatterntests:>Evaluate[(InducedDecomposition[leftupindicesnopattern,{h,u}]/.rulesprojected/.Scalar->ProtectMyScalar)];
temp=(tempdown/.rulepert)/.FixProtectScalar;
tobecan=FullToInducedDerivative[temp,CD,cd];
oncecan=ContractMetricOutSideProjector[Expand@tobecan,cd];
rhs=ToCanonical[oncecan,UseMetricOnVBundle->None];

IndexRule[((First@#)/.hp_HoldPattern:>First@hp),rhs])&/@Listlhsrhsnoreplace;

RulesCovDs
]
]
]
];

ProjectionAndBackgroundRuleQ[rule_,h_?InducedMetricQ,n_]:=(PerturbationOrder[RemovePatterns[First@rule]]===0)&&(Projector[h][Last@rule]===0||OrthogonalToVectorQ[n][Last@rule]);

SplitPerturbations[expr_,ListPairs_List,h_?InducedMetricQ]:=Module[{res,restemp,resfinal,restemp2,res0,res1,res2,res3,res4,res4bis,res6bis,res6ter,res6quater,res4ter,res4quater,res5,res6,res7,res8,rules,rulesprojectedandbackground},
With[{g=First@InducedFrom@h,u=Last@InducedFrom@h},
With[{CD=CovDOfMetric[g],cd=CovDOfMetric[h]},

(* This is a proposal to deal with the normal vector automatically *)
(* It remains to see if this design is interesting or if the rules for the normal vector should be placed in ListPairs by the user *)
(* Deactivated for the moment *)

rulesprojectedandbackground=Select[ListPairs,ProjectionAndBackgroundRuleQ[#,h,u]&];

rules=Flatten@(RulesCovDsOfTensor[expr,#,rulesprojectedandbackground,h]&/@ListPairs);
If[$DebugInfoQ,Print["rules ",rules]];

(* TODO This part needs commenting. A lot of commenting. *)
res=AbsoluteTiming[
res0=FullToInducedDerivativeAndCDDown[org[((expr/.$RulesVanishingBackgroundFields[h])//ProjectorToMetric//GaussCodazzi[#,h]&)/.Projector[h]->ProjectWith[h]],CD,cd]//ProjectorToMetric;

 If[$DebugInfoQ,Print["Stage 0", res0]];

(* 7/12/14: Michael Kopp Found a bug below. I have replaced ToCanonical by ToCanonical@NoScalar and this seems to correct for the bug. The problem is the treatment of Scalars by ToCanonical. In some very strange configurations this did not work. This is very strange though and we should warn JMM about it. *)
(* 8/12/14 In fact the solution is to make sure that we expand before using PutScalar. So that there is noever any sum inside a Scalar Head. Just to be sure, I add both solutions.*)
res1 =ToCanonical[NoScalar[(res0/.rules)/.ih_?InertHeadQ[ex_]:>ih[ToCanonical[ex]]],UseMetricOnVBundle->None];
If[$DebugInfoQ,Print["Stage 1  ", res1," \n ",res1/.ListPairs ]];

res2=((res1//.ListPairs)/.$BackgroundFieldRule)/.Projector[h][exp___]:>Projector[h][Expand@NoScalar[exp]];
If[$DebugInfoQ,Print["Stage 2  ",res2]];

res3=(res2/.ProtectMyScalar[exp_]:>IndicesDown[ToCanonical[NoScalar@exp]](*Expand[NoScalar@exp]*));
If[$DebugInfoQ,Print["Stage 3  ",res3]];

res4=NoScalar[res3/.Projector[h]->ProjectWith[h]](*//ContractMetric//ToCanonical*);
If[$DebugInfoQ,Print["Ready for canonicalization."]];

res4bis=SameDummies@ContractMetric@res4;
(*If[$DebugInfoQ,Print["Stage 4 bis"]];*)
If[$DebugInfoQ,Print["Number of terms to be canonicalized before Expand : ",Length[res4bis],"  ",LeafCount[res4bis],". Head is ",Head[res4bis]];];

res4quater=res4bis;

res6=((*res5*)res4quater//MetricToProjector[#,h]&//CommuteCDSafe[#,cd]&);
If[$DebugInfoQ,Print["Stage 6"]];

If[FlatSpaceBool[SpaceType[h]],res6quater=res6;,
res6bis=FixedPoint[SameDummies@Expand[#]&,(res6/.ConstantsOfStructureRules[h])];
If[$DebugInfoQ,Print["Canonicalisation of ",Length[res6bis]," terms starting"];];

res6ter=SameDummies@ToCanonical@ContractMetric[res6bis];

res6quater=FixedPoint[SameDummies@Expand[#]&,(res6ter/.ConstantsDecompositionRules[h])];
If[$DebugInfoQ,Print["Canonicalisation of ",Length[res6quater]," terms starting"];];
];

res7=SameDummies@ToCanonical@ContractMetric[res6quater];
If[$DebugInfoQ,Print["Riemann replaced and simplified. "(*,res7*)]];

res7

];

If[$DebugInfoQ,Print["Constants of structure opened and simplified. "]];

CheckSTFTensors[Last@res,h,{g,u,SymbolJoin[Det,g],Riemann[cd],Ricci[cd],RicciScalar[cd],h,CS[h],nt[h],av[h],\[ScriptK][h],delta,epsilon[h],epsilon[g]}];

If[$DefInfoQ,Print["The Splitting of ",Take[expr,1]," +... was performed in ",First@res ," seconds."];];

res8=(*ToCanonical@ContractMetric@*)NoScalar[Last@res];

(* Brute force contractions of induced derivatives with indiced covD.*)
resfinal=PostProcess[h][res8];

(* The default in xTensor is to not commute the induced Covariant derivatives. So we set it back. Indeed, in some cases I have noteiced it can conflict with xPert and Conformal... I am not sure why. Anyway, It is better to avoird automatic commutation of cd in general except in 'SplitPerturbations' since we want maximal simplification.*)

(*If[$SortCovDAutomatic,Off[Unset::norep];Block[{Print},SortCovDsStop[cd]];On[Unset::norep];];*)

resfinal
]
]
];

IndicesDownOnLieDerivatives[expression_]:=expression/.expr_?DefProjectedTensorQ[LI[n_],LI[q_?(#>=1&)],inds__] :>IndicesDown[expr[LI[n],LI[q],inds]];

PostProcess[h_][expr_]:=If[AnisotropyBool@SpaceType[h]||BianchiBool@SpaceType[h],
ScreenDollarIndices@collect@SameDummies@Expand[IndicesDownOnLieDerivatives[expr]//.$Rulecdh[h]],
ScreenDollarIndices@collect[expr]
]

SplitPerturbations[expr_,h_?InducedMetricQ]:=SplitPerturbations[expr,{},h]
SetNumberOfArguments[SplitPerturbations,{2,3}]
Protect[SplitPerturbations];


ToxPandFromRules[expr_,RulesList_List,h_?InducedMetricQ,n_?IntegerQ]:=
With[{g=First@InducedFrom@h},
If[Not[DefTensorQ[ConformalMetricName[g,a[h]]]]&&SpaceType[h]=!="Minkowski",
If[$DefInfoQ,
Print["** Warning: The conformally related metric ",ConformalMetricName[g,a[h]],"  was not previously defined **"];
Print["** We call DefConformalMetric in order to define it. **"];
];
DefConformalMetric[g,a[h]];
];
SplitPerturbations[ExpandPerturbation@Perturbed[Conformal[g,ConformalMetricName[g,If[SpaceType[h]=!="Minkowski",a[h],1]]][expr],n],RulesList,h]
];

ToxPandFromRules[expr_,h_?InducedMetricQ,n_]:=ToxPandFromRules[expr,{},h,n];

SetNumberOfArguments[ToxPandFromRules,{3,4}]
Protect[ToxPandFromRules]


(* 0.4.0 version of ToxPand*)
(*
ToxPand[expr_,dg_,uf_,duf_,h_?InducedMetricQ,gauge_?GaugeQ,order_?IntegerQ,tiltvalue_:"NotTilted"]:=Module[{ruleslist},
With[{g=First@InducedFrom@h},
ruleslist=Join[SplitMetric[g,dg,h,gauge],SplitMatter[uf,duf,-1,h,gauge,order,tiltvalue]];
If[Not[DefTensorQ[Evaluate[ConformalMetricName[g,a[h]]]]]&&SpaceType[h]=!="Minkowski",
If[$DefInfoQ,
Print["** Warning: The conformally related metric ",ConformalMetricName[g,a[h]],"  was not previously defined **"];
Print["** We call DefConformalMetric in order to define it. **"];
];
DefConformalMetric[g,a[h]];
];
SplitPerturbations[ExpandPerturbation@Perturbed[Conformal[g,ConformalMetricName[g,If[SpaceType[h]=!="Minkowski",a[h],1]]][expr],order],ruleslist,h]
]
];
*)

ToxPand[expr_,dg_,uf_,duf_,h_?InducedMetricQ,gauge_?GaugeQ,order_?IntegerQ,tiltvalue_:"NotTilted"] := ToxPand[expr,dg,{{uf,duf}},h,gauge,order,tiltvalue]

ToxPand[expr_,dg_,ufduflist_List,h_?InducedMetricQ,gauge_?GaugeQ,order_?IntegerQ,tiltvalue_:"NotTilted"]:=Module[{ruleslist},
With[{g=First@InducedFrom@h},
ruleslist=Flatten@Join[SplitMetric[g,dg,h,gauge],Map[SplitMatter[#[[1]],#[[2]],-1,h,gauge,order,tiltvalue]&,ufduflist]  ];
If[Not[DefTensorQ[Evaluate[ConformalMetricName[g,a[h]]]]]&&SpaceType[h]=!="Minkowski",
If[$DefInfoQ,
Print["** Warning: The conformally related metric ",ConformalMetricName[g,a[h]],"  was not previously defined **"];
Print["** We call DefConformalMetric in order to define it. **"];
];
DefConformalMetric[g,a[h]];
];
SplitPerturbations[ExpandPerturbation@Perturbed[Conformal[g,ConformalMetricName[g,If[SpaceType[h]=!="Minkowski",a[h],1]]][expr],order],ruleslist,h]
]
];

SetNumberOfArguments[ToxPand,{6,8}];
Protect[ToxPand];


ExtractComponents[expr_,h_?InducedMetricQ,proj_List,ListIndsToContract_List]:=
Catch@With[{frees=List@@FindFreeIndices[expr],u=Last@InducedFrom@h},
If[Length[proj]=!=Length[frees],Throw@Message[xPanding::error,"Number of projectors not equal to number of free indices."]];
If[Sort[ListIndsToContract]=!=frees,Throw@Message[xPanding::error, "List of indices on which the projection is taken is not equal to the list of free indices"]];
(* frees is the list of free indices in the expression*)

If[Length[Select[proj,And[#=!=u,#=!=h,#=!="Space",#=!="Time"]&]]!=0,Throw@Message[ExtractComponents::invalidprojector,proj]];

(* Then it performs the contraction according to what is written in proj___*)
(* NOte that we also apply /.$RulesVanishingBackgroundFields[h]. In principle this is unnecessary because the expression should have beend computed with SplitPerturbations, so all quantities which vanish on the background are already removed. Just ine case we perform this extra list of replacements, in case the user had prepared the expression in a different manner.*)
PostProcess[h][NoScalar@org@Module[{dummy},
Fold[(dummy=Function[{ind},If[UpIndexQ[ind],UpIndex,DownIndex][DummyIn[VBundleOfIndex[ind]]]][#2[[1]]];Switch[#2[[2]],h,h[-dummy,#2[[1]]],"Space",h[-dummy,#2[[1]]],u,If[UpIndexQ[dummy],-1*If[$ConformalTime,1,a[h][]],1*If[$ConformalTime,1,1/a[h][]]]*u[-dummy],"Time",If[UpIndexQ[dummy],-1*If[$ConformalTime,1,a[h][]],1*If[$ConformalTime,1,1/a[h][]]]*u[-dummy]]ReplaceIndex[#1,#2[[1]]->dummy])&,(expr/.$RulesVanishingBackgroundFields[h]),Transpose[{ListIndsToContract,proj}]]]]]

ExtractComponents[expr_,h_?InducedMetricQ,proj_List]:=With[{frees=List@@FindFreeIndices[expr]},ExtractComponents[expr,h,proj,frees]]

ExtractComponents[expr_,h_?InducedMetricQ]:=VisualizeTensor[expr,h](*expr*)/;Length@IndicesOf[Free][expr]===2||Length@IndicesOf[Free][expr]===1

(*Even when there is nothing to do, we remove all quantities which should cancel on the background. Just in case the expression provided was not obtained from SplitPerturbations, where such cancellation have been made*)
ExtractComponents[expr_,h_?InducedMetricQ]:=(expr/.$RulesVanishingBackgroundFields[h])/;Length@IndicesOf[Free][expr]>2||Length@IndicesOf[Free][expr]===0

SetNumberOfArguments[ExtractComponents,{2,4}];
Protect[ExtractComponents];


VisualizeTensor[expr_,h_?InducedMetricQ]:=With[{u=Last@InducedFrom@h},Grid[{{Null,u,h},{u,ExtractComponents[expr,h,{"Time","Time"}]//$PrePrint,ExtractComponents[expr,h,{"Time","Space"}]//$PrePrint},{h,ExtractComponents[expr,h,{"Space","Time"}]//$PrePrint,ExtractComponents[expr,h,{"Space","Space"}]//$PrePrint}},Frame->All]]/;Length@IndicesOf[Free][expr]===2

VisualizeTensor[expr_,h_?InducedMetricQ]:=With[{u=Last@InducedFrom@h},Grid[{{u,ExtractComponents[expr,h,{"Time"}]//$PrePrint},{h,ExtractComponents[expr,h,{"Space"}]//$PrePrint}},Frame->All]]/;Length@IndicesOf[Free][expr]===1

VisualizeTensor[expr_,h_?InducedMetricQ]:=expr/;Not[Length@IndicesOf[Free][expr]===2]&&Not[Length@IndicesOf[Free][expr]===1]

SetNumberOfArguments[VisualizeTensor,2];
Protect[VisualizeTensor];


ExtractOrder[expr_,n_]:=SeriesCoefficient[expr,{$PerturbationParameter,0,n}]

SetNumberOfArguments[ExtractOrder,2];
Protect[ExtractOrder];


SplitGaugeChange[expr_,ListPairs_List,my\[Xi]_,h_?InducedMetricQ,ord_?IntegerQ]:=Catch@Module[{temp,temp1,rulebackgroundfield,res,oldgauge,ReplaceConfFactorByUnperturbed,DummyConf,confpert,q,n},
With[{conffactor=a[h],g=First@InducedFrom@h},
If[Not[ScalarTensorQ[conffactor]]&&conffactor=!=1,Throw[Message[ToxPand::invalidconffactor]]];

With[{CD=CovDOfMetric[g],M=ManifoldOfCovD@CovDOfMetric[g]},
Off[TagUnset::"norep"];

(* We temporarily authorize the perturbation of the conformal factor *)
If[conffactor=!=1,
conffactor/:Perturbation[conffactor[LI[0],LI[q_],indices___],n_]=.;
conffactor/:Perturbation[conffactor[LI[0],LI[q_],indices___]]=.;
confpert=Perturbed[(Conformal[g,SymbolJoin[g,conffactor,2]][expr]),ord];,
confpert=Perturbed[expr,ord]
];
If[$DebugInfoQ,Print["confpert ",confpert]];

rulebackgroundfield=If[BackgroundFieldMethod,{CD[_][my\[Xi][ind___]]:>0},{}];

temp1=GaugeChange[confpert,my\[Xi]];
If[$DebugInfoQ,Print["temp1 ",temp1]];

(* Then we perform the GaugeChange using the function GaugeChange of xPert *)
temp=PutScalar@Expand[LieDToCovD[ExpandPerturbation[temp1],CD]/.rulebackgroundfield];
If[$DebugInfoQ,Print["temp ",temp]];

(*Then we say again that the conformal factor should not be perturbed*)
If[conffactor=!=1,
Perturbation[conffactor[LI[0],LI[q_],indices___],n_]^:=0/;n>=1;
Perturbation[conffactor[LI[0],LI[q_],indices___]]^:=0;
];

(* And we split the result with all the rules the user wants to use, including those for the vector field used in the conformal transformation*)
res=SplitPerturbations[(temp),ListPairs,h];
On[TagUnset::"norep"];
res
]
]
];
SplitGaugeChange[expr_,my\[Xi]_,h_?InducedMetricQ,ord_]:=SplitGaugeChange[expr,{},my\[Xi],h,ord];

SetNumberOfArguments[SplitGaugeChange,{4,5}];
Protect[SplitGaugeChange];


SplitFieldsAndGaugeChange[expr_,g_?MetricQ,dg_,uf_,duf_,h_?InducedMetricQ,order_?IntegerQ,tiltvalue_:"NotTilted"]:=
Module[{ord,n,res,ruleslist,RulesGaugeChange},
With[{M=ManifoldOfCovD@CovDOfMetric@g,cd=CovDOfMetric@h,u=Last@InducedFrom@h},
Block[{ind=DummyIn@Tangent@M},
If[Not[DefTensorQ[\[Xi][h]]],
DefTensor[Evaluate[\[Xi][h]][LI[ord],ind],M];
];
RulesGaugeChange={\[Xi][h][LI[n_],ind_]:>Evaluate[T[h][LI[n],LI[0]]u[ind]+Lv[h][LI[n],LI[0],ind]+cd[ind][Ls[h][LI[n],LI[0]]]]};
ruleslist=Join[RulesGaugeChange,SplitMetric[g,dg,h,"AnyGauge"],SplitMatter[uf,duf,-1,h,"AnyGauge",order,tiltvalue]];
SplitGaugeChange[expr,ruleslist,\[Xi][h],h,order]
]]]

SetNumberOfArguments[SplitFieldsAndGaugeChange,{7,8}];
Protect[SplitFieldsAndGaugeChange];


RulesVelocitySpatial[h_?InducedMetricQ,uf_,duf_,NormVectorSquare_,gauge_?GaugeQ,order_?IntegerQ,TiltedBool_:False]:=Module[{n,q},With[{M=ManifoldOfCovD@CovDOfMetric@h,cd=CovDOfMetric@h,u=Last@InducedFrom[h]},With[{ind1=DummyIn@Tangent@M,ind2=DummyIn@Tangent@M},
Flatten@Join[
Flatten@Join[
If[gauge==="ScalarFieldComovingGauge",{\[CurlyPhi][LI[n_?(#>=1&)],LI[q_]]:>0},{}],
If[gauge==="IsoDensityGauge",{\[Rho][uf][LI[n_?(#>=1&)],LI[q_]]:>0},{}],
If[TiltedBool,{BuildRule[Evaluate[{uf[ind1],(Sqrt[Scalar[Vspat[h,uf][LI[0],LI[0],ind2]Vspat[h,uf][LI[0],LI[0],-ind2]]-NormVectorSquare])u[ind1]+Vspat[h,uf][LI[0],LI[0],ind1]}]]},
{BuildRule[Evaluate[{uf[ind1],u[ind1]}]]}],

{BuildRule@Evaluate[{duf[LI[1],ind1],
Evaluate[V0[h,uf][LI[1],LI[0]]u[ind1]+cd[ind1][If[gauge==="FluidComovingGauge",-Bs[h][LI[1],LI[0]] ,Vs[h,uf][LI[1],LI[0]] ]]+ BV1 Vv[h,uf][LI[1],LI[0],ind1]]}]}],
Table[
BuildRule@Evaluate[{duf[LI[i],ind1],
Evaluate[V0[h,uf][LI[i],LI[0]]u[ind1]+cd[ind1][If[gauge==="FluidComovingGauge",-Bs[h][LI[i],LI[0]] ,Vs[h,uf][LI[i],LI[0]] ]]+ Vv[h,uf][LI[i],LI[0],ind1]]}],{i,2,order}]
]
]
]
]


AllRulesFromProjectedRule[vector_,NormVectorSquare_,RulesVspat_,h_?InducedMetricQ,order_?IntegerQ]:=
Module[{Rulebackgroundvelocity,RuleBoost,LocalRulesVspat,RuleBoostUpton,RulesBoostUptoOrder},
With[{g=First@InducedFrom@h,u=Last@InducedFrom@h,Boost=V0[h,vector],Vspatial=Vspat[h,vector]},
With[{M=ManifoldOfCovD@CovDOfMetric[g]},
With[{i1=DummyIn@Tangent@M,i2=DummyIn@Tangent@M},

LocalRulesVspat=RulesVspat;(* This is because I have problems with n being replaced automatically in RulesVspat. Honestly I am not sure of why this happens. But this is necessary.*)

RuleBoost[0]:={};

RuleBoostUpton[m_]:=Flatten@Table[ToCanonical@ContractMetric[RuleBoost[i]],{i,0,m}];

RuleBoost[n_?(#>=1&)]:=MakeRule[Evaluate[{
Boost[LI[n],LI[0]],ToCanonical@ContractMetric@PutScalar@Expand[(Boost[LI[n],LI[0]]/.First@IndexSolve[org[ExpandPerturbation@Perturbation[g[-i1,-i2]vector[i1]vector[i2],n]/.LocalRulesVspat]==0,Boost[LI[n],LI[0]]])/.RuleBoostUpton[n-1]/.LocalRulesVspat]/.$BackgroundFieldRule}]];

Join[RulesVspat/.RuleBoostUpton[order]]

]]]];


SplitMatter[uf_,duf_,normvector_,h_?InducedMetricQ,gauge_?GaugeQ,order_?IntegerQ,tiltvalue_:"NotTilted"]:=(
If[Not[DefTensorQ[\[Rho][uf]]]||Not[DefTensorQ[duf]],Print["** Warning: The perturbed velocity, or the fields required to parameterize the splitting of matter fields perturbations splitting  were not previously defined **"];
Print["** DefMatterFields is called to build the perturbation of the vector field and the projected fields needed for future splitting **"];
DefMatterFields[uf,duf,h];
];
Switch[tiltvalue,"Tilted",
AllRulesFromProjectedRule[uf,normvector,RulesVelocitySpatial[h,uf,duf,normvector,gauge,order,True],h,order],
"NotTilted",AllRulesFromProjectedRule[uf,normvector,RulesVelocitySpatial[h,uf,duf,normvector,gauge,order,False],h,order]]);


SetNumberOfArguments[SplitMatter,{6,7}]
Protect[SplitMatter]


\[ScriptCapitalN]0[h_?InducedMetricQ]:=SymbolJoin[\[ScriptCapitalN]0,h];
\[ScriptCapitalN]i[h_?InducedMetricQ]:=SymbolJoin[\[ScriptCapitalN]i,h];
\[ScriptCapitalN][h_?InducedMetricQ]:=SymbolJoin[\[ScriptCapitalN],h];
d\[ScriptCapitalN][h_?InducedMetricQ]:=SymbolJoin[d\[ScriptCapitalN],h];


$ListNormalFields[h_?InducedMetricQ]:={{\[ScriptCapitalN]0[h][],"\[ScriptCapitalN]0"}}


DefNormalFields[h_?InducedMetricQ,PerturbParameter_:\[Epsilon]]:=Module[{ord},

With[{n=Last@InducedFrom@h,
M=ManifoldOfCovD@CovDOfMetric@First@InducedFrom@h,NN=\[ScriptCapitalN][h],dn=d\[ScriptCapitalN][h],Ni=\[ScriptCapitalN]i[h],N0=\[ScriptCapitalN]0[h]},

Block[{\[Mu]},
{\[Mu]}=GetIndicesOfVBundle[Tangent@M,1];

(DefProjectedTensor[#[[1]],h,TensorProperties->{},SpaceTimesOfDefinition->{"Background","Perturbed"},PrintAs->#[[2]]])&/@$ListNormalFields[h];

If[Not[DefTensorQ[Ni]],
DefTensor[Ni[\[Mu]],M];,
If[$DefInfoQ,Print["** Warning: Tensor ",Ni , "is already defined. It cannot be redefined without undefining it. **"];];
];

If[Not[DefTensorQ[NN]],
DefTensor[NN[-\[Mu]],M];,
If[$DefInfoQ,Print["** Warning: Tensor ",NN , "is already defined. It cannot be redefined without undefining it. **"];];
];

If[Not[DefTensorQ[dn]],
DefTensorPerturbation[dn[LI[ord],-\[Mu]],NN[-\[Mu]],M],PrintAs->"\[Delta]"<>PrintAs[n];,
If[$DefInfoQ,Print["** Warning: Tensor ",dn , "is already defined. It cannot be redefined without undefining it. **"];];
];
(* Note here, and that's very important, that the perturbation of the normal vector is defined with a down idnex. That way it is always proportionnal to the time form Subscript[d\[Eta], \[Mu]].*)

If[Not@DefinedPerturbationParameter[$PerturbationParameter],$PerturbationParameter=PerturbParameter;DefinedPerturbationParameter[$PerturbationParameter]=True];

dn[\[Mu]_?AIndexQ]:=dn[LI[0],\[Mu]];
Ni[\[Mu]_]:=(-Projector[h][N0[]NN[\[Mu]]]/.Projector[h]->ProjectWith[h])//ProjectorToMetric;

]
]
]

SetNumberOfArguments[DefNormalFields,{1,2}]
Protect[DefNormalFields];


RulesNormalFields[h_?InducedMetricQ,order_?IntegerQ]:=Module[{},

With[{M=ManifoldOfCovD@CovDOfMetric@First@InducedFrom@h,cd=CovDOfMetric@h,n=Last@InducedFrom@h,dn=d\[ScriptCapitalN][h],N0=\[ScriptCapitalN]0[h]},

With[{ind1=DummyIn@Tangent@M,ind2=DummyIn@Tangent@M},

Flatten@Join[
{BuildRule[Evaluate[{dn[ind1],n[ind1]}]]},
{BuildRule[Evaluate[{N0[LI[0],LI[0]],1}]]},
(* \[ScriptCapitalN]0 is the lapse function in ADM formalism. So the background lapse is 1.*)

(* The perturbation of Subscript[n, \[Mu]] has no spatial components since it is proportionnal to Subscript[d\[Eta], \[Mu]] *)
(* It is THIS PROPERTY which says that this vector, timelike, is the one normal to constant time hypersurfaces *)

Table[
BuildRule@Evaluate[{dn[LI[i],ind1],
Evaluate[N0[LI[i],LI[0]]n[ind1]]}],{i,1,order}]

(* We have to remember that dn is the perturbation of the normal vector when indices are down. These rules are defined with both positions of indices, so the down position is the correct one, and the up position means raising with the background metric.*)

]
]
]
]


RulesForLapse[RulesNspat_,h_?InducedMetricQ,order_?IntegerQ]:=Module[{RuleBoost,LocalRulesNspat,RuleBoostUpton},

With[{g=First@InducedFrom@h,vector=\[ScriptCapitalN][h],Boost=\[ScriptCapitalN]0[h]},

With[{M=ManifoldOfCovD@CovDOfMetric@g},

With[{ind1=DummyIn@Tangent@M,ind2=DummyIn@Tangent@M},

LocalRulesNspat=RulesNspat;
(* Even though one might think that this definition is useless and that I could use RulesNspat directly. This creates a conflict further that I could never understand.*)

RuleBoost[0]:={BuildRule[Evaluate[{Boost[LI[0],LI[0]],1}]]};

RuleBoostUpton[m_]:=Flatten@Table[ToCanonical@ContractMetric[RuleBoost[i]],{i,0,m}];

(* OLD Implementation. *)
(*Now that IndexSolve is corrected to deal correctly with Label indices in version 1.0.5 of xAct, this implementation is obsolete.*)
(*RuleBoost[p_?(#\[GreaterEqual]1&)]:=(
Boost[LI[p],LI[0]]\[RuleDelayed]Evaluate[ToCanonical@ContractMetric@PutScalar[(X[]/.First@IndexSolve[org[ExpandPerturbation@Perturbation[g[ind1,ind2]vector[-ind1]vector[-ind2],p]/.LocalRulesNspat/.Boost[LI[p],LI[0]]\[Rule]X[]]\[Equal]0,X[]])/.RuleBoostUpton[p-1]/.LocalRulesNspat]/.$BackgroundFieldRule]);*)

RuleBoost[p_?(#>=1&)]:=
MakeRule[Evaluate[{Boost[LI[p],LI[0]],ToCanonical@ContractMetric@PutScalar@Expand[(Boost[LI[p],LI[0]]/.First@IndexSolve[org[ExpandPerturbation@Perturbation[g[ind1,ind2]vector[-ind1]vector[-ind2],p]/.LocalRulesNspat]==0,Boost[LI[p],LI[0]]])/.RuleBoostUpton[p-1]/.LocalRulesNspat]/.$BackgroundFieldRule}]];

RuleBoostUpton[order]
]
]
]
]


RulesForPerturbedNormal[RulesNspat_,h_?InducedMetricQ,order_?IntegerQ]:=Join[RulesNspat/.RulesForLapse[RulesNormalFields[h,order],h,order]];



SplitNormalVector[h_?InducedMetricQ,order_?IntegerQ]:=(
If[Not[DefTensorQ[\[ScriptCapitalN]0[h]]]||Not[DefTensorQ[d\[ScriptCapitalN][h]]||Not[DefTensorQ[\[ScriptCapitalN]i[h]]]],Print["** Warning: The perturbed normal vector, or the fields required to parameterize its splitting were not previously defined **"];
Print["** DefNormalFields is called to build the perturbation of the normal vector and the projected fields **"];
DefNormalFields[h];
];
Join[RulesForPerturbedNormal[RulesNormalFields[h,order],h,order],RulesForLapse[RulesNormalFields[h,order],h,order]]);


SetNumberOfArguments[SplitNormalVector,{2}]
Protect[SplitNormalVector]


On[RuleDelayed::rhs];


End[]
EndPackage[]