(*
Code by Ofek Birnholtz, ofek.birnholtz@mail.huji.ac.il
Cite as
	Birnholtz, Ofek & Hadar, Shahar 2014
	Radiation-reaction of gravitation in general dimension
*)

DefConsts={MS[r]-> CS CS2 (hd^2)-2CS2 hd(hd+1)(w^2)(r^2)+(((hd+1))^2)(w^4)(r^4),α -> -CS CS2 hd (1+hd) ,β[r]-> ((CS CS2)/r) (CS hd-(1+hd)(w^2)(r^2)),β'[r]-> D[((CS CS2)/r) (CS hd-(1+hd)(w^2)(r^2)),r],γ->((CS CS2)/(r^2)) (2(1+hd)(w^2)(r^2)-CS hd)}
DefCS2={CS2->CS-(hd+1)};
DefCS={CS->L(L+hd)};
DefTrick={(α γ-((β[r])^2))->-(((CS^2) (CS2^2) MS[r])/(r^2))}
Print["test:"]
TestTrick=FullSimplify[FullSimplify[(α γ-((β[r])^2))+(((CS^2) (CS2^2) MS[r])/(r^2))/.DefConsts]/.DefCS2]
Print["Homogenous L Scalar:"]
ActionLHomo=FullSimplify[((r^(1+hd))/(MS[r]))(α (∂_r)ϕ[r] (∂_r)Cϕ[r] + β[r] (ϕ[r] (∂_r)Cϕ[r]+Cϕ[r](∂_r)ϕ[r])+γ ϕ[r]Cϕ[r])/.DefCS2]
Print["Source L Scalar:"]
SourceL=(r^(1+hd))(Cϕ[r]J[r]+CJ[r]ϕ[r])
Print["L:"]
ActionL=ActionLHomo+SourceL
ActionByParts=FullSimplify[((∂_r)(((r^(1+hd))/(MS[r]))α (∂_r)ϕ[r])Cϕ[r] +( (∂_r)(((r^(1+hd))/(MS[r]))β[r] ϕ[r]) Cϕ[r]+ ((r^(1+hd))/(MS[r]))β[r] Cϕ[r](∂_r)ϕ[r])+((r^(1+hd))/(MS[r]))γ ϕ[r]Cϕ[r])/.DefCS2];
DefPiFields=FullSimplify[Solve[{pi[r]==D[ActionL,(Cϕ)'[r]],Cpi[r]==D[ActionL,ϕ'[r]]},{ϕ'[r],(Cϕ)'[r]}][[1]]]
Print["H:"]
HamiltonianH=FullSimplify[FullSimplify[pi[r](Cϕ)'[r]+Cpi[r]ϕ'[r]-ActionL/.DefPiFields]/.DefTrick]
DefNewFields={pi[r]->η[r] ϕ2[r],Cpi[r]->η[r] Cϕ2[r],ϕ[r]->-(1/(η[r]))π2[r],Cϕ[r]->-(1/(η[r])) Cπ2[r]}
expη=((hd-1)/2);
DefEta={η[r]->-CS2 (r^(expη)),η'[r]->(∂_r)(-CS2 (r^(expη)))(*((hd-1)/(2 r))η[r]*)}
Print["H2:"]
HamiltonianH2=Simplify[Simplify[Simplify[HamiltonianH+D[η[r],r](ϕ[r]Cϕ2[r]+Cϕ[r]ϕ2[r])/.DefNewFields]/.DefEta]/.DefEta]
Print["test:"]
FullSimplify[D[D[HamiltonianH2,π2[r]],Cϕ2[r]]-D[D[HamiltonianH2,Cπ2[r]],ϕ2[r]]]
DefPi2Fields=FullSimplify[Solve[{(ϕ2)'[r]==D[HamiltonianH2,Cπ2[r]],(Cϕ2)'[r]==D[HamiltonianH2,π2[r]]},{π2[r],Cπ2[r]}][[1]]]
Print["L2:"]
ActionL2=FullSimplify[π2[r](Cϕ2)'[r]+Cπ2[r](ϕ2)'[r]-HamiltonianH2/.DefPi2Fields]
Print["test:"]
FullSimplify[D[D[ActionL2,(ϕ2)'[r]],Cϕ2[r]]-D[D[ActionL2,(Cϕ2)'[r]],ϕ2[r]]]
fixL2term=D[-((α )/(2 (CS^2)))(Cϕ2[r]ϕ2[r]+ϕ2[r]Cϕ2[r]) (((β[r])/α)-FullSimplify[((η'[r])/(η[r]))/.DefEta]) ,r];
Print["L2 fixed:"]
ActionL2fixed=FullSimplify[ActionL2-fixL2term]
Print["L2 fixed no consts:"]
ActionL2fixedOpen=FullSimplify[FullSimplify[ActionL2fixed/.DefConsts]/.DefCS2]
Print["Kinetic term:"]
Simplify[D[ActionL2fixed,(ϕ2)'[r]]/.DefConsts]
Simplify[D[ActionL2fixed,(Cϕ2)'[r]]/.DefConsts]
KineticTerm=Simplify[D[D[ActionL2fixed,(ϕ2)'[r]],(Cϕ2)'[r]]/.DefConsts]
Print["Potential term:"]
PotentialTerm=Simplify[Simplify[D[D[ActionL2fixed,ϕ2[r]],Cϕ2[r]]/.DefConsts]/.DefCS2]
PotentialTerm2=-((CS2 hd(hd+1))/CS)((1/(r^2))(((((hd^2)-1))/4)+ CS)-(w^2))
FullSimplify[FreeTerm-FreeTerm2/.DefCS2]
Print["Source term:"]
SourceTerm=Simplify[Simplify[D[ActionL2fixed,CJ[r]]/.DefConsts]/.DefCS2]
Print["Source term Conjugate:"]
SourceTermC=Simplify[Simplify[D[ActionL2fixed,J[r]]/.DefConsts]/.DefCS2]
Print["test:"]
FullSimplify[SourceTerm-SourceTermC/.{CJ[r]->J[r],Cϕ2[r]->ϕ2[r],(Cϕ2)'[r]->(ϕ2)'[r]}]
fixL2termSource=D[-((hd (1+hd))/CS)(r^((hd+3)/2)) (CJ[r]ϕ2[r]+J[r]Cϕ2[r]),r];
Print["Action fixed with Source:"]
ActionL2fixedSource=Simplify[Simplify[ActionL2fixed-fixL2termSource/.DefConsts]/.DefCS2]
Print["Kinetic term fixed:"]
KineticTermFixed=Simplify[D[ActionL2fixedSource,(Cϕ2)'[r]]/.DefConsts]
Simplify[D[ActionL2fixedSource,(ϕ2)'[r]]/.DefConsts]
Print["Potential term fixed:"]
PotentialTermFixed=Simplify[Simplify[D[D[ActionL2fixedSource,ϕ2[r]],Cϕ2[r]]/.DefConsts]/.DefCS2]
Print["Field-Source Coupling:"]
CouplingTerm=FullSimplify[FullSimplify[ActionL2fixedSource-Cϕ2[r]ϕ2[r]*PotentialTermFixed-(Cϕ2)'[r]KineticTermFixed/.DefConsts]/.DefCS2]
Print["Equation with Source:"]
EqWithSource=Simplify[Simplify[D[ActionL2fixedSource,Cϕ2[r]]/.DefConsts]/.DefCS2]

Print["Zerilli to Master form:"]
Simplify[Simplify[Simplify[(((ϕ2)'[r])^2)+(((CS+(((hd^2)-1)/4))/(r^2))-(w^2))((ϕ2[r])^2)-(∂_r)((r^(1+hd+2 L)) ((∂_r)((((A[r]))^2)/2)+((1+hd+2 L)/(2r))((A[r])^2)))/.{ϕ2[r]-> (r^(L+((hd+1)/2)))A[r],(ϕ2)'[r]->(∂_r)((r^(L+((hd+1)/2)))A[r])}]/.DefCS2]/.DefCS];
DefNoSources={J[r]->0,CJ[r]->0,J'[r]->0,CJ'[r]->0};
fixtermEqCA=D[-((hd(hd+1)CS2)/CS) (r^(1+hd+2 L)) CA[r]((∂_r)A[r]+((1+hd+2 L)/(2r))A[r]),r];
Simplify[Simplify[Simplify[Simplify[ActionL2fixedSource-fixtermEqCA/.{ϕ2[r]-> (r^(L+((hd+1)/2)))A[r],(ϕ2)'[r]->(∂_r)((r^(L+((hd+1)/2)))A[r]),Cϕ2[r]-> (r^(L+((hd+1)/2)))CA[r],(Cϕ2)'[r]->(∂_r)((r^(L+((hd+1)/2)))CA[r])}]/.DefCS2]/.DefCS]/.DefNoSources]
ActionL2Canonical=FullSimplify[FullSimplify[FullSimplify[ActionL2fixedSource-fixtermEqCA/.{ϕ2[r]-> (r^(L+((hd+1)/2)))A[r],(ϕ2)'[r]->(∂_r)((r^(L+((hd+1)/2)))A[r]),Cϕ2[r]-> (r^(L+((hd+1)/2)))CA[r],(Cϕ2)'[r]->(∂_r)((r^(L+((hd+1)/2)))CA[r])}]/.DefCS2]/.DefCS]
Print["Equation with Source:"]
EqWithSourceCanonical=Simplify[Simplify[D[ActionL2Canonical,CA[r]]/.DefConsts]/.DefCS2]
EqWithSourceCanonical2=((hd (1+hd) CS2)/CS) (r^(2L+hd+1))((w^2) A[r]+A''[r]+((2L+hd+1)/r) A'[r])+((r^(L+hd+3))/CS)(((hd+1)(w^2)-((hd(CS-2(hd+1)))/(r^2))) J[r]+((hd (1+hd))/r) J'[r])
Print["test:"]
FullSimplify[FullSimplify[EqWithSourceCanonical2-EqWithSourceCanonical/.DefCS2]/.DefCS]
SourceCanonical=Simplify[EqWithSourceCanonical2/.{A[r]->0,A'[r]->0,A''[r]->0}]

J[r]=Ttt[r]+((2I w)/(r^(hd+1)))D[(hd+1)(((hd+1)(w^2)(r^2)-hd CS2)/(MS[r]))(r^(hd+3))Ttr[r],r]+((2(((hd+1))^2)I (w^3)(r^3))/(MS[r]))Ttr[r]-(((hd+1))^2)D[(1/(r^(hd+1)))D[(((w^2)(r^(hd+5)))/(MS[r]))Trr[r],r],r]-(((hd+1))^2)D[(((w^2)(r^3))/(((MS[r]))^2))(CS CS2(hd-1)(hd^2)-2hd CS2 (w^2)(r^2)+(hd+3)(w^4)(r^4))Trr[r],r]+(1/(((MS[r]))^2))((CS^2)(CS2^2)(hd^3)-CS CS2 (hd^2)(hd+1)(2CS+(hd+1)((hd^2)+hd-4))(w^2)(r^2)+CS2 hd(((hd+1))^2)(CS+2hd(hd+1)(hd+3))(w^4)(r^4)-(((hd+1))^4)(hd+2)(hd+3)(w^6)(r^6))Trr[r]+((((hd+1))^3)/(r^(hd+2)))D[(((w^2)(r^(hd+7)))/(MS[r]))Ts[r],r]+(((hd+1)CS2 (r^2)(CS hd-(hd+1)(w^2)(r^2)))/(MS[r]))Ts[r]
SourceSubT=FullSimplify[((r^(L+hd+3))/CS)(((hd+1)(w^2)-((hd(CS-2(hd+1)))/(r^2))) J[r]+((hd (1+hd))/r) J'[r])/DefConsts]