(* Package Requirement : 
 The package GRTensors (RG&TC-Code) of Sotirios Bonanos should be loaded before loading SurfaceCharges
*)

(* Current Version v1.01 *)

SurfaceCharges-Code

(* BeginPackage["SurfaceCharges`"];*)

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

Unprotect[metric];

LiediffU::usage = "LiediffU[vector2_,coord_,vector_] : Lie differential of vectors 
Input: two vector fields vector and vector2 and the list of coordinates coord 
Output: the commutator [vector,vector2]";

LiediffU[vector2_List,coord_List,vector_List]:=Table[
Sum[D[vector2[[indmu]],coord[[indrho]]]vector[[indrho]]-D[vector[[indmu]],coord[[indrho]]]vector2[[indrho]],
{indrho, Length[coord]}],{indmu, Length[vector2]}] /; Dimensions[vector]== Dimensions[vector2]&&Dimensions[vector]== Dimensions[coord];

Liediffd::usage = "Liediffd[form_,coord_,xiU_] : Lie differential of forms 
Input: the Lie vector field xiU, a one-form form and the list of coordinates coord 
Output: the Lie derivative of form with respect to xiU";

Liediffd[form_List,coord_List,xiU_List] := Table[
Sum[D[form[[indmu]],coord[[indrho]]]xiU[[indrho]]+D[xiU[[indrho]],coord[[indmu]]] form[[indrho]],{indrho,Length[form]}],{indmu, Length[coord]}] /; Dimensions[form]== Dimensions[coord]&&Dimensions[xiU]== Dimensions[coord];

Liediffdd::usage ="Liediffdd[gdd_,coord_,xiU_] : Lie differential of tensor with two covariant indices 
Input: the Lie vector field xiU, a tensor gdd and the list of coordinates coord 
Output: the Lie derivative of gdd with respect to xiU";

Liediffdd[gdd_List,coord_List,xiU_List] :=Module[{temp00,indrho,indmu,indnu,dimmetric},
If[Length[Dimensions[gdd]]==2,,Return["Invalid arguments"];];
If[Dimensions[gdd][[1]]== Dimensions[coord]&&Dimensions[gdd][[2]]== Dimensions[coord]&&Dimensions[xiU]== Dimensions[coord],,Return["Incompatible dimensions among the arguments"];];
dimmetric = Length[coord];  Table[
temp00 = 0;
For[indrho=1,indrho≤dimmetric, indrho++,
temp00 =temp00 + D[gdd[[indmu,indnu]],coord[[indrho]]]xiU[[indrho]]+D[xiU[[indrho]],coord[[indmu]]] gdd[[indrho,indnu]]+D[xiU[[indrho]],coord[[indnu]]] gdd[[indmu,indrho]] ;
]; temp00,{indmu,dimmetric},{indnu,dimmetric}] ]

Commutator::usage = "Commutator[v1_,v2_] : Commutator of v1 and v2 in the coordinates 'coord' that should be specified.";

Commutator[v1_,v2_] := LiediffU[v2,coord,v1]/;(ValueQ[coord]&& Dimensions[coord]==Dimensions[v1]&& Dimensions[coord]==Dimensions[v2])

Cotton::usage ="Cotton[signatureg_] : Cotton tensor of a 3d metric
Needed : RGTensors executed on that metric
Input: Signature (+1 or -1) of the metric
Output: CotUd and Cotdd and the epsilon tensors : epsddd, epsUUU 
Convention : epsUUU[[1,2,3]] = + 1 / Sqrt[Abs[detg]]";

Cotton[signatureg_Integer] :=Module[{g2dd,g2ddjd}, 
If[ValueQ[gdd] && ValueQ[Rdd] && ValueQ[R] &&ValueQ[detg],,Print["RGtensors should be executed prior to this procedure."]; Return[];];
If[Dimensions[gdd] =!= {3,3},Print["The Cotton tensor exists only in 3 dimensions. Incoherent dimensions detected."]; Return[]; ];
If[ValueQ[epsUUU] || ValueQ[epsddd] || ValueQ[CotUd] || ValueQ[Cotdd] ,Print["Variables epsUUU, epsddd, CotUd, Cotdd already assigned. They will be updated."];];
If[signatureg =!=  1 && signatureg =!=  -1 ,Return["Incoherent signature detected"];, Print["Signature of the metric: ",signatureg," ."];]; 
g2dd = Rdd-1/4 gdd R;
epsUUU = 1/ Sqrt[signatureg detg]  {{{0,0,0},{0,0,1},{0,-1,0}},{{0,0,-1},{0,0,0},{1,0,0}},{{0,1,0},{-1,0,0},{0,0,0}}};
epsddd = signatureg Sqrt[signatureg detg]  {{{0,0,0},{0,0,1},{0,-1,0}},{{0,0,-1},{0,0,0},{1,0,0}},{{0,1,0},{-1,0,0},{0,0,0}}};
g2ddjd = covD[g2dd];
CotUd = multiDot[epsUUU,g2ddjd,{2,3},{3,1}];
Cotdd = Lower[CotUd,1];
Print["epsUUU,epsddd,CotUd,Cotdd have been computed."];
]

GenerateEpsilon::usage ="GenerateEpsilon[dim_] : Generate the totally antisymmetric numerically invariant pseudo-tensor Epsilon
Input: dimension dim > 0
Output: Epsilon
Convention: Epsilon[[1,...,dim]] = +1";

GenerateEpsilon[dim_Integer] := Module[{String0,String1,perm00,k,i},
If[dim<1, Return["Invalid arguments."];];
Epsilon=.;
Epsilon = ConstantArray[0,Table[dim,{k,1,dim}]];
perm00 = Permutations[Table[k,{k,1,dim}]];
For[i=1,i≤Length[perm00],i++,
String0="";
For[j=1,j<dim,j++,String0 = String0 <> ToString[perm00[[i]][[j]]] <> ","  ];
String0 = String0 <> ToString[perm00 [[i]][[dim]]];
String1 = "Epsilon[[" <> String0 <> "]] = Signature[" <> ToString[perm00] <>"[[" <>ToString[i] <> "]]];";
ToExpression[String1];
];Epsilon
];

ChargesEinstein::usage = "ChargesEinstein[xiU_,hdd_,signatureg_,{index1_,index2_}] compute the surface charge in Einstein gravity
Need = metric and coord should be defined and RGtensors should have been executed prior to this procedure
Input = A Killing vector or asymptotic Killing vector 'xiU' 
        A metric perturbation 'hdd'
	    The signature of the metric 'signatureg' (1 or -1)
        A surface of codimension 2 specified by the two indices '{index1,index2}' in the coordinates coord that are fixed
Output = The components {index1,index2} of the surface charge associated with the vector xiU under the perturbation
Extra output : ChargeSuppl : supplementary term with respect to the Iyer-Wald formula.
Convention : The Einstein Lagrangian has the form 1/(16π)Sqrt[signatureg g] R
";

ChargesEinstein[xiU_List,hdd_List,signatureg_Integer,{index1_Integer,index2_Integer}] :=Module[{hUd,hUU,deth,xid,xiUjd,xiUjU,xidjU,dethjd,dethjU,hUUjd,hUUjU,Integrand1UU,Integrand2UU,KillingUd,IntegrandEUU,Charge},
If[ValueQ[gdd] && NameQ["Raise"] && NameQ["Lower"]&&NameQ["covD" ]&&ValueQ[gUU] &&ValueQ[detg]   ,, Return["RGtensors should be executed prior to this procedure."];];
If[ValueQ[coord] && ValueQ[metric]   ,, Return["metric and coord not defined."];];
If[signatureg =!=  1 && signatureg =!=  -1 ,Return["Incoherent signature detected"];, Print["Signature of the metric: ",signatureg," ."];];
If[Dimensions[hdd] =!=  Dimensions[metric],Return["Incompatible dimensions detected"]; ];
If[Dimensions[coord] =!= Dimensions[xiU],Return["Incompatible dimensions detected"];];
If[1≤index1 ≤Dimensions[coord][[1]] &&1≤index1 ≤Dimensions[coord][[1]], Print["The surface charge is inversely proportional to Newton constant and should be integrated on the ",coord[[index1]]," , ",coord[[index2]]," = constant surface"];
,Return["Invalid arguments for the choice of n-2 surface"];];
hUd = Raise[hdd,1];hUU = Raise[hUd,2];deth = multiDot[gUU,hdd,{1,1},{2,2}];
xid = Lower[xiU,1];xiUjd = covD[xiU,{1}];xiUjU = Raise[xiUjd,2];xidjU = Lower[xiUjU,1];
dethjd = covD[deth];dethjU = Raise[dethjd,1];hUUjd = covD[hUU,{1,2}];hUUjU = Raise[hUUjd,3];
Integrand1UU = Outer[Times,xiU,Contract[hUUjd,{2,3}]-dethjU];
Integrand2UU = -Transpose[multiDot[xid,hUUjU,{1,2}],{2,1}]-1/2 deth Transpose[xiUjU,{2,1}]+ 1/2 multiDot[hUU,xiUjd-Transpose[xidjU,{2,1}],{2,2}];
KillingUd = xiUjd+Transpose[xidjU,{2,1}];
IntegrandEUU = -1/2 multiDot[hUU,KillingUd,{2,2}];
ChargeSuppl =  1/16/Pi  Sqrt[signatureg detg]  (IntegrandEUU[[index1,index2]]-IntegrandEUU[[index2,index1]]);
Charge = 1/16/Pi  Sqrt[signatureg detg]  (Integrand1UU[[index1,index2]]-Integrand1UU[[index2,index1]]+Integrand2UU[[index1,index2]]-Integrand2UU[[index2,index1]])
];

ChargesEinsteinMaxwellScalar::usage = "ChargesEinsteinMaxwellScalar[{xiU_,λ_},{Ad_,ϕ_},{hdd_,ad_,varϕ_},signatureg_,{index1_,index2_},{kinϕ_,kinA_,ϕUneval_}] compute the surface charge in Einstein gravity coupled to Maxwell and one scalar field
Need = metric and coord should be defined and RGtensors should have been executed prior to this procedure
Input = A couple of (asymptotic) symmetry parameters : 'xiU' and 'λ' 
        The background gauge field 'Ad' and scalar fields 'ϕ'
        The variation of the metric 'hdd', gauge field 'ad' and scalar field 'varϕ'
	    The signature of the metric 'signatureg' (1 or -1)
        A surface of codimension 2 specified by the two indices '{index1,index2}' in the coordinates coord that are fixed
	    The normalisation constants 'kinA[ϕUneval]', 'kinϕ[ϕUneval]' in the Lagrangian as a function of the special variable 'ϕUneval' 
Output = The components {index1,index2} of the surface charge associated with the symmetry (xiU,λ) under the perturbation
Extra output : ChargeSuppl : supplementary term with respect to the Iyer-Wald formula.
Convention : The Lagrangian has the form
        L = \!\(\*FractionBox[SqrtBox[\(signatureg\\ g\)], \(16  π\)]\)(R - (1/4)kinA[ϕ] (F^(μν))(Subscript[F, μν]) - (1/2)kinϕ[ϕ](∂_μ)ϕ(∂^μ)ϕ)
";

ChargesEinsteinMaxwellScalar[{xiU_List,λ_},{Ad_List,ϕ_},{hdd_List,ad_List,varϕ_},signatureg_Integer,{index1_Integer,index2_Integer},{kinϕ_,kinA_,ϕUneval_Symbol}] :=Module[{hUd,hUU,deth,xid,xiUjd,xiUjU,xidjU,dethjd,dethjU,hUUjd,hUUjU,Integrand1UU,Integrand2UU,KillingUd,IntegrandEUU,ChargeEinstein,dim,Fdd,FdU,FUU,Integrand3UU,temp,xiA,LieAU,aU,fdd,fUU,Integrand4UU,Integrand5UU,ϕjd,ϕjU,ContribMaxwell,ContribScal,DkEval,kinAEval,kinϕEval},
If[ValueQ[gdd] && NameQ["Raise"] && NameQ["Lower"]&&NameQ["covD" ]&&ValueQ[gUU] &&ValueQ[detg]  &&NameQ["multiDot"] ,, Return["RGtensors should be executed prior to this procedure."];];
If[ValueQ[coord] && ValueQ[metric]   ,, Return["metric and coord not defined."];];
If[signatureg =!=  1 && signatureg =!=  -1 ,Return["Incoherent signature detected"];, Print["Signature of the metric: ",signatureg," ."];];
If[Dimensions[hdd] =!=  Dimensions[metric],Return["Incompatible dimensions detected"]; ];
If[Dimensions[coord] =!= Dimensions[xiU] && Dimensions[coord] =!= Dimensions[Ad]&&Dimensions[coord] =!= Dimensions[ad],Return["Incompatible dimensions detected"];];
If[Dimensions[kinϕ] =!={}&& Dimensions[kinA] =!={},Return["Incoherent kinetic terms detected"];];
If[1≤index1 ≤Dimensions[coord][[1]] &&1≤index1 ≤Dimensions[coord][[1]], Print["The surface charge is inversely proportional to Newton constant and should be integrated on the ",coord[[index1]]," , ",coord[[index2]]," = constant surface"];
,Return["Invalid arguments for the choice of n-2 surface"];];
dim = Dimensions[coord][[1]];
hUd = Raise[hdd,1];hUU = Raise[hUd,2];deth = multiDot[gUU,hdd,{1,1},{2,2}];
xid = Lower[xiU,1];xiUjd = covD[xiU,{1}];xiUjU = Raise[xiUjd,2];xidjU = Lower[xiUjU,1];
dethjd = covD[deth];dethjU = Raise[dethjd,1];hUUjd = covD[hUU,{1,2}];hUUjU = Raise[hUUjd,3];
Integrand1UU = Outer[Times,xiU,Contract[hUUjd,{2,3}]-dethjU];
Integrand2UU = -Transpose[multiDot[xid,hUUjU,{1,2}],{2,1}]-1/2 deth Transpose[xiUjU,{2,1}]+ 1/2 multiDot[hUU,xiUjd-Transpose[xidjU,{2,1}],{2,2}];
KillingUd = xiUjd+Transpose[xidjU,{2,1}];
IntegrandEUU = -1/2 multiDot[hUU,KillingUd,{2,2}];
ChargeEinstein = 1/16/Pi  Sqrt[signatureg detg]  (Integrand1UU[[index1,index2]]-Integrand1UU[[index2,index1]]+Integrand2UU[[index1,index2]]-Integrand2UU[[index2,index1]]);
(* Maxwell Part *)
Fdd = Table[Table[D[Ad[[i]],coord[[j]]]-D[Ad[[j]],coord[[i]]] ,{i,1, dim}],{j,1,dim}];
FdU = Raise[Fdd,2];FUU = Raise[FdU,1];
DkEval = D[kinA,ϕUneval] /. {ϕUneval -> ϕ};
kinAEval = kinA/. {ϕUneval -> ϕ};
kinϕEval = kinϕ /. {ϕUneval -> ϕ};
Integrand3UU= kinAEval Outer[Times,xiU,multiDot[FUU,ad,{1,1}]];
temp = multiDot[hUU,FdU,{2,1}];
xiA = xiU . Ad + λ;
LieAU = Raise[multiDot[xiU,Fdd,{1,1}]+Table[D[xiA,coord[[i]]] ,{i,1, dim}],1];
fdd = Table[Table[D[ad[[i]],coord[[j]]]-D[ad[[j]],coord[[i]]] ,{i,1, dim}],{j,1,dim}];
fUU =  Raise[fdd,1,2];
aU = Raise[ad,1];
Integrand4UU= kinAEval((-1/4 deth xiA -1/2(xiU . ad) )FUU +  xiA temp -1/2 xiA fUU) - 1/2 xiA DkEval  varϕ  FUU;
Integrand5UU=-1/2 kinAEval Outer[Times,aU,LieAU] ;
ContribMaxwell=1/(16 Pi)Sqrt[signatureg detg] (Integrand3UU[[index1,index2]]-Integrand3UU[[index2,index1]]+ Integrand4UU[[index1,index2]]- Integrand4UU[[index2,index1]]+Integrand5UU[[index1,index2]]-Integrand5UU[[index2,index1]]);
(* Scalar Part *)
ϕjd = Table[D[ϕ,coord[[i]]] ,{i,1, dim}];ϕjU = Raise[ϕjd,1];
ContribScal=-kinϕEval /(16 Pi)Sqrt[signatureg detg]varϕ(xiU[[index1]] ϕjU[[index2]]-xiU[[index2]] ϕjU[[index1]]);
(* Result*)
ChargeSuppl =  1/16/Pi  Sqrt[signatureg detg] (IntegrandEUU[[index1,index2]]-IntegrandEUU[[index2,index1]]+Integrand5UU[[index1,index2]]-Integrand5UU[[index2,index1]]);
Charge = ChargeEinstein + ContribMaxwell+  ContribScal
];

ChargesTMG::usage = "ChargesTMG[xiU_,hdd_,mu_,{signatureg_,signepsUUU_},{index1_,index2_}] compute the surface charge in Topologically massive gravity
Need = metric and coord should be defined and RGtensors[...] should have been executed prior to this procedure
Input = A Killing vector or asymptotic Killing vector 'xiU' 
        A metric perturbation 'hdd'
	    The signature of the metric 'signatureg' (1 or -1) and the convention for the epsilon tensor with upper indices 
               eps^{123} = 'signepsUUU' where 123 is the order of the coordinates 'coord'
        mu : Gravitational Chern-Simons coupling
        A surface of codimension 2 specified by the two indices '{index1,index2}' in the coordinates coord that are fixed
Output = The components {index1,index2} of the surface charge associated with the vector xiU under the perturbation
Extra output : ChargeSuppl : supplementary term with respect to the covariant phase space formula.
Convention : The Lagrangian has the form (1/(16π) Sqrt[signatureg g] (R + 1/mu Topological Chern-Simons). 
             The equations of motion are G_{\mu\nu}+1/mu C_{\mu\nu} = 0
Origin : The formulae are taken from 0808.1911. The charges -('Charges' - 'ChargeSuppl') are equal to the charges defined in 0706.0263.
";

ChargesTMG[xiU_List,hdd_List,mu_,{signatureg_Integer,signepsUUU_Integer},{index1_Integer,index2_Integer}] :=Module[{hUd,hUU,deth,xid,xidjd,xiUjd,dethjd,dethjU,hUUjd,hUUjU,xidjU,epsUUU,etaU,xitotU,xitotd,xitotUjd,xitotUjU,IntegrandEUU ,result1,KillingUd,hddjd,hddjdjd,dethjdjd,hddjdjU,Boxhdd,linRcdd,linRccontr,handR,linEindd,hEUd,
linEinUd,linEincontr,temp1UUU,xidualUU,temp1dd,temp1,temp2d,KillingUdjd,temp1d,resultCS,Charge,Integrand1UU ,Integrand2UU,ContribCS1UU,ContribCS2UU ,ContribCS3UU},
If[ValueQ[gdd] && NameQ["Raise"] && NameQ["Lower"]&&NameQ["covD" ]&&ValueQ[gUU] &&ValueQ[detg]  &&NameQ["multiDot"]&&NameQ["Contract"]&&ValueQ[EUd]  &&ValueQ[R],, Return["RGtensors should be executed prior to this procedure."];];
If[ValueQ[coord] && ValueQ[metric]   ,, Return["metric and coord not defined."];];
If[signatureg =!=  1 && signatureg =!=  -1 ,Return["Incoherent signature detected"];, Print["Signature of the metric: ",signatureg," ."];];
If[signepsUUU =!=  1 && signepsUUU =!=  -1 ,Return["Incoherent convention for epsilon detected"];, Print["Convention for epsilon^{123}: ",signepsUUU," ."];];
If[Dimensions[hdd] =!=  Dimensions[metric],Return["Incompatible dimensions detected"]; ];
If[Dimensions[coord] =!= Dimensions[xiU],Return["Incompatible dimensions detected"];];
If[1≤index1 ≤Dimensions[coord][[1]] &&1≤index1 ≤Dimensions[coord][[1]], Print["The surface charge is inversely proportional to Newton constant and should be integrated on the ",coord[[index1]]," , ",coord[[index2]]," = constant surface"];
,Return["Invalid arguments for the choice of n-2 surface"];];
hUd = Raise[hdd,1];hUU = Raise[hUd,2];deth = Contract[hUd,{1,2}];
xid = Lower[xiU,1];xiUjd = covD[xiU,{1}];xidjd = covD[xid];xidjU = Raise[xidjd,2];
dethjd = covD[deth];dethjU = Raise[dethjd,1];
hUUjd = covD[hUU,{1,2}];hUUjU = Raise[hUUjd,3];
(* Total vector xitotU got a contribution from TMG *)
epsUUU =signepsUUU/Sqrt[signatureg detg] {{{0,0,0},{0,0,1},{0,-1,0}},{{0,0,-1},{0,0,0},{1,0,0}},{{0,1,0},{-1,0,0},{0,0,0}}};
etaU = multiDot[epsUUU,xidjd,{2,2},{3,1}]/2;
xitotU = xiU + etaU/mu;xitotd = Lower[xitotU,1];
xitotUjd = covD[xitotU,{1}];xitotUjU = Raise[xitotUjd,2];
(* Einstein charges but for the total vector xitotU *)
Integrand1UU = Outer[Times,xitotU,Contract[hUUjd,{2,3}]-dethjU];
Integrand2UU = -Transpose[multiDot[xitotd,hUUjU,{1,2}],{2,1}]-1/2 deth Transpose[xitotUjU,{2,1}]+  multiDot[hUU,xitotUjd,{2,2}];
KillingUd = xiUjd+Transpose[xidjU,{2,1}];
IntegrandEUU = -1/2 multiDot[hUU,KillingUd,{2,2}];
result1 =  Integrand1UU[[index1,index2]]-Integrand1UU[[index2,index1]]+Integrand2UU[[index1,index2]]-Integrand2UU[[index2,index1]]+IntegrandEUU[[index1,index2]]-IntegrandEUU[[index2,index1]];
(* Contibutions coming only from TMG *)
hddjd = covD[hdd]; hddjdjd = covD[hddjd];
hddjdjU = Raise[hddjdjd,4]; temp1dd = Contract[hddjdjU ,{1,4}];
Boxhdd = Contract[hddjdjU ,{3,4}];dethjdjd = covD[dethjd];
linRcdd =( temp1dd +Transpose[temp1dd,{2,1}]-Boxhdd-dethjdjd)/2;
linRccontr = multiDot[gUU,linRcdd,{1,1},{2,2}];
handR = multiDot[hUU,Rdd,{1,1},{2,2}];
linEindd = linRcdd -R  hdd /2+(handR -linRccontr)gdd / 2;
hEUd = multiDot[hUd,EUd,{2,1}];linEinUd = Raise[linEindd,1]-hEUd;
linEincontr = Contract[linEinUd,{1,2}];temp1UUU = multiDot[epsUUU,linEinUd,{3,2}];
xidualUU = multiDot[epsUUU,xid,{3,1}];
ContribCS1UU = multiDot[temp1UUU,xid,{3,1}]-1/2 linEincontr xidualUU;
temp1 = Contract[hEUd,{1,2}];temp2d =  multiDot[xid,EUd,{1,1}]+1/2 R xid ;
ContribCS2UU = 1/2 temp1 xidualUU +1/4 deth multiDot[epsUUU,temp2d,{3,1}];
(* Extra contribution of the Chern-Simons term (see 0808.1911 ) *)
KillingUdjd = covD[KillingUd,{1}];temp1d = Contract[KillingUdjd,{1,3}];temp2d = Contract[KillingUdjd,{1,2}];
ContribCS3UU = 1/4 deth multiDot[epsUUU,temp1d-temp2d,{3,1}] ;(* TO CHANGE *)
resultCS = ContribCS1UU[[index1,index2]]-ContribCS1UU[[index2,index1]]+ContribCS2UU[[index1,index2]]-ContribCS2UU[[index2,index1]]+(ContribCS3UU[[index1,index2]]-ContribCS3UU[[index2,index1]]);
ChargeSuppl =  1/16/Pi  Sqrt[signatureg detg] (IntegrandEUU[[index1,index2]]-IntegrandEUU[[index2,index1]] -  (ContribCS3UU[[index1,index2]]-ContribCS3UU[[index2,index1]])/mu);
Charge = 1/16/Pi  Sqrt[signatureg detg] (result1 - resultCS/mu)
];

ChargespForm::usage = "ChargespForm[{xiU_,λ_},{Addotsd_,ϕ_},{hdd_,addotsd_,varϕ_},signatureg_,{index1_,index2_},{kinA_,ϕUneval_}] compute the surface charge of the p-form action
Need = metric and coord should be defined and RGtensors should have been executed prior to this procedure
Input = A couple of (asymptotic) symmetry parameters : 'xiU' and 'λ' 
        The background gauge field 'Addotsd' and scalar fields 'ϕ'
        The variation of the metric 'hdd', gauge field 'addotsd' and scalar field 'varϕ'
	    The signature of the metric 'signatureg' (1 or -1)
        A surface of codimension 2 specified by the two indices '{index1,index2}' in the coordinates coord that are fixed
	    The normalisation constants 'kinA[ϕUneval]' in the Lagrangian as a function of the special variable 'ϕUneval' 
Output = The components {index1,index2} of the surface charge associated with the symmetry (xiU,λ) under the perturbation
Extra output : ChargeSuppl : supplementary term with respect to the Iyer-Wald formula.
Convention : The Lagrangian has the form
        L = \!\(\*FractionBox[SqrtBox[\(signatureg\\ g\)], \(16  π\)]\)( - \!\(\*FractionBox[\(1\), \(2\\ \(\((p + 1)\)!\)\)]\)kinA[ϕ] (F^((Subscript[μ, 1]) ..  (Subscript[μ, p + 1])))(Subscript[F, (Subscript[μ, 1]) ...  (Subscript[μ, p + 1])]) )
";

ChargespForm[{xiU_List,λ_},{Addotsd_List,ϕ_},{hdd_List,addotsd_List,varϕ_},signatureg_Integer,{index1_Integer,index2_Integer},{kinA_,ϕUneval_Symbol}] :=Module[{hUd,hUU,deth,p,dim,perm,permp,dAddotsd,Fddotsd,FUdotsU,Integrand3UU,temp,xiAddotd,xiaddotd,FxiAUU,FxiaUU,hFxiAUU,hFUdotsU,xiAhHUU,preLie ,LieAddotsd ,LieAUdotsU,daddotsd,fddotsd,fUdotsU,fxiAUU,aUdotsd,Integrand4UU,Integrand5UU, ContribpForm,DkEval,kinAEval},
If[Dimensions[hdd] =!=  Dimensions[gdd],Print["Incompatible dimensions detected"];Return[]; ];
If[Dimensions[coord] =!= Dimensions[xiU],Print["Incompatible dimensions detected"];Return[];];
If[signatureg =!=  1 && signatureg =!=  -1 ,Print["Incoherent signature detected"];Return[];, Print["Signature of the metric: ",signatureg," ."];];
If[IntegerQ[index1] &&  IntegerQ[index2],
If[1≤index1 ≤Dimensions[coord][[1]] &&1≤index1 ≤Dimensions[coord][[1]], Print["The final answer is inversely proportional to Newton constant and should be integrated on the ",coord[[index1]]," , ",coord[[index2]]," = constant surface"];
,Print["Incompatible dimensions detected in the indices"];Return[];
]
,Print["Non-integer arguments for the choice of n-2 surface"];Return[];];
If[Length[Dimensions[Addotsd]]<2,Print["This procedure deals only with p-forms of rank at least 2"];Return[];];

dim = Dimensions[coord][[1]];
p = Length[Dimensions[Addotsd]];
perm = Table[RotateLeft[Table[i,{i,1,p+1}],k],{k,1,p+1}];
permp = Table[RotateLeft[Table[i,{i,1,p}],k],{k,1,p}];
hUd = Raise[hdd,1];hUU = Raise[hUd,2];deth = multiDot[gUU,hdd,{1,1},{2,2}];
 
dAddotsd = Table[D[Addotsd,coord[[j]] ],{j,1,dim}];
Fddotsd =Sum[Signature[perm[[i]]]  Transpose[dAddotsd,perm[[i]]],{i,1,p+1}];
FUdotsU = Raise[Fddotsd,Table[i,{i,1,p+1}]];(* Raise accepts lists as second argument *)
DkEval = D[kinA,ϕUneval] /. {ϕUneval -> ϕ};
kinAEval = kinA/. {ϕUneval -> ϕ};
temp = Outer[Times,FUdotsU,addotsd];(* multiDot does not accept lists as third argument so other implementation *)
For[i=0,i<p,i++,temp = Contract[temp,{2,p-i+2}]];
Integrand3UU=-  2/p kinAEval Outer[Times,xiU,temp ];
xiAddotd =  multiDot[xiU ,Addotsd,{1,1} ]+ λ ;
xiaddotd =  multiDot[xiU ,addotsd,{1,1} ];

temp = Outer[Times,FUdotsU,xiAddotd];
For[i=0,i<p-1,i++,temp = Contract[temp,{3,p-i+2}]];
FxiAUU = temp;
(* FxiAUU = multiDot[FUdotsU,xiAddotd,Table[{i+2,i},{i,1,p-1}] ];*)
temp = Outer[Times,FUdotsU,xiaddotd];
For[i=0,i<p-1,i++,temp = Contract[temp,{3,p-i+2}]];
FxiaUU = temp;
(*FxiaUU = multiDot[FUdotsU,xiaddotd,Table[{i+2,i},{i,1,p-1}] ];*)
hFxiAUU = multiDot[hUd,FxiAUU,{2,1}];
hFUdotsU = multiDot[hUd,FUdotsU,{2,3}];
temp = Table[i,{i,1,p+1}];temp[[1]]=3;temp[[2]] = 1;temp[[3]]=2;
hFUdotsU = Transpose[hFUdotsU ,temp];

temp = Outer[Times,hFUdotsU,xiAddotd];
For[i=0,i<p-1,i++,temp = Contract[temp,{3,p-i+2}]];
xiAhHUU = temp;
(* xiAhHUU = multiDot[xiAddotd,hFUdotsU,Table[{i,i+2},{i,1,p-1}] ];*)
preLie = Table[D[xiAddotd,coord[[i]]],{i,1,dim}]; 

LieAddotsd = multiDot[xiU,Fddotsd,{1,1}]+Sum[ Signature[permp[[i]]] Transpose[preLie,permp[[i]]],{i,1,p}];
LieAUdotsU =Raise[LieAddotsd ,Table[i,{i,1,p}]];
daddotsd = Table[D[addotsd,coord[[j]] ],{j,1,dim}];
fddotsd =Sum[Signature[perm[[i]]]  Transpose[daddotsd,perm[[i]]],{i,1,p+1}];
fUdotsU = Raise[fddotsd,Table[i,{i,1,p+1}]];

temp = Outer[Times,fUdotsU,xiAddotd];
For[i=0,i<p-1,i++,temp = Contract[temp,{3,p-i+2}]];
fxiAUU = temp;
(* fxiaUU = multiDot[fUdotsU,xiAddotd,Table[{i+2,i},{i,1,p-1}] ];*)
Integrand4UU= kinAEval(-1/2 deth  FxiAUU -FxiaUU + (p-1) xiAhHUU+2  hFxiAUU -fxiAUU) - DkEval  varϕ FxiAUU;

aUdotsd = Raise[addotsd,1];
temp = Outer[Times,LieAUdotsU,aUdotsd];
For[i=0,i<p-1,i++,temp = Contract[temp,{2,p- i+2}]];
Integrand5UU= kinAEval temp;
 ContribpForm=1/(2 16 Pi (p-1)!)Sqrt[signatureg detg] (Integrand3UU[[index1,index2]]-Integrand3UU[[index2,index1]]+ Integrand4UU[[index1,index2]]- Integrand4UU[[index2,index1]]+Integrand5UU[[index1,index2]]-Integrand5UU[[index2,index1]]);
(* Result*)
ChargeSuppl = 1/(2 16 Pi (p-1)!)Sqrt[signatureg detg] (Integrand5UU[[index1,index2]]-Integrand5UU[[index2,index1]]);
Charge = ContribpForm
];

KomarIntegral::usage = "KomarIntegral[xiU_,signatureg_,{index1_,index2_}] compute the Komar Integral in Einstein gravity
Need = metric and coord should be defined and RGtensors should have been executed prior to this procedure
Input = A Killing vector or asymptotic Killing vector 'xiU' 
	    The signature of the metric 'signatureg' (1 or -1)
        A surface of codimension 2 specified by the two indices '{index1,index2}' in the coordinates coord that are fixed
Output = The components {index1,index2} of the Komar integral associated with the vector xiU
Convention : The Komar integral has the form 1/(16π) Sqrt[signatureg g] ((D^α)(xiU^β)-(D^β)(xiU^α))((d^(n - 2))x)_(αβ)
";

KomarIntegral[xiU_List,signatureg_Integer,{index1_Integer,index2_Integer}] :=Module[{xiUjd,xiUjU,Charge},
If[ValueQ[gdd] && NameQ["Raise"] && NameQ["Lower"]&&NameQ["covD" ]&&ValueQ[gUU] &&ValueQ[detg]   ,, Return["RGtensors should be executed prior to this procedure."];];
If[ValueQ[coord] && ValueQ[metric]   ,, Return["metric and coord not defined."];];
If[signatureg =!=  1 && signatureg =!=  -1 ,Return["Incoherent signature detected"];, Print["Signature of the metric: ",signatureg," ."];];
If[Dimensions[coord] =!= Dimensions[xiU],Return["Incompatible dimensions detected"];];
If[1≤index1 ≤Dimensions[coord][[1]] &&1≤index1 ≤Dimensions[coord][[1]], Print["The surface charge is inversely proportional to Newton constant and should be integrated on the ",coord[[index1]]," , ",coord[[index2]]," = constant surface"];
,Return["Invalid arguments for the choice of n-2 surface"];];
xiUjd = covD[xiU,{1}];xiUjU = Raise[xiUjd,2];
Charge = 1/16/Pi  Sqrt[signatureg detg]  (xiUjU[[index1,index2]]-xiUjU[[index2,index1]])
];

(* Functions used in Asymptotic Analyses  *)

SeriesWithNTerms::usage = "SeriesWithNTerms[expr_,{coordAsympt_,LimitAsympt_},NumberTerms_]  : Develop in series up to a fixed number of terms 
Input: the expression to be developed expr which can also be an array
           The coordinate to expand onto and the value of the coordinate at the limit  
           The number of terms wished
Output: the series expansion";

SeriesWithNTerms[expr_,{coordAsympt_,LimitAsympt_},NumberTerms_Integer] := 
Module[{temp,result},
If[NumberTerms === Null || coordAsympt === Null || LimitAsympt === Null,Print["Missing arguments"]; Return[expr]];
If[ListQ[expr],SeriesWithNTerms[#,{coordAsympt,LimitAsympt},NumberTerms]&/@ expr,
temp = Series[expr,{coordAsympt,LimitAsympt,-1000}];
If[Length[temp] ≥ 6,If[Part[temp,1] === coordAsympt && Part[temp,2] === LimitAsympt,result = Series[expr,{coordAsympt,LimitAsympt,Ceiling[temp[[4]]/temp[[6]]]+NumberTerms-1}], result=expr],result = expr];result]];

CottonAsympt::usage ="CottonAsympt[signatureg_,{coordAsympt_,LimitAsympt_}, order_] : Asymptotic form of the Cotton tensor of a 3d metric
Needed : RGTensors executed on that metric
Input: Signature (+1 or -1) of the metric,
       The coordinate to expand onto, the value of the coordinate at the limit and the number of terms wished in the final expansion 
Output: CotUd and Cotdd and the epsilon tensors : epsddd, epsUUU 
Convention : epsUUU[[1,2,3]] = + 1 / Sqrt[Abs[detg]]";

CottonAsympt[signatureg_Integer,{coordAsympt_,LimitAsympt_}, order_Integer] :=Module[{g2dd,g2ddjd},
If[ValueQ[gdd] && ValueQ[Rdd] && ValueQ[R] &&ValueQ[detg],,Print["RGtensors should be executed prior to this procedure."]; Return[];];
If[Dimensions[gdd] =!= {3,3},Print["The Cotton tensor exists only in 3 dimensions. Incoherent dimensions detected."]; Return[]; ];
If[ValueQ[epsUUU] || ValueQ[epsddd] || ValueQ[CotUd] || ValueQ[Cotdd] ,Print["Variables epsUUU, epsddd, CotUd, Cotdd already assigned. They will be updated."];];
If[signatureg =!=  1 && signatureg =!=  -1 ,Return["Incoherent signature detected"];, Print["Signature of the metric: ",signatureg," ."];];
g2dd = SeriesWithNTerms[Rdd-1/4 gdd R ,{coordAsympt,LimitAsympt},order];
epsUUU = 1/ Sqrt[signatureg detg] {{{0,0,0},{0,0,1},{0,-1,0}},{{0,0,-1},{0,0,0},{1,0,0}},{{0,1,0},{-1,0,0},{0,0,0}}};
epsddd = signatureg Sqrt[signatureg detg] {{{0,0,0},{0,0,1},{0,-1,0}},{{0,0,-1},{0,0,0},{1,0,0}},{{0,1,0},{-1,0,0},{0,0,0}}};
g2ddjd =SeriesWithNTerms[ covD[g2dd],{coordAsympt,LimitAsympt},order];
CotUd =SeriesWithNTerms[ multiDot[epsUUU,g2ddjd,{2,3},{3,1}],{coordAsympt,LimitAsympt},order];
Cotdd =SeriesWithNTerms[Lower[CotUd,1],{coordAsympt,LimitAsympt},order];
Print["epsUUU,epsddd,CotUd,Cotdd have been computed."];
]

ChargesEinsteinAsympt::usage = "ChargesEinsteinAsympt[xiU_,hdd_,signatureg_,{index1_,index2_},{coordAsympt_,LimitAsympt_}, order_] compute the surface charge asymptotically in Einstein gravity
Need = metric and coord should be defined and RGtensors should have been executed prior to this procedure
Input = A Killing vector or asymptotic Killing vector 'xiU' 
        A metric perturbation 'hdd'
	    The signature of the metric 'signatureg' (1 or -1)
        A surface of codimension 2 specified by the two indices '{index1,index2}' in the coordinates coord that are fixed
        The coordinate to expand onto 'coordAsympt', the value of the coordinate at the limit 'LimitAsympt' and the number of terms wished in the final expansion 'order'
Output = The components {index1,index2} of the surface charge associated with the vector xiU under the perturbation
Extra output : ChargeSuppl : supplementary term with respect to the Iyer-Wald formula.
Convention : The Einstein Lagrangian has the form 1/(16π)Sqrt[signatureg g] R
";

ChargesEinsteinAsympt[xiU_List,hdd_List,signatureg_Integer,{index1_Integer,index2_Integer},{coordAsympt_,LimitAsympt_}, order_Integer] :=
ChargesEinstein[xiU,SeriesWithNTerms[hdd,{coordAsympt,LimitAsympt},order],signatureg,{index1,index2}];

ChargesEinsteinMaxwellScalarAsympt::usage = "ChargesEinsteinMaxwellScalarAsympt[{xiU_,λ_},{Ad_,ϕ_},{hdd_,ad_,varϕ_},signatureg_,{index1_,index2_},{kinϕ_,kinA_,ϕUneval_},{coordAsympt_,LimitAsympt_}, order_] compute asymptotically the surface charge in Einstein gravity coupled to Maxwell and one scalar field
Need = metric and coord should be defined and RGtensors should have been executed prior to this procedure
Input = A couple of (asymptotic) symmetry parameters : 'xiU' and 'λ' 
        The background gauge field 'Ad' and scalar fields 'ϕ'
        The variation of the metric 'hdd', gauge field 'ad' and scalar field 'varϕ'
	    The signature of the metric 'signatureg' (1 or -1)
        A surface of codimension 2 specified by the two indices '{index1,index2}' in the coordinates coord that are fixed
	    The normalisation constants 'kinA[ϕUneval]', 'kinϕ[ϕUneval]' in the Lagrangian as a function of the special variable 'ϕUneval'
        The coordinate to expand onto 'coordAsympt', the value of the coordinate at the limit 'LimitAsympt' and the number of terms wished in the final expansion 'order' 
Output = The components {index1,index2} of the surface charge associated with the symmetry (xiU,λ) under the perturbation
Extra output : ChargeSuppl : supplementary term with respect to the Iyer-Wald formula.
Convention : The Lagrangian has the form
        L = \!\(\*FractionBox[SqrtBox[\(signatureg\\ g\)], \(16  π\)]\)(R - (1/4)kinA[ϕ] (F^(μν))(Subscript[F, μν]) - (1/2)kinϕ[ϕ](∂_μ)ϕ(∂^μ)ϕ)
";

ChargesEinsteinMaxwellScalarAsympt[{xiU_List,λ_},{Ad_List,ϕ_},{hdd_List,ad_List,varϕ_},signatureg_Integer,{index1_Integer,index2_Integer},{kinϕ_,kinA_,ϕUneval_Symbol},{coordAsympt_,LimitAsympt_}, order_Integer] :=ChargesEinsteinMaxwellScalar[{xiU,λ},{SeriesWithNTerms[Ad,{coordAsympt,LimitAsympt},order],SeriesWithNTerms[ϕ,{coordAsympt,LimitAsympt},order]},{SeriesWithNTerms[hdd,{coordAsympt,LimitAsympt},order],SeriesWithNTerms[ad,{coordAsympt,LimitAsympt},order],SeriesWithNTerms[varϕ,{coordAsympt,LimitAsympt},order]},signatureg,{index1,index2},{kinϕ,kinA,ϕUneval}] ;

ChargesTMGAsympt::usage = "ChargesTMGAsympt[xiU_,hdd_,mu_,{signatureg_,signepsUUU_},{index1_,index2_},{coordAsympt_,LimitAsympt_}, order_] compute asymptotically the surface charge in Topologically massive gravity
Need = metric and coord should be defined and RGtensors[...] should have been executed prior to this procedure
Input = A Killing vector or asymptotic Killing vector 'xiU' 
        A metric perturbation 'hdd'
	    The signature of the metric 'signatureg' (1 or -1) and the convention for the epsilon tensor with upper indices 
               eps^{123} = 'signepsUUU' where 123 is the order of the coordinates 'coord'
        mu : Gravitational Chern-Simons coupling
        A surface of codimension 2 specified by the two indices '{index1,index2}' in the coordinates coord that are fixed
        The coordinate to expand onto 'coordAsympt', the value of the coordinate at the limit 'LimitAsympt' and the number of terms wished in the final expansion 'order' 
Output = The components {index1,index2} of the surface charge associated with the vector xiU under the perturbation
Extra output : ChargeSuppl : supplementary term with respect to the covariant phase space formula.
Convention : The Lagrangian has the form (1/(16π) Sqrt[signatureg g] (R + 1/mu Topological Chern-Simons). 
             The equations of motion are G_{\mu\nu}+1/mu C_{\mu\nu} = 0
Origin : The formulae are taken from 0808.1911. The charges -('Charges' - 'ChargeSuppl') are equal to the charges defined in 0706.0263.
";

ChargesTMGAsympt[xiU_List,hdd_List,mu_,{signatureg_Integer,signepsUUU_Integer},{index1_Integer,index2_Integer},{coordAsympt_,LimitAsympt_}, order_Integer] :=
ChargesTMG[xiU,SeriesWithNTerms[hdd,{coordAsympt,LimitAsympt},order],mu,{signatureg,signepsUUU},{index1,index2}];

ChargespFormAsympt::usage = "ChargespFormAsympt[{xiU_,λ_},{Addotsd_,ϕ_},{hdd_,addotsd_,varϕ_},signatureg_,{index1_,index2_},{kinA_,ϕUneval_},{coordAsympt_,LimitAsympt_}, order_] compute the asymptotic form of the surface charge of the p-form action
Need = metric and coord should be defined and RGtensors should have been executed prior to this procedure
Input = A couple of (asymptotic) symmetry parameters : 'xiU' and 'λ' 
        The background gauge field 'Addotsd' and scalar fields 'ϕ'
        The variation of the metric 'hdd', gauge field 'addotsd' and scalar field 'varϕ'
	    The signature of the metric 'signatureg' (1 or -1)
        A surface of codimension 2 specified by the two indices '{index1,index2}' in the coordinates coord that are fixed
	    The normalisation constants 'kinA[ϕUneval]' in the Lagrangian as a function of the special variable 'ϕUneval'
        The coordinate to expand onto 'coordAsympt', the value of the coordinate at the limit 'LimitAsympt' and the number of terms wished in the final expansion 'order'  
Output = The components {index1,index2} of the surface charge associated with the symmetry (xiU,λ) under the perturbation
Extra output : ChargeSuppl : supplementary term with respect to the Iyer-Wald formula.
Convention : The Lagrangian has the form
        L = \!\(\*FractionBox[SqrtBox[\(signatureg\\ g\)], \(16  π\)]\)( - \!\(\*FractionBox[\(1\), \(2\\ \(\((p + 1)\)!\)\)]\)kinA[ϕ] (F^((Subscript[μ, 1]) ..  (Subscript[μ, p + 1])))(Subscript[F, (Subscript[μ, 1]) ...  (Subscript[μ, p + 1])]) )
";

ChargespFormAsympt[{xiU_List,λ_},{Addotsd_List,ϕ_},{hdd_List,addotsd_List,varϕ_},signatureg_Integer,{index1_Integer,index2_Integer},{kinA_,ϕUneval_Symbol},{coordAsympt_,LimitAsympt_}, order_Integer] :=ChargespForm[{xiU,λ},{SeriesWithNTerms[Addotsd,{coordAsympt,LimitAsympt},order],SeriesWithNTerms[ϕ,{coordAsympt,LimitAsympt},order]},{SeriesWithNTerms[hdd,{coordAsympt,LimitAsympt},order],SeriesWithNTerms[addotsd,{coordAsympt,LimitAsympt},order],SeriesWithNTerms[varϕ,{coordAsympt,LimitAsympt},order]},signatureg,{index1,index2},{kinA,ϕUneval}];

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

(* EndPackage[]; *)