(*^

::[paletteColors = 128; 
	fontset = title, "New York", 24, L3, center, bold, nohscroll;
	fontset = subtitle, "New York", 18, L2, center, bold, nohscroll;
	fontset = subsubtitle, "New York", 14, L2, center, bold, nohscroll;
	fontset = section, "New York", 14, L2, bold, nohscroll, grayBox;
	fontset = subsection, "New York", 12, L2, bold, nohscroll, blackBox;
	fontset = subsubsection, "New York", 10, L2, bold, nohscroll, whiteBox;
	fontset = text, "New York", 12, L2, nohscroll;
	fontset = smalltext, "New York", 10, L2, nohscroll;
	fontset = input, "Courier", 12, L2, bold, nowordwrap;
	fontset = output, "Courier", 12, L2, nowordwrap;
	fontset = message, "Courier", 12, L2, R65535, nowordwrap;
	fontset = print, "Courier", 12, L2, nowordwrap;
	fontset = info, "Courier", 12, L2, nowordwrap;
	fontset = postscript, "Courier", 12, L2, nowordwrap;
	fontset = name, "Geneva", 10, L2, italic, B65535, nowordwrap, nohscroll;
	fontset = header, "Times", 10, L2;
	fontset = footer, "Times", 12, L2, center;
	fontset = help, "Geneva", 10, L2, nohscroll;
	fontset = clipboard, "New York", 12, L2;
	fontset = completions, "New York", 12, L2, nowordwrap;
	fontset = network, "Courier", 10, L2, nowordwrap;
	fontset = graphlabel, "Courier", 12, L2, nowordwrap;
	fontset = special1, "New York", 12, L2, nowordwrap;
	fontset = special2, "New York", 12, L2, center, nowordwrap;
	fontset = special3, "New York", 12, L2, right, nowordwrap;
	fontset = special4, "New York", 12, L2, nowordwrap;
	fontset = special5, "New York", 12, L2, nowordwrap;]
:[font = input; initialization; ]
*)
(* Package written by 
      Pekka Janhunen
      Finnish Meteorological Institute
      Geophysics Dept.
      pjanhune@finsun.csc.fi  *)

BeginPackage["EinsteinTensor`"]

EinsteinTensor::usage = "EinsteinTensor[g,x] with g a nxn-matrix
  (the metric with lower indices) and x n-vector (the coordinates)
  gives the Einstein tensor (a nxn-matrix) with lower indices."

Begin["`Private`"]

EinsteinTensor[metric_,x_]:=
  Block[ {Dim,Metric, PreChristoffel, Christoffel, Riemann,
          PreRiemann, Ricci, CurvatureScalar,
          sigma, mu, nu, alpha, beta, gamma},
          Dim = Length[x];
          Metric = Simplify[Inverse[metric]];
            (* Metric with upper indices *)
          PreChristoffel =
            Table[ D[metric[[gamma,alpha]],x[[beta]]]
                 + D[metric[[beta,gamma]],x[[alpha]]]
           		    - D[metric[[alpha,beta]],x[[gamma]]],
           	   	 {gamma,Dim}, {alpha,Dim}, {beta,Dim} ];
           	 (* The "lower index part" of Christoffel symbols *)
          PreChristoffel = Simplify[PreChristoffel];
          Christoffel = (1/2) Metric . PreChristoffel;
             (* The full Christoffel symbols *)
          Christoffel = Simplify[Christoffel];
          PreRiemann = 
             Table[ D[Christoffel[[sigma,alpha,nu]],x[[mu]]]
                    + Sum[Christoffel[[gamma,alpha,nu]]
                            Christoffel[[sigma,gamma,mu]],
                          {gamma,Dim} ],
                    {sigma,Dim}, {alpha,Dim}, {mu,Dim}, {nu,Dim} ];
           	(* PreRiemann has to be antisymmetrized to yield
           	   Riemann tensor: *)
          Riemann = Table[ PreRiemann[[sigma,alpha,mu,nu]]
                         - PreRiemann[[sigma,alpha,nu,mu]],
                           {sigma,Dim}, {alpha,Dim},
                           {mu,Dim}, {nu,Dim} ];
          Ricci = Table[ Sum[Riemann[[sigma,alpha,sigma,beta]],
                             {sigma,Dim}],
                         {alpha,Dim}, {beta,Dim} ];
          CurvatureScalar = Sum[ Metric[[alpha,beta]]
                                 Ricci[[alpha,beta]],
                                 {alpha,Dim}, {beta,Dim} ];
          (* Return Einstein tensor: *)
          Ricci - (1/2) CurvatureScalar metric ]

End[]

EndPackage[]

Print[{EinsteinTensor}]
(*
;[s]
2:0,1;2289,0;2290,-1;
2:1,14,10,Courier,1,12,0,0,0;1,12,9,Courier,1,10,0,0,0;
^*)