(* ::Package:: *)

BeginPackage["grt`"];

grt::usage=
"
grt: a package to compute the Curvature tensor, Ricci tensor, 
Einstein tensor, Christoffel symbols, etc. for arbitrary metrics.

(c) Pascal M. Vaudrevange (pascal@vaudrevange.com)
released under GPL
Version 22.Feb.2010


First, load a metric with 

	LoadMetric[coords, metric]

where coords is a list of coordinates of length ndim and metric is 
the ndim \[Times] ndim metric tensor. Note that the first coordinate has index 1. 
For example, Minkowski space would be

  coords={t,x,y,z}
  metric={{1,0,0,0},{0,-1,0,0},{0,0,-1},{0,0,0,-1}}

All objects are indexed the following way:

	A[\!\(\*SuperscriptBox[\"\[Mu]\", \"-\"]\),\!\(\*SubscriptBox[\"\[Nu]\", \"-\"]\)] = \!\(\*SubscriptBox[SuperscriptBox[\"A\", \"\[Mu]\"], \"\[Nu]\"]\)
	A[\!\(\*SubscriptBox[\"\[Mu]\", \"-\"]\),\!\(\*SubscriptBox[\"\[Nu]\", \"-\"]\)] = \!\(\*SubscriptBox[\"A\", \"\[Mu]\[Nu]\"]\)
	A[\!\(\*SuperscriptBox[\"\[Mu]\", \"-\"]\),\!\(\*SuperscriptBox[\"\[Nu]\", \"-\"]\)] = \!\(\*SuperscriptBox[\"A\", \"\[Mu]\[Nu]\"]\)
	etc.

where \[Mu],\[Nu] are any numbers between 1, ndim.
The indices of all tensor objects MUST indicate whether it is a lower 
or upper index. For example, A[1,1] is illegal. To test for the position 
of an index, use the supplied functions Uppper[] and Lower[].

For the curious: The upper/lower index is defined by 

	SuperMinus[\[Mu]_Integer]={\[Mu],'u'};
	SubMinus[\[Mu]_Integer]={\[Mu],'l'};

For a list of defined symbols, enter 

	?grt`*"


LoadMetric::usage=
"Loads the metric and coordinates into the grt package."

g::usage=
"Metric tensor g[\!\(\*SubscriptBox[\"a\", \"-\"]\),\!\(\*SubscriptBox[\"a\", \"-\"]\)]=\!\(\*SubscriptBox[\"g\", \"ab\"]\)"

\[CapitalGamma]::usage=
"Christoffel Symbol 

	\[CapitalGamma][\!\(\*SuperscriptBox[\"a\", \"-\"]\),\!\(\*SubscriptBox[\"m\", \"-\"]\),\!\(\*SubscriptBox[\"n\", \"-\"]\)] = \!\(\*SubscriptBox[SuperscriptBox[\"\[CapitalGamma]\", \"a\"], \"mn\"]\) = \!\(\*FractionBox[\"1\", \"2\"]\)\!\(\*SuperscriptBox[\"g\", \"ab\"]\) (\!\(\*SubscriptBox[\"g\", 
RowBox[{\"bm\", \",\", \"n\"}]]\)-\!\(\*SubscriptBox[\"g\", 
RowBox[{\"gn\", \",\", \"m\"}]]\)+\!\(\*SubscriptBox[\"g\", 
RowBox[{\"mn\", \",\", \"b\"}]]\))

where \!\(\*SubscriptBox[\"f\", 
RowBox[{\",\", \"n\"}]]\) is the partial derivative \!\(\*SubscriptBox[\"\[PartialD]\", \"n\"]\)f.
"

R::usage=
"Curvature (Riemann) tensor 

	R[\!\(\*SuperscriptBox[\"b\", \"-\"]\),\!\(\*SubscriptBox[\"m\", \"-\"]\),\!\(\*SubscriptBox[\"s\", \"-\"]\),\!\(\*SubscriptBox[\"q\", \"-\"]\)]= \!\(\*SubscriptBox[SuperscriptBox[\"R\", \"b\"], \"msq\"]\) = \!\(\*SubscriptBox[SuperscriptBox[\"\[CapitalGamma]\", \"b\"], 
RowBox[{\"mq\", \",\", \"s\"}]]\)-\!\(\*SubscriptBox[SuperscriptBox[\"\[CapitalGamma]\", \"b\"], 
RowBox[{\"ms\", \",\", \"q\"}]]\)+\!\(\*SubscriptBox[SuperscriptBox[\"\[CapitalGamma]\", \"b\"], \"ns\"]\)\!\(\*SubscriptBox[SuperscriptBox[\"\[CapitalGamma]\", \"n\"], \"mq\"]\)-\!\(\*SubscriptBox[SuperscriptBox[\"\[CapitalGamma]\", \"b\"], \"nq\"]\)\!\(\*SubscriptBox[SuperscriptBox[\"\[CapitalGamma]\", \"n\"], \"ms\"]\),

or

Ricci Tensor

	R[\!\(\*SubscriptBox[\"m\", \"-\"]\),\!\(\*SubscriptBox[\"q\", \"-\"]\)]=\!\(\*SubscriptBox[\"R\", \"mq\"]\)= \!\(\*SubscriptBox[SuperscriptBox[\"R\", \"s\"], \"msq\"]\)

where \!\(\*SubscriptBox[\"f\", 
RowBox[{\",\", \"n\"}]]\) is the partial derivative \!\(\*SubscriptBox[\"\[PartialD]\", \"n\"]\)f.
"

RicciScalar::usage=
"RicciScalar=\!\(\*SubscriptBox[SuperscriptBox[\"R\", \"m\"], \"m\"]\)"

G::usage=
"Einstein tensor G[\!\(\*SubscriptBox[\"a\", \"-\"]\),\!\(\*SubscriptBox[\"b\", \"-\"]\)]=\!\(\*SubscriptBox[\"G\", \"ab\"]\)=\!\(\*SubscriptBox[\"R\", \"ab\"]\)-\!\(\*FractionBox[\"1\", \"2\"]\) R \!\(\*SubscriptBox[\"g\", \"ab\"]\)"

Detg::usage=
"Determinant Detg=|g| (can be negative)"

ValidP::usage=
"ValidP[\[Mu]] tests if \[Mu] is a valid upper or lower index,automatically threads over lists."

Upper::usage=
"Upper[\[Mu]] tests if \[Mu] is an upper index, automatically threads over lists.

Examples:
	Upper[\!\(\*SuperscriptBox[\"1\", \"-\"]\)]=True
	Upper[\!\(\*SubscriptBox[\"1\", \"-\"]\)]=False"

Lower::usage=
"Lower[\[Mu]] tests if \[Mu] is a lower index, automatically threads over lists.

Examples:
	Lower[\!\(\*SubscriptBox[\"1\", \"-\"]\)]=True
	Lower[\!\(\*SuperscriptBox[\"1\", \"-\"]\)]=False"

Quabla::usage=
"Quabla[f]=\[Square]f= \!\(\*SuperscriptBox[\"g\", \"ab\"]\) \!\(\*SubscriptBox[\"f\", 
RowBox[{\",\", 
RowBox[{\"a\", \";\", \"b\"}]}]]\) is the covariant box operator acting on scalar 
functions, where \!\(\*SubscriptBox[\"f\", 
RowBox[{\",\", 
RowBox[{\"a\", \";\", \"b\"}]}]]\)=\!\(\*SubscriptBox[\"D\", 
RowBox[{\"b\", \" \"}]]\)\!\(\*SubscriptBox[\"\[PartialD]\", \"a\"]\)f"

Show\[CapitalGamma]::usage=
"Show\[CapitalGamma] prints the non-zero components of \!\(\*SubscriptBox[SuperscriptBox[\"\[CapitalGamma]\", \"a\"], \"bc\"]\)"

ShowCurvatureTensor::usage=
"ShowCurvatureTensor prints the non-zero compontents of the Curvature (Riemann) tensor \!\(\*SubscriptBox[SuperscriptBox[\"R\", \"a\"], \"bcd\"]\)."

ShowRicciTensor::usage=
"ShowRicciTensor prints the non-zero compontents of the Ricci tensor \!\(\*SubscriptBox[\"R\", \"ab\"]\)."

ShowEinsteinTensor::usage=
"ShowEinstenTensor prints the non-zero components of the Einstein Tensor \!\(\*SubscriptBox[\"G\", \"ab\"]\)."


Begin["`Private`"];
Print["grt: a package to compute Christoffel symbols, Curvature tensor, Ricci Tensor, 
Einstein tensor etc. for arbitrary metrics.

(c) Pascal M. Vaudrevange (pascal@vaudrevange.com)
released under GPL
Version 22. Feb. 2010"];


LoadMetric[InCoord_,InMetric_]:=Module[{},
				coord=InCoord;
				metric=InMetric;
				Print["\!\(\*SubscriptBox[\"g\", \"\[Mu]\[Nu]\"]\)=",metric//MatrixForm];
				Print["coordinates = ", coord];
				ndim=Length[coord];
				If[UnsameQ[Dimensions[metric],{4,4}], 
					Print["Metric does not have the correct dimensions"];
				];
				inversemetric:=inversemetric=Simplify[Inverse[metric]];
Connection:=Connection=Module[{a,b,m,n},
			      Simplify[
				       Table[
					     1/2 Sum[
						     g[SuperMinus[a],SuperMinus[b] ](D[g[SubMinus[b],SubMinus[m]],coord[[n]]] +D[g[SubMinus[b],SubMinus[n]],coord[[m]]]- D[g[SubMinus[m],SubMinus[n]],coord[[b]]]),
						     {b,1,ndim}
						     ],
					     {a,1,ndim},
					     {m,1,ndim},
					     {n,1,ndim}
					     ]
				       ]
			      ];
CurvatureTensor:=CurvatureTensor=Module[{b,m,s,q,n},
		Simplify[
			Table[
				D[Connection[[b,m,q]],coord[[s]]]-
				D[Connection[[b,m,s]],coord[[q]]]+
				Sum[
					Connection[[b,n,s]] Connection[[n,m,q]] - Connection[[b,n,q]] Connection[[n,m,s]],
					{n,1,ndim}
				],
			{b,1,ndim},
			{s,1,ndim},
			{q,1,ndim},
			{m,1,ndim}
			]
		]
	];

RicciTensor:=RicciTensor=Module[{b, m, q, n}, 
            Simplify[
	    Table[
            Sum[
	    D[Connection[[b, m, q]], coord[[b]]] - D[Connection[[b, m, b]], coord[[q]]] + Sum[
	    Connection[[b, n, b]] Connection[[n, m, q]] - Connection[[b, n, q]] Connection[[n, m, b]], {n, 1,ndim}],
	    {b, 1, ndim}],
	    {q, 1, ndim}, {m, 1, ndim}]]];

EinsteinTensor:=EinsteinTensor=Simplify[RicciTensor-1/2 RicciScalar  * metric];
]


SuperMinus[\[Mu]_Integer]={\[Mu],"u"};


SubMinus[\[Mu]_Integer]={\[Mu],"l"};


ValidP[\[Mu]_]:=Module[{res,i},
		       res=True;
		       Switch[ Depth[\[Mu]],
		       2,res=(Length[\[Mu]]==2),
                       3,Do[
			    res=(res && ValidP[\[Mu][[i]]]),
			 {i,1,Length[\[Mu]]}
			 ],
		       _,res=False;Print["Error in ValidP"]
		       ];
		       res
		]


Upper[\[Mu]_]:=Module[{res,i},
		      res=True;
		      Switch[Depth[\[Mu]],
		      2,res=(\[Mu][[2]]=="u"),
                      3,Do[
			   res=(res && Upper[\[Mu][[i]]]),
			{i,1,Length[\[Mu]]}
			],
		      _, res=False; Print["Error in Upper"]
		      ];
		      res
	      ]


Lower[\[Mu]_]:=Module[{res,i},
			res=True;
			Switch[Depth[\[Mu]],
				2,res=(\[Mu][[2]]=="l"),
				3,Do[
					res=(res && Lower[\[Mu][[i]]]),
					{i,1,Length[\[Mu]]}
					],
				_, res=False; Print["Error in Lower"]
			];
			res
		]


g[\[Mu]_,\[Nu]_]:=Module[{res},
			If[(Length[\[Mu]]==  2 ) && ( Length[\[Nu]]==2),
                Which[
					Lower[{\[Mu],\[Nu]}], res=metric[[\[Mu][[1]],\[Nu][[1]]]],
					Upper[{\[Mu],\[Nu]}], res=inversemetric[[\[Mu][[1]],\[Nu][[1]]]],
					True,If[SameQ[\[Mu][[1]],\[Nu][[1]]],res=1,res=0]
			],
			Print["Error in \!\(\*SubscriptBox[\"g\", \"\[Mu]\[Nu]\"]\): \[Mu]=", \[Mu], ", \[Nu]=",\[Nu]];
			];
		res
	];


\[CapitalGamma][\[Mu]_,\[Nu]_,\[Lambda]_]:=Module[{\[Sigma],\[Rho],\[Alpha]},
				If[ValidP[{\[Mu],\[Nu],\[Lambda]}],
					Sum[
						g[\[Mu],SubMinus[\[Alpha]]] g[\[Nu],SuperMinus[\[Rho]]] g[\[Lambda],SuperMinus[\[Sigma]]] Connection[[\[Alpha], \[Rho], \[Sigma]]],
						{\[Sigma],1,ndim},
						{\[Rho],1,ndim},
						{\[Alpha],1,ndim}
					],
					Print["Error in \[CapitalGamma]"];
				]
			];


R[\[Mu]_,\[Nu]_,\[Rho]_,\[Sigma]_]:=Module[{\[Alpha],\[Beta],\[Gamma],\[Delta]},
			If[ValidP[{\[Mu],\[Nu],\[Rho],\[Sigma]}],
				Sum[
					g[\[Mu],SubMinus[\[Alpha]]] g[\[Nu],SubMinus[\[Beta]]] g[\[Rho],SuperMinus[\[Gamma]]] g[\[Sigma],SuperMinus[\[Delta]]]  CurvatureTensor[[\[Alpha],\[Beta],\[Gamma],\[Delta]]],
					{\[Alpha],1,ndim},
					{\[Beta],1,ndim},
					{\[Gamma],1,ndim},
					{\[Delta],1,ndim}
				],
				Print["Error in \!\(\*SubscriptBox[\"R\", \"\[Mu]\[Nu]\[Rho]\[Sigma]\"]\)"];
			]
		];



R[\[Mu]_,\[Nu]_]:=Module[{\[Alpha],\[Beta]},
			If[ValidP[{\[Mu],\[Nu]}],
				Sum[
					g[\[Mu],SuperMinus[\[Alpha]]] g[\[Nu],SuperMinus[\[Beta]] ] RicciTensor[[\[Alpha],\[Beta]]],
					{\[Alpha],1,ndim},
					{\[Beta],1,ndim}
				],
				Print["Error in Ricci tensor \!\(\*SubscriptBox[\"R\", \"\[Mu]\[Nu]\"]\)"];
			]
		];


RicciScalar:=Simplify[Tr[inversemetric.RicciTensor]];


G[\[Mu]_,\[Nu]_]:=Module[{\[Alpha],\[Beta]},
			If[ValidP[{\[Mu],\[Nu]}],
				Sum[
					g[\[Mu],SuperMinus[\[Alpha]]] g[\[Nu],SuperMinus[\[Beta]] ] EinsteinTensor[[\[Alpha],\[Beta]]],
					{\[Alpha],1,ndim},
					{\[Beta],1,ndim}
				],
				Print["Error in Einstein tensor \!\(\*SubscriptBox[\"G\", \"\[Mu]\[Nu]\"]\)"];
			]
		];


Detg:=Module[{},
		Det[metric]
];


NonZero[object_]:=Module[{},Print["NonZero: General Printing routine: Implement me!"]];


Quabla[f_]:=Module[{a,b,n},
			Simplify[Sum[
				g[SuperMinus[a],SuperMinus[b]] (D[f,coord[[a]],coord[[b]]]-
							Sum[
								Connection[[n,a,b]] D[f,coord[[n]]],
								{n,1,ndim}
							]
						),
				{b,1,ndim},
				{a,1,ndim}
			]
			]
		];


Show\[CapitalGamma]:=Module[{a,b,c},
			Table[
				If[
					UnsameQ[Connection[[a,b,c]],0],
					Print[Subscript[Superscript["\[CapitalGamma]",a-1],b-1,c-1],"=",Connection[[a,b,c]]];
				],
				{a,1,ndim},
				{b,1,ndim},
				{c,1,ndim}
			];
		];


ShowCurvatureTensor:=Module[{a,b,c,d},
			Table[
				If[
					UnsameQ[CurvatureTensor[[a,b,c,d]],0],
					Print[Subscript[Superscript["R",a-1],b-1,c-1,d-1],"=",CurvatureTensor[[a,b,c,d]]];
				],
				{a,1,ndim},
				{b,1,ndim},
				{c,1,ndim},
				{d,1,ndim}
			];
		];


ShowRicciTensor:=Module[{a,b},
			Table[
				If[
					UnsameQ[RicciTensor[[a,b]],0],
					Print[Subscript["R",a-1,b-1],"=",RicciTensor[[a,b]]];
				],
				{a,1,ndim},
				{b,1,ndim}
			];
	];


ShowEinsteinTensor:=Module[{a,b},
			Table[
				If[
					UnsameQ[EinsteinTensor[[a,b]],0],
					Print[Subscript["G",a-1,b-1],"=",EinsteinTensor[[a,b]]];
				],
				{a,1,ndim},
				{b,1,ndim}
			];
	];


End[];
EndPackage[];
