(* The template code for computing the Curvature tensors and the Geodesic equations. It can be copied/pasted into a Wolfram Cloud notebook *)
(* wolframscript -file U02_ComputerLab_CurvatureWithMathematicaTemplate_Mathematica.m > out.m 
   Then copy/paste all the output in out.m into a notbook cell, and evaluate. You will see the results printed nicely. *)

Clear[coord, metric, inversemetric, affine, riemann, lriemann, uriemann, ricci, scalar, einstein, weyl, geodesic, R, G, \[Tau], i, j, k, l, s];
Clear[r, \[Theta], \[Phi], t, \[Chi], a, m]; 

(*----------------------------------------------------------------------------*)

(* This is what you need to set: *)

coord  = {t, \[Chi], \[Theta], \[Phi]};
n      = Length[coord];
metric = {
   {-1, 0      , 0                    ,                                    0},
   { 0, a[t]^2 , 0                    ,                                    0},
   { 0, 0      , a[t]^2 Sin[\[Chi]]^2 ,                                    0},
   { 0, 0      , 0                    , a[t]^2 Sin[\[Chi]]^2 Sin[\[Theta]]^2}};
(*----------------------------------------------------------------------------*)
inversemetric = FullSimplify[Inverse[metric]];
Print["(*-----------------------------------------------------------------------*)"];
Print["(*The Manifold has dimension   n= ", n,"\nCoordinate system:             ", coord," *)"];
Print["(*-----------------------------------------------------------------------*)"];
Print["(*\!\(\*SubscriptBox[\(g\), \(\[Mu]\[Nu]\)]\)= *)\n"  ,       metric // MatrixForm];
Print["(*\!\(\*SuperscriptBox[\(g\), \(\[Mu]\[Nu]\)]\)= *)\n",inversemetric // MatrixForm];
Print["(*g  = *)\n", Det[metric] // FullSimplify];
affine := affine = FullSimplify[Table[
    (1/2)*Sum[
        (inversemetric[[i, s]])*
        (     D[metric[[s, j]], coord[[k]] ] + 
              D[metric[[s, k]], coord[[j]] ] - 
              D[metric[[j, k]], coord[[s]] ]),
        {s, 1, n}],
    {i, 1, n}, {j, 1, n}, {k, 1, n}] ];
(*The non zero Christoffel symbols are computed and selected below: *)
listaffine := Table[
    If[
        UnsameQ[affine[[i, j, k]], 0],
        {Subscript[Superscript[\[CapitalGamma], i], j, k], affine[[i, j, k]]}
      ] ,
    {i, 1, n}, {j, 1, n}, {k, 1, j}];
Print["(*-----------------------------------------------------------------------*)"];
Print["(* Christoffel Symbols: *)"];
Print[TableForm[Partition[DeleteCases[Flatten[listaffine], Null], 2], TableSpacing -> {2, 2}]];
riemann := riemann = FullSimplify[Table[
    D[  affine[[i, l, j]], coord[[k]]          ] - 
    D[  affine[[i, k, j]], coord[[l]]          ] +      
    Sum[affine[[i, k, s]] affine[[s, l, j]] -    
        affine[[i, l, s]] affine[[s, k, j]] ,
        {s, 1, n}
       ],
    {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n} ] ];
listriemann := Table[
    If[
        UnsameQ[riemann[[i, j, k, l]], 0],
        {Subscript[Superscript[R, i], j, k, l], riemann[[i, j, k, l]]}
      ] ,
    {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, k - 1}];
Print["(*-----------------------------------------------------------------------*)"];
Print["(* Riemann Tensor: *)"];
Print[TableForm[Partition[DeleteCases[Flatten[listriemann], Null], 2], TableSpacing -> {2, 2}]];
lriemann := lriemann = FullSimplify[Table[
    Sum[metric[[i, ii]] riemann[[ii, j, k, l]], {ii, 1, n}],
    {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}]];
listlriemann := Table[
    If[
        UnsameQ[lriemann[[i, j, k, l]], 0],
        {Subscript    [R, i, j, k, l] , lriemann[[i, j, k, l]]}
      ] , {i, 1, n}, {j, 1, i - 1}, {k, 1, n}, {l, 1, k - 1}];
Print["(*-----------------------------------------------------------------------*)"];
Print["(* Contravariant Riemann Tensor: *)"];
Print[TableForm[
   Partition[DeleteCases[Flatten[listlriemann], Null], 2], TableSpacing -> {2, 2}]];
uriemann := uriemann = FullSimplify[Table[
    Sum[
        inversemetric[[j, jj]] inversemetric[[k, kk]] inversemetric[[l,ll]] riemann[[i, jj, kk, ll]], 
        {jj, 1, n}, {kk, 1, n}, {ll, 1,n}
       ],
    {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}]];
listuriemann := Table[
    If[
        UnsameQ[uriemann[[i, j, k, l]],0], 
        {Superscript[Superscript[Superscript[Superscript[R, i], j], k], l], 
         uriemann[[i, j, k, l]]}
      ] , {i, 1, n}, {j, 1, i - 1}, {k, 1, n}, {l, 1, k - 1}];
Print["(*-----------------------------------------------------------------------*)"];
Print["(* Covariant Riemann Tensor: *)"];
Print[TableForm[
   Partition[DeleteCases[Flatten[listuriemann], Null], 2], TableSpacing -> {2, 2}]];
r2 = FullSimplify[
   Sum[lriemann[[i, j, k, l]] uriemann[[i, j, k, l]], {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}]];
Print["(*-----------------------------------------------------------------------*)"];
Print["(*\!\(\*SuperscriptBox[\(R\), \(2\)]\): *)\n", r2];
ricci := ricci = FullSimplify[Table[
    Sum[
        riemann[[i, j, i, l]],
        {i, 1, n}
       ], {j, 1, n}, {l, 1, n}] ];
listricci := Table[
    If[
        UnsameQ[             ricci[[j, l]], 0],
        {Subscript[R, j, l], ricci[[j, l]]}
      ] , {j, 1, n}, {l, 1, j}];
Print["(*-----------------------------------------------------------------------*)"];
Print["(* Ricci Tensor: *)"];
Print[TableForm[Partition[DeleteCases[Flatten[listricci], Null], 2],TableSpacing -> {2, 2}]];
scalar = FullSimplify[
    Sum[inversemetric[[i, j]] ricci[[i, j]], {i, 1, n}, {j, 1, n}] 
                     ];
Print["(*-----------------------------------------------------------------------*)"];
Print["(* Curvature Scalar: *)"];
Print[scalar];
einstein := einstein = FullSimplify[ricci - (1/2) scalar*metric];
listeinstein := Table[
    If[
        UnsameQ[             einstein[[j, l]], 0],
        {Subscript[G, j, l], einstein[[j, l]]}
      ] , 
    {j, 1, n}, {l, 1, j}];
Print["(*-----------------------------------------------------------------------*)"];
Print["(* Einstein Tensor: *)"];
Print[TableForm[
   Partition[DeleteCases[Flatten[listeinstein], Null], 2], TableSpacing -> {2, 2}]];
weyl := weyl = FullSimplify[Table[
     If[n > 3,
      lriemann[[i, j, k, l]]
       - 1/(n - 2)           (metric[[i, k]] ricci[[l, j]] - 
                              metric[[i, l]] ricci[[k, j]] - 
                              metric[[j, k]] ricci[[l, i]] + 
                              metric[[j, l]] ricci[[k, i]])
       + 1/((n - 1) (n - 2)) (metric[[i, k]] metric[[l, j]] - 
                              metric[[i, l]] metric[[k, j]]) scalar
      (*else, if n<= 3 return 0:*), 0],
     {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}]];
listweyl := Table[
   If[
    UnsameQ[weyl[[i, j, k, l]], 0],
    {Subscript[\[GothicCapitalC], i, j, k, l], weyl[[i, j, k, l]]}
    ] , {i, 1, n}, {j, 1, i - 1}, {k, 1, n}, {l, 1, k - 1}];
Print["(*-----------------------------------------------------------------------*)"];
Print["(* Weyl Tensor: *)"];
Print[TableForm[Partition[DeleteCases[Flatten[listweyl], Null], 2], TableSpacing -> {2, 2}]];
geodesic := 
  geodesic = 
   Simplify[
    Table[-Sum[affine[[i, j, k]] u[j] u[k], {j, 1, n}, {k, 1, n}], {i, 1, n}]];
subst = Table[u[i] -> Subscript[coord[[i]], \[Tau]], {i, 1, n}];
nlistgeodesic := 
  Table[{Subscript[coord[[i]], \[Tau]\[Tau]]," \"+\" ", -geodesic[[i]] /. subst,"  \"=  0\"   "}, {i, 1, n}];
Print["(*-----------------------------------------------------------------------*)"];
Print["(* Geodesic Equations: *)"];
Print[TableForm[nlistgeodesic, TableSpacing -> {2}]];
