BeginPackage["EinsteinVariation`"];

MetricFromDS::usage = "MetricFromDS[ds, coords] - generates g_{uv} from a 
line element."

CC::usage = "CC[gll,guu,X,l,u,v] - returns Christoffel^l_{uv} element for the input metric gll (=g_{uv} and inverse guu = g^{uv}) and coordinate system X."

Riemann::usage = "Riemann[gll,guu,X,s,p,u,v] - returns R^s_{puv} element for input metric gll (=g_{uv} and inverse guu = g^{uv}) and coordinate system X."

GetFieldEquations::usage = "GetFieldEquations[DS,ab,X, NonKill] - returns field equations for the list of functions 'ab', the line element DS, coordinates X, and dependent coordinates NonKill."

GetSchwarzschildgll::usage = "GetSchwarzschildgll[M,r,T] - returns a matrix representation of the Schwarzshild metric (g_{uv}) with coordinates printed as r, T, and mass parameter M."

GetKerrgll::usage = "GetKerrgll[M,a,r,T] - returns a matrix representation of the Kerr metric (g_{uv}) with coordinates printed as r and T, mass parameter M and spin a."

GetRicci::usage = "GetRicci[gll,X] - returns R_{uv} for a given metric gll = g_{uv} and coordinates X."

GetRiemann::usage = "GetRiemann[gll,X] - returns R^u_{vab} given metric gll = g_{uv} and coordinates X."

GetCCull::usage="GetCCull[gll,X] - returns Christoffel G^i_{jk} given metric gll = g_{uv} and coordinates X."

Sym::usage = "Sym[T...., rank] - returns the fully symmetrized tensor T(.....) with of the given rank."

TF::usage = "TF[T(.....), rank, dim] - returns the completely tracefree (T_{aa} not T_a^a) tensor of the given rank and dimension given a symmetrized tensor T(.....)."

TFScalar::usage = "TFScalar[T(.....), rank, dim, gll] - returns the completely tracefree (T^a_a not T_{aa}) tensor of the given rank and dimension starting from a symmetrized T(.....) and using the metric gll (g_{uv} in matrix form)."

SymTF::usage = "SymTF[T...., rank, dim] - returns the totally symmetric, tracefree (T_{aa}) part of T...., a tensor of the given rank and dimension."

SymTFScalar::usage="SymTFScalar[T....., rank, dim, gll] - returns the totally symmetric, tracefree (T^a_a) part of T....., a tensor of the given rank and dimension, using gll (=g_{uv} in matrix form)."

GetRicciS::usage = "GetRicciS[gll,X] - returns the Ricci scalar R for a given metric gll = g_{uv} and coordinates X."

Begin["`Private`"]

Needs["Calculus`VariationalMethods`"];

GetRicciS[gll_,X_] := Module[{Rll, guu, Rout},
	Rll = Simplify[GetRicci[gll,X]];
	guu = Simplify[Inverse[gll]];
	Rout = Sum[Rll[[aa,bb]] guu[[aa,bb]],{aa,1,Length[X]},{bb,1,Length[X]}];
	Rout
]

SymTF[M_, rank_, dim_] := TF[Sym[M, rank], rank, dim]

SymTFScalar[M_, rank_, dim_, gll_] := TFScalar[Sym[M, rank], rank, dim,gll]


TFScalar[M_, rank_, dim_, 
	 gll_] := Module[{ret, index, Metricsll, Metricsuu, Mpart,
    tindex, dtindex, sindex, tpart, guu},
    Metricsll = gll;
    guu = Simplify[Inverse[gll]];
    Metricsuu = guu;
    ret = aa[0, rank] M;
    For[index = 1, index <= Floor[rank/2], index++,
          dtindex = Table[k[i], {i, 1, 2 index}];
          tindex = Table[k[2 index + i], {i, 1, rank - 2 index}];
          sindex = Table[{j[2*i - 1], j[2*i]}, {i, 1, index}];
          (* Print["Delta indices = ", dtindex];
            Print["Tensor indices = ", tindex];
	   Print["Summation indices = ", sindex];*)
          tpart = Hold[Apply[Table, Table[If[i == 0, Extract[Metricsll, dtindex] Extract[M, Flatten[Join[tindex, sindex]]] Extract[Metricsuu, Flatten[sindex]], {k[i], 1, dim}], {i, 0, rank}]]];
          Mpart = 
            Apply[Sum, Table[If[i == 0, ReleaseHold[tpart], {j[i], 1, dim}], {i, 0, 2*index}]];
          Metricsll = Outer[Times, Metricsll, gll];
          Metricsuu = Outer[Times, Metricsuu, guu];
          (*Print["Mp = ", Sym[Mpart, rank]];*)
          ret = ret + aa[index, rank] Sym[Mpart, rank];
          ] ret //. Null -> 1
    ]


TF[M_, rank_, 
    dim_] := Module[{ret, index, Deltas, Delta, Mpart, tindex,
       dtindex, sindex, tpart},
    Delta = Table[KroneckerDelta[a, b], {a, 1, dim}, {b, 1, dim}];
    Deltas = Delta;
    ret = aa[0, rank] M;
    (* Print[ret]; *)
    For[index = 1, index ² Floor[rank/2], index++,
          dtindex = Table[k[i], {i, 1, 2 index}];
          tindex = Table[k[2 index + i], {i, 1, rank - 2 index}];
          sindex = Table[{j[i], j[i]}, {i, 1, index}];
          (* Print["Delta indices = ", dtindex];
            Print["Tensor indices = ", tindex];
            Print["Summation indices = ", sindex]; *)
          tpart = Hold[Apply[Table, Table[If[i == 0, 
				 Extract[Deltas, dtindex ]Extract[M, Flatten[Join[tindex, sindex]]], {k[i], 1, dim}], {i, 0, rank }] ]];
          Mpart = Apply[Sum, Table[If[i == 0, 
          ReleaseHold[tpart], {j[i], 1, dim}], {i, 0, index}]];
          Deltas = Outer[Times, Deltas, Delta];
          (*  Print["Mp = ", Sym[Mpart, rank]]; *)
	ret = ret + aa[index, rank] Sym[Mpart, rank];      
          ] ret //. Null -> 1
    ]

Sym[M_, rank_] := Module[{perms, ret},
    perms = Permutations[Table[j, {j, 1, rank}]];
    ret = (1/rank!)Sum[Transpose[M, perms[[j]]], {j, 1, Length[perms]}];
    ret
    ]

aa[n_, ll_] := (-1)^n (ll!(2 ll - 2 n - 1)!!)/((ll - 2n)!(2 ll - 1)!! (2 n)!!)

MetricFromDS[ds_, coords_] := Module[{gllout, dsl, dcoords, cterms, i, j},
    dcoords = Table[ToExpression["d" <> ToString[coords[[j]]]], {j, 1, Length[coords]}];
    (* dsl = Collect[Expand[ds], dcoords]; *)
    dsl = Expand[ds];				     
    cterms = Table[{dcoords[[i]] dcoords[[j]] -> ToExpression[ToString[dcoords[[i]]] <> ToString[dcoords[[j]]]], dcoords[[j]] dcoords[[i]] -> ToExpression[ToString[dcoords[[i]]] <> ToString[dcoords[[j]]]]}, {i, 1, Length[dcoords]}, {j, i, Length[dcoords]}];
    cterms = Flatten[cterms];
    dsl = dsl //. cterms;

    gllout = Table[0, {i, 1, Length[coords]}, {j, 1, Length[coords]}];
    For[i = 1, i <= Length[coords], i++,
      For[j = 1, j <= Length[coords], j++,
            mterm = Coefficient[dsl,  
              ToExpression[ToString[dcoords[[i]]] <>ToString[dcoords[[j]]]]];
          If[mterm == 0,
             mterm =  Coefficient[dsl,  ToExpression[ToString[dcoords[[j]]] <> ToString[dcoords[[i]]]]];
            ];
          If[i ­ j,
            gllout[[i, j]] = mterm/2 , gllout[[i, j]] = mterm
            ];
          ];
      ];
			       
    gllout
    ]

CC[gll_,guu_,X_, l_, u_, v_] := (1/2)Sum[guu[[l, ppJJ]](D[gll[[ppJJ, u]], X[[v]]] + D[gll[[ppJJ, v]], X[[u]]] - D[gll[[u, v]], X[[ppJJ]]]), {ppJJ, 1, Length[X]}]

GetCCull[gll_, X_] := Module[{guu},
			     guu = Simplify[Inverse[gll]];
			     Table[Simplify[CC[gll,guu,X,llJJ, uuJJ, vvJJ]], {llJJ, 1, Length[X]}, {uuJJ, 1, Length[X]}, {vvJJ, 1, Length[X]}]
]

Riemann[gll_,guu_, X_, ss_, pp_, uu_, vv_] := -D[CC[gll,guu,X,ss, uu, pp], X[[vv]]] + D[CC[gll,guu,X,ss, vv, pp], X[[uu]]] - Sum[CC[gll,guu,X,aaJJ, uu, pp] CC[gll,guu,X,ss, aaJJ, vv] - CC[gll,guu,X,aaJJ, vv, pp] CC[gll,guu,X,ss, aaJJ, uu], {aaJJ, 1, Length[X]}] 

GetFieldEquations[DS_,ab_,X_, NonKill_] := Module[{gll, guu, CCull,Rulll, Rill, R,SR, dg, retfields},
    gll = FullSimplify[MetricFromDS[DS,X]];
    guu = FullSimplify[Inverse[gll]];
    CCull = Table[Simplify[CC[gll,guu,X,llJJ, uuJJ, vvJJ]], {llJJ, 1, Length[X]}, {uuJJ, 1, Length[X]}, {vvJJ, 1, Length[X]}];
    Rulll = Table[Riemann[gll,guu,X,ii, jj, kk, ll], {ii, 1, Length[X]}, {jj, 1, Length[X]}, {kk, 1, Length[X]}, {ll, 1, Length[X]}];
    Rill = Table[Sum[Rulll[[aaJJ,bbJJ,aaJJ,ccJJ]],{aaJJ,1,Length[X]}],{bbJJ,1,Length[X]},{ccJJ,1,Length[X]}];
    R = Sum[Rill[[aaJJ,bbJJ]] guu[[bbJJ,aaJJ]],{aaJJ,1,Length[X]},{bbJJ,1,Length[X]}];
    dg = FullSimplify[Sqrt[-Det[gll]]];
    SR = FullSimplify[dg R]; 
    retfields = Table[VariationalD[SR,ab[[jjJJ]],NonKill],{jjJJ,1,Length[ab]}];
    {retfields, R, SR, Rill}
]

GetRicci[gll_,X_] := Module[{ guu, CCull,Rulll, Rill},
    guu = FullSimplify[Inverse[gll]];
    CCull = Table[Simplify[CC[gll,guu,X,llJJ, uuJJ, vvJJ]], {llJJ, 1, Length[X]}, {uuJJ, 1, Length[X]}, {vvJJ, 1, Length[X]}];
    Rulll = Table[Riemann[gll,guu,X,ii, jj, kk, ll], {ii, 1, Length[X]}, {jj, 1, Length[X]}, {kk, 1, Length[X]}, {ll, 1, Length[X]}];
    Rill = Table[Sum[Rulll[[aaJJ,bbJJ,aaJJ,ccJJ]],{aaJJ,1,Length[X]}],{bbJJ,1,Length[X]},{ccJJ,1,Length[X]}];
    Rill 
]

GetRiemann[gll_,X_] := Module[{ guu, CCull,Rulll, Rill},
    guu = FullSimplify[Inverse[gll]];
    CCull = Table[Simplify[CC[gll,guu,X,llJJ, uuJJ, vvJJ]], {llJJ, 1, Length[X]}, {uuJJ, 1, Length[X]}, {vvJJ, 1, Length[X]}];
    Rulll = Table[Riemann[gll,guu,X,ii, jj, kk, ll], {ii, 1, Length[X]}, {jj, 1, Length[X]}, {kk, 1, Length[X]}, {ll, 1, Length[X]}];
    Rulll
]


GetSchwarzschildgll[M_,r_,T_] := {{-1 + 2 M/r,0,0,0},{0,1/(1 - 2 M/r),0,0},{0,0,r^2,0},{0,0,0,r^2 Sin[T]^2}}

GetKerrgll[M_,a_,r_,T_] := FullSimplify[{{-1 + 2 M r/(a^2 Cos[T]^2 + r^2),0,0,-2 a M r Sin[T]^2/(r^2 + a^2 Cos[T]^2)},
			      {0, (a^2 Cos[T]^2 + r^2)/(a^2 - 2 M r + r^2),0,0},
			      {0,0,a^2 Cos[T]^2 + r^2,0},
			      {-2 a M r Sin[T]^2/(r^2 + a^2 Cos[T]^2),0,0,(Sin[T]^2 ( (a^2 + r^2)^2 - a^2 (a^2 - 2 M r + r^2) Sin[T]^2))/(a^2 Cos[T]^2 + r^2)}}]


End[ ]
EndPackage[ ]