-- metric tensor gdd = zero(4,4) gdd[1,1] = -exp(2 Phi(r)) gdd[2,2] = exp(2 Lambda(r)) gdd[3,3] = r^2 gdd[4,4] = r^2 sin(theta)^2 X = (t,r,theta,phi) -- for computing gradients -- Step 1. Calculate guu guu = inv(gdd) -- Step 2. Calculate connection coefficients ("Gravitation" by MTW p. 210) -- -- Gamma = 1/2 (g + g - g ) -- abc ab,c ac,b bc,a gddd = d(gdd,X) GAMDDD = 1/2 (gddd + transpose(gddd,2,3) - transpose(gddd,2,3,1,2)) -- transpose bc,a to (,a)bc GAMUDD = dot(guu,GAMDDD) -- raise first index -- Step 3. Calculate Riemann tensor (MTW p. 219) -- -- a a a a u a u -- R = Gamma - Gamma + Gamma Gamma - Gamma Gamma -- bcd bd,c bc,d uc bd ud bc GAMUDDD = d(GAMUDD,X) GAMGAM = dot(transpose(GAMUDD,2,3),GAMUDD) RUDDD = transpose(GAMUDDD,3,4) - GAMUDDD + transpose(GAMGAM,2,3) - transpose(GAMGAM,2,3,3,4) -- Step 4. Calculate Ricci tensor (MTW p. 343) -- -- a -- R = R -- uv uav RDD = contract(RUDDD,1,3) -- Step 5. Calculate Ricci scalar (MTW p. 343) -- -- u -- R = R -- u R = contract(dot(guu,RDD)) -- Step 6. Calculate Einstein tensor (MTW p. 343) -- -- G = R - 1/2 g R -- uv uv uv GDD = RDD - 1/2 gdd R -- Check GDD Gtt = exp(2 Phi(r)) d(r (1 - exp(-2 Lambda(r))),r) / r^2 Grr = (1 - exp(2 Lambda(r))) / r^2 + 2 d(Phi(r),r) / r Gthetatheta = r^2 exp(-2 Lambda(r)) ( d(d(Phi(r),r),r) + d(Phi(r),r)^2 + d(Phi(r),r) / r - d(Phi(r),r) d(Lambda(r),r) - d(Lambda(r),r) / r) Gphiphi = Gthetatheta sin(theta)^2 G = zero(4,4) G[1,1] = Gtt G[2,2] = Grr G[3,3] = Gthetatheta G[4,4] = Gphiphi "Verify Einstein tensor" check(GDD == G) "ok" "Non-zero components of Einstein tensor" Gtt Grr Gthetatheta Gphiphi Run