-- gdd is the metric tensor gdd = ((-exp(2*Phi(r)), 0, 0, 0), ( 0, exp(2*Lambda(r)), 0, 0), ( 0, 0, r^2, 0), ( 0, 0, 0, r^2*sin(theta)^2)) -- Note: "dd" stands for two "down" indices, "uu" stands for two "up" indices. -- X is for computing gradients. X = (t,r,theta,phi) -- Step 1: Calculate guu. guu = inv(gdd) -- Step 2: Calculate the connection coefficients. ("Gravitation" by MTW p. 210) -- -- Gamma = 1/2 (g + g - g ) -- abc ab,c ac,b bc,a -- -- Note: The comma means gradient which increases the rank of gdd by 1. gddd = d(gdd,X) -- Transpose indices to match abc. GAMDDD = 1/2 * (gddd + -- indices are already in correct order transpose(gddd,2,3) - -- transpose c and b transpose(gddd,2,3,1,2)) -- transpose c and a, then b and a -- Raise first index. -- -- a au -- Gamma = g Gamma -- bc ubc -- -- Note: Sum over index u means contraction. GAMUDD = dot(guu,GAMDDD) -- Step 3. Calculate the Riemann tensor. ("Gravitation" by MTW p. 219) -- -- a is alpha -- b is beta -- c is gamma -- d is delta -- u is mu -- -- a a a a u a u -- R = Gamma - Gamma + Gamma Gamma - Gamma Gamma -- bcd bd,c bc,d uc bd ud bc -- -- Do the gradient once and save in a temporary variable. T1 = d(GAMUDD,X) -- T2 is the product Gamma Gamma contracted over u. T2 = dot(transpose(GAMUDD,2,3),GAMUDD) -- Now put it all together. Transpose indices to match abcd. RUDDD = transpose(T1,3,4) - -- transpose d and c T1 + -- already in correct order transpose(T2,2,3) - -- transpose c and b transpose(T2,2,3,3,4) -- transpose d and b, then d and c -- Step 4: Calculate the Ricci tensor. ("Gravitation" by MTW p. 343) -- -- a -- R = R -- uv uav -- -- Contract over "a" (1st and 3rd indices). RDD = contract(RUDDD,1,3) -- Step 5: Calculate the Ricci scalar. ("Gravitation" by MTW p. 343) -- -- uv -- R = g R -- uv R = contract(dot(guu,transpose(RDD))) -- Step 6: Finally, calculate the Einstein tensor. ("Gravitation" by 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 "Non-zero components of Einstein tensor" Gtt Grr Gthetatheta Gphiphi T = ((Gtt, 0, 0, 0), ( 0, Grr, 0, 0), ( 0, 0, Gthetatheta, 0), ( 0, 0, 0, Gphiphi)) "Verify Einstein tensor (1=ok)" GDD == T