(** Mathematica code for computing          **
 ** connections and covariant derivatives   **
 ** Courtesy of Prof. Francesco Bullo, UIUC **)

BeginPackage["MechSys`"]

LeviCivita::usage = "LeviCivita[M,x] computes 
the Christoffel symbols of the Levi-Civita 
connection for the mechanical system with inertia
matrix M with respect to the coordinates x.";

CovariantDer::usage = "CovariantDer[X,Y,Nabla,x]
computes the covariant derivative of Y along X
with respect to a connection Nabla in coordinates x 
(that is, Nabla_XY).";

Begin["Private`"]

LeviCivita[M_, x_] := 
 Module[{Minv=Inverse[M],i,j,k,h,N=Length[x]},
    Table[ Sum[ Minv[[h,k]] ( D[M[[h,j]], x[[i]]] +
                              D[M[[i,h]], x[[j]]] - 
                              D[M[[i,j]], x[[h]]] ) / 2,
                {h,N}], 
           {k,N}, {j,N}, {i,N} ]];

CovariantDer[X_, Y_, Nabla_, x_] := 
 Module[{i,j,k,N=Length[x]},
     Table[ Sum[ D[Y[[i]],x[[j]]] X[[j]] +
                 Sum[ Nabla[[i,j,k]] X[[j]] Y[[k]], 
                      {k,N}],
                 {j,N}], 
            {i,N} ] ];

End[] EndPackage[]

(** To use the package in Mathematica, first enter  **
 ** Needs["MechSys'"]                               **
 ** then you can use the functions                  **)