(* ::Package:: *)

BeginPackage["GR`"];

Chr::usage =
 	"Christoffel Symbols of the first kind";

Einstein::usage =
	"Einstein tenosr G"

Imetric::usage =
	"Inverse of the metric tensor (fully covariant)"

Riemann::usage =
	"Riemann tensor (1-3)"

Ricci::usage =
	"Ricci tensor"

RicciScalar::usage =
	"Ricci scalar of curvature"

TensorC::usage =
	"Tensor Contraction"


Derivative/:MakeBoxes[Derivative[\[Alpha]__][f1_][vars__Symbol],TraditionalForm]:=Module[{bb,dd,sp},

	MakeBoxes[dd,_]^=If[Length[{\[Alpha]}]==1,"\[DifferentialD]","\[PartialD]"];
	MakeBoxes[sp,_]^="\[ThinSpace]";
	
	bb/:MakeBoxes[bb[x__],_]:=RowBox[Map[ToBoxes[#]&,{x}]];

	TemplateBox[{ToBoxes[bb[dd^Plus[\[Alpha]],f1]],ToBoxes[Apply[bb,Riffle[Map[bb[dd,#]&,Select[({vars}^{\[Alpha]}),(#=!=1&)]],sp]]],
		ToBoxes[Derivative[\[Alpha]][f1][vars]]},"ShortFraction",DisplayFunction:>(FractionBox[#1,#2]&),
			InterpretationFunction:>(#3&),Tooltip->Automatic]];


Begin["`Private`"];

Print["Tensor-Manipulation package for General Relativity"];


Imetric[metric_]:=Inverse[metric];


Chr[metric_,y_]:=Module[{l1,l2,l3,l4,Ig,dim = Length[y],Christoffel},

	Ig = Imetric[metric];	

	Christoffel = 1/2 Table[Sum[Ig[[l1,l4]](D[metric[[l2,l4]],y[[l3]]]
		+D[metric[[l4,l3]],y[[l2]]]-D[metric[[l2,l3]],y[[l4]]]),{l4,1,dim}],{l1,1,dim},{l2,1,dim},{l3,1,dim}];

	Christoffel
]


Riemann[metric_,y_]:=Module[{l1,l2,l3,l4,l5,dim = Length[y],Christoffel,Riem},
		
	Christoffel = Chr[metric,y];
	
Riem = Table[D[Christoffel[[l1]][[l2,l4]],y[[l3]]],{l1,1,dim},{l2,1,dim},{l3,1,dim},{l4,1,dim}]
			-Table[D[Christoffel[[l1]][[l2,l3]],y[[l4]]],{l1,1,dim},{l2,1,dim},{l3,1,dim},{l4,1,dim}]
			-Table[Sum[Christoffel[[l1]][[l5,l4]]Christoffel[[l5]][[l2,l3]],{l5,1,dim}],{l1,1,dim},{l2,1,dim},{l3,1,dim},{l4,1,dim}]
			+Table[Sum[Christoffel[[l1]][[l5,l3]]Christoffel[[l5]][[l2,l4]],{l5,1,dim}],{l1,1,dim},{l2,1,dim},{l3,1,dim},{l4,1,dim}];	
			
			Riem
]


Ricci[metric_,y_]:=Module[{l1,l2,l3,dim = Length[y],RicciTensor,Riem},
		
	Riem = Riemann[metric,y];
	
	RicciTensor = Table[Sum[Riem[[l3,l1,l3,l2]],{l3,dim}],{l1,dim},{l2,dim}];

	RicciTensor
]


RicciScalar[metric_,y_]:=Module[{l1,l2,Ig,dim = Length[y],RicciScal,RicciTensor},

	Ig = Imetric[metric];

	RicciTensor = Ricci[metric,y];

	RicciScal = Sum[Ig[[l1,l2]]RicciTensor[[l1,l2]],{l1,1,dim},{l2,1,dim}];

	RicciScal
]


Einstein[metric_,y_]:=Module[{l1,l2,dim = Length[y],RicciScal,RicciTensor,G,\[Epsilon]},

	RicciTensor = Ricci[metric,y];
		RicciScal = RicciScalar[metric,y];

	G = Table[RicciTensor[[l1,l2]]-1/2metric[[l1,l2]]RicciScal,{l1,dim},{l2,dim}];

	G
]


TensorC[T1_,T2_]:=Module[{Tensor},

	Tensor = Inner[Times,T1,T2];
	
	Return[FullSimplify[Tensor]];
]


End[];

EndPackage[];
