(*^

::[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", 10, L2, nohscroll;
	fontset = smalltext, "New York", 9, L2, nohscroll;
	fontset = input, "Courier", 10, L2, bold, nowordwrap;
	fontset = output, "Courier", 10, L2, nowordwrap;
	fontset = message, "Courier", 10, L2, R65535, nowordwrap;
	fontset = print, "Courier", 10, L2, nowordwrap;
	fontset = info, "Courier", 10, 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; startGroup; ]
<<einsteintensor.m
:[font = print; inactive; endGroup; ]
{EinsteinTensor}
:[font = input; startGroup; ]
(* Spherical coordinates: *)
x = {t,r,theta,phi}
:[font = output; inactive; output; endGroup; ]
{t, r, theta, phi}
;[o]
{t, r, theta, phi}
:[font = input; startGroup; ]
(* The Schwarzschild metric: *)
(metric =
  DiagonalMatrix[{-Exp[2 ny[r]], Exp[2 lambda[r]], r^2,
                  (r Sin[theta])^2}]) //MatrixForm
:[font = output; inactive; output; endGroup; ]
MatrixForm[{{-E^(2*ny[r]), 0, 0, 0}, {0, E^(2*lambda[r]), 0, 0}, 
   {0, 0, r^2, 0}, {0, 0, 0, r^2*Sin[theta]^2}}]
;[o]
  2 ny[r]
-E              0               0               0
 
                 2 lambda[r]
0               E               0               0
 
                                 2
0               0               r               0
 
                                                 2           2
0               0               0               r  Sin[theta]
:[font = input; startGroup; ]
(* The following result is found in all textbooks on General
   Relativity. When the differential equations Einstein==0
   are solved, the Schwarzschild solution emerges. This cannot
   (yet) be done using Mathematica's DSolve (one gets the message
   DSolve::NotYet!). However, the solution of the equation is
   relatively easy to do by hand, the toughest part is to 
   calculate the Einstein tensor.*)
(Einstein = 
   Inverse[metric].Simplify[EinsteinTensor[metric,x]])
:[font = output; inactive; output; endGroup; ]
{{-((-1 + E^(2*lambda[r]) + 2*r*Derivative[1][lambda][r])/
      (E^(2*lambda[r])*r^2)), 0, 0, 0}, 
  {0, -((-1 + E^(2*lambda[r]) - 2*r*Derivative[1][ny][r])/
      (E^(2*lambda[r])*r^2)), 0, 0}, 
  {0, 0, (-Derivative[1][lambda][r] + Derivative[1][ny][r] - 
      r*Derivative[1][lambda][r]*Derivative[1][ny][r] + 
      r*Derivative[1][ny][r]^2 + r*Derivative[2][ny][r])/
    (E^(2*lambda[r])*r), 0}, {0, 0, 0, 
   (-Derivative[1][lambda][r] + Derivative[1][ny][r] - 
      r*Derivative[1][lambda][r]*Derivative[1][ny][r] + 
      r*Derivative[1][ny][r]^2 + r*Derivative[2][ny][r])/
    (E^(2*lambda[r])*r)}}
;[o]
          2 lambda[r]
    -1 + E            + 2 r lambda'[r]
{{-(----------------------------------), 0, 0, 0}, 
              2 lambda[r]  2
             E            r
 
              2 lambda[r]
        -1 + E            - 2 r ny'[r]
  {0, -(------------------------------), 0, 0}, 
                2 lambda[r]  2
               E            r
 
  {0, 0, 
 
                                                         2
    -lambda'[r] + ny'[r] - r lambda'[r] ny'[r] + r ny'[r]  + r ny''[r]
    ------------------------------------------------------------------, 
                               2 lambda[r]
                              E            r
 
   0}, {0, 0, 0, 
 
                                                         2
    -lambda'[r] + ny'[r] - r lambda'[r] ny'[r] + r ny'[r]  + r ny''[r]
    ------------------------------------------------------------------}}
                               2 lambda[r]
                              E            r
:[font = input; startGroup; ]
(* Confirm that the Einstein tensor of the Schwarzschild
   solution vanishes (i.e. that the Schwarzschild solution
   satisfies Einstein's vacuum field equations): *)
EinsteinTensor[DiagonalMatrix[{-(1-2M/r),1/(1-2M/r),
                                r^2, (r Sin[theta])^2}],
               {t,r,theta,phi}] // Simplify
:[font = output; inactive; output; endGroup; ]
{{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}
;[o]
{{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}
:[font = input; startGroup; ]
(* Check that the Einstein tensor of flat spacetime vanishes: *)
EinsteinTensor[DiagonalMatrix[{-1,1,1,1}],{x1,x2,x3,x4}]
:[font = output; inactive; output; endGroup; ]
{{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}
;[o]
{{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}
:[font = input; ]
(* Then, consider the Kerr-Newman metric (a charged,
   rotating black hole): *)
delta:= r^2 - 2M r + a^2 + Q^2;
rho2:= r^2 + a^2 Cos[theta]^2;
ds2 = -(delta/rho2)(d[0]-a Sin[theta]^2 d[3])^2 +
      (Sin[theta]^2/rho2)((r^2+a^2)d[3]-a d[0])^2 +
      rho2 d[1]^2/delta + rho2 d[2]^2;
:[font = input; startGroup; ]
ds2 = Collect[ds2,Array[d,4,0]]
:[font = output; inactive; output; endGroup; ]
(-(Q^2/(r^2 + a^2*Cos[theta]^2)) - a^2/(r^2 + a^2*Cos[theta]^2) + 
     (2*M*r)/(r^2 + a^2*Cos[theta]^2) - r^2/(r^2 + a^2*Cos[theta]^2) + 
     (a^2*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2))*d[0]^2 + 
  (r^2/(Q^2 + a^2 - 2*M*r + r^2) + 
     (a^2*Cos[theta]^2)/(Q^2 + a^2 - 2*M*r + r^2))*d[1]^2 + 
  (r^2 + a^2*Cos[theta]^2)*d[2]^2 + 
  ((2*Q^2*a*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2) - 
     (4*M*a*r*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2))*d[0]*d[3] + 
  ((a^4*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2) + 
     (2*a^2*r^2*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2) + 
     (r^4*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2) - 
     (Q^2*a^2*Sin[theta]^4)/(r^2 + a^2*Cos[theta]^2) - 
     (a^4*Sin[theta]^4)/(r^2 + a^2*Cos[theta]^2) + 
     (2*M*a^2*r*Sin[theta]^4)/(r^2 + a^2*Cos[theta]^2) - 
     (a^2*r^2*Sin[theta]^4)/(r^2 + a^2*Cos[theta]^2))*d[3]^2
;[o]
            2                      2
           Q                      a                    2 M r
(-(-------------------) - ------------------- + ------------------- - 
    2    2           2     2    2           2    2    2           2
   r  + a  Cos[theta]     r  + a  Cos[theta]    r  + a  Cos[theta]
 
              2               2           2
             r               a  Sin[theta]          2
     ------------------- + -------------------) d[0]  + 
      2    2           2    2    2           2
     r  + a  Cos[theta]    r  + a  Cos[theta]
 
             2                2           2
            r                a  Cos[theta]          2
  (-------------------- + --------------------) d[1]  + 
    2    2            2    2    2            2
   Q  + a  - 2 M r + r    Q  + a  - 2 M r + r
 
    2    2           2      2
  (r  + a  Cos[theta] ) d[2]  + 
 
      2             2                      2
   2 Q  a Sin[theta]     4 M a r Sin[theta]
  (------------------- - -------------------) d[0] d[3] + 
    2    2           2    2    2           2
   r  + a  Cos[theta]    r  + a  Cos[theta]
 
      4           2         2  2           2      4           2
     a  Sin[theta]       2 a  r  Sin[theta]      r  Sin[theta]
  (------------------- + ------------------- + ------------------- - 
    2    2           2    2    2           2    2    2           2
   r  + a  Cos[theta]    r  + a  Cos[theta]    r  + a  Cos[theta]
 
       2  2           4       4           4           2             4
      Q  a  Sin[theta]       a  Sin[theta]       2 M a  r Sin[theta]
     ------------------- - ------------------- + -------------------- - 
      2    2           2    2    2           2    2    2           2
     r  + a  Cos[theta]    r  + a  Cos[theta]    r  + a  Cos[theta]
 
       2  2           4
      a  r  Sin[theta]        2
     -------------------) d[3]
      2    2           2
     r  + a  Cos[theta]
:[font = input; startGroup; ]
KerrNewman = 
  Simplify[Table[Coefficient[ds2,d[alpha]d[beta]],
                 {alpha,0,3}, {beta,0,3} ] ]
:[font = output; inactive; output; endGroup; ]
{{-((Q^2 + a^2 - 2*M*r + r^2 - a^2*Sin[theta]^2)/
      (r^2 + a^2*Cos[theta]^2)), 0, 0, 
   (2*a*(Q^2 - 2*M*r)*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2)}, 
  {0, (r^2 + a^2*Cos[theta]^2)/(Q^2 + a^2 - 2*M*r + r^2), 0, 0}, 
  {0, 0, r^2 + a^2*Cos[theta]^2, 0}, 
  {(2*a*(Q^2 - 2*M*r)*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2), 0, 0, 
   -((Sin[theta]^2*(-a^4 - 2*a^2*r^2 - r^4 + Q^2*a^2*Sin[theta]^2 + 
          a^4*Sin[theta]^2 - 2*M*a^2*r*Sin[theta]^2 + 
          a^2*r^2*Sin[theta]^2))/(r^2 + a^2*Cos[theta]^2))}}
;[o]
     2    2            2    2           2
    Q  + a  - 2 M r + r  - a  Sin[theta]
{{-(-------------------------------------), 0, 0, 
              2    2           2
             r  + a  Cos[theta]
 
         2                    2        2    2           2
   2 a (Q  - 2 M r) Sin[theta]        r  + a  Cos[theta]
   ----------------------------}, {0, --------------------, 0, 0}, 
        2    2           2             2    2            2
       r  + a  Cos[theta]             Q  + a  - 2 M r + r
 
          2    2           2
  {0, 0, r  + a  Cos[theta] , 0}, 
 
         2                    2
   2 a (Q  - 2 M r) Sin[theta]
  {----------------------------, 0, 0, 
        2    2           2
       r  + a  Cos[theta]
 
                2    4      2  2    4    2  2           2
   -((Sin[theta]  (-a  - 2 a  r  - r  + Q  a  Sin[theta]  + 
 
           4           2        2             2    2  2           2
          a  Sin[theta]  - 2 M a  r Sin[theta]  + a  r  Sin[theta] )) / 
 
        2    2           2
      (r  + a  Cos[theta] ))}}
:[font = input; startGroup; ]
(* Consider the limiting cases the Kerr-Newman metric:
   The Kerr metric is a non-charged Kerr-Newman: *)
Kerr = Simplify[Limit[KerrNewman,Q->0]]
:[font = output; inactive; output; endGroup; ]
{{-((a^2 - 2*M*r + r^2 - a^2*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2)), 0, 
   0, (-4*M*a*r*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2)}, 
  {0, (r^2 + a^2*Cos[theta]^2)/(a^2 - 2*M*r + r^2), 0, 0}, 
  {0, 0, r^2 + a^2*Cos[theta]^2, 0}, 
  {(-4*M*a*r*Sin[theta]^2)/(r^2 + a^2*Cos[theta]^2), 0, 0, 
   -((Sin[theta]^2*(-a^4 - 2*a^2*r^2 - r^4 + a^4*Sin[theta]^2 - 
          2*M*a^2*r*Sin[theta]^2 + a^2*r^2*Sin[theta]^2))/
      (r^2 + a^2*Cos[theta]^2))}}
;[o]
     2            2    2           2                            2
    a  - 2 M r + r  - a  Sin[theta]          -4 M a r Sin[theta]
{{-(--------------------------------), 0, 0, --------------------}, 
           2    2           2                 2    2           2
          r  + a  Cos[theta]                 r  + a  Cos[theta]
 
       2    2           2
      r  + a  Cos[theta]                  2    2           2
  {0, -------------------, 0, 0}, {0, 0, r  + a  Cos[theta] , 0}, 
         2            2
        a  - 2 M r + r
 
                      2
   -4 M a r Sin[theta]
  {--------------------, 0, 0, 
    2    2           2
   r  + a  Cos[theta]
 
                2    4      2  2    4    4           2
   -((Sin[theta]  (-a  - 2 a  r  - r  + a  Sin[theta]  - 
 
               2             2    2  2           2
          2 M a  r Sin[theta]  + a  r  Sin[theta] )) / 
 
        2    2           2
      (r  + a  Cos[theta] ))}}
:[font = input; startGroup; ]
(* And the Schwarzschild metric is the same as nonrotating
   Kerr: *)
Schwarzschild = Simplify[Limit[Kerr, a->0]]
:[font = output; inactive; output; endGroup; ]
{{-1 + (2*M)/r, 0, 0, 0}, {0, r/(-2*M + r), 0, 0}, {0, 0, r^2, 0}, 
  {0, 0, 0, r^2*Sin[theta]^2}}
;[o]
       2 M                   r                     2
{{-1 + ---, 0, 0, 0}, {0, --------, 0, 0}, {0, 0, r , 0}, 
        r                 -2 M + r
 
             2           2
  {0, 0, 0, r  Sin[theta] }}
:[font = input; startGroup; ]
MaxMemoryUsed[]
:[font = output; inactive; output; endGroup; ]
1531008
;[o]
1531008
:[font = input; ]
(* Warning: This calculation takes very long time and requires
   plenty of memory. *)
Einstein = EinsteinTensor[KerrNewman,{t,r,theta,phi}] ;
:[font = input; startGroup; ]
MaxMemoryUsed[]
:[font = output; inactive; output; endGroup; ]
4950592
;[o]
4950592
:[font = input; startGroup; ]
(* The expression Einstein should be zero! It is far too long
   to be simplified using Simplify! Even computing numerical
   values takes too much time on my computer. By adding 
   extra Simplify -calls in the intermediate stages of
   EinsteinTensor[], I succeeded to verify that EinsteinTensor
   vanishes for the Kerr metric (not documented in this file)*)
Einstein//Short
:[font = output; inactive; output; endGroup; ]
-interrupted-
;[o]
              r <<3>>
{{<<4>>}, {0, ------- + <<20>>, <<1>>, 0}, <<1>>, {<<4>>}}
              2 <<1>>
:[font = input; ]
(* the following calculation was interrupted by user.. *)
Einstein /. {a->1.2,Q->0.001,M->15.34, r->4.5345,theta->0.34}
^*)