(**************************************************************************
*                                                                         *
*                               Ricci                                     *
*        A Mathematica package for doing tensor calculations              *
*                      in differential geometry                           *
*                                                                         *
*                            Version 1.51                                 *
*                                                                         *
*                                                                         *
*                           By John M. Lee                                *
*    assisted by Dale Lear, John Roth, Lee Nave, and Larry Peterson       *
*                                                                         *
*                Copyright (c) 1992 - 2004 John M. Lee                    *
*                        All rights reserved                              *
*                                                                         *
*         Development of this software was supported in part              *
*                by NSF grants DMS-9101832, DMS-9404107                   *
*                                                                         *
*                                                                         *
* This software package and its accompanying documentation are provided   *
* as is, without guarantee of support or maintenance.  The copyright      *
* holder makes no express or implied warranty of any kind with respect    *
* to this software, including implied warranties of merchantability or    *
* fitness for a particular purpose, and is not liable for any damages     *
* resulting in any way from its use.                                      *
*                                                                         *
* Everyone is granted permission to copy, modify and redistribute this    *
* software package and its accompanying documentation, provided that:     *
*  1. All copies contain this notice in the main program file and in the  *
*     supporting documentation.                                           *
*  2. All modified copies carry a prominent notice stating who made the   *
*     last modification and the date of such modification.                *
*  3. No charge is made for this software or works derived from it, with  *
*     the exception of a distribution fee to cover the cost of materials  *
*     and/or transmission.                                                *
*                                                                         *
**************************************************************************)

BeginPackage["Ricci`"];

Print[" -- Ricci Version 1.51 (June 28, 2004) --"];
Print["    Copyright 1992 - 2004 John M. Lee"];
Print["    Problem reports or suggestions to:"];
Print["      lee@math.washington.edu"];

If [ $VersionNumber >= 3.0,
     Print[" -- For formatted output in a Version 3, 4, or 5 notebook:"];
     Print["    Cell menu"];
     Print["      -> Default Output Format Type"];
     Print["         -> OutputForm"]
    ];


(* This function is used to append our messages to the system messages
   for system functions that are modified by Ricci.  If the system
   message doesn't have head String, then ignore it.  Otherwise, see
   if our message is already there.  If so, leave it alone; if not,
   append our message to the existing one. *)

Begin["`Private`"];

AppendMessage[ sysmsg_, string_ ] :=
  Which[
    Head[sysmsg] =!= String,
      string,
    StringMatchQ[ sysmsg, Evaluate[ "*" <> string <> "*" ] ],
      sysmsg,
    True,
      sysmsg <> "\n\n" <> string];

End[ (* "Ricci`Private`" *) ];

(*****************************************************************)

(** Initialize usage messages.  These messages define all the
    user-accessible symbols, defined in the "Ricci`" context. **)

AbsorbMetrics::usage = "AbsorbMetrics[expr] simplifies expr by
eliminating any metrics that are contracted with other tensors, and
using them to raise or lower indices. AbsorbMetrics[x,n] or
AbsorbMetrics[x,{n1,n2,...}] applies AbsorbMetrics only to term n or
to terms n1,n2,... of x.
\nOption:
\n* Mode -> All or NoOneDims. NoOneDims means that AbsorbMetrics
should not absorb one-dimensional metrics (unless they are paired with
other metrics). Default is All.";

Alt::usage = "Alt[x] is the alternating part of the tensor expression
x.";

Alternating::usage = "Alternating is a value for the Symmetries option
of DefineTensor.";

Any::usage = "Any can be used as a value for the Bundle option of the
DefineTensor command. Bundle -> Any means that indices from any
bundle can be inserted.";

Bar::usage =
"The internal form for a barred index is L[i[Bar]] or U[i[Bar]]. In
input form, these can be abbreviated LB[i] and UB[i].
\nThe
internal form for the conjugate of a tensor is Tensor[name[Bar],{...},{...}].
In input form, this is typed Conjugate[name] [...] [...].";

Basis::usage = "Basis is the name used for generic basis vectors and
covectors for any bundle. Basis vectors are generated automatically
by BasisExpand. For example, if index i is associated with bundle b,
then Basis[L[i]] and Basis[U[i]] represent contravariant and covariant
basis elements for b (i.e., basis elements for b and its dual),
respectively; Basis[LB[i]] and Basis[UB[i]] represent basis elements
for the conjugate bundle Conjugate[b].
\nIf you insert indices into
Basis vectors, you get Kronecker delta functions, metrics, or zero as
appropriate.";

BasisExpand::usage = "BasisExpand[x] converts x to a sum of component
expressions multiplied by basis vectors and covectors, Basis[L[i]] and
Basis[U[j]]. BasisExpand[x,n] or BasisExpand[x,{n1,n2,...}] applies
BasisExpand only to term n or terms n1,n2,... of x.";

BasisGather::usage = "BasisGather[x,tensor] attempts to recognize the
basis expression for \"tensor\" in x, and replaces it by \"tensor\".
BasisGather[x,{tensor1,tensor2,...}] does the same thing for several
tensors at once.";

BasisName::usage = "BasisName -> name is an option for DefineBundle.
It specifies the name to be used in the input and output forms of the
basis vectors for the bundle being defined. Default = Basis.";

BianchiRules::usage = "BianchiRules[i,j,k] converts tensors with
Symmetries -> RiemannSymmetries containing the indices i,j,k in order
to two-term sums, using the first and second Bianchi identities.";

Bundle::usage = "Bundle -> name is a DefineTensor option, specifying
the name of the bundle the tensor is to be associated with. This can
be a list of bundle names, meaning that the tensor is associated with
the direct sum of the bundles in the list.
\nThe function Bundle[i]
returns the name of the bundle associated with the index i.";

BundleDummy::usage =
"BundleDummy[bundle] returns the symbol that is used for
computer-generated dummy indices associated with the bundle. This is the last
index name in the list BundleIndices[bundle]. For most bundles,
new dummy indices will be of the form \"kn\", where k=BundleDummy[bundle] and
n is an integer. For one-dimensional bundles, n is omitted.";

BundleIndices::usage =
"BundleIndices[bundle] returns a list of the index names currently associated
with bundle.";

BundleQ::usage = "BundleQ[x] returns True if x is the name of a bundle,
and False otherwise.";

Bundles::usage = "Bundles[x] returns a list of the bundles
associated with index positions in the tensor expression x; each
entry in the list is a sublist giving the allowed bundles for the
corresponding index position.";

Co::usage = "Co is an abbreviation for Covariant.";

CoBasisName::usage = "CoBasisName -> name is an option for
DefineBundle. It specifies the name to be used in input and output
form for the basis covectors for the bundle being defined. Default is
the name given in the BasisName option, or Basis if BasisName is not
specified.";

CollectConstants::usage = "CollectConstants[x] groups together terms
in the tensor expression x having the same tensor factors but
different constant factors. CollectConstants[x,n] or
CollectConstants[x,{n1,n2,...}] applies CollectConstants only to term
n or to terms n1,n2,... of x.";

CommuteCovD::usage = "CommuteCovD[ x, L[i], L[j] ] changes all
adjacent occurrences of indices L[i], L[j] after the \";\" to L[j],
L[i] by adding appropriate curvature and torsion terms.";

CommutingFrame::usage = "CommutingFrame is an option for DefineBundle,
which is meaningful only for tangent bundles and their subbundles. If
True, it means that the default frame for that bundle is always
assumed to consist of commuting vector fields.";

CompatibilityRule::usage = "CompatibilityRule is a rule that
transforms derivatives of metric components into connection
coefficients, using the fact that the default connection is compatible
with the metric.";

Complex::usage = Ricci`Private`AppendMessage[Complex::usage, "Ricci
recognizes Complex as a value for the Type option of DefineBundle,
DefineConstant, DefineMathFunction and DefineTensor."];

Con::usage = "Con is an abbreviation for Contravariant.";

Conjugate::usage = Ricci`Private`AppendMessage[Conjugate::usage, "The
Ricci package modifies Conjugate to handle tensor expressions and
indices. The behavior of a tensor, constant, index, mathematical
function, or bundle under conjugation is determined by the Type option
when the object is defined."];

Conn::usage = "Conn is the generic connection form for the default
connection on any bundle. It is defined as a product tensor of rank
{2,1}; when the first two indices are inserted, it represents the
matrix of connection 1-forms associated with the default basis. When
all three indices are inserted, it represents the Christoffel symbols
of the default connection relative to the default basis.";

Connection::usage = "Connection -> cn is an option for some of the
Ricci differentiation functions. It specifies that covariant
derivatives are to be taken with respect to the connection cn instead
of the default connection. Ordinarily, cn will be an expression of
the form \"Conn + diff\", where diff is a 3-tensor expression
representing the difference tensor between the default connection and
cn.";

ConnToMetricRule::usage = "ConnToMetricRule is a rule that causes
components of the default connection Conn to be converted to
expressions involving directional derivatives of the metric. This
rule applies only to fully-indexed connection components whose indices
all refer to a single Riemannian tangent bundle with the option
CommutingFrame -> True.";

ConstantFactor::usage = "ConstantFactor[x] returns the product
of all the constant factors in x, which should be a monomial.";

ConstantQ::usage = "ConstantQ[x] returns True if there are no explicit
tensors in x, and False otherwise.";

ContractedBianchiRules::usage = "ContractedBianchiRule is a rule that
simplifies contracted covariant derivatives of the Riemannian and
Ricci curvature tensors using contracted versions of the second
Bianchi identity.";

Contravariant::usage = "Contravariant is a value for the Variance option
of DefineTensor, and may be abbreviated Con. If an index slot is
Contravariant, indices in that slot are upper by default.";

CorrectAllVariances::usage = "CorrectAllVariances[x] changes the
variance (upper to lower or lower to upper) of indices in x whose
variances are not correct for their positions, by inserting
appropriate metric coefficients.
\nOption:
\n* Mode -> All or OneDims. OneDims means that CorrectAllVariances
should correct variances only of one-dimensional indices and indices
that appear inside differential operators such as Del or Extd.
Default is All.";

CovD::usage = "If x is a component expression (no unfilled index
slots), then x[ L[i], L[j] ] or CovD[ x, {L[i],L[j]} ] is the
component of the covariant derivative of x in the L[i],L[j]
directions. In output form, covariant derivatives of a tensor are
represented by indices following a semicolon."<>
"\nOptions (for the second format only):
\n* Connection -> cn specifies that covariant derivatives are to be
taken with respect to the connection cn instead of the default
connection.
\n* Metric -> g: a symmetric 2-tensor expression, indicating that
the covariant derivatives are to be taken with respect to the
Levi-Civita connection of g instead of the default metric for x's bundle
(which is assumed to be Riemannian).";

CovDExpand::usage = "CovDExpand[x] converts all covariant derivatives
in x to ordinary directional derivatives (represented as
Del[Basis[L[i]],...]) and connection coefficients. CovDExpand[x,n] or
CovDExpand[x,{n1,n2,...}] applies CovDExpand only to term n or to terms
n1,n2,... of x.";

CovDSimplify::usage = "CovDSimplify[x] attempts to simplify x as much
as possible by ordering all dummy indices, including those that occur
after the \";\" in component expressions. CovDSimplify[x,n] or
CovDSimplify[x,{n1,n2,...}] applies CovDSimplify only to term n or to
terms n1,n2,... of x. CovDSimplify first calls TensorSimplify, then
OrderCovD, then TensorSimplify again.";

Covariant::usage = "Covariant is a tensor Variance. If an index slot
is Covariant, indices in that slot are lower by default. May be
abbreviated Co.";

Curv::usage = "Curv is the generic curvature tensor associated with
the default connection on any bundle. It is generated when covariant
derivatives are commuted, and when SecondStructureRule is applied to
derivatives of connection forms. Curv is a product tensor with rank
{2,2}, which is Alternating in its last two indices. Inserting the
first two indices yields the matrix of curvature 2-forms associated
with the default basis. Inserting all four indices yields the
coefficients of the curvature tensor.";

Curvature::usage = "Curvature[cn] is the curvature tensor associated
to the connection cn. Ordinarily, cn will be an expression of the
form \"Conn + diff\", where diff is a 3-tensor expression representing
the difference tensor between the default connection and cn.";

CurvToConnRule::usage = "CurvToConnRule is a rule that converts
components of curvature tensors to connection coefficients and their
directional derivatives.";

Declare::usage = "Declare[name,options] can be used to change certain
options for a previously-defined bundle, constant, index, tensor, or
mathematical function. Declare[{name1,name2,...},options] changes
options for several names at once. For constants and math functions,
only the Type option is allowed. For indices, the only allowable
option is TeXFormat. For bundles, the allowable options are
FlatConnection, ParallelFrame, OrthonormalFrame, PositiveDefinite,
CommutingFrame, and TorsionFree. For tensors, the allowable options
are Type, TeXFormat, and Bundle.";

DeclareBundle::usage = Declare::usage;

DeclareConstant::usage = Declare::usage;

DeclareIndex::usage = Declare::usage;

DeclareTensor::usage = Declare::usage;

DeclareMathFunction::usage = Declare::usage;

DefineBundle::usage =
"DefineBundle[ name, dim, metric, {indices} ] defines a bundle.
\n* name: A symbol that uniquely identifies the bundle.
\n* dim: The dimension of the bundle. A positive integer or
symbolic constant.
\n* metric: A name for the bundle's metric.
\n* indices: A list of index names to be associated with the
bundle."<>
"\nOptions:
\n* Type -> Complex or Real. Default is Real.
\n* TangentBundle -> bundle. The name of the underlying tangent
bundle for the bundle being defined. Default is $DefaultTangentBundle
if defined, otherwise this bundle itself (and its conjugate if the
bundle is complex)."<> "\n* FlatConnection -> True or False. Default
is False.
\n* ParallelFrame -> True or False. Default is False.
\n* OrthonormalFrame -> True or False. Default is False.
\n* CommutingFrame -> True or False. Default is False.
\n* PositiveDefinite -> True or False. Default is True.
\n* TorsionFree -> True or False. Default is True.
\n* BasisName -> name. A name for the default basis for this bundle.
\n* CoBasisName -> name. A name for the default covariant basis for
this bundle.
\n* MetricType -> Riemannian or Normal. Default is Normal. If Riemannian
is specified, there are additional options for curvature
conventions.
\n* Quiet -> True or False. Default is False, or the value of $Quiet if set.";

DefineConstant::usage = "DefineConstant[symbol] defines symbol to be a
constant with respect to differentiation in all space variables.
\nOptions: 
\n* Type -> type. This can be a single keyword or a list of
keywords, chosen from among Real, Complex, Imaginary, Positive,
Negative, NonPositive, NonNegative, Integer, Even, or Odd. Default is
Real.
\n* Quiet -> True or False. Default is False, or the value of
$Quiet if set.";

DefineIndex::usage = "DefineIndex[ {i,j,k}, bundle ]
causes the index names i,j,k to be associated with bundle.
\n* The first argument must be a symbol or list of symbols.
\n* The second argument must be a the name of a bundle.
\nOptions:
\n* TeXFormat -> \"texformat\". Default
is the index name itself.
\n* Quiet -> True or False. Default is False, or the value of $Quiet if set.";

DefineMathFunction::usage = "DefineMathFunction[f] declares f to be a
scalar-valued function of one real or complex variable that will be
used in tensor expressions.
\nOptions:
\n* Type -> type. This can be a single keyword or a list of keywords,
chosen from among Real, Complex, Imaginary, Automatic, Positive,
Negative, NonPositive, or NonNegative. Default is Real.
\n* Quiet -> True or False. Default is False, or the value of
$Quiet if set.";

DefineRelation::usage = "DefineRelation[ tensor, expr,
cond ] defines a relation of the form
\n	tensor := expr /; cond
\n associated with the tensor\'s name. The \"cond\" is optional. The
first argument \"tensor\" represents a tensor name with or without
indices, and can be written in Ricci's input form. If index names used
in \"tensor\" are associated with bundles, then the relation will only
be applied when indices from those bundles appear in those positions.
The relation will substitute x (suitably modified) in place of
any expression that matches tensor after insertion of indices,
covariant differentiation, conjugation, or raising or lowering of certain
indices."<>
"\nOption:
\nQuiet -> True or False. Default is False, or the value of $Quiet if set.";

DefineRule::usage = "DefineRule[ rulename, lhs, rhs, cond ] defines a
rule named \"rulename\" of the form lhs :> rhs /; cond. The \"cond\"
is optional. The rule is appended to previous rules defined with the
same name, unless the option NewRule -> True is specified. If index
names used on the left-hand side are associated with bundles, then the
rule is applied only when indices from those bundles appear in
those positions. If lhs is a single tensor with or without indices,
then the rule will substitute rhs (suitably modified) in place of any
expression that matches lhs after insertion of indices, covariant
differentiation, conjugation, or raising or lowering of certain
indices."<>
"\nOptions:
\nQuiet -> True or False. Default is False or the value of $Quiet if
set.
\nNewRule -> True or False. Default is False.";

DefineTensor::usage =
"DefineTensor[ name, rank ] defines a tensor.
\n* name: A symbol that uniquely identifies the tensor.
\n* rank: Rank of the tensor. For ordinary tensors, this
must be a non-negative integer. For product tensors, this
is a list of non-negative integers." <>
"\nOptions:
\n* Symmetries -> sym:
Symmetries of the tensor.
Default is NoSymmetries.
\n* Type -> type. For rank-0 tensors, this can be a single keyword or
a list of keywords, chosen from among Complex, Real, Imaginary,
Positive, Negative, NonPositive, or NonNegative. For higher-rank
tensors, it can be Real, Complex, or Imaginary. The default is always
Real.
\n* TeXFormat -> \"texformat\":    Default is the tensor's name." <>
"\n* Bundle -> bundle:  The bundle with
which the tensor is associated. Can be a list, meaning
the tensor is associated with the direct sum of all the
bundles in the list. Can be Any, to indicate that the tensor
will accept indices from any bundle.
Default is $DefaultTangentBundle if
set; otherwise Bundle must be specified.
\n* Variance -> Covariant or Contravariant:  Default is Covariant.
This can be a list whose length is equal to the total rank
of the tensor, in which case each entry in the list specifies
the variance of the corresponding index slot. May be
abbreviated Co and Con.
\n* Quiet -> True or False:  Default is False, or the value
of $Quiet if set."

DefineTensorSymmetries::usage =
"DefineTensorSymmetries[name,{perm1,sgn1,...,permN,sgnN}] defines a
TensorSymmetry object that can be used in DefineTensor. The permi
must be lists that are nontrivial permutations of {1,2,...,d} and the
sgni must be constants (ordinarily plus or minus 1).";

Del::usage = "Del[x] is the total covariant derivative of
the tensor expression x. If x is a k-tensor, then Del[x] is a
k+1-tensor. The last index of Del[x] is the one generated by
differentiation.
\nDel[v,x] is the covariant derivative of x in
the direction v; v must be a 1-tensor expression.
\nOptions:
\n* Connection -> cn specifies that the covariant derivative is to be
taken with respect to the connection cn instead of the default
connection.
\n* Metric -> g: a symmetric 2-tensor expression, indicating that
the covariant derivative is to be taken with respect to the
Levi-Civita connection of g instead of the default metric for x's bundle
(which is assumed to be Riemannian).";

Det::usage = Ricci`Private`AppendMessage[ Det::usage,
"If x is a 2-tensor expression, Ricci interprets Det[x] as the determinant
of x, using the metric to convert x to a {Contravariant,Covariant} tensor
if necessary." ];

Dimension::usage =
"Dimension[bundle] returns the bundle\'s dimension.";

Div::usage = "Div[x] is the divergence of the tensor expression x,
which is the covariant derivative of x contracted on its last two
indices. If x is a k-tensor, Div[x] is a (k-1)-tensor. Div is the
formal adjoint of -Del.
\nOptions:
\n* Connection -> cn specifies that covariant derivatives are to be
taken with respect to the connection cn instead of the default
connection.
\n* Metric -> g: a symmetric 2-tensor expression, indicating that the
divergence is to be taken with respect to the Levi-Civita connection
of g instead of the default metric for x's bundle (which is assumed to
be Riemannian).";

DivGrad::usage = "DivGrad is a value for the global variable
$LaplacianConvention.";

Dot::usage = Ricci`Private`AppendMessage[Dot::usage,
"When x and y are tensor expressions, x.y represents the contraction of
the last index of x with the first index of y."];

Even::usage = "Even is a value for the Type option of DefineConstant.";

ERROR::usage = "If you attempt to plug in the wrong number of
indices in a tensor expression, Ricci will return ERROR[expression].";

Expand::usage = Ricci`Private`AppendMessage[Expand::usage,
"If you apply the Mathematica function Expand to an expression that
includes tensors with indices, Ricci automatically converts it to
TensorExpand."];

Extd::usage = "Extd[x] is the exterior
derivative of x, which must be an alternating
covariant tensor expression. In output form, Extd[x] prints as d[x].";

ExtdStar::usage = "ExtdStar[x] is the adjoint of the operator Extd
applied to x, which must be an alternating covariant tensor
expression.
\nOption:
\n* Metric -> g: a symmetric 2-tensor expression,
representing a metric to be used in place of the default
metric for x's bundle (which is assumed to be Riemannian).";

FactorConstants::usage = "FactorConstants[x] applies the Mathematica
function Factor to the constant factor in each term of the tensor
expression x. FactorConstants[x,n] or FactorConstants[x,{n1,n2,...}]
applies FactorConstants only to term n or to terms n1,n2,... of x.";

FirstBianchiRule::usage = "FirstBianchiRule is a rule that attempts to
turn sums containing two Riemannian curvature tensors into a single
term, using the first Bianchi identity.";

FirstStructureRule::usage = "FirstStructureRule is a rule that implements
the first structure equation for exterior derivatives of basis covectors.";

FirstUp::usage = "FirstUp is a value for the RiemannConvention option
of DefineBundle.";

FlatConnection::usage = "FlatConnection is a DefineBundle option.
\"FlatConnection -> True\" means the default connection has zero
curvature. The default is False.
\nThe function FlatConnection[bundle] returns True or False.";

FormQ::usage =
"FormQ[x] returns True if x is a Covariant alternating tensor
expression, and False otherwise."

Grad::usage =
"Grad[x] is the gradient of the tensor expression x. It is the same
as Del[x], except that the last index is Contravariant instead of
Covariant. If x is a function (0-tensor), then Grad[x] is a vector
field.";

Hermitian::usage =
"Hermitian is a value for the Symmetry option of DefineTensor. It is
valid only for 2-tensors.
A Hermitian
tensor is actually a Real, Symmetric 2-tensor associated with a bundle and
its conjugate, with the additional property that h[L[i],L[j]] and
h[LB[i],LB[j]] are both zero if i and j are associated with the same
bundle.";

HodgeInner::usage =
"HodgeInner[x,y] represents the Hodge inner product of the alternating
tensors x and y. In output form, HodgeInner[x,y]
appears as <<x, y>>.";

HodgeNorm::usage = "HodgeNorm[x] is the norm of the alternating
tensor expression x, with respect to the Hodge inner product. It
is automatically converted by Ricci to Sqrt[HodgeInner[x,Conjugate[x]]].";

Im::usage = Ricci`Private`AppendMessage[Im::usage,
"The Ricci package modifies Im to handle tensor expressions. Im[x]
is converted to (x - Conjugate[x])/(2I)."];

Imaginary::usage = "Imaginary is a value for the Type option of
DefineConstant, DefineTensor, and DefineMathFunction.";

IndexOrderedQ::usage = "IndexOrderedQ[{indices}] returns True if
the indices are ordered correctly according to Ricci's index ordering
rules:  first by name, then by altitude (lower before upper), and
False otherwise.";

IndexQ::usage = "IndexQ[i] returns True if i is an index name, and
False otherwise.";

Inner::usage = Ricci`Private`AppendMessage[Inner::usage, "If x and y are
tensor expressions of the same rank, Ricci interprets Inner[x,y] as the
inner product of x and y. In output form, Inner[x,y] appears as <x, y>."];

InsertIndices::usage = "If x is a tensor expression of rank k, then x[
i1,...,ik ] or InsertIndices[ x, {i1,...,ik} ] causes the indices
i1,...,ik to be inserted into the k index slots in order. If x is a
scalar (rank 0) expression, then InsertIndices[ x, {} ] or x[]
converts it to a component expression by generating dummy indices as
necessary, while InsertIndices[ x, {i1,...,ik} ] is automatically
replaced by CovD[ x, {i1,...,ik} ].";

Int::usage = "If x and y are alternating tensor expressions with
Rank[x] <= Rank[y], Int[x,y] represents the generalized interior
product of x into y. When x is a vector field, Int[x,y] is interior
multiplication of x into y (with a numerical factor depending on
$WedgeConvention). In general, Int[x,_] is the adjoint (with respect
to the Hodge inner product) of wedging with x on the left.";

Integer::usage = Ricci`Private`AppendMessage[Integer::usage, "Ricci
recognizes Integer as a value for the Type option of
DefineConstant."];

Inverse::usage = Ricci`Private`AppendMessage[Inverse::usage, "If x is a
2-tensor expression, Ricci interprets Inverse[x] as the inverse of x."];

Kronecker::usage =
"Kronecker[L[i],U[j]] is the Kronecker delta tensor.";

L::usage = "L[i] represents a lower index i.";

LaplaceBeltrami::usage = "LaplaceBeltrami[x] is the
Laplace-Beltrami operator applied to the differential form x.
It is automatically replaced by Extd[ExtdStar[x]] + ExtdStar[Extd[x]].
\nOption:
\n* Metric -> g: a symmetric 2-tensor expression, representing a
metric to be used in place of the default metric for x's bundle
(which is assumed to be Riemannian).";

Laplacian::usage = "Laplacian[x] is the covariant Laplacian of the
tensor expression x. It is automatically replaced by Div[Grad[x]] if
$LaplacianConvention = DivGrad (the default), and by -Div[Grad[x]] if
$LaplacianConvention = PositiveSpectrum.
\nOptions:
\n* Metric -> g: a symmetric 2-tensor expression, representing a
metric to be used in place of the default metric for the underlying
tangent bundle of x.
\n* Connection -> cn specifies that covariant derivatives are to be
taken with respect to the connection cn instead of the default
connection.";

LB::usage = "LB[i] represents a lower barred index i.";

LeviCivitaConnection::usage = "LeviCivitaConnection[g] represents the
Levi-Civita connection (as a tensor of rank {2,1}, i.e. a matrix of
one-forms) of the arbitrary metric g. When indices are inserted, the
components of the connection are computed in terms of the background
connection of g's bundle (assumed to be Riemannian) and covariant
derivatives of g.";

Lie::usage = "Lie[v,x] is the Lie derivative of the tensor expression
x in the direction v; v must be a vector field (a contravariant 1-tensor
expression).";

LieRule::usage = "LieRule transforms Lie derivatives of differential
forms to expressions involving exterior derivatives and Int.";

LowerAllIndices::usage = "LowerAllIndices[x] lowers all of the
indices in x by inserting appropriate metrics with raised indices.
LowerAllIndices[x,n] or LowerAllIndices[x,{n1,n2,...}] applies
LowerAllIndices only to term n or to terms n1,n2,... of x.";

Method::usage = Ricci`Private`AppendMessage[ Method::usage, "Method is an
option for the Ricci simplification command OrderDummy, which specifies
how hard the command should work to simplify the expression. The
allowable values are 0, 1, and 2."];

Metric::usage = "Metric[bundle] returns the bundle\'s metric.
\nMetric -> metric is an option for Del, Div, Grad, ExtdStar, CovD,
Laplacian, and LaplaceBeltrami.";

MetricQ::usage =
"MetricQ[x] returns True if x is a metric with or without indices,
and False otherwise.";

MetricType::usage = "MetricType is a DefineBundle option,
which specifies whether the bundle's metric has special properties.
Allowable values are MetricType -> Normal (no special
properties), and MetricType ->
Riemannian (for a Riemannian tangent bundle). The
default is Normal. If Riemannian is specified, the following additional
options may also be given:
\n* RiemannTensor -> name: a name to
be given to the Riemannian curvature tensor
for this bundle. Default is Rm.
\n* RicciTensor -> name: a name to be given to the Ricci tensor
for this bundle. Default is Rc.
\n* ScalarCurv -> name: a name to be given to the scalar curvature function
for this bundle. Default is Sc.
\n* RiemannConvention -> SecondUp or FirstUp: determines the sign convention
used for the Riemannian curvature tensor for this bundle. Default is SecondUp,
or the value of the global variable $RiemannConvention if it has been set.";

Negative::usage = Ricci`Private`AppendMessage[Negative::usage, "Ricci
recognizes Negative as a value for the Type option of DefineConstant,
DefineTensor, DefineMathFunction."];

NewDummy::usage = "NewDummy[x] converts all dummy indices occurring in
the tensor expression x to computer-generated dummy indices.
NewDummy[x,n] or NewDummy[x,{n1,n2,...}] applies NewDummy only to term
n or to terms n1,n2,... of x. For most bundles, the dummy names are
of the form \"kn\", where k=BundleDummy[bundle] and n is an integer.
For one-dimensional bundles, only the dummy name k itself is
generated.";

NewRule::usage = "NewRule is a DefineRule option. NewRule -> True
specifies that previous rules defined with the same name are to be
erased before defining the new rule. The default is False, which
means that this rule is to be appended to the list of already-existing
rules by the same name.";

NonNegative::usage = Ricci`Private`AppendMessage[NonNegative::usage, "Ricci
recognizes NonNegative as a value for the Type option of DefineConstant,
DefineTensor, DefineMathFunction."];

NonPositive::usage = Ricci`Private`AppendMessage[NonPositive::usage, 
"NonPositive is a value for the Type option of
DefineConstant, DefineTensor, DefineMathFunction."];

NoOneDims::usage = "NoOneDims is a value for the Mode option of
AbsorbMetrics.";

Norm::usage =  Ricci`Private`AppendMessage[Norm::usage,"Norm[x]
is the norm of the tensor expression x. It
is automatically converted by Ricci to Sqrt[Inner[x,Conjugate[x]]]."];

Normal::usage = Ricci`Private`AppendMessage[Normal::usage, 
"Normal is a value for the MetricType option of
DefineBundle."];

NoSymmetries::usage = "NoSymmetries is a value for the
Symmetries option of DefineTensor.";

Odd::usage = "Odd is a value for the Type option of DefineConstant.";

OneDims::usage = "OneDims is an AbsorbMetrics option. OneDims -> False
means that AbsorbMetrics should not absorb one-dimensional metrics
(unless they are contracted with other metrics). Default is True.";

OrderCovD::usage = "OrderCovD[x] orders all of the indices appearing
after \";\" in the tensor expression x, by adding appropriate
curvature and torsion terms. OrderCovD[x,n] or
OrderCovD[x,{n1,n2,...}] applies OrderCovD only to term n or to terms
n1,n2,... of x. ";

OrderDummy::usage = "OrderDummy[x] attempts to put the dummy indices
occurring in the tensor expression x in a \"canonical form\".
OrderDummy[x,n] or OrderDummy[x,{n1,n2,...}] applies OrderDummy only
to term n or to terms n1,n2,... of x. Dummy index pairs are ordered
so that the lower member appears first whenever possible. OrderDummy
tries various rearrangements of dummy index names, and chooses the
lexically smallest version of the expression that results." <>
"\nOption:
\n* Method -> n specifies how hard OrderDummy should work to find the
best possible version of the expression. The default is Method -> 1,
which means that dummy index names will be interchanged in pairs only.
Method -> 2 causes OrderDummy to try all possible permutations of the
dummy index names. Method -> 0 means don't try interchanging names at
all.";

OrthonormalFrame::usage = "OrthonormalFrame is an option for
DefineBundle. If True, it specifies that the metric coefficients are
always assumed to be constants in the default basis. For
one-dimensional bundles, the metric coefficient is always taken to be
1. Default is False";

ParallelFrame::usage = "ParallelFrame is an option for DefineBundle.
If True, it specifies that the default basis is always assumed to be
parallel. Thus the connection and curvature forms will always be
zero, and the bundle will be flat. Default is False.";

Plus::usage = Ricci`Private`AppendMessage[Plus::usage, "Ricci
transforms expressions of the form (a + b)[L[k],...] into
InsertIndices[a+b,{L[k],...}]."];

Positive::usage = Ricci`Private`AppendMessage[Positive::usage, "Ricci
recognizes Positive as a value for the Type option of DefineConstant,
DefineTensor, DefineMathFunction."];

PositiveDefinite::usage = "PositiveDefinite is a DefineBundle option.
PositiveDefinite -> False means that the bundle's metric is not assumed
to be positive definite. Default is True.";

PositiveInteger::usage = "PositiveInteger is an obsolete value for the
Type option of DefineConstant. Use Type -> {Positive,Integer}
instead.";

PositiveInteger := (Message[PositiveInteger::usage];{Positive,Integer});

PositiveReal::usage = "PositiveReal is an obsolete value for the Type
option of DefineConstant, DefineTensor, and DefineMathFunction. Use
Type -> {Positive,Real} instead.";

PositiveReal := (Message[PositiveReal::usage];{Positive,Real});

PositiveSpectrum::usage = "PositiveSpectrum is a value for the global
variable $LaplacianConvention.";

Power::usage = Ricci`Private`AppendMessage[Power::usage, "If x is a
tensor expression and p is a positive integer, Ricci interprets x^p as
the p-th symmetric power of x with itself. \nRicci transforms
expressions of the form (x^p)[L[i],...] into
InsertIndices[x^p,{L[i],...}] or CovD[x^p,{L[i],...}]." <> "\nRicci
causes a power of a product such as (a b)^p to be expanded into a
product of powers a^p * b^p, provided a and b do not contain any
indices; if they do contain indices, then (a b)^p is transformed to
Summation[a b]^p, to prevent the expression from being expanded to a^p
b^p. Summation is not printed in output form."];

PowerSimplify::usage = "PowerSimplify[x] attempts to simplify negative and
nonintegral powers that appear in the tensor expression x, by
expanding and collecting constants in the base and the exponent of
each such power. PowerSimplify[x,n] or PowerSimplify[x,{n1,n2,...}] applies
PowerSimplify only to term n or to terms n1,n2,... of x.";

ProductExpand::usage = "ProductExpand[x] expands out symmetric
products and wedge products of 1-tensors that occur in x, and rewrites
them in terms of tensor products.";

Quiet::usage = "Quiet is an option for some of the defining and
undefining commands in the Ricci package. The option Quiet -> True will
silence the usual messages printed by these commands. The default is
the value of the global variable $Quiet, which is initially False."

Rank::usage = "Rank[x] returns the rank of the tensor expression x.";

Re::usage = Ricci`Private`AppendMessage[Re::usage,
"The Ricci package modifies Re to handle tensor expressions.
Re[x] is converted to (x + Conjugate[x])/2."];

Real::usage = Ricci`Private`AppendMessage[Real::usage, "Ricci
recognizes Real as a value for the Type option of DefineBundle,
DefineConstant, DefineTensor, and DefineMathFunction."];

RenameDummy::usage = "RenameDummy[x] changes the names of dummy
indices in x to standard names. RenameDummy[x,n] or
RenameDummy[x,{n1,n2,...}] applies RenameDummy only to term n or to
terms n1,n2,... of x. RenameDummy chooses names in alphabetical order
from the list of index names associated with the appropriate bundle,
skipping those names that already appear in x as free indices. When
the list of index names is exhausted, computer-generated names of the
form \"kn\" are used, where k is the last index name in the list and n
is an integer. For one-dimensional bundles, n is omitted.";

Ricci::usage = "Ricci is a Mathematica package for doing symbolic
computations that arise in differential geometry. To load it, type
\"<<Ricci.m\". To find out more information about a Ricci command,
object, or function, type \"?name\".";

RicciSave::usage = "RicciSave[\"filename\"] causes all of your current Ricci
definitions to be saved into the file \"filename\". The previous
contents of \"filename\" are erased. The definitions can be read back
in by the command <<filename.";

RicciTensor::usage = "RicciTensor -> name is a DefineBundle option,
which specifies a name to be used for the Ricci tensor for this
bundle. The default is Rc. This option takes effect only when
MetricType -> Riemannian has been specified.
\nRicciTensor[g] represents the Ricci curvature of the arbitrary
metric g.";

RiemannConvention::usage = "RiemannConvention -> SecondUp or FirstUp
is a DefineBundle option, which specifies the index convention of the
Riemannian curvature tensor. The default is SecondUp, or the value of
the global variable $RiemannConvention if it has been set. This
option takes effect only when MetricType -> Riemannian has been
specified.";

Riemannian::usage =
"MetricType -> Riemannian is a DefineBundle option for defining a
Riemannian tangent bundle. If this option
is specified, the following additional
options may be given:
\n* RiemannTensor -> name: the name for the Riemannian curvature
tensor for this bundle. Default is Rm.
\n* RicciTensor -> name: the name for the Ricci curvature
for this bundle. Default is Rc.
\n* ScalarCurv -> name: the name for the scalar
curvature for this bundle. Default is Sc.
\n* RiemannConvention -> FirstUp or SecondUp: controls the index convention
of the Riemannian curvature tensor. Default is SecondUp.";

RiemannSymmetries::usage = "RiemannSymmetries is a value for the
Symmetries option of DefineTensor. Tensors having RiemannSymmetries
must be of degree 4. The symmetries are those of the Riemann
curvature tensor, namely skew-symmetric in the first two indices,
skew-symmetric in the last two indices, symmetric when the first two
and last two are interchanged. The (first) Bianchi identity can be
explicitly applied to any tensor with these symmetries by means of
FirstBianchiRule or BianchiRules.";

RiemannTensor::usage = "RiemannTensor -> name is a DefineBundle
option, which specifies a name to be used for the Riemannian curvature
tensor for this bundle. The default is Rm. This option takes effect
only when MetricType -> Riemannian has been specified.
\nRiemannTensor[g] represents the Riemannian curvature tensor of the
arbitrary metric g.";

Same::usage = "Bundle -> Same is a special DefineTensor option for
product tensors. It means that indices from any bundle can be
inserted in the first set of index slots, but indices in the remaining
slots must be from the same bundle or its conjugate.";

ScalarCurv::usage = "ScalarCurv -> name is a DefineBundle option,
which specifies a name to be used for the scalar curvature for this
bundle. The default is Sc. This option takes effect only when
MetricType -> Riemannian has been specified.
\nScalarCurv[g] represents the scalar curvature of the arbitrary
metric g.";

ScalarQ::usage = "ScalarQ[x] is True if x is a rank 0 tensor
expression with no free indices, and False otherwise.";

SecondBianchiRule::usage = "SecondBianchiRule is a rule that attempts
to turn sums containing two differentiated Riemannian curvature
tensors into a single term using the second Bianchi identity.";

SecondStructureRule::usage = "SecondStructureRule is a rule that
implements the second structure equation for exterior derivatives of
the generic connection forms.";

SecondUp::usage = "SecondUp is a value for the RiemannConvention option
of DefineBundle.";

SimplifyConstants::usage = "SimplifyConstants[x] applies the
Mathematica function Simplify to the constant factor in each term of
the tensor expression x. SimplifyConstants[x,n] or
SimplifyConstants[x,{n1,n2,...}] applies SimplifyConstants only to
term n or to terms n1,n2,... of x.";

Skew::usage = "Skew is a synonym for Alternating.";

SkewHermitian::usage = "SkewHermitian is a value for the Symmetry
option of DefineTensor. A SkewHermitian tensor is actually a Real,
Alternating tensor associated with a bundle and its conjugate, with
the additional property that h[L[i],L[j]] and h[LB[i],LB[j]] are both
zero if i and j are associated with the same bundle.";

SkewQ::usage = "SkewQ[x] is True if x is an alternating or
skew-Hermitian tensor expression, and False otherwise.";

StructureRules::usage = "StructureRules is the union of
CompatibilityRule, FirstStructureRule, and SecondStructureRule.";

Summation::usage = "If a and b contain indices, then Ricci transforms (a
b)^p internally to Summation[a b]^p, to prevent the expression from
being expanded to a^p b^p. Summation is not printed in output
form.";

SuperSimplify::usage = "SuperSimplify[x] attempts to put the tensor
expression x into a canonical form, so that two expressions that are
equal will usually be identical. SuperSimplify[x,n] or
SuperSimplify[x,{n1,n2,...}] applies SuperSimplify only to term n or
to terms n1,n2,... of x. SuperSimplify works exactly the same way as
TensorSimplify, except that it calls OrderDummy with the option Method
-> 2, so that all possible permutations of dummy index names will be
tried. This will generally be much slower for expressions having more
than 4 dummy index pairs.";

Sym::usage = "Sym[x] represents the symmetrization of the tensor
expression x.";

Symmetric::usage = "Symmetric is a value for the Symmetries option
of DefineTensor.";

SymmetricProduct::usage = "If x, y, and z are tensor expressions,
SymmetricProduct[x,y,z] or x * y * z or x y z represents their symmetric
product. Symmetric products are represented internally by ordinary
multiplication. In output form, Ricci inserts explicit asterisks
whenever non-scalar tensor expressions are multiplied together, to
remind the user that multiplication is being interpreted as symmetric
product.";

SymmetricQ::usage = "SymmetricQ[x] is True if x is a symmetric or
Hermitian tensor expression, and False otherwise.";

Symmetries::usage = "Symmetries -> sym is a DefineTensor option.
Allowed values are NoSymmetries, Symmetric, Alternating, Hermitian,
SkewHermitian, RiemannSymmetries, or a symmetry defined with
DefineTensorSymmetries. Skew is a synonym for Alternating. Hermitian
and SkewHermitian can be used only with 2-tensors. For product
tensors, sym can be a list of symmetries, one for each rank in the
rank list.";

TangentBundle::usage = "TangentBundle is a DefineBundle option. If
bundle x is defined with the option \"TangentBundle -> {y,z}\", this
means that the underlying tangent bundle of x is the direct sum of y
and z. The form \"TangentBundle -> y\" may be used when there is only
one bundle. If $DefaultTangentBundle is defined, it will be used as
the default when TangentBundle is not specified; otherwise by default
the tangent bundle of x is assumed to be x itself.
\nThe function
TangentBundle[x] returns the tangent bundle list for the bundle x.";

Tensor::usage = "Internally, Ricci represents tensors in the form
Tensor[name,{i,j,...},{k,l,...}], where i,j,... are the tensor indices
and k,l,... are indices resulting from covariant differentiation. In
input form, an unindexed tensor is represented just by typing its
name. Indices are inserted by typing them in brackets after the
tensor name. Once all index slots are full, indices resulting from
covariant differentiation can be typed in a second set of brackets.
Tensors are created by DefineTensor and removed by UndefineTensor.";

TensorCancel::usage = "TensorCancel[x] attempts to simplify each term
of x by canceling common factors, even when the factors have
different names for their dummy indices. TensorCancel[x,n] or
TensorCancel[x,{n1,n2,...}] applies TensorCancel only to term n or to
terms n1,n2,... of x.";

TensorData::usage = "TensorData[name] is a list containing data for
the tensor \"name\", used internally by Ricci.";

TensorExpand::usage = "TensorExpand[x] expands products and
positive integral powers in x, just as Expand does, but maintains
correct dummy index conventions and does not expand constant factors.
TensorExpand[x,n] or TensorExpand[x,{n1,n2,...}] applies TensorExpand
only to term n or to terms n1,n2,... of x.";

TensorFactor::usage = "TensorFactor[expr] returns the product of all the
non-constant factors in expr, which should be a monomial.";

TensorMetricQ::usage = "TensorMetricQ[tensorname] returns True if
tensorname is the metric of some bundle, and False otherwise.";

TensorProduct::usage =
"TensorProduct[x,y,z] or TProd[x,y,z] or x ~TProd~ y ~TProd~ z
represents the tensor product of x, y, and z.
In output form, tensor products appear as
\n   x (X) y (X) z.";

TensorQ::usage = "TensorQ[name] returns True if name is the name of a
tensor, and False otherwise.";

TensorRankList::usage = "TensorRankList stores the list of ranks for a
product tensor. Used internally by Ricci.";

TensorSimplify::usage = "TensorSimplify[x] attempts to put the tensor
expression x into a canonical form, so that two expressions that are
equal will usually be identical. TensorSimplify[x,n] or
TensorSimplify[x,{n1,n2,...}] applies TensorSimplify only to term n or
to terms n1,n2,... of x. TensorSimplify expands products and positive
integer powers, uses metrics to raise and lower indices, tries to
rename all dummy indices in a canonical order, and collects all terms
containing the same tensor factors but different constant factors. It
does not reorder indices after the \";\" (use OrderCovD or
CovDSimplify to do that). TensorSimplify calls CorrectVariances,
TensorExpand, AbsorbMetrics, PowerSimplify, RenameDummy, OrderDummy, and
CollectConstants.";

TensorSymmetry::usage = "A tensor symmetry is an object of the form
TensorSymmetry[name,d,{perm1,sgn1,...,permN,sgnN}], where d is a positive
integer, permi is a non-trivial permutation of {1,2,...,d}, and sgni is
plus or minus 1. TensorSymmetry objects can be defined with
DefineTensorSymmetries and used in the
DefineTensor command.";

TeXFormat::usage = "TeXFormat -> \"string\" is an option for DefineTensor
and DefineIndex. It specifies the way the tensor or index name will
appear in TeXForm. The default is the tensor's or index's name.";

Times::usage = Ricci`Private`AppendMessage[Times::usage,
"Ricci uses ordinary multiplication
to represent symmetric products of tensors.
Ricci transforms an
expression of the form (a * b)[i]
into InsertIndices[a*b,{i}] whenever \"i\" is
one or more indices (L[j] or U[j]).
\nIn output form, Ricci modifies Mathematica's usual
ordering of factors: constants are printed first, then fully-indexed tensors,
then other tensor expressions."];

Tor::usage = "Tor is the name for the generic torsion tensor for the
default connection in any bundle. It is a product tensor of rank
{2,1}; inserting the first two indices yields the torsion 1-forms
associated with the default basis. Inserting all three indices yields
the components of the torsion tensor.";

TorsionFree::usage = "TorsionFree is a DefineBundle option.
\"TorsionFree -> False\" means the default connection has nonvanishing
torsion. The default is True.
\nThe function TorsionFree[bundle]
returns True or False.";

TotalRank::usage = "TotalRank[tensorname] returns the total rank of
the tensor tensorname.";

TProd::usage = "TProd is an abbreviation for TensorProduct.";

Tr::usage = Ricci`Private`AppendMessage[ Tr::usage,
"If x is a 2-tensor expression, Ricci interprets
Tr[x] as the trace of the x with respect to the metric of x's bundle." ];

Transpose::usage = Ricci`Private`AppendMessage[ Transpose::usage, "If x
is a tensor expression, Ricci interprets Transpose[x] as x with its
index positions reversed." ];

Type::usage = "Type is an option for DefineBundle, DefineConstant,
DefineTensor, and DefineMathFunction.";

U::usage = "U[i] represents an upper index i.";

UB::usage = "UB[i] represents an upper barred index i.";

UndefineBundle::usage = "UndefineBundle[bundle] clears the definition
of bundle. It also clears the definition of the bundle's metric and
indices.
\nOption:
\n* Quiet -> True or False. Default is False, or the value of $Quiet if set.";

UndefineConstant::usage =
"UndefineConstant[symbol] removes symbol's definition as a constant.
\nOption:
\n* Quiet -> True or False. Default is False, or the value of $Quiet if set.";

UndefineIndex::usage =
"UndefineIndex[i] or UndefineIndex[{i,j,k}] removes the association
of indices with their bundles.
\nOption:
\n* Quiet -> True or False. Default is False, or the value of $Quiet if set.";

UndefineRelation::usage =
"UndefineRelation[tensor] deletes the relation previously defined for tensor.
The tensor must exactly match the first argument of the corresponding call
to DefineRelation.
\nOptions:
\n* Quiet -> True or False. Default is False, or the value of $Quiet if set.";

UndefineTensor::usage = "UndefineTensor[tensorname] clears the definition
of tensorname.
\nOption:
\n* Quiet -> True or False. Default is False, or the value of $Quiet
if set.";

UndefineTensorSymmetries::usage = "UndefineTensorSymmetries[name] deletes
the TensorSymmetry object created by DefineTensorSymmetries[name,...].";

UnderlyingTangentBundle::usage = "UnderlyingTangentBundle[x] returns a
list of bundles representing the underlying tangent bundle of the
expression. The tangent bundle is the direct sum of the bundles in
the list. It is assumed that all tensors used in a given expression
have the same underlying tangent bundle.";

Variance::usage = Ricci`Private`AppendMessage[ Variance::usage,
"Variance -> Covariant or Contravariant is a
DefineTensor option. May be a list whose length is equal to the total
rank of the tensor, in which each entry specifies the variance of the
corresponding index slot. May be abbreviated Co and Con. Default is
Covariant.
\nIf x is a tensor expression, Variance[x] returns a list of variances of
x, one for each index slot." ];

VectorFieldQ::usage = "VectorFieldQ[x] returns True if x is a
Contravariant 1-tensor expression, and False otherwise.";

Wedge::usage = "If x, y, and z are alternating tensor expressions,
then x ~Wedge~ y ~Wedge~ z or Wedge[x,y,z] represents their wedge
product. The interpretation of Wedge in terms of tensor products is
determined by the global variable $WedgeConvention. In output form,
wedge products print as x ^ y ^ z.";

(************************* Global variables ***********************)

$DefaultTangentBundle::usage = "The global variable
$DefaultTangentBundle can be set by the user to a bundle or list of
bundles. It is used by DefineBundle as the default value for the
TangentBundle option, and by DefineTensor as the default value for the
Bundle option. By default, it is set to the first bundle the user
defines (or the direct sum of this bundle and its conjugate if the bundle
is complex).";

$LaplacianConvention::usage = "$LaplacianConvention determines which sign
convention is used for the covariant Laplacian on functions and tensors.
$LaplacianConvention = DivGrad means that Laplacian[x] = Div[Grad[x]],
while $LaplacianConvention = PositiveSpectrum means that Laplacian[x] =
-Div[Grad[x]]. Default is DivGrad.";

$MathFunctions::usage = "$MathFunctions is a list of names that have
been defined as scalar mathematical functions for use by Ricci.";

$Quiet::usage = "The global variable $Quiet is used by the Ricci package
to determine whether the defining and undefining commands report on
what they are doing. Setting $Quiet=True will silence these chatty
commands. The default is False. It can be overridden for a
particular command by specifying the Quiet option as part of the
command call.";

$RiemannConvention::usage = "The global variable $RiemannConvention
can be set by the user to specify a default value for the
RiemannConvention option of DefineBundle and MakeRiemannian. Default
is SecondUp.";

$TensorFormatting::usage = "The global variable $TensorFormatting can
be set to True or False by the user to turn on or off Ricci's special
output formatting of tensors and indices. Default is True.";

$TensorTeXFormatting::usage = "The global variable
$TensorTeXFormatting can be set to True or False by the user to turn
on or off Ricci's special formatting of tensors in TeXForm. Default is
True.";

$WedgeConvention::usage = "The global variable $WedgeConvention can be
set by the user to determine the interpretation of wedge products.
The allowed values are Alt and Det. The default is Alt.";

(********************* Error messages *************************)

General::invalid = "\"`1`\" is not a valid `2`.";

Index::error = "Incorrect number of indices provided for `1`.";

Bundle::error = "No bundle specified or implied for `1`.";

General::oneterm = "`1` has only one term.";

General::termspec = "`1` is not a valid term number or list
of term numbers.";

(************ Initialize default values for user parameters *************)

$Quiet=False;
$DefaultTangentBundle = Null;
$WedgeConvention = Alt;
$RiemannConvention = SecondUp;
$TensorFormatting = True;
$TensorTeXFormatting = True;
$MathFunctions = {};
$LaplacianConvention = DivGrad;

(******************** Switch to private context ************************)

Begin["`Private`"];

Unprotect[TensorQ];
SetAttributes[ { TensorQ,Tensor,TensorData,TensorMetricQ,TotalRank,rank }, 
 HoldFirst];
SetAttributes[ ERROR, HoldAll ];
Protect[TensorQ];

(********************** General utility routines **************************)

(** MakeLinear can be used to make a function behave linearly with respect
    to sums and constant multiples. **)

MakeLinear[f_Symbol] := (
  HoldPattern[ f[(c_?ConstantQ) * v_] ] := c f[v];
  HoldPattern[ f[x_Plus] ] := f /@ x;
  );

(* PartAt[ exp, pos ] behaves like Part[exp,pos], except that "pos" is
   a list.  *)

PartAt[ exp_, {pos___Integer} ] := Part[exp,pos];
HeldPartAt[ exp_, {pos___Integer} ] := HeldPart[exp,pos];

(* Occurrences[exp,pattern] returns a list of the actual occurrences
   of "pattern" in exp. *)

Occurrences[ exp_, pattern_ ] := PartAt[ exp, # ]& /@ Position[exp,pattern];

(* ContainsQ[list1,list2] is True if list2 is a subset of list1, 
   False otherwise. *)

ContainsQ[ a_, b_ ] := And @@ (MemberQ[a,#]& /@ b);

(** optionsOK is an internal subroutine to check the validity of optional
   arguments.  optionsOK[command, options] returns True if the "options"
   represents a valid list of options for "command", and prints a message
   and returns False if not. **)

optionsOK[subr_,opts___] := Module[{badoptions},
 badoptions = Select[{opts},
   !(Head[#]===Rule && MemberQ[optionNames[Options[subr]],#[[1]]])&];
 If[badoptions==={},
   (*then OK*)    Return[True],
   (*else error*) Message[
                    Evaluate[ToExpression[ToString[subr]<>"::invalid"]],
                    badoptions[[1]],
                    ToString[subr]<>" option"];
                  Return[False]
   ]];

(** optionNames takes a list of rules and gives back the list of
    left-hand sides.  For example:
    optionNames[{a->b,c->d}] = {a,c}.
**)

optionNames[optionlist_] := #[[1]]& /@ optionlist;

(* newsymQ checks that a symbol hasn't been used before as a tensor or bundle
   name, and that it hasn't been assigned a value. *)

SetAttributes[newsymQ,HoldFirst];

newsymQ[Conjugate[x_]] := newsymQ[x];
newsymQ[x_] :=
  (!ValueQ[x] && Head[x]===Symbol && Attributes[x]==={} &&
   !BundleQ[x] && !IndexQ[x]);

tfQ[True|False] := True;
tfQ[_] := False;

(************************** RicciSave **********************************)

RicciSave[] := Message[RicciSave::argx,"RicciSave",0];
RicciSave[_,extras__] := Message[RicciSave::argx,"RicciSave",Length[{extras}]+1];

RicciSave[file_Symbol] := RicciSave[ToString[file]];

RicciSave[file_String] := Module[{tform,protectlist},
  Put @@ Append[ deflist["Ricci`$"], file ];
  tform = $TensorFormatting;
  $TensorFormatting = False;
  protectlist = Unprotect[Evaluate[namelist[$Context]]];
  PutAppend @@ Append[ deflist[$Context], file ];
  PutAppend[ "DefineMathFunction[$MathFunctions,Quiet->True]"//OutputForm,
             file ];
  Unprotect[ Curv,Tor,Conn,Basis,Kronecker];
  PutAppend[ "Unprotect[Curv,Tor,Conn,Basis,Kronecker]"//OutputForm, file];
  PutAppend[ Definition[Curv,Tor,Conn,Basis,Kronecker], file ];
  PutAppend[ "Protect[Curv,Tor,Conn,Basis,Kronecker]"//OutputForm, file];
  PutAppend[ "Protect["//OutputForm, file];
  PutAppend[ protectlist, file];
  PutAppend[ "]"//OutputForm, file];
  $TensorFormatting = tform;
  Protect[ Curv,Tor,Conn,Basis,Kronecker];
  Protect[ Evaluate[protectlist]];
  ];

RicciSave[x_] := Message[RicciSave::invalid, x, "file name"];

namelist[prefix_] := Names[ StringJoin[prefix, "*"] ];

deflist[prefix_] :=
  Select[ Definition /@ namelist[prefix],
          ToString[#] =!= "Null" &];

(*************************** Bundle *****************************)

SetAttributes[{TensorMetricQ,DefineBundle},HoldFirst];

(* MetricQ[exp] gives True iff exp is a metric 2-tensor with or without
   indices.  If exp is a tensor, this is translated into 
   TensorMetricQ, which is a HoldFirst function so that it
   can be associated with the name of the tensor.  *)

MetricQ[Tensor[g_,{_,_},{}]] := TensorMetricQ[g];
MetricQ[Tensor[g_,{},{}]] := TensorMetricQ[g];
MetricQ[_] = False;

(* Initialize TensorMetricQ[anything] to be False. 
   DefineBundle will set this to True for any tensors that are metrics.  *)

TensorMetricQ[_] = False;

(* BundleQ[name] = True iff name is a bundle.  Initially false for any 
   name, it is set to True by DefineBundle. *)

BundleQ[_] = False;

(* Tell the various bundle information functions how to deal with conjugate
   bundles. *)
BundleQ[Conjugate[b_]] := BundleQ[b];
Dimension[Conjugate[b_]] := Dimension[b];
Metric[Conjugate[b_]] := Metric[b];
TorsionFree[Conjugate[b_]] := TorsionFree[b];
FlatConnection[Conjugate[b_]] := FlatConnection[b];
ParallelFrame[Conjugate[b_]] := ParallelFrame[b];
OrthonormalFrame[Conjugate[b_]] := OrthonormalFrame[b];
CommutingFrame[Conjugate[b_]] := CommutingFrame[b];
Type[Conjugate[b_]] := Type[b];
BundleIndices[Conjugate[b_]] := Conjugate /@ BundleIndices[b];
TangentBundle[Conjugate[b_]] := TangentBundle[b];
PositiveDefinite[Conjugate[b_]] := PositiveDefinite[b];

Dimension[Any] := Plus @@ Dimension /@ Flatten[{$DefaultTangentBundle}];

(* These are needed by UnderlyingTangentBundle. *)
TangentBundle[Any] := $DefaultTangentBundle;
TangentBundle[Null] := {Null};

(* The following are needed when computing conjugates of bundles *)

Unprotect[{Null,Automatic}];
Any/:       Conjugate[Any] := Any;
Null/:      Conjugate[Null] := Null;
Automatic/: Conjugate[Automatic] := Automatic;
Protect[{Null,Automatic}];

(* Initialize BundleIndices[any bundle] to an empty list *)
BundleIndices[_] := {};

(* BundleDummy returns the last index in BundleIndices *)
BundleDummy[bun_] := Last @ BundleIndices @ bun;

(*********** argument-checking subroutines for DefineBundle **********)

dimQ[d_] := ConstantQ[d] && (!NumberQ[d] || (IntegerQ[d] && d > 0));

decompQ[Automatic] = True;
decompQ[Null] = True;
decompQ[x_List] := And @@ ((BundleQ[#] || newsymQ[#])& /@ x);
decompQ[x_] := decompQ[{x}];

indQ[{i__},name_]:= And @@ (indQ[#,name]&) /@ {i}; 
indQ[i_Symbol,name_] := newsymQ[i] || (IndexQ[i] && Bundle[i]===name);
indQ[_,_] := False;

mtypeQ[t_] := t===Normal || t===Riemannian;

rmconvQ[FirstUp] = True;
rmconvQ[SecondUp] = True;
rmconvQ[_] = False;

(* Define the default names for Riemannian curvature tensors.
   These should not be evaluated until used, so that the names
   are defined in the user's context. *)

defaultRm := ToExpression["Rm"];
defaultRc := ToExpression["Rc"];
defaultSc := ToExpression["Sc"];

(***************** DefineBundle ***********************************)

Options[DefineBundle] :=
  {Type -> Real,
   TangentBundle -> Automatic,
   FlatConnection -> False,
   TorsionFree -> True,
   Quiet -> $Quiet,
   RiemannTensor -> Null,
   RicciTensor -> Null,
   ScalarCurv -> Null,
   RiemannConvention -> $RiemannConvention,
   MetricType -> Normal,
   BasisName -> Automatic,
   CoBasisName -> Automatic,
   ParallelFrame -> False,
   OrthonormalFrame -> False,
   CommutingFrame -> False,
   PositiveDefinite -> True};

DefineBundle[]      := Message[DefineBundle::argm,"DefineBundle",0,4];
DefineBundle[_]     := Message[DefineBundle::argmu,"DefineBundle",4];
DefineBundle[_,_]   := Message[DefineBundle::argm,"DefineBundle",2,4];
DefineBundle[_,_,_] := Message[DefineBundle::argm,"DefineBundle",3,4];

DefineBundle[name_, dim_, metric_, indices_, opts___] := 
  Module[{type,compat,tdecomp,flat,tfree,quiet,rm,rc,sc,rmc,mtype,
          riemanntensor,riccitensor,scalarcurv,bname,cobname,pframe,
          oframe,cframe,pdef},
  If [!optionsOK[DefineBundle,opts], Return[]];
  {type,tdecomp,flat,tfree,quiet,rm,rc,sc,
   rmc,mtype,bname,cobname,pframe,oframe,cframe,pdef} =
    {Type,TangentBundle,FlatConnection,TorsionFree,Quiet,
     RiemannTensor,RicciTensor,ScalarCurv,RiemannConvention,
     MetricType,BasisName,CoBasisName,ParallelFrame,
     OrthonormalFrame,CommutingFrame,PositiveDefinite} /.
    {opts} /. Options[DefineBundle];
  Which[
   !newsymQ[name], 
     Message[DefineBundle::invalid,HoldForm[name],"bundle name"],
   !newsymQ[metric],
     Message[DefineBundle::invalid,metric,"metric name"],
   !dimQ[dim],
     Message[DefineBundle::invalid,dim,"dimension"],
   !indQ[indices,name],
     Message[DefineBundle::invalid,indices,"list of indices"],
   dim===1 && !Length[Flatten[{indices}]]===1,
     Print["ERROR: a one-dimensional bundle can have only one index name"],
   !MemberQ[{Real,Complex}, type],
     Message[DefineBundle::invalid,type,"Type"],
   !decompQ[tdecomp],
     Message[DefineBundle::invalid,tdecomp,"tangent bundle"],
   !tfQ[flat],
     Message[DefineBundle::invalid,flat,"value for FlatConnection"],
   !tfQ[tfree],
     Message[DefineBundle::invalid,tfree,"value for TorsionFree"],
   !mtypeQ[mtype],
     Message[DefineBundle::invalid,mtype,"MetricType"],
   !rmconvQ[rmc],
     Message[DefineBundle::invalid,rmc,"RiemannConvention"],
   rm =!= Null && !newsymQ[Evaluate[rm]],
     Message[DefineBundle::invalid,rm,"RiemannTensor name"],
   mtype === Riemannian && rm === Null && !newsymQ[Evaluate[defaultRm]],
     Message[DefineBundle::invalid,defaultRm,"RiemannTensor name"],
   rc =!= Null && !newsymQ[Evaluate[rc]],
     Message[DefineBundle::invalid,rc,"RicciTensor name"],
   mtype === Riemannian && rc === Null && !newsymQ[Evaluate[defaultRc]],
     Message[DefineBundle::invalid,defaultRc,"RicciTensor name"],
   sc =!= Null && !newsymQ[Evaluate[sc]],
     Message[DefineBundle::invalid,sc,"ScalarCurv name"],
   mtype === Riemannian && sc === Null && !newsymQ[Evaluate[defaultSc]],
     Message[DefineBundle::invalid,defaultSc,"ScalarCurv name"],
   bname =!= Automatic && !newsymQ[Evaluate[bname]],
     Message[DefineBundle::invalid,bname,"BasisName"],
   cobname =!= Automatic && !newsymQ[Evaluate[bname]],
     Message[DefineBundle::invalid,cobname,"CoBasisName"],
   !tfQ[pframe],
     Message[DefineBundle::invalid,pframe,"value for ParallelFrame"],
   !tfQ[oframe],
     Message[DefineBundle::invalid,oframe,"value for OrthonormalFrame"],
   !tfQ[cframe],
     Message[DefineBundle::invalid,cframe,"value for CommutingFrame"],
   !tfQ[pdef],
     Message[DefineBundle::invalid,pdef,"value for PositiveDefinite"],
   (* Now all the arguments are OK.  Let's define the bundle. *)
   True,
     name/: BundleQ[name] = True;
     If[ type===Real, (*then*) (name/: Conjugate[name] = name)];
     If[ newsymQ[dim],
         DefineConstant[dim,Type->{Positive,Integer}] ];
     name/: Dimension[name] = dim;
     (* If $DefaultTangentBundle hasn't been set, then set it as best
        we can from the information available. *)
     Which[ $DefaultTangentBundle =!= Null,
              Null, (* do nothing *)
            tdecomp =!= Automatic,
              $DefaultTangentBundle = tdecomp,
            True,
              $DefaultTangentBundle = Union[{name,Conjugate[name]}] ];
     (* Flatten appears below to make sure we get a list, even if
        the user just specified a single bundle *)
     name/: TangentBundle[name] = Flatten[{tdecomp} /.
                                  Automatic -> $DefaultTangentBundle ];
     name/: Type[name] = type;
     (* Associate the indices with the bundle *)
     DefineIndex[indices,name,Quiet->quiet];
     (* Define the metric as a symmetric or hermitian 2-tensor *)
     DefineTensor[metric,2,Type->Real,
                  Symmetries->If[type===Real, Symmetric, Hermitian],
                  Bundle->name, Quiet->quiet ];
     name/: Metric[name] = metric;
     (* Tensors that are metrics have additional properties. *)
     defBunMetric[metric,name];
     If[ bname =!= Automatic || cobname =!= Automatic,
         (*then*) defineBasisNames[bname,cobname,name,quiet]];
     name/: MetricType[name] = mtype;
     (* If the user requests, make this bundle Riemannian.  If we use
        default names for the curvature tensors, create them in the
        user's context so they'll be saved with the user's stuff. *)
     If[ mtype===Riemannian,
           (*then*)
           riemanntensor = If[ rm=!=Null, rm, defaultRm];
           riccitensor   = If[ rc=!=Null, rc, defaultRc];
           scalarcurv    = If[ sc=!=Null, sc, defaultSc];
           makeRiemannian[ name, riemanntensor, riccitensor, scalarcurv, 
                           rmc, quiet ] ]; 
     (* Now protect the bundle name against accidental modification *)
     Protect[name];
     DeclareBundle[name, TorsionFree -> tfree,     
                         FlatConnection -> flat,
                         ParallelFrame -> pframe,  
                         OrthonormalFrame -> oframe,
                         CommutingFrame -> cframe, 
                         PositiveDefinite -> pdef];
     (* Tell the user what we've done. *)
     If[!quiet,
       Print["Bundle ",name," defined."];
       Print["   Metric = ",Metric[name],
                 "   Dimension = ",Dimension[name],
                 "   Indices = ",BundleIndices[name]];
       Print["   Bundle Type = ",Type[name],
         "   Metric Type = ",mtype];
       Print["   Tangent Bundle = ", TangentBundle[name]];
       Print["   Connection is ",If[flat, "flat, ", ""],
                If[!tfree, "not ", ""],"torsion free."];
       ];
     Return[name]
  ]];

(* The following function defines the extra rules needed by metric 
   tensors. *)

defBunMetric[Tensor[metric_,{},{}],bundlename_] := (
  Unprotect[metric];
  metric/: TensorMetricQ[metric] = True;
  (** Define the rule
              j                j
      metric    := Kronecker
            i               i      **)
  metric/: HoldPattern[Tensor[metric,{L[i_],U[j_]},{}]]
           := Tensor[Kronecker,{L[i],U[j]},{}];
  metric/: HoldPattern[Tensor[metric,{U[i_],L[j_]},{}]]
           := Tensor[Kronecker,{U[i],L[j]},{}];
  (** Always assume the metric is parallel. **)
  metric/: HoldPattern[Tensor[metric,{i_,j_},{__}]] := 0;
  (** If the frame is orthonormal, derivatives of the metric vanish **)
  If[ OrthonormalFrame[bundlename],
      metric/: HoldPattern[Tensor[metric,{_,_},{__}]] := 0;
      (** For 1-dimensional bundles with orthonormal frames, 
          the metric coefficients are all 1 **)
      If[ Dimension[bundlename] === 1,
          metric/: HoldPattern[Tensor[metric,{_,_},{}]] := 1;
        ]
    ];
  Protect[metric];
  );

defineBasisNames[bname_,cobname_,bundle_,quiet_] := (
  Unprotect[{Basis,bname}];
  bname/: HoldPattern[bname[L[i_]]] := 
            Tensor[Basis,{L[i]},{}];
  Basis/: HoldPattern[TensorFormat[Basis,{L[i_]},{}]] := 
            TensorFormat[bname,{L[i]},{}] /;
            MemberQ[{bundle,Conjugate[bundle]},Bundle[i]];
  Basis/: HoldPattern[TensorTeXFormat[Basis,{L[i_]},{}]] := 
            TensorTeXFormat[bname,{L[i]},{}] /;
            MemberQ[{bundle,Conjugate[bundle]},Bundle[i]];
  If[ cobname === Automatic,
      (*then*) bname/: HoldPattern[bname[U[i_]]] := 
                         Tensor[Basis,{U[i]},{}];
               Basis/: HoldPattern[TensorFormat[Basis,{U[i_]},{}]] := 
                         TensorFormat[bname,{U[i]},{}] /;
                         MemberQ[{bundle,Conjugate[bundle]},Bundle[i]];
               Basis/: HoldPattern[TensorTeXFormat[Basis,{U[i_]},{}]] := 
                         TensorTeXFormat[bname,{U[i]},{}] /;
                         MemberQ[{bundle,Conjugate[bundle]},Bundle[i]],
      (*else*) Unprotect[cobname];
               cobname/: HoldPattern[cobname[U[i_]]] :=
                         Tensor[Basis,{U[i]},{}];
               Basis/: HoldPattern[TensorFormat[Basis,{U[i_]},{}]] := 
                         TensorFormat[cobname,{U[i]},{}] /;
                         MemberQ[{bundle,Conjugate[bundle]},Bundle[i]];
               Basis/: HoldPattern[TensorTeXFormat[Basis,{U[i_]},{}]] := 
                         TensorTeXFormat[cobname,{U[i]},{}] /;
                         MemberQ[{bundle,Conjugate[bundle]},Bundle[i]];
               Protect[cobname];
    ];
  Protect[{Basis,bname}];
  );

(******************* DeclareBundle ****************************)

Options[DeclareBundle] := {FlatConnection   -> Null,
                           ParallelFrame    -> Null,
                           OrthonormalFrame -> Null,
                           PositiveDefinite -> Null,
                           CommutingFrame   -> Null,
                           TorsionFree      -> Null};

DeclareBundle[ b_, opts___ ] := 
  Module[ {flat,tfree,pframe,oframe,cframe,pdef},
  {flat,tfree,pframe,oframe,cframe,pdef} =
    {FlatConnection,TorsionFree,ParallelFrame,
     OrthonormalFrame,CommutingFrame,PositiveDefinite} /.
    {opts} /. Options[DeclareBundle];
  Which[
   flat =!= Null && !tfQ[flat],
     Message[DefineBundle::invalid,flat,"value for FlatConnection"],
   tfree =!= Null && !tfQ[tfree],
     Message[DefineBundle::invalid,tfree,"value for TorsionFree"],
   pframe =!= Null && !tfQ[pframe],
     Message[DefineBundle::invalid,pframe,"value for ParallelFrame"],
   oframe =!= Null && !tfQ[oframe],
     Message[DefineBundle::invalid,oframe,"value for OrthonormalFrame"],
   cframe =!= Null && !tfQ[cframe],
     Message[DefineBundle::invalid,cframe,"value for CommutingFrame"],
   pdef =!= Null && !tfQ[pdef],
     Message[DefineBundle::invalid,pdef,"value for PositiveDefinite"],
   (* Now all the arguments are OK.  Let's change the options. *)
   True,
     Unprotect[b];
     If[ tfree  =!= Null,  b/: TorsionFree[b] = tfree ];
     If[ pframe ||
         flat   =!= Null,  b/: FlatConnection[b] = flat ];
     If[ pframe =!= Null,  b/: ParallelFrame[b] = pframe ];
     If[ oframe =!= Null,  b/: OrthonormalFrame[b] = oframe ];
     If[ cframe =!= Null,  b/: CommutingFrame[b] = cframe ];
     If[ pdef   =!= Null,  b/: PositiveDefinite[b] = pdef ];
     Protect[b];
     Return[b]
  ] ];
  
(******************* UndefineBundle ****************************)

Attributes[UndefineBundle] = {Listable};

Options[UndefineBundle] := {Quiet -> $Quiet};

UndefineBundle[bundle_,opts___] := Module[{quiet},
  If[!optionsOK[UndefineBundle,opts],Return[]];
  quiet = Quiet /. {opts} /. Options[UndefineBundle];
  If[!BundleQ[bundle],
     Message[UndefineBundle::invalid,bundle,"bundle name"];
     Return[]
    ];
  UndefineIndex[BundleIndices[bundle],Quiet->quiet];
  UndefineTensor[Metric[bundle],Quiet->quiet];
  If[ MetricType[bundle] === Riemannian,
      UndefineTensor[RiemannTensor[bundle]];
      UndefineTensor[RicciTensor[bundle]];
      UndefineTensor[ScalarCurv[bundle]] ];
  If[ !quiet , Print["Undefining bundle ",bundle]];
  Unprotect[bundle];
  ClearAll[bundle];
  ];

(*************************** Constant *****************************)

(********************* ConstantQ ***********************************)

(* This function simply checks whether there are any tensors
   in the expression.  The reason for not using the Constant
   attribute to distinguish constants is that a symbol that
   is a constant in the "space" variables (i.e. with respect to
   covariant differentiation) may be a parameter which should
   act as a variable when applying D or Dt. *)

ConstantQ[x_] := FreeQ[x,Tensor[_,_,_]];

(******************** DefineConstant *******************************)

SetAttributes[DefineConstant,{Listable}];

Options[DefineConstant] := {Type -> Real,
                            Quiet -> $Quiet};

contypeQ[t_List] := 
  (t==={Complex}) ||
  (t==={Imaginary}) ||
  (ContainsQ[ {Real,Positive,Negative,NonPositive,
               NonNegative,Integer,Even,Odd}, t ] &&
   !ContainsQ[ t, {Positive,Negative} ] &&
   !ContainsQ[ t, {Positive,NonPositive} ] &&
   !ContainsQ[ t, {Negative,NonNegative} ] &&
   !ContainsQ[ t, {Odd,Even} ]);
contypeQ[t_] := contypeQ[{t}];

DefineConstant[c_, opts___] :=
  Module[{type,quiet},
    If[!optionsOK[DefineConstant,opts],Return[]];
    {type,quiet} = {Type,Quiet} /. {opts} /. Options[DefineConstant];
    Which[
      !newsymQ[c],   
        Message[DefineConstant::invalid,HoldForm[c],"constant name"],
      !contypeQ[type],
        Message[DefineConstant::invalid,type,"constant type"],
      True,
        DeclareConstant[c,Type->type];
        c/: c[i:(_L|_U)...] := InsertIndices[c,{i}];
        If[ !quiet, 
            Print["Constant ",c," defined.  Conjugate[",c,"] = ",
                   Conjugate[c],"."]
          ];
    ];
    Return[c];
  ];

(********************* UndefineConstant *******************************)

SetAttributes[UndefineConstant, {Listable,HoldFirst}];

Options[UndefineConstant] := {Quiet -> $Quiet};

UndefineConstant[c_,opts___] := Module[{quiet},
  If[!optionsOK[UndefineConstant,opts],Return[]];
  quiet = Quiet /. {opts} /. Options[UndefineBundle];
  Which[
    !optionsOK[UndefineConstant,opts],  
               Return[],
    !(Head[c]===Symbol && ConstantQ[c]),
               Message[UndefineConstant::invalid,HoldForm[c],"constant"],
    True,      Clear[c];
               If[!quiet,Print["Constant ",c," undefined."]]
    ]
  ];

(********************* Declare, DeclareConstant *************************)

Attributes[Declare] = {Listable};

Options[Declare] := Union[ Options[DeclareTensor],
                           Options[DeclareBundle],
                           Options[DeclareMathFunction],
                           Options[DeclareIndex],
                           Options[DeclareConstant] ];

Declare[x_, opts___] := Module[{sub},
  Switch[x,
    Tensor[t_,{},{}],       sub = DeclareTensor,
    b_Symbol?BundleQ,       sub = DeclareBundle,
    i_Symbol?IndexQ,        sub = DeclareIndex,
    f_Symbol?mathFunctionQ, sub = DeclareMathFunction,
    c_Symbol?ConstantQ,     sub = DeclareConstant,
    _,                      Message[Declare::invalid,x,
                              "tensor, bundle, index, math function, or constant"];
                            Return[];
    ];
  If[ !optionsOK[sub,opts], Return[]];
  sub[x,opts]
  ];

Options[DeclareConstant] := {Type -> Null};

DeclareConstant[ c_, opts___ ] := Module[{type},
  If[MemberQ[optionNames[{opts}], Type],
    type = Flatten[{Type /. {opts}}];
    If[ !contypeQ[type], 
        (*then*) Message[Declare::invalid, type, "constant type"];
                 Return[]
      ];
    If[ ValueQ[Sign[c]],        c/: Sign[c]        =.];
    If[ ValueQ[Positive[c]],    c/: Positive[c]    =.];
    If[ ValueQ[NonNegative[c]], c/: NonNegative[c] =.];
    If[ ValueQ[NonPositive[c]], c/: NonPositive[c] =.];
    If[ ValueQ[Negative[c]],    c/: Negative[c]    =.];
    If[ ValueQ[Conjugate[c]],   c/: Conjugate[c]   =.];
    If[ EvenQ[c],               c/: EvenQ[c]       =.];
    If[ OddQ[c],                c/: OddQ[c]        =.];
    If[ IntegerQ[c],            c/: IntegerQ[c]    =.];
    If[ MemberQ[type, Real|Positive|Negative|
                      NonNegative|NonPositive|
                      Integer|Even|Odd],       c/: Conjugate[c]   = c ];
    If[ MemberQ[type, Imaginary],              c/: Conjugate[c]   = -c ];
    If[ MemberQ[type, Positive],               c/: Sign[c]        = 1;
                                               c/: Positive[c]    = True;
                                               c/: NonNegative[c] = True ];
    If[ MemberQ[type, NonNegative],            c/: NonNegative[c] = True;
                                               c/: Negative[c]    = False ];
    If[ MemberQ[type, NonPositive],            c/: Positive[c]    = False;
                                               c/: NonPositive[c] = True ];
    If[ MemberQ[type, Negative],               c/: Sign[c]        = -1;
                                               c/: Positive[c]    = False;
                                               c/: NonNegative[c] = False ];
    If[ MemberQ[type, Integer|Even|Odd],       c/: IntegerQ[c]    = True ];
    If[ MemberQ[type, Even],                   c/: EvenQ[c] = True ];
    If[ MemberQ[type, Odd],                    c/: OddQ[c]  = True ];
  ];
  Return[c] ];


(*************************** DefineRelation *****************************)

(************ DefineRelation, UndefineRelation, & DefineRule ****************)

(* These commands are both implemented by the internal subroutine defineR.
   The exact function to perform is indicated by the argument to defineR:

   rhs = Null:       UndefineRelation
   rulename = Null:  DefineRelation
   rulename != Null: DefineRule
*)

SetAttributes[{DefineRelation,UndefineRelation,DefineRule},HoldAll];

Options[DefineRelation] := {Quiet -> $Quiet};

Options[UndefineRelation] := {Quiet -> $Quiet};

Options[DefineRule] := {Quiet   -> $Quiet,
                        NewRule -> False};

DefineRelation[]    := Message[DefineRelation::argm,"DefineRelation",0,2];
DefineRelation[_]   := Message[DefineRelation::argmu,"DefineRelation",2];
DefineRule[]        := Message[DefineRule::argm, "DefineRule",0,3];
DefineRule[_]       := Message[DefineRule::argmu,"DefineRule",3];
DefineRule[_,_]     := Message[DefineRule::argm,"DefineRule",2,3];
UndefineRelation[]  := Message[UndefineRelation::argm,"UndefineRelation",0,1];

UndefineRelation[t_,opts___] := 
  defineR[Hold[t],Hold[Null],Hold[Null],Hold[Null],UndefineRelation,opts];

DefineRelation[lhs_,rhs_,opts___Rule] := 
  defineR[Hold[lhs],Hold[rhs],Hold[Null],Hold[Null],DefineRelation,opts];

DefineRelation[lhs_, rhs_, cond_, opts___] := 
  defineR[Hold[lhs],Hold[rhs],Hold[cond],Hold[Null],DefineRelation,opts];

DefineRule[rulename_,lhs_,rhs_,opts___Rule] := 
  defineR[Hold[lhs],Hold[rhs],Hold[Null],Hold[rulename],DefineRule,opts];

DefineRule[rulename_, lhs_, rhs_, cond_, opts___] := 
  defineR[Hold[lhs],Hold[rhs],Hold[cond],Hold[rulename],DefineRule,opts];

defineR[Hold[lhs_],Hold[rhs_],Hold[cond_],Hold[rulename_],
        fcnname_,opts___] := Module[{newrule, quiet, holdlhs,
                                    singletensorQ, lhsdums, lhsfrees,
                                    lhsones, lhsonesloc, newlhsones,
                                    newlhsonenames, newlhsonepatterns,
                                    lhsdumnames, pairlhsdumnames,
                                    lhsdumrules, lhsfreerules,
                                    rhsfreerules, lhsindexlist,
                                    bundles, lhsnames, bunpairs,
                                    holdcond, complexnames},

  If[!optionsOK[fcnname,opts], Return[]];

  {newrule,quiet} = {NewRule,Quiet} /. {opts} /. Options[fcnname];

  (* Check that rulename is valid--either a new symbol or a list of rules *)

  If[!MatchQ[Hold[rulename],Hold[_Symbol]] ||
     (ValueQ[rulename] && !MatchQ[rulename,_List]),
     Message[Ricci::invalid, rulename, "rule name"];
     Return[]];

  (* Convert input tensors in lhs to internal form; if any are in
     internal form to begin with, first put them back in input form to
     avoid infinite loops. *)

  holdlhs = Hold[lhs] //. HoldPattern[Tensor[t_,{i___},{j___}]] :> t[i][j]
                      /. {t_Symbol?TensorQ :> Tensor[t,{},{}]}
                      //. HoldPattern[Inverse[x_][i__]] :> 
                            Tensor[Inverse[x],{i},{}]
                      //. HoldPattern[Conjugate[Tensor[Inverse[x_],{},{}]]] :>
                            Evaluate[Tensor[Inverse[Conjugate[x]],{},{}]]
                      //. HoldPattern[Conjugate[Tensor[t_[Bar],{},{}]]] :>
                            Tensor[t,{},{}]
                      //. HoldPattern[Conjugate[Tensor[t_,i_,j_]]] :>
                            Tensor[t[Bar],{},{}]
                      //.{HoldPattern[Tensor[t_,{i___},{j__}][(k___)?IndexQ]] :> 
                            Tensor[t,{i},{j,k}],
                          (Tensor[t_,{i___},{}])?slotsfullQ [(j___)?IndexQ] :> 
                            Tensor[t,{i},{j}],
                          HoldPattern[Tensor[t_,{i___},{}][(j___)?IndexQ]] :> 
                            Tensor[t,{i,j},{}]}
                      /. barrule;

  (* Require that lhs is a single tensor (with or without indices)
     for DefineRelation & UndefineRelation *)

  singletensorQ = MatchQ[holdlhs,(Hold[Tensor[_,_,_]])];

  If[ (fcnname === DefineRelation || fcnname === UndefineRelation) &&
    !singletensorQ,
       (*then*) Print["ERROR: first argument to ",fcnname,
                      " must be a single tensor"];
                Return[]];

  (* OK, arguments look fine.  Let's proceed. *)

  (* NewRule -> True means erase previously existing rules by the same name *)

  If[newrule, rulename = .];

  (* BUILD LEFT-HAND SIDE.  Begin transformations of indices in lhs.
     We separate lhs indices into three groups: dummy indices, free
     indices, and one-dimensional indices.  First get separate lists
     of the three categories. *)

  lhsdums = GetDummyIndices[ holdlhs ];
  lhsfrees = GetFreeIndices[ holdlhs ];
  lhsones = Select[ lhsfrees, Dimension[Bundle[#]] === 1 & ];
  lhsfrees = Complement[ lhsfrees, lhsones ];

  (* Deal with the one-dimensional indices first.  Find exactly where
     they occur. *)

  lhsonesloc = Position[ holdlhs, Alternatives @@ lhsones ];
  lhsones = PartAt[holdlhs,#]& /@ lhsonesloc;

  (* Give new names to the one-dimensional indices, since when the
     relation is invoked, different occurrences of the index may
     appear at different altitudes. *)

  newlhsones = newdummyname /@ lhsones;
  newlhsonenames = First /@ newlhsones /. x_[Bar] :> x; 
  newlhsonepatterns = Pattern[#,Blank[]]& /@ newlhsonenames;

  (* Replace all the one-dimensional indices on lhs by their new names. *)

  holdlhs = Fold[ ReplacePart[ #1, #2[[1]], #2[[2]] ]&,
                  holdlhs,
                  Transpose[{newlhsonepatterns,lhsonesloc}] ];

  (* Next deal with dummy pairs on lhs.  We must create a rule that
     will recognize either "L[i],U[i]" or "U[i],L[i]".  This is done
     by using a pattern of the form "i_,pairi_" on the lhs, and then
     checking that pairi===Pair[i] before applying the rule.  The new
     list "lhsdums" will be of the form {i, j, ... } *)

  lhsdumnames = First /@ lhsdums /. x_[Bar] :> x;

  (* For lhsdumnames = {i, j, ... } created above, create a new list of
     names of the form newdums = {pairi, pairj, ... }. *)

  pairlhsdumnames = newpairname /@ lhsdumnames;

  (* dumpairs = { {i,pairi}, {j,pairj}, ... }. *)

  dumpairs = Transpose[{lhsdumnames,pairlhsdumnames}];

  (* "lhsdumrules" is a list of rules that change L[i], U[i] to i_,
     pairi_.  The condition built by buildcond will make sure that
     these indices form a matched pair before applying the rule. *)

  lhsdumrules = Flatten[ Replace[ #, {a_,b_} :> 
                            {(L[a]|L[a[Bar]]) :> Pattern[#,Blank[]]& [a],
                             (U[a]|U[a[Bar]]) :> Pattern[#,Blank[]]& [b]} ]& /@
                         dumpairs ];

  (* Now apply the rules to lhs and change dummies to patterns. *)

  holdlhs = holdlhs /. lhsdumrules;

  (* Finally deal with free indices on lhs.  "lhsfreerules" is a list
     of rules that will convert free lhs indices of the form L[i],
     U[i], L[i[Bar]], or U[i[Bar]] to i_.  This is how we get the rule
     to be applied regardless of whether the free indices are up or
     down. *)

  lhsfreerules = Cases[ lhsfrees /. { 
    x:((L[i_Symbol])|(U[i_Symbol])|(L[i_Symbol[Bar]])|(U[i_Symbol[Bar]])) :> 
      (x :> Pattern[#,Blank[]]& [i])},
    __RuleDelayed];

  (* Apply the rules. *)

  holdlhs = holdlhs /. lhsfreerules;

  (* Finally, if the left-hand side is a single tensor with or without
     indices, insert dummy variables ind1___ and ind2___ into the
     left-hand side as place holders for additional indices that might
     be inserted by the user. *)

  holdlhs = Replace[holdlhs, Hold[Tensor[t_,{i___},{j___}]] :> 
                             Hold[Tensor[t,{i,ind1___},{j,ind2___}]]];

  (* BUILD THE RHS.  First we replace all dummy indices on rhs by new
     dummy names, to prevent conflicts with index names on the LHS or
     those used when the rule is invoked; we also hold the RHS to
     prevent evaluation, and convert LB[_] and UB[_] to L[_[Bar]] and
     U[_[Bar]]. *)

  holdrhs = NewDummy[Hold[rhs] /. barrule];

  (* Now build a set of rules for transforming the rhs free indices. *)

  rhsfreerules = Cases[ lhsfrees /. { 
    x:((L[i_Symbol])|(U[i_Symbol])|(L[i_Symbol[Bar]])|(U[i_Symbol[Bar]])) :> 
      (x :> i)},
    __RuleDelayed ];

  (* Apply the rules. *)

  holdrhs = holdrhs /. rhsfreerules;

  (* If there were one-dimensional indices on the lhs, insert metric
     factors into rhs to compensate if the user puts them at different
     altitudes. *)

  holdrhs = Fold[ insertmetricfactor, 
                  holdrhs, 
                  Transpose[{lhsones,newlhsonenames}] ] ;

  (* Lastly, we enclose rhs in a "Module" so that dummy index names
     will be newly created every time the relation is invoked; this is
     to avoid conflict with dummy names such as "c2" in the user's
     call.  Begin by getting a list of dummy names to pass to
     "buildmodule". *)

  allindices = GetIndices[holdrhs];
  rhsdumnames = Union [ First /@ Intersection[
            (Pair /@ Cases[allindices,U[i_]/;Dimension[Bundle[i]]=!=1]), 
                     Cases[allindices,L[i_]/;Dimension[Bundle[i]]=!=1]
          ] /. x_[Bar] :> x ];
  holdrhs = buildmodule[holdrhs,rhsdumnames];

  (* BUILD THE CONDITION.  "bundles" is just a list of the bundle
     names associated with all the indices in lhs.  This will be Null
     if there's no bundle.  Similarly, "lhsnames" is the list of names
     of all indices that appear in lhs.  *)

  lhsindexlist = Flatten[Join[lhsfrees,lhsdums,newlhsones]];
  bundles = Bundle /@ Lower /@ lhsindexlist;
  lhsnames = (First /@ lhsindexlist) /. x_[Bar] :> x;

  (* bunpairs = { {i, buni}, {j, bunj}, ... }.  (Discard any entries
     for which the bundle name is Null.)  This will be used to
     construct the condition for applying the rule: the rule will get
     applied only if Bundle[i] = buni, etc.  (This has to be modified
     somewhat for complex indices.  See "buildcond" below.) *)

  bunpairs = Union[Cases[Transpose[{lhsnames,bundles}],
                         {_Symbol,_?(#=!=Null&)}]];

  (* Now build the condition for applying the rule. *)  

  holdcond = buildcond[ Hold[cond], bunpairs, dumpairs ];
  holdcond = holdcond /. rhsfreerules;

  (* DEFINE THE RELATION.  "defOne" defines one rule or relation, with
     appropriate side conditions. *)

  defOne[ holdlhs, holdrhs, holdcond, Hold[rulename]];

  (* See if we have to define the conjugate relation.  We do this only
     if (a) lhs is a single tensor, and (b) either the tensor name or
     one or more indices are complex (cxQ[holdlhs] == True). *)

  (* "complexnames" is a list of lhs index pattern names that are
     associated with complex bundles.  These
     are candidates for conjugation when the conjugate rule is
     generated. *)

  complexnames = Select[ Join[lhsnames, pairlhsdumnames],
                           Type[Bundle[#]]===Complex&
                       ];

  If[MatchQ[holdlhs,Hold[Tensor[_,_,_]]] && cxQ[holdlhs],
   (*then*) defOne[holdlhs /. conjTensorRule /. conjIndexRule,
                   holdrhs /. (conjrule /@ complexnames) /.
                              Hold[x_] :> Hold[Conjugate[x]] /; 
                                !(Hold[x]===Hold[Null]),
                   holdcond /. conjTensorRule /. conjIndexRule /. 
                               conjBundleRule,
                   Hold[rulename]];
    ];

  (* Tell the user what we've done. *)

  If[!quiet, Which[
    Hold[rhs]===Hold[Null],  
      Print["Relation undefined."],
    Hold[rulename]===Hold[Null],
      Print["Relation defined."],
    True,
      Print["Rule defined."]
    ]];
  ];

(************** buildcond ***********************************************)

(* This builds a condition for deciding whether to apply the rule.
   This consists of three parts: 
   (1) the user's condition, if specified; 
   (2) the index/bundle pairs, to make sure the user's indices are associated
       with the correct bundles;
   (3) the dummy pairs {i,pairi}, to make sure they form a paired set. 

   Checking the index/bundle pairs is complicated by the presence of complex
   indices.  Basically, if a is an index associated with the complex bundle E, 
   and the user specifies L[a] in lhs, we want the rule to be applied whenever 
   the supplied index is a lower E index or an upper Conjugate[E] index.  The 
   way to do this is to check that the LOWERED version of the index is associated
   with E.
*)

buildcond[Hold[Null], {}, {}] := Hold[Null];

buildcond[Hold[Null], {}, {dumpairs__}] :=
  Hold[And @@ (PairQ @@ # &) /@ {dumpairs}];

buildcond[Hold[Null], {bunpairs__}, {}] := 
  Hold[And @@ (Bundle @ Lower[ #[[1]] ] === #[[2]] &) /@ {bunpairs}];

buildcond[Hold[Null], {bunpairs__}, {dumpairs__}] :=
  Hold[And @@ (Bundle @ Lower[ #[[1]] ] === #[[2]] &) /@ {bunpairs} &&
       And @@ (PairQ @@ # &) /@ {dumpairs}];

buildcond[Hold[cond_], {}, {}] := Hold[cond] /. barrule;

buildcond[Hold[cond_], {}, {dumpairs__}] :=
  Hold[cond && 
       And @@ (PairQ @@ # &) /@ {dumpairs}] /. barrule;

buildcond[Hold[cond_], {bunpairs__}, {}] := 
  Hold[cond && 
       And @@ (Bundle @ Lower[ #[[1]] ] === #[[2]] &) /@ {bunpairs}] /. barrule;

buildcond[Hold[cond_], {bunpairs__}, {dumpairs__}] :=
  Hold[cond && 
       And @@ (Bundle @ Lower[ #[[1]] ] === #[[2]] &) /@ {bunpairs} &&
       And @@ (PairQ @@ # &) /@ {dumpairs}] /. barrule;

(*********************** defOne **************************************)

(* For DefineRelation and UndefineRelation, the case of a barred tensor
   on the lhs has to be handled separately, since we have to extract
   the tensor name as a tag.
*)

(* If rhs is Null, then undefine relation. *)

defOne[Hold[Tensor[t_[Bar], i_, j_]],Hold[Null],_,_] := (
  Unprotect[t];
  t/: Tensor[t[Bar], i, j]=.;
  Protect[t]
  );

defOne[Hold[Tensor[t_, i_, j_]],Hold[Null],_,_] := (
  Unprotect[t];
  t/: Tensor[t, i, j]=.;
  Protect[t]
  );

(* If rulename is Null, then define relation. *)

(* Case 1: no condition. *)

defOne[Hold[Tensor[t_[Bar], i_, j_]], Hold[rhs_], Hold[Null], 
    Hold[Null] ] := (
  Unprotect[t];
  t/: Tensor[t[Bar], i, j] := evalRelation[rhs,{ind1},{ind2}];
  Protect[t]
  );

defOne[Hold[Tensor[t_, i_, j_]], Hold[rhs_], Hold[Null], 
    Hold[Null] ] := (
  Unprotect[t];
  t/: Tensor[t, i, j] := evalRelation[rhs,{ind1},{ind2}];
  Protect[t]
  );

(* Case 2: condition to check. *)

defOne[Hold[Tensor[t_, i_, j_]], Hold[rhs_], Hold[cond_], Hold[Null] ] := (
  Unprotect[t];
  t/: Tensor[t, i, j] := evalRelation[rhs,{ind1},{ind2}] /; cond;
  Protect[t]
  );

defOne[Hold[Tensor[t_[Bar], i_, j_]], Hold[rhs_], Hold[cond_], Hold[Null]
    ] := (
  Unprotect[t];
  t/: Tensor[t[Bar], i, j] := evalRelation[rhs,{ind1},{ind2}] /; cond;
  Protect[t]
  );

(* If rulename is not Null, then defining a rule. *)

(* Case A: lhs is a single tensor.  Same 2 cases as for DefineRelation. *)

defOne[Hold[Tensor[t_,i_,j_]], Hold[rhs_], Hold[Null], Hold[rulename_] ] := 
  rulename = AppendRule[rulename, 
             HoldPattern[Tensor[t, i, j]] :> evalRelation[rhs,{ind1},{ind2}]];

defOne[Hold[Tensor[t_,i_,j_]], Hold[rhs_], Hold[cond_], Hold[rulename_] ] := 
  rulename = AppendRule[rulename,
    HoldPattern[Tensor[t, i, j]] :> evalRelation[rhs,{ind1},{ind2}] /; cond];

(* Case B: general left-hand side. Again two cases. *)

defOne[Hold[lhs_], Hold[rhs_], Hold[Null], Hold[rulename_] ] := 
  rulename = AppendRule[rulename, HoldPattern[lhs] :> NewDummy[rhs]];

defOne[Hold[lhs_], Hold[rhs_], Hold[cond_], Hold[rulename_] ] := 
  rulename = AppendRule[rulename, HoldPattern[lhs] :> NewDummy[rhs] /; cond];

defOne[___] := Print["ERROR: Internal error in DefineRule."];

(*************** evalRelation ******************************************)

(* Called when the right-hand side is evaluated, to decide whether 
   indices need to be inserted or not. *)

evalRelation[rhs_,{i___},{j___}] := 
  If[{i,j}==={}, 
        (*then*)  NewDummy[rhs],
        (*else*)  CovD[InsertIndices[NewDummy[rhs],{i}],{j}]
    ];

(****************** Misc. subroutines for defineR **********************)

newdummyname[ L[i_[Bar]] ] := L[ Unique[ToString[i]] [Bar] ];
newdummyname[ U[i_[Bar]] ] := U[ Unique[ToString[i]] [Bar] ];
newdummyname[ L[i_] ]      := L[ Unique[ToString[i]] ];
newdummyname[ U[i_] ]      := U[ Unique[ToString[i]] ];

insertmetricfactor[ exp_, {oldindex_,newindex_} ] := 
  Replace[ exp, Hold[y_] :> 
    Hold[ y * Metric[Bundle[oldindex]] [Pair[oldindex], newindex] ] ];

buildmodule[ Hold[exp_], {} ] := Hold[exp];

buildmodule[ Hold[exp_], dummies_ ] :=
    Hold[ Module[ dummies, exp ] ];    

conjTensorRule = {
   Tensor[t_[Bar],{i___},{j___}] :> 
     Tensor[t, {i}, {j}];
   Tensor[t_,{i___},{j___}] :> 
     Tensor[t[Bar], {i}, {j}] /; (Type[t])==={Complex}};

conjIndexRule = {
   L[i_[Bar]] :> L[i],
   U[i_[Bar]] :> U[i],
   L[i_] :> L[i[Bar]] /; Type[Bundle[i]]===Complex,
   U[i_] :> U[i[Bar]] /; Type[Bundle[i]]===Complex};

conjBundleRule = {
   Conjugate[b_Symbol?BundleQ] :> b,
   b_Symbol?BundleQ      :> Conjugate[b]};

conjrule[i_] := (i -> Conjugate[i]);

pname[HoldPattern[x_Pattern[Bar]]] := x[[1]];
pname[HoldPattern[_[x_Pattern[Bar]]]] := x[[1]];
pname[HoldPattern[x_Pattern]] := x[[1]];
pname[HoldPattern[_[x_Pattern]]] := x[[1]];
pname[_] := Null;

cxQ[Hold[Tensor[t_,{i___},{j___}]]] := 
  Type[t]==={Complex} ||
  Or @@ (Type[#]===Complex&) /@ 
    Bundle /@ pname /@ {i,j};

AppendRule[{r___}, rule_] := {r,rule};
AppendRule[r_, rule_] := {rule};

barrule = {HoldPattern[LB[x_]] :> L[x[Bar]], HoldPattern[UB[x_]] :> U[x[Bar]]};

newpairname[x_] := Module[{name,spellon,spell1on},
  spellon  = !MatchQ[ General::spell,  _$Off ];
  spell1on = !MatchQ[ General::spell1, _$Off ];
  Off[General::spell,General::spell1];
  name = ToExpression[ "pair" <> ToString[x] ];
  If[ spellon,  On[General::spell]];
  If[ spell1on, On[General::spell1]];
  Return[name]
  ];

(*************************** Deriv *****************************)

(************ Del[exp] = total covariant derivative of exp *****************)

Options[Del] := {Metric -> Automatic,
                 Conn   -> Conn};

Del[_?ConstantQ, ___Rule] := 0;
Del[f_Plus, con___Rule] := (Del[#,con]&) /@ f;

(* Leibniz rule for function multiples only *)

Del[(f_ /; Rank[f]===0) g_, con___Rule] := 
  f Del[g,con] + TensorProduct[g, Extd[f]];

HoldPattern[Del[Tensor[g_,{},{}]]] := 0 /; MetricQ[g];

HoldPattern[Del[g_Tensor,Metric->g_]] := 0;

HoldPattern[Del[___, Tensor[Basis,{i_},{}]]] := 0 /; ParallelFrame[Bundle[i]];

HoldPattern[Del[Summation[f_],con___Rule]] := Del[f,con] ;

HoldPattern[Del[(f_ /; Rank[f]===0), con___Rule]] := Extd[f];

Del[f_[g_],con___Rule] := f'[NewDummy[g]] * Del[g,con] /;
  mathFunctionQ[f];

(*********** Del[v,exp] = covariant derivative in direction v *************)

HoldPattern[Del[f_ * v_, x_, con___Rule]] := f Del[v,x,con] /; Rank[f]===0;
HoldPattern[Del[v_ + w_, x_, con___Rule]] := Del[v,x,con] + Del[w,x,con];
Del[0, _, ___Rule] := 0;
Del[v_, x_Plus, con___Rule] := Del[v,#,con]& /@ x;
Del[v_, f_ * g_, con___Rule] := Del[v,f,con] g + f Del[v,g,con];

Del[v_, f_^p_, con___Rule] := 
  NewDummy[p] NewDummy[f]^(p-1) Del[v,f,con] + 
  Log[NewDummy[f]] f^NewDummy[p] Del[v,p,con];

HoldPattern[Del[v_, Wedge[x_,y__], con___Rule]] := 
  Wedge[Del[v,x,con],y] + Wedge[x,Del[v,Wedge[y],con]];

HoldPattern[Del[v_, TensorProduct[x_,y__], con___Rule]] := 
  TensorProduct[Del[v,x,con],y] + TensorProduct[x,Del[v,TensorProduct[y],con]];

HoldPattern[Del[v_, Sym[x_], con___Rule]] := Sym[Del[v,x,con]];
HoldPattern[Del[v_, Alt[x_], con___Rule]] := Alt[Del[v,x,con]];

HoldPattern[Del[_, Tensor[g_,{},{}]]] := 0 /; MetricQ[g];

HoldPattern[Del[_, Tensor[Kronecker, {_,_}, {}], ___]] := 0;

HoldPattern[Del[_, 
            Tensor[(g_ /; MetricQ[g] && 
                          OrthonormalFrame[ Bundles[g][[1,1]] ]),
                   {_,_},{}]]] := 0;

HoldPattern[Del[v_, Tensor[g_?MetricQ, {U[i_],U[j_]}, {}]]] := 
  - (g[U[i],U[#1]] g[U[j],U[#2]] Del[v, g[L[#1],L[#2]]]&) [ 
    NewBundleSymbol[Conjugate[Bundle[i]]], 
    NewBundleSymbol[Conjugate[Bundle[j]]] ];

HoldPattern[Del[v_, Tensor[Inverse[x_],{i_,j_},{}]]] :=
  - Plus @@ Flatten [ Outer @@ Join[
      {InsertIndices[NewDummy[Inverse[x]], {i,Pair[#1]}] *
       InsertIndices[NewDummy[Inverse[x]], {Pair[#2],j}] *
       Del[v, InsertIndices[x, {#1,#2}] ] & },
      NewIndex /@ Transpose[{Bundles[x],Variance[x]}] ]];

HoldPattern[Del[v_, Summation[f_],con___Rule]] := Del[v,f,con] ;

(* If basis elements commute, derivatives of functions commute too. *)

HoldPattern[ Del[ Tensor[Basis,{L[i_]},{}],
              Del[ Tensor[Basis,{L[j_]},{}], 
                   (f_ /; Rank[f]===0) ] ] ] :=
  Del[ Tensor[Basis,{L[j]},{}], Del[ Tensor[Basis,{L[i]},{}], f ]] /;
  Bundle[i]===Bundle[j] && CommutingFrame[Bundle[i]] &&
  !IndexOrderedQ[{L[i],L[j]}];

(* For the contraction operators Inner, Int, & Dot, Del[v,_] satisfies
   the product rule unless a non-default connection is specified.*)

HoldPattern[Del[v_, Inner[x_,y_]]] := 
  Inner[Del[v,x],y] + Inner[x,Del[v,y]];

HoldPattern[Del[v_, Int[x_,y_]]] := 
  Int[Del[v,x],y] + Int[x,Del[v,y]];
HoldPattern[Del[v_, Dot[x_,y__]]] :=
  Dot[ Del[v,x], y] + Dot[ x, Del[v, Dot[y] ] ];

HoldPattern[Del[v_, Inverse[x_], con___Rule ]] :=
  - Dot[ Dot[ NewDummy[Inverse[x]], 
              Del[v,x,con]], 
         NewDummy[Inverse[x]] ];

HoldPattern[Del[v_, Det[x_], con___Rule ]] :=
  Det[x] Tr[ Inverse[x] . Del[v,x,con] ] /; rank[x]===2;

Del[v_, c_?ConstantQ, con___Rule] := 0 /; Head[c] =!= Rule;

Del[v_, f_[g_], con___Rule] := f'[NewDummy[g]] * Del[v,g,con] /;
  mathFunctionQ[f];

Del/: Conjugate[HoldPattern[Del[x_, y___, con_Rule]]] := 
      Del[Conjugate[x],Conjugate[y],con];
Del/: Conjugate[HoldPattern[Del[x__]]] := Conjugate /@ Del[x];

(******** Insertion of indices **********)

(* Here's where the actual computation of Del is done. *)

(* First handle a special case: derivatives of Conn are left
   in the form Del[ Basis[_], Conn[_,_,_] ] to avoid "covariant
   derivatives" of Conn. *)

Del/: TCompute[ HoldPattern[ Del[ Tensor[Conn,{i_,j_,k_},{}] ] ], 
                {l_} ] :=
  Del[ Basis[l], Tensor[Conn,{i,j,k},{}] ];

(* Now the general case.  The main complication here arises when the
   tensor being differentiated has free indices, for then its total
   covariant derivative involves some connection coefficients.
   Because of this, any other function that computes covariant
   derivatives of arbitrary tensor expressions must call this one.  
   Also, this is where we compute components of covariant derivatives
   with respect to non-default connections.  *)

Del/: TCompute[ HoldPattern[Del[x_, con___Rule]], {j___,i_}] := 
  Module[{n,diff,freeindices,connection},
  freeindices = GetFreeIndices[x];
  Switch[ {con},
    {},                connection = Conn,
    {Connection -> _}, connection = con[[2]],
    {Metric -> _},     connection = LeviCivitaConnection[ con[[2]] ],
    _,                 Message[Del::invalid, con, "Del option"];
                       Return[ERROR[Del[x,con][j,i]]] ];
  diff = connection - Conn;
  CovD[ InsertIndices[x,{j}], {i} ] +
  (* Add non-default connection terms associated to previous indices *)
  If[ diff===0,
      0,
      Sum[
        Plus @@ Flatten [ Outer @@ Join[
        {InsertIndices[ x, Join[Take[{j},n-1],{#},Drop[{j},n]] ] *
           If[ Variance[x][[n]]===Covariant,
               (*then covariant index*) 
                 - InsertIndices[ diff, { {j}[[n]], Pair[#], i } ],
               (*else contravariant*)
                   InsertIndices[ diff, { Pair[#], {j}[[n]], i } ]
             ] & },
        {NewIndex [ {Bundles[x][[n]], Variance[{j}[[n]]]} ]} ] ],
      {n,Length[{j}]}]
    ] +
  (* Finally, add connection terms associated to free indices. *)
  Sum[ 
    InsertIndices[ x /. freeindices[[n]] -> #, {j} ] *
    If[IndexAltitude[#]===Upper,
      (*then contravariant index*)
        - InsertIndices[Conn, { Pair[#], freeindices[[n]], i } ],
      (*else covariant*)
        + InsertIndices[Conn, { freeindices[[n]], Pair[#], i } ]
      ] & @ 
    dumdex[ freeindices[[n]] ],
    { n, Length[freeindices] }
    ]
  ];

(* Now Del[v,x].  First do the simpler case when v is a basis vector. *)

Del/: TCompute[HoldPattern[Del[Tensor[Basis,{k_},{}], x_, con___Rule]], {j___}] := 
  InsertIndices[ Del[x,con], {j,k} ];

(* Now the general case. *)

Del/: TCompute[ HoldPattern[ Del[v_, x_, con___Rule] ], {j___}] := 
  Module[{vindices,i},
  If[!(Rank[v]===1),
    Print["ERROR: ",v," is not a 1-tensor"];
    Return[ERROR[Del[v,x,con][j]]]
    ];
  vindices = NewLowerIndex /@ Flatten[Bundles[v]];
  Return[ Sum[ InsertIndices[ v, {Pair[ vindices[[i]] ]} ] * 
               InsertIndices[ Del[x,con], {j,vindices[[i]]} ],
               {i,Length[vindices]}]
        ]
  ];

(*** Grad ******************************************************)

Options[Grad] := {Metric -> Automatic,
                 Conn   -> Conn};

Grad[f_Plus, opt___] := (Grad[#,opt]&) /@ f;
Grad[_?ConstantQ, opts___] := 0;

(* Leibniz rule for function multiples only *)

Grad[(f_ /; Rank[f]===0) g_, opt___] := 
  f Grad[g,opt] + TensorProduct[g, Grad[f,opt]];

HoldPattern[Grad[f_ ^ p_, opt___]] :=
  NewDummy[p] NewDummy[f]^(p-1) Grad[f,opt] + 
  Log[NewDummy[f]] f^NewDummy[p] Grad[p,opt] /; 
  Rank[f]===0 && Rank[p]===0;

HoldPattern[Grad[Tensor[g_,{},{}]]] := 0 /; MetricQ[g];

HoldPattern[Grad[g_Tensor,Metric->g_]] := 0;

HoldPattern[Grad[Summation[f_],opt___]] := Grad[f,opt] ;

HoldPattern[Grad[Tensor[Basis,{i_},{}]]] := 0 /; ParallelFrame[Bundle[i]];

Grad[f_[g_],opt___] := f'[NewDummy[g]] * Grad[g,opt] /;
  mathFunctionQ[f];

Grad[f_?ConstantQ] := 0;

Grad/: Conjugate[HoldPattern[Grad[x_, y___, opt_Rule]]] := 
       Grad[Conjugate[x],Conjugate[y],opt];
Grad/: Conjugate[HoldPattern[Grad[x__]]] := Conjugate /@ Grad[x];

Grad/: TCompute[ HoldPattern[Grad[x_, opt___]], {j___,i_}] := 
  Module[{met,con},
    {met, con} = {Metric, Conn} /. {opt} /. Options[Grad];
    If[ met === Automatic,
        (*then*) Return[ Del[x,opt] [j,i] ],
        (*else*) Return[ (Del[x,opt] . Inverse[met]) [j,i] ]
      ]
    ];

(*** Div ******************************************************)

HoldPattern[Div[f_Plus,opts___Rule]] := (Div[#,opts]&) /@ f;
HoldPattern[Div[(c_?ConstantQ) f_, opts___Rule]] := c Div[f,opts];

HoldPattern[Div[ (f_ /; Rank[f]===0) * x_, opts___Rule ]] :=
  f Div[x,opts] + 
    If[ Last[Variance[x]] === Contravariant,
        (*then*) x . Extd[f],
        (*else*) x . Grad[f, opts] ];

HoldPattern[Div[Summation[x_],opts___Rule]] := Div[x,opts];

HoldPattern[Div[t_,opts___Rule]] := 0 /; Rank[t]===0;

Options[Div] := {Metric -> Automatic,
                 Connection -> Conn};

Div/: TCompute[ HoldPattern[Div[t_,opts___Rule]], {i___}] := 
  Module[{bun,met,con},
  bun = UnderlyingTangentBundle[t];
  Which[
    MemberQ[bun,Null],
      Message[Bundle::error,Div[t,opts]];
      Return[ERROR[Div[t,opts][i]]],
    !optionsOK[Div, opts],
      Return[ERROR[Div[t,opts][i]]],
    {opts} === {},
      Return[Plus @@ ( Del[t,opts] [i, #, Pair[#] ] &) /@ NewLowerIndex /@ bun],
    True,
      {met,con} = {Metric,Connection} /. {opts} /. Options[Div];
      If[ con === Conn,
          con = LeviCivitaConnection[met]];
      If[ Last[Variance[t]] === Contravariant,
          (*then*) Return[ Plus @@ 
                             ( (Del[t,Connection->con] [i, Pair[#], # ] &)
                               /@ NewLowerIndex /@ bun ) ],
          (*else*) Return[ Plus @@ ( InsertIndices[
                         If[ met===Automatic,
                             Metric[Bundle[#[[1]]]],
                             Inverse[ met ]],
                         {Pair[ #[[1]] ],Pair[ #[[2]] ]}] *
                       Del[t,Connection->con] [i, #[[1]], #[[2]] ] &)
                         /@ (Table[NewLowerIndex[#],{2}]&) /@ bun ]
          ]    
    ];
  ];

Div/:  Conjugate[HoldPattern[Div[x_,opts___Rule]]] := 
           Div[Conjugate[x],opts];

(****************** Laplacian *****************************************)

Laplacian[x_,opts___] := 
  Switch[ $LaplacianConvention,
          DivGrad,            Div[ Grad[ x, opts ], opts],
          PositiveSpectrum, - Div[ Grad[ x, opts ], opts],
          _,                Message[ $LaplacianConvention::invalid, 
                                     $LaplacianConvention, 
                                     "$LaplacianConvention" ];
                            ERROR[ Laplacian[x,opts] ]
        ];

(****************** LaplaceBeltrami **************************************)

LaplaceBeltrami[x_,opts___] := 
  Extd[ ExtdStar[ x, opts ] ] + ExtdStar[ Extd[x], opts];

(************************ Extd ************************************)

MakeLinear[Extd];

HoldPattern[Extd[_Extd]] = 0;

HoldPattern[Extd[f_ * g_]] :=
  Wedge[Extd[f], g] + f Extd[g] /; Rank[f] === 0;

HoldPattern[Extd[Wedge[f_,g__]]] := 
  Wedge[Extd[f],g] + (-1)^Rank[f] Wedge[f,Extd[Wedge[g]]];

HoldPattern[Extd[f_ ^ p_]] :=
  NewDummy[p] NewDummy[f]^(p-1) Extd[f] + 
  Log[NewDummy[f]] f^NewDummy[p] Extd[p] /; 
  Rank[f]===0 && Rank[p]===0;

HoldPattern[Extd[Summation[f_]]] := Extd[f] ;

HoldPattern[Extd[Tensor[Kronecker,{_,_},{}]]] := 0;

HoldPattern[Extd[Tensor[(g_ /; MetricQ[g] && 
                           OrthonormalFrame[ Bundles[g][[1,1]] ]),
                    {_,_},{}]]] := 0;

HoldPattern[Extd[Tensor[Inverse[x_],{i_,j_},{}]]] :=
  - Plus @@ Flatten [ Outer @@ Join[
      {InsertIndices[NewDummy[Inverse[x]], {i,Pair[#1]}] *
       InsertIndices[NewDummy[Inverse[x]], {Pair[#2],j}] *
       Extd[ InsertIndices[x, {#1,#2}] ] & },
      NewIndex /@ Transpose[{Bundles[x],Variance[x]}] ]];

HoldPattern[Extd[Tensor[g_?MetricQ,{U[i_],U[j_]},{}]]] :=
  - (g[U[i],U[#1]] g[U[j],U[#2]] Extd[g[L[#1],L[#2]]]&) [ 
    NewBundleSymbol[Conjugate[Bundle[i]]], 
    NewBundleSymbol[Conjugate[Bundle[j]]] ];

HoldPattern[Extd[f_?ConstantQ]] := 0;

HoldPattern[Extd[f_[g_]]] := f'[NewDummy[g]] * Extd[g] /;
  mathFunctionQ[f];

(** Here's where we insert indices. **)

Extd/: TCompute[ HoldPattern[Extd[f_]], {i__} ] := 
  Module[{rk,tanbun,dummyindices,m,n,k},
  rk = rank[f];
  If[!FormQ[f],
    Print["ERROR: ",f," is not a differential form"];
    Return[ERROR[Extd[f][i]]]
    ];
  If[rk===0,
    (*then*) 
      (* Here Del[f] must be kept Unevaluated, or it will be
         turned back into Extd[f] by the rules for Del. *)
      TCompute[ Unevaluated[Del[f]], {i} ],
    (*else*) 
      tanbun = UnderlyingTangentBundle[f];
      (* Here's the actual computation:
         First the numerical factor (depending on WedgeConvention *)
      Signature[{i}] * WedgeFactor[{1,rk}] / (rk+1) * (
      (* Next the antisymmetrized derivatives *)
      Sum[ 
        ((-1)^rk * Signature[#] * InsertIndices[Del[f],#] &) 
           @ RotateLeft[{i},n],
        {n,rk+1}] +
      (* Finally the torsion terms, if there are any free indices *)
      If[And @@ TorsionFree /@ tanbun,
        (*then*) 
          0,
        (*else*)
          dummyindices = NewLowerIndex /@ tanbun;
          Sum[
            ((1/2) * Signature[#] *
              InsertIndices[f, Join[ {dummyindices[[k]]}, Take[#,rk-1] ] ] * 
              InsertIndices[Tor, Join[ {Pair[dummyindices[[k]]]}, 
                                       Drop[#,rk-1] ] ]&
            )
              @ (Join[ RotateLeft[ Drop[#,-1], m ], Take[#,-1] ] & @
                RotateLeft[{i},n]),
            {k,Length[dummyindices]},
            {m,rk},
            {n,rk+1}
            ]
        ]
      ) 
    ]
  ];

Extd/: Conjugate[HoldPattern[Extd[x_]]] := Extd[Conjugate[x]];

(****************************** ExtdStar ***********************************)

HoldPattern[ExtdStar[f_Plus,opts___Rule]] := (ExtdStar[#,opts]&) /@ f;
HoldPattern[ExtdStar[(c_?ConstantQ) f_, opts___Rule]] := c ExtdStar[f,opts];

HoldPattern[ExtdStar[Summation[x_],opts___Rule]] := ExtdStar[x,opts];

HoldPattern[ExtdStar[t_,opts___Rule]] := 0 /; rank[t]===0;

HoldPattern[ExtdStar[f_*t_]] := f ExtdStar[t] - Int[Extd[f], t] /; rank[f]===0;
HoldPattern[ExtdStar[ExtdStar[f_,opts___Rule],opts___Rule]] := 0;

Options[ExtdStar] := {Metric -> Automatic};

ExtdStar/: TCompute[ HoldPattern[ExtdStar[f_,opt___Rule]], {i___} ] := 
  Module[{rk,tanbun,dummyindices,m,n,k,met,con},
  If[ !FormQ[f],
    Print["ERROR: ",f," is not a differential form"];
    Return[ERROR[ExtdStar[f][i]]]];
  rk = Rank[f];
  tanbun = UnderlyingTangentBundle[f];
  If[ MemberQ[bun,Null],
      Message[Bundle::error,ExtdStar[t,opt]];
      Return[ERROR[ExtdStar[t,opt][i]]]];
  If[ !optionsOK[ExtdStar, opt],
      Return[ERROR[ExtdStar[t,opt][i]]]];
  (* Here's the computation:
     First the numerical factor, depending on WedgeConvention *)
  (-1)^(rk-1) * WedgeFactor[{rk-1,1}] * 
   IntFactor[rk,rk] / IntFactor[rk-1,rk-1] *
   (* Next the derivative terms *)
   (InsertIndices[ - Div[f,opt], {i}] +
   (* Finally the torsion terms, if there are free indices *)
    If[rk === 1 || (And @@ TorsionFree /@ tanbun),
      (*then*)
        0,
      (*else*)
        dummyindices = Table[NewLowerIndex /@ tanbun, {2}];
        - (1/2) * 
        Sum[ 
         (Signature[#] *
          InsertIndices[f, 
            Join[ {dummyindices[[1,m]],dummyindices[[2,n]]},
                  Take[#,rk-2] ] ] *
          InsertIndices[Tor,
            Join[ Drop[#,rk-2], 
                  Pair /@ {dummyindices[[1,m]],dummyindices[[2,m]]} ]] &)
          @ RotateLeft[{i},k],
          {m,Length[dummyindices[[1]]]},
          {n,Length[dummyindices[[2]]]},
          {k,rk-1}]
      ]
    )
  ];

ExtdStar/:  Conjugate[HoldPattern[ExtdStar[x_]]] := ExtdStar[Conjugate[x]];

(****************** Generic structure equations ******************)

CompatibilityRule := {
  HoldPattern[Extd[Tensor[g_?MetricQ,{L[i_],L[j_]},{}]]] :>
    Tensor[Conn,{L[i],L[j]},{}] + Tensor[Conn,{L[j],L[i]},{}],
  HoldPattern[Del[v_,Tensor[g_?MetricQ,{L[i_],L[j_]},{}]]] :>
    Tensor[Conn,{L[i],L[j]},{}] . v +
    Tensor[Conn,{L[j],L[i]},{}] . v,
  HoldPattern[Extd[Tensor[g_?MetricQ,{U[i_],U[j_]},{}]]] :> 
    - Tensor[Conn,{U[i],U[j]},{}] - Tensor[Conn,{U[j],U[i]},{}],
  HoldPattern[Del[v_,Tensor[g_?MetricQ,{U[i_],U[j_]},{}]]] :>
    - Tensor[Conn,{U[i],U[j]},{}] . v -
      Tensor[Conn,{U[j],U[i]},{}] . v  
  };

FirstStructureRule := {HoldPattern[Extd[Tensor[Basis,{U[i_]},{}]]] :>
  (WedgeFactor[2]/2) Tensor[Tor,{U[i]},{}] +
  (Tensor[Basis,{Pair[#]},{}] ~Wedge~ Tensor[Conn,{#,U[i]},{}] &)
   @ NewLowerIndex[Bundle[i]] };

SecondStructureRule := {
  HoldPattern[ Extd[ Tensor[Conn,{L[i_],U[j_]},{}] ] ] :>
    (WedgeFactor[2]/2) Tensor[Curv,{L[i],U[j]},{}] +
    (Tensor[Conn,{L[i],Pair[#]},{}] ~Wedge~ Tensor[Conn,{#,U[j]},{}] &)
      @ NewLowerIndex[Bundle[i]],
  HoldPattern[ Del[ Tensor[Basis,{L[k_]},{}],
                Tensor[Conn,{L[i_],U[j_],L[l_]},{}] ] *
           Wedge[ a___, 
                  Tensor[Basis,{U[k_]},{}],
                  b___,
                  Tensor[Basis,{U[l_]},{}],
                  c___ ] ] :>
    ( (1/2) Curv[L[i],U[j],L[k],L[l]] +
    (Conn[L[i],U[#],L[k]] Conn[L[#],U[j],L[l]] &) @
        NewBundleSymbol[Bundle[i]] -
      Plus @@ (
        (Conn[L[i],U[j],L[#]] Conn[L[k],U[#],L[l]] +
          (1/2) Conn[L[i],U[j],L[#]] Tor[U[#],L[k],L[l]] &) /@
        NewBundleSymbol /@
        TangentBundle[Bundle[i]]
        ) ) *
    Wedge[ a,Tensor[Basis,{U[k]},{}],b,Tensor[Basis,{U[l]},{}],c],
  HoldPattern[ Del[ Tensor[Basis,{L[k_]},{}],
                Tensor[Conn,{L[i_],U[j_],L[l_]},{}] ] *
           Wedge[ a___, 
                  Tensor[Basis,{U[l_]},{}],
                  b___,
                  Tensor[Basis,{U[k_]},{}],
                  c___ ] ] :>
    ( (1/2) Curv[L[i],U[j],L[k],L[l]] +
    (Conn[L[i],U[#],L[k]] Conn[L[#],U[j],L[l]] &) @
        NewBundleSymbol[Bundle[i]] -
      Plus @@ (
        (Conn[L[i],U[j],L[#]] Conn[L[k],U[#],L[l]] +
          (1/2) Conn[L[i],U[j],L[#]] Tor[U[#],L[k],L[l]] &) /@
        NewBundleSymbol /@
        TangentBundle[Bundle[i]]
        ) ) *
    Wedge[ a,Tensor[Basis,{U[l]},{}],b,Tensor[Basis,{U[k]},{}],c],
  HoldPattern[ a_. * Del[ Tensor[Basis,{L[k_]},{}],
                      Tensor[Conn,{L[i_],U[j_],L[l_]},{}] ] +
           b_. * Del[ Tensor[Basis,{L[l_]},{}],
                      Tensor[Conn,{L[i_],U[j_],L[k_]},{}] ] ] :>
    b * (
      Plus @@ (
      (Conn[L[i], U[j], L[#]]*Tor[U[#], L[k], L[l]] +
       Conn[L[i], U[j], L[#]]*Conn[L[k], U[#], L[l]] - 
       Conn[L[i], U[j], L[#]]*Conn[L[l], U[#], L[k]] & ) /@
         NewBundleSymbol /@
         TangentBundle[Bundle[i]] ) + 
      (Conn[L[i], U[#], L[l]]*Conn[L[#], U[j], L[k]] - 
       Conn[L[i], U[#], L[k]]*Conn[L[#], U[j], L[l]] & ) @
         NewBundleSymbol[Bundle[i]] - 
      Curv[L[i], U[j], L[k], L[l]]
      ) /; a+b===0 
  };

StructureRules := 
  Join[CompatibilityRule,FirstStructureRule,SecondStructureRule];

(******************** Curvature for arbitrary connections *******)

Curvature/: TCompute[ HoldPattern[Curvature[cn_]],
                      {a_,b_,i_,j_}] := Module[{diff},
  diff = cn - Conn;
  TensorExpand[
  Plus @@ Flatten[ {
   InsertIndices[ Curv, {a,b,i,j} ] +
   CovD[ InsertIndices[ diff, {a, b, j} ], {i}] - 
     CovD[ InsertIndices[ diff, {a, b, i} ], {j}],
   (InsertIndices[ diff, {a, U[#], j} ] * 
       InsertIndices[ diff, {L[#], b, i} ] -
    InsertIndices[ diff, {a, U[#], i} ] * 
       InsertIndices[ diff, {L[#], b, j} ] &) /@ 
         NewBundleSymbol /@ Bundles[diff][[1]],
   (InsertIndices[ diff, {a, b, L[#]} ] * 
       Tor[U[#], i, j] &) /@
         NewBundleSymbol /@ TangentBundle[Bundle[a]]}
  ]]
  ];

Curvature/: TCompute[ HoldPattern[Curvature[conn_]], 
                      {i_,j_}] := 
  TensorExpand[
    Plus @@ Flatten[ Outer @@ Join[ {
      InsertIndices[ Curvature[conn], {i,j,L[#1],L[#2]} ] * 
        (1/WedgeFactor[2]) *
        Wedge[ Basis[U[#1]], Basis[U[#2]] ]
      &},
      Table[ NewBundleSymbol /@ UnderlyingTangentBundle[conn[L[i],L[j]]], {2} ]
  ]]];

Curvature/: TCompute[ HoldPattern[Curvature[h_]], {i___} ] := (
  Message[ Index::error, "Curvature" ];
  Return[ERROR[ InsertIndices[ Curvature[h], {i} ]]];
  ) /; Length[{i}] =!= 2 && Length[{i}] =!= 4;

(********************* Lie **************************************)

Lie[v_ + w_, t_] := Lie[v,t] + Lie[w,t];
Lie[(c_?ConstantQ) * v_, t_] := c Lie[v,t];
Lie[f_ * v_, w_] := f Lie[v,w] - Inner[Extd[f],w] * v /; 
  VectorFieldQ[v] && VectorFieldQ[w];
Lie[0, _] := 0;
Lie[v_, t_Plus] := Lie[v,#]& /@ t;
Lie[v_, f_ * t_] := Lie[v,f] * t + f * Lie[v,t];

Lie[v_, f_^p_] := 
  NewDummy[p] NewDummy[f]^(p-1) Lie[v,f] + 
  Log[NewDummy[f]] f^NewDummy[p] Lie[v,p];

HoldPattern[Lie[v_, Wedge[x_,y__]]] := 
  Wedge[Lie[v,x],y] + Wedge[x,Lie[v,Wedge[y]]];

HoldPattern[Lie[v_, TensorProduct[x_,y__]]] :=
  TensorProduct[Lie[v,x],y] + TensorProduct[x,Lie[v,TensorProduct[y]]];

HoldPattern[Lie[v_, Inverse[x_]]] := 
  - Dot[ Dot[ NewDummy[Inverse[x]], 
              Lie[v,x]], 
         NewDummy[Inverse[x]] ]; 

HoldPattern[Lie[v_, Sym[x_]]] := Sym[Lie[v,x]];
HoldPattern[Lie[v_, Alt[x_]]] := Alt[Lie[v,x]];

HoldPattern[Lie[v_, f_]] := Del[v,f] /; Rank[f]===0 && VectorFieldQ[v];
HoldPattern[Lie[v_, w_]] := - Lie[w, v] /; 
  VectorFieldQ[v] && VectorFieldQ[w] && !OrderedQ[{v,w}];
HoldPattern[Lie[v_, v_]] := 0;

HoldPattern[ Lie[ Tensor[Basis,{L[i_]},{}], Tensor[Basis,{L[j_]},{}] ] ] :=
  0 /; Bundle[i]===Bundle[j] && CommutingFrame[Bundle[i]];

Lie/: Conjugate[HoldPattern[Lie[x__]]] := Conjugate /@ Lie[x];

Lie/: TCompute[ HoldPattern[Lie[v_, t_]], {i___} ] := 
  Module[{vindices,tindices,variances,m,n},
  If[!VectorFieldQ[v],
    Print["ERROR: ",v," is not a vector field."];
    Return[ ERROR[ Lie[v,t][i] ] ] ];
  vindices = NewLowerIndex /@ UnderlyingTangentBundle[t];
  tindices = NewLowerIndex /@ UnderlyingTangentBundle[t];
  variances = Variance[t];
  InsertIndices[Del[v,t],{i}] +
    Sum[ 
      If[variances[[n]]===Covariant,
        (*then covariant index*)
          InsertIndices[t, Join[ Take[{i},n-1], {vindices[[m]]}, Drop[{i},n] ] ] *
          ( InsertIndices[ Unevaluated[Del[v]], 
                           {Pair[vindices[[m]]], {i}[[n]] } ] +
            Plus @@ ( InsertIndices[ v, {Pair[#]} ] * 
                      InsertIndices[ Tor, {Pair[vindices[[m]]], 
                                           #, 
                                           {i}[[n]]} ] & /@ tindices)
          ),
        (*else contravariant index*)
          InsertIndices[t, Join[ Take[{i},n-1], {Pair[vindices[[m]]]}, Drop[{i},n] ] ] *
          (- InsertIndices[ Unevaluated[Del[v]],
                            { {i}[[n]], vindices[[m]] } ] -
            Plus @@ ( InsertIndices[ v, {Pair[#]} ] * 
                      InsertIndices[ Tor, { {i}[[n]], 
                                           #, 
                                           vindices[[m]]} ] & /@ tindices)
          )
        ],
      {m, Length[vindices]},
      {n, Length[{i}]}
      ]
  ];

LieRule = 
  {HoldPattern[Lie[v_,t_]] :> 
      Int[v, Extd[t]] + Extd[Int[v,t]] /; FormQ[t]
  };

(****************** CovD *************************************************)

(* Covariant differentiation with no indices inserted is a no-op. *)

CovD[x_,{},opts___] := x;

(* First check that the arguments are OK. *)

Options[CovD] := {Metric -> Automatic,
                  Connection -> Conn};

CovD[x_, {i__}, opts___] := Which[
  !IndexQ[i],
    Print["ERROR: invalid index ",i," given to CovD"];
    Return[ERROR[CovD[x,{i},opts]]],
  !optionsOK[CovD,opts], 
    Return[ERROR[CovD[x,{i},opts]]],
  Head[rank[x]]=!=Integer,
    Print["ERROR: invalid tensor expression ",x];
    Return[ERROR[CovD[x,{i},opts]]],
  rank[x]=!=0,    
    Print["ERROR: Covariant derivative of invalid scalar expression ",x];
    Return[ERROR[CovD[x,{i},opts]]],
  True,           
    cd[x,{i},opts]
  ];

cd[x_,{i_,j__},opts___] := cd[ cd[x,{i},opts], {j}, opts];

cd[a_Plus,{i_},opts___] := cd[#,{i},opts]& /@ a;

cd[t_,{},opts___] := t;

cd[f_ * g_, {i_}, opts___] := 
  f*cd[g,{i},opts] + g*cd[f,{i},opts];

cd[f_^p_, {i_}, opts___] :=
  NewDummy[p] NewDummy[f]^(p-1) cd[f,{i},opts] + 
  Log[NewDummy[f]] f^NewDummy[p] cd[p,{i},opts];

cd[HoldPattern[Det[h_]], {i_}, opts___] :=
  InsertIndices[ Del[ Basis[i], Det[h] ], {} ];

cd[HoldPattern[Summation[f_]], {i_}] := 
  cd[f,{i}];

cd[ HoldPattern[ Tensor[Conn,{i_,j_,k_},{}] ], {l_} ] :=
  expandedCovD[ Tensor[Conn,{i,j,k},{}], {l} ];

(* The following special case has to be here because 
   Del[Basis[...],Conn[...]] is considered a component expression. *)

cd[ HoldPattern[ Del[ Tensor[Basis,{i_},{}], x_ ] ], {l_} ] :=
  expandedCovD[ Del[Tensor[Basis,{i},{}],x], {l} ];

cd[ HoldPattern[Tensor[t__]], {k_} ] := 
  TCompute[Tensor[t],{k}];

cd[f_[g_],{i_},opt___] := f'[NewDummy[g]] * cd[g,{i},opt] /;
  mathFunctionQ[f];

cd[Tensor[t_,{a___},{b___}],{k_},Connection->con_] := 
  InsertIndices[ Del[ Nest[Del,t,Length[{b}]], Connection->con], 
                 {a,b,k} ];

cd[Tensor[t_,{a___},{b___}],{L[k_]},Metric->metric_] := 
  cd[Tensor[t,{a},{b}],{L[k]},Connection->LeviCivitaConnection[metric]];

cd[Tensor[t_,{a___},{b___}],{U[k_]},Metric->metric_] := 
  (InsertIndices[Inverse[metric], {U[k],#}] *
    cd[Tensor[t,{a},{b}],{Pair[#]},Connection->LeviCivitaConnection[metric]]& @
    dumdex[ U[k] ]
    );

cd[c_?ConstantQ,{i_},opts___] := 0;

(* Covariant derivative of anything else produces an error message. *)

cd[x_,{i__},opts___] := (
  Print["ERROR: Covariant derivative of invalid component expression "<>x];
  ERROR[ CovD[x,{i},opts] ]
  );

(*************************** subroutines **************************)

(** various subroutines used by the differentiation functions *)

(* dumdex[index] returns a new or upper index for Bundle[index]. *)

dumdex[L[i_]] :=
  NewLowerIndex[Bundle[i]];

dumdex[U[i_]] :=
  NewUpperIndex[Bundle[i]];

(* tandumdex[index] gives a list of new indices for the tangent
   decomposition. *)

tandumdex[L[i_]] := 
  NewLowerIndex /@ TangentBundle[Bundle[i]];

tandumdex[U[i_]] :=
  NewUpperIndex /@ TangentBundle[Bundle[i]];

(********************* Format *********************************)

SetAttributes[{TensorFormat,IndexFormat,IndexTeXFormat,TensorTeXFormat,
               BracketFormat,TimesFormat,InputTensorFormat,
               SymmetricProductFormat,TensorPowerFormat,
               holdlist,indexformQ,SubscriptedNameFormat
               },
              HoldAll];

(** All of the following formatting routines have to be very careful
    not to evaluate their arguments.  That's why there are so many
    HoldForm's and all these functions have the HoldAll attribute. **)

(**** Decide which format to use for tensors. ****)

Tensor/: HoldPattern[Format[Tensor[t_,{i___},{j___}]]] :=
           TensorFormat[t,{i},{j}] /; 
           indexformQ[i,j] && $TensorFormatting;

Tensor/: HoldPattern[Format[Tensor[t_,{i___},{j___}],TeXForm]] :=
           TensorTeXFormat[t,{i},{j}] /; 
           indexformQ[i,j] && $TensorTeXFormatting;

Tensor/: HoldPattern[Format[Tensor[t_,{i___},{j___}],InputForm]] :=
           InputTensorFormat[t,{i},{j}] /; 
           indexformQ[i,j] && $TensorFormatting;

(** indexformQ[i,j,...] does the same thing as IndexQ, but without
    evaluating its arguments. **)

indexformQ[] := True;
indexformQ[L[_Symbol]] := True;
indexformQ[U[_Symbol]] := True;
indexformQ[L[_Symbol[Bar]]] := True;
indexformQ[U[_Symbol[Bar]]] := True;
indexformQ[_] := False;
indexformQ[i_,j__] := indexformQ[i] && indexformQ[j];

(************************* Tensor Formatting *******************)

(* Inverse tensors with indices have to be handled specially. *)

HoldPattern[TensorFormat[Inverse[Tensor[t_,{},{}]],{i_,j_},{}]] := 
  SequenceForm[ "(", Format[Inverse[Tensor[t,{},{}]]], ")", 
                IndexFormat[{i,j},{}] ];

HoldPattern[TensorFormat[Inverse[x_],{i_,j_},{}]] := 
  SequenceForm[ "(", Format[Inverse[HoldForm[x]]], ")", 
                IndexFormat[{i,j},{}] ];

TensorFormat[t_,{i__},{}] := 
  SequenceForm[ TensorFormat[t], 
                IndexFormat[{i},{}] ];

TensorFormat[t_,{i___},{j__}] := 
  SequenceForm[ TensorFormat[t], 
                IndexFormat[{i},{j}] ];

TensorFormat[t_[Bar],{},{}] := 
  BarFormat[TensorFormat[t]];

TensorFormat[t_[Bar]] := 
  BarFormat[TensorFormat[t]];

TensorFormat[t_,{},{}] := TensorFormat[t];

TensorFormat[t_] := HoldForm[t];

(*********** TeX formatting of tensors ***************************)

HoldPattern[TensorTeXFormat[Inverse[Tensor[t_,{},{}]],{i_,j_},{}]] := 
  SequenceForm[ "\\tensor{(", Format[Inverse[Tensor[t,{},{}]],TeXForm], ")}{", 
                IndexTeXFormat[{i,j},{}], "}" ];

HoldPattern[TensorTeXFormat[Inverse[x_],{i_,j_},{}]] := 
  SequenceForm[ "\\tensor{(", Format[Inverse[x],TeXForm], ")}{", 
                IndexTeXFormat[{i,j},{}], "}" ];

TensorTeXFormat[t_,{},{}] := TensorTeXFormat[t];

TensorTeXFormat[t_,{i___},{j___}] :=
  SequenceForm[ "\\tensor{", TensorTeXFormat[t], "}{",
         IndexTeXFormat[{i},{j}], "}" ];

TensorTeXFormat[t_[Bar]] :=
  SequenceForm["\\overline{",TensorTeXFormat[t],"}"];

TensorTeXFormat[t_] := (TeXFormat /. TensorData[t] /. 
  {"" -> HoldForm[t]}) /; TensorQ[t];

TensorTeXFormat[t_] := Format[HoldForm[t],TeXForm];

(************************ Index Formatting ***************************)

IndexFormat[{L[i_[Bar]]},{}] :=
  SequenceForm[ ColumnForm[{"_",HoldForm[i]},Left,Below], " " ];

IndexFormat[{U[i_[Bar]]},{}] :=
  SequenceForm[ ColumnForm[{"_",HoldForm[i]," "},Left,Above], " " ];

IndexFormat[{L[i_]},{}] := Subscript[SequenceForm[HoldForm[i]," "]];

IndexFormat[{U[i_]},{}] := Superscript[SequenceForm[HoldForm[i]," "]];

IndexFormat[{},{}] := "";

IndexFormat[{i_,j__},{}] := 
  SequenceForm[  IndexFormat[{i},{}],IndexFormat[{j},{}] ];

IndexFormat[{i___},{j__}] := 
  SequenceForm[ IndexFormat[{i},{}],Subscript[";"],IndexFormat[{j},{}] ];

(*********************** Input format for tensors **********************)

InputTensorFormat[Inverse[Tensor[t_,{},{}]],{i_,j_},{}] := 
  Format[SequenceForm[ "Inverse[",
                       HoldForm[t],
                       "]",
                       BracketFormat[{i,j},{}] ],
         OutputForm];

InputTensorFormat