(* :Name: Calculus`Tensor2Analysis` *)

(* :Title: Three-Dimensional Vector and Tensor Analysis *)

(* :Author: 
    Initial Package 
       VectorAnalysis - Stephen Wolfram, Bruce K. Sawhill, Jerry B. Keiper, 
    Actual Version 
      Tensor2Analysis - Andrei Constantinescu, Alexander Korsunsky *) 

(* :Summary:
This package performs standard vector and second rank tensor differential 
operations in  orthogonal coordinate systems.  Coordinate transformations are \

also provided.
*)  

(* :Context: Calculus`Tensor2Analysis` *)

(* :Package Version: 3. *)

(* :Copyright: Copyright 1990-1999,  Wolfram Research, Inc.
*)

(* :History:

    VectorAnalysis
	Version 1.0 by Stephen Wolfram, 1988.
	Revised by Bruce K. Sawhill, August 1990.
	Extensively revised by Jerry B. Keiper, November 1990.
	Version 2.1 by M. Trott and J. Keiper, October 1994:
		revised default coordinate symbols.

    Tensor2Analysis
    Version 3. by A.Constantinescu, A Korsunsky August 2002 
          extensily revised computation of scale factors 
          and introduced operators for second order tensors 
	  as presented in the book

          Elasticity with Mathematica (r)
          by A. Constantinescu and A.M. Korsunsky
          published by Cambridge University Press, 2007
          (ISBN-13: 9780521842013) 
    
*)

(* :Keywords:
	coordinates, orthogonal coordinates, curvilinear coordinates,
	div, curl, grad, Laplacian, biharmonic, Jacobian, coordinate
	systems, coordinate transformations, Cartesian, cylindrical,
	spherical, parabolic cylindrical, paraboloidal, elliptic
	cylindrical, prolate spheroidal, oblate spheroidal, bipolar,
	bispherical, toroidal, conical, confocal ellipsoidal, confocal
	paraboloidal.
*)
	
(* :Source:

Version 1:

	Murray R. Spiegel, Schaum's Mathematical Handbook of 
	Formulas and Tables, McGraw-Hill, New York, 1968, pp. 116-130.

	Philip M. Morse and Herman Feshbach, Methods of Theoretical
	Physics, McGraw-Hill, New York, 1953.

	George B. Arfken, Mathematical Methods for Physicists, 3rd ed.,
	Academic Press, New York, 1985.
	
Version 3:

    Malvern, L.E. - Introduction to the Mechanics of Continous Media 
    Prentice-Hall 1969 
    
    Soos, E. and Teodosiu, C. - Calcul tensorial cu aplicatii in mecanica 
    solider, Editura Stiintifica si Enciclopedica, Bucuresti, 1983 	

    A. Constantinescu and A.M. Korsunsky-Elasticity with Mathematica (r)
          Cambridge University Press, 2007
          

*)

(* :Mathematica Version: 2.0 *)

(* :Limitation:
	The branches associated with the signs of the square roots in the
	coordinate transformations between Cartesian coordinates and
	confocal ellipsoidal or confocal paraboloidal cannot be dealt
	with in the manner used in this package.  In particular the
	transformation from the confocal ellipsoidal system always gives
	a point in the octant with Xx, Yy, Zz > 0.  The transformation from
	the confocal paraboloidal coordinate system always gives a point
	in the quadrant with Xx, Yy > 0.
	
	Only right-handed coordinate systems are implemented.

	Parameters and coordinates cannot take on complex values in some
	of the coordinate systems.

	Only three-dimensional coordinate systems are implemented.
	
	Only basic operators for vectors and second-rank tensors are implemented. 
*)

BeginPackage["Calculus`Tensor2Analysis`"]

Unprotect[Cartesian, Cylindrical, Spherical, ParabolicCylindrical,
    Paraboloidal, EllipticCylindrical, ProlateSpheroidal, OblateSpheroidal,
    Bipolar, Bispherical, Toroidal, Conical, ConfocalEllipsoidal,
    ConfocalParaboloidal];

Unprotect[DotProduct, CrossProduct, ScalarTripleProduct, SetCoordinates,
    Coordinates, Parameters, CoordinateRanges, CoordinatesToCartesian,
    CoordinatesFromCartesian, FieldFromCartesian, FieldToCartesian, 
ScaleFactors, ArcLengthFactor,
    JacobianDeterminant, JacobianMatrix, Hessian3Tensor, PreGrad, PostGrad, 
    Grad, PreDiv, PostDiv, Div, Curl, Inc, Nabla, PreNabla, PostNabla, Laplacian,
    Biharmonic, CoordinateSystem,Twirl,GDot,GTr, $VecQ];

Unprotect[Eeta, Llambda, Mmu, Nnu, Pphi, Rr, Ttheta, Uu, Vv, Xx, Xxi, Yy, Zz]

DotProduct::usage = "DotProduct[v1, v2] gives the dot product (sometimes
called inner product or scalar product) of the two vectors v1, v2 in three
space in the default coordinate system.  DotProduct[v1, v2, coordsys] gives
the dot product of v1 and v2 in the coordinate system coordsys."

CrossProduct::usage = "CrossProduct[v1, v2] gives the cross product
(sometimes called vector product) of the two vectors v1, v2 in three space
in the default coordinate system.  CrossProduct[v1, v2, coordsys] gives
the cross product of v1 and v2 in the coordinate system coordsys."

ScalarTripleProduct::usage = "ScalarTripleProduct[v1, v2, v3] gives the \
scalar
triple product of the three vectors v1, v2, and v3 in three space in the
default coordinate system.  ScalarTripleProduct[v1, v2, v3, coordsys] gives
the scalar triple product of v1, v2, and v3 in the coordinate system \
coordsys."

SetCoordinates::usage = "SetCoordinates[coordsys] sets the default
coordinate system to be coordsys with default variables.
SetCoordinates[coordsys[c1, c2, c3]] sets the default coordinate system to be
coordsys with variables c1, c2, and c3.  Certain coordinate systems
have parameters associated with them, and these parameters can be set by
specifying the full description of the coordinate system.  For example,
SetCoordinates[Conical[lambda, mu, nu]] sets only the variables of the
conical coordinate system, but SetCoordinates[Conical[lambda, mu, nu, a, b]]
sets both the variables and the parameters of the conical coordinate
system. * * * * * * * Defined options for the coordinate systems are:
	Cartesian, Cylindrical, Spherical, 
	ParabolicCylindrical, Paraboloidal, EllipticCylindrical, 
	ProlateSpheroidal, OblateSpheroidal, Bipolar,
	Bispherical, Toroidal, Conical, ConfocalEllipsoidal, 
	ConfocaParaboloidal."

CoordinateSystem::usage = "CoordinateSystem gives the name of the default
coordinate system."

Coordinates::usage = "Coordinates[ ] gives a list of the default
coordinate variables in the default coordinate system.
Coordinates[coordsys] gives a list of the default coordinate variables
in the coordinate system coordsys."

Parameters::usage = "Parameters[ ] gives a list of the default parameters
of the default coordinate system.  Parameters[coordsys] gives a list of the
default parameters in the coordinate system coordsys.  (Many of the
coordinate systems do not have parameters.)"

CoordinateRanges::usage = "CoordinateRanges[ ] gives the intervals over which
each of the coordinates in the default coordinate system may range.
CoordinateRanges[coordsys] gives the ranges for each of the coordinates
in the coordinate system coordsys."

ParameterRanges::usage = "ParameterRanges[ ] gives the intervals over which
each the of the parameters (if any) in the default coordinate system may \
range.
ParameterRanges[coordsys] gives the ranges for each of the parameters in
the coordinate system coordsys."

CoordinatesToCartesian::usage =
"CoordinatesToCartesian[pt] gives the Cartesian coordinates
of the point pt given in the default coordinate system.
CoordinatesToCartesian[pt, coordsys] gives the Cartesian
coordinates of the point given in the coordinate system coordsys."

CoordinatesFromCartesian::usage =
"CoordinatesFromCartesian[pt] gives the coordinates in the
default coordinate system of the point pt given in Cartesian coordinates.
CoordinatesFromCartesian[pt, coordsys] gives the coordinates
in the coordinate system coordsys of the point given in Cartesian \
coordinates."

ScaleFactors::usage = "ScaleFactors[pt] gives a list of the scale factors
at the point pt in the default coordinate system.  ScaleFactors[pt, coordsys]
gives a list of the scale factors at the point pt in the coordinate system
coordsys.  If pt is not given, the default coordinate variables are
used."

ArcLengthFactor::usage = "ArcLengthFactor[pt, t] gives the derivative of
the arc length of a curve given by pt with respect to the parameter t in
the default coordinate system.  ArcLengthFactor[pt, t, coordsys] gives the
derivative of the arc length of a curve given by pt with respect to the
parameter t in the coordinate system coordsys. If pt is not given, the
default coordinate variables are used."

FieldFromCartesian::usage =
"FieldFromCartesian[f] gives the components of the scalar, vector or \
second-rank 
tensor field f in the default coordinate system.
FieldFromCartesian[f, coordsys] gives the components of the scalar, vector or \
second-rank 
tensor field f in the coordinate system coordsys."

FieldToCartesian::usage =
"FieldFromCartesian[f] gives the components of the scalar, vector or \
second-rank 
tensor field f in the Cartesian coordinates system if actual components are \
in the 
default coordinate system.
FieldFromCartesian[f, coordsys] gives the components of the scalar, vector or \
second-rank 
tensor field f in the Cartesian coordinates system if actual components are \
in the 
coordinate system coordsys."

JacobianDeterminant::usage = "JacobianDeterminant[pt] gives the determinate
of the Jacobian matrix of the transformation from the default coordinate
system to the Cartesian coordinate system at the point pt.
JacobianDeterminant[pt, coordsys] gives the determinate of the Jacobian
matrix of the transformation from the coordinate system coordsys to the
Cartesian coordinate system at the point pt.  If pt is not given, the
default coordinate variables are used."

JacobianMatrix::usage = "JacobianMatrix[pt] gives the Jacobian matrix of
the transformation from the default coordinate system to the Cartesian
coordinate system at the point pt.  JacobianMatrix[pt, coordsys] gives the
Jacobian matrix of the transformation from the coordinate system coordsys
to the Cartesian coordinate system at the point pt.  If pt is not given,
the default coordinate variables are used."

Hessian3Tensor::usage = "Hessian3Tensor[pt] gives the second derivative of
the transformation from the default coordinate system to the Cartesian
coordinate system at the point pt. Its a third order tensor.  
Hessian3Tensor[pt, coordsys] gives the
second derivative of the transformation from the coordinate system coordsys
to the Cartesian coordinate system at the point pt.  If pt is not given,
the default coordinate variables are used."

PreNabla::usage = "PreNabla[f] is a differential operator which in cartesian 
coordinates applies as follows: 
PreNabla_a f_bc = [f_bc,a] e_a x e_b x e_c"

PostNabla::usage = "PostNabla[f] is a differential operator which in cartesian 
coordinates applies as follows: 
f_ab PostNabla_c = [f_ab,c] e_a x e_b x e_c"

Nabla::usage = "Nabla[f] is a differential operator which is equivalent to PreNabla 
and in cartesian coordinates applies as follows: 
Nabla_a f_bc = [f_bc,a] e_a x e_b x e_c"

IntegrateGrad::usage = "IntegrateGrad[nf_vector] is an integral 
operator which finds a field f such that Nabla[f] = nf"


PreGrad::usage = "PreGrad[f] gives the gradient of the scalar, vector, or second \
rank tensor valued function f in the default coordinate system. 
PreGrad[f, coordsys] gives the gradient of f in the coordinate system coordsys.
PreGrad in cartesian coordinates operates on the following indices:
PreGrad_a f_bc = [f_bc,a] e_a x e_b x e_c"

PostGrad::usage = "PostGrad[f] gives the gradient of the scalar, vector, or second \
rank tensor valued function f in the default coordinate system. 
PostGrad[f, coordsys] gives the gradient of f in the coordinate system coordsys.
PostGrad in cartesian coordinates operates on the following indices:
f_ab PostGrad_c = [f_ab,c] e_a x e_b x e_c"

Grad::usage = "Grad[f] gives the gradient of the scalar, vector, or second \
rank tensor valued function f in the default coordinate system. 
Grad[f, coordsys] gives the gradient of f in the coordinate system coordsys.
Grad in cartesian coordinates operates on the following indices:
f_ab Grad_c = [f_ab,c] e_a x e_b x e_c
and is equivalent to PostGrad[f]"

PreDiv::usage = "PreDiv[f] gives the divergence of the vector- or second-rank 
tensor-valued function f in the default coordinate system. 
PreDiv[f, coordsys] gives the divergence of f
in the coordinate system coordsys.
PreDiv[f] in cartesian coordinates operates on the following indices:
PreDiv[f_ab]= f_ab,a"

PostDiv::usage = "PostDiv[f] gives the divergence of the vector- or second-rank 
tensor-valued function f in the default coordinate system. 
PostDiv[f, coordsys] gives the divergence of f
in the coordinate system coordsys.
PostDiv[f] in cartesian coordinates operates on the following indices:
PostDiv[f_ab]= f_ab,b"

Div::usage = "Div[f] gives the divergence of the vector- or second-rank 
tensor-valued function f in the default coordinate system. 
Div[f, coordsys] gives the divergence of f
in the coordinate system coordsys.
Div[f] in cartesian coordinates operates on the following indices:
Div[f_ab]= f_ab,b
and is equivalent to PostDiv[f]"

Curl::usage = "Curl[f] gives the curl of the vector- or second-rank \
tensor-valued function f 
in the default coordinate system. Curl[f, coordsys] gives the curl of f in \
the
coordinate system coordsys."

Laplacian::usage = "Laplacian[f] gives the Laplacian of the scalar,
vector or second-rank tensor-valued function f in the default coordinate \
system.
Laplacian[f, coordsys] gives the Laplacian of f in the coordinate
system coordsys."

Biharmonic::usage = "Biharmonic[f] gives the Laplacian of the Laplacian
of the scalar vector or second-rank tensor-valued function f in the default \
coordinate system.
Biharmonic[f, coordsys] gives the biharmonic of f in the coordinate system
coordsys."

Inc::usage = "Inc[t] computes the incompatibility tensor of the tensor t in \
the default coordinate 
system. Inc[t, coordsys] computes the incompatibility in the coordinate \
system cordsys."

Cartesian::usage = "Cartesian represents the Cartesian coordinate system
with default variables Xx, Yy, and Zz.  Cartesian[x, y, z] represents the
Cartesian coordinate system with variables x, y, and z."

Cylindrical::usage = "Cylindrical represents the cylindrical coordinate \
system
with default variables Rr, Ttheta, and Zz.  Cylindrical[r, theta, z]
represents the cylindrical coordinate system with variables r, theta, and
z."

Spherical::usage = "Spherical represents the spherical coordinate system
with default variables Rr, Ttheta, and Pphi.  Spherical[r, theta, phi]
represents the spherical coordinate system with variables r, theta, and
phi."

ParabolicCylindrical::usage = "ParabolicCylindrical represents the parabolic
cylindrical coordinate system with default variables Uu, Vv, and Zz.
ParabolicCylindrical[u, v, z] represents the parabolic cylindrical coordinate
system with variables u, v, and z."
	(* Parabolic cylindrical coordinates: Holding either of the first
	two variables constant while the third variable is held constant
	produces parabolas facing in opposite directions; the third
	variable specifies distance along the axis of common focus. *)
	
Paraboloidal::usage = "Paraboloidal represents the paraboloidal coordinate
system with default variables Uu, Vv, and Pphi.  Paraboloidal[u, v, phi]
represents the paraboloidal coordinate system with variables u, v, and
phi."
	(* Paraboloidal coordinates: Holding either of the first two
	variables constant while the third variable is held constant
	produces opposite facing parabolas; the third parameter
	describes rotations about their common bisectors. *)

EllipticCylindrical::usage = "EllipticCylindrical represents the elliptic
cylindrical coordinate system with default variables Uu, Vv, and Zz and
default parameter value 1.  EllipticCylindrical[u, v, z] represents the
elliptic cylindrical coordinate system with variables u, v, and z and
default parameter value 1.  EllipticCylindrical[u, v, z, a] represents the
elliptic cylindrical coordinate system with variables u, v, and z and
parameter a."
	(* Elliptic cylindrical coordinates: The parameter a is one-half
	the distance between the two foci.  Holding the first variable
	constant produces confocal ellipses, holding the second variable
	constant produces confocal hyperbolas, and the third variable
	specifies distance along the axis of common focus. *)
	
ProlateSpheroidal::usage = "ProlateSpheroidal represents the prolate \
spheroidal
coordinate system with default variables Xxi, Eeta, and Pphi and
default parameter value 1.  ProlateSpheroidal[xi, eta, phi] represents the
prolate spheroidal coordinate system with variables xi, eta, and phi and
default parameter value 1.  ProlateSpheroidal[xi, eta, phi, a] represents the
prolate spheroidal coordinate system with variables xi, eta, and phi and
parameter a." 
	(* Prolate spheroidal coordinates: These coordinates
	are obtained by rotating elliptic cylindrical coordinates 
	about the axis joining the two foci, which are separated by
	a distance 2*a.  The third variable parameterizes the rotation. *)

OblateSpheroidal::usage = "OblateSpheroidal represents the oblate spheroidal
coordinate system with default variables Xxi, Eeta, and Pphi and 
default parameter value 1.  OblateSpheroidal[xi, eta, phi] represents the
oblate spheroidal coordinate system with variables xi, eta, and phi and
default parameter value 1.  OblateSpheroidal[xi, eta, phi, a] represents the
oblate spheroidal coordinate system with variables xi, eta, and phi and
parameter a."
	(* Oblate spheroidal coordinates: These coordinates are
    	obtained by rotating elliptic cylindrical coordinates about an
    	axis perpendicular to the axis joining the two foci, which are
    	separated by a distance 2*a.  The third variable parameterizes the
    	rotation. *)
         
Bipolar::usage = "Bipolar represents the bipolar coordinate system with 
default variables Uu, Vv, and Zz and default parameter value 1.
Bipolar[u, v, z] represents the bipolar coordinate system with variables
u, v, and z and default parameter value 1.  Bipolar[u, v, z, a] represents \
the
bipolar coordinate system with variables u, v, and z, and parameter a."
	(* Bipolar coordinates: These coordinates are built around two
	foci, separated by a distance 2*a.  Holding the first variable
	fixed produces circles that pass through both foci.  Holding
	the second variable fixed produces degenerate ellipses about one
	of the foci. The third variable parameterizes distance along the
	axis of common foci. *)

Bispherical::usage = "Bispherical represents the bispherical coordinate \
system
with default variables Uu, Vv, and Pphi and default parameter value 1.
Bispherical[u, v, phi] represents the bispherical coordinate system with
variables u, v, and phi and default parameter value 1.
Bispherical[u, v, phi, a] represents the bispherical coordinate system with
variables u, v, and phi and parameter a."
	(* Bispherical coordinates: These coordinates are related
	to bipolar coordinates except that the third coordinate
	measures azimuthal angle rather than Zz. *)

Toroidal::usage = "Toroidal represents the toroidal coordinate system with
default variables Uu, Vv, and Pphi and default parameter value 1.
Toroidal[u, v, phi] represents the toroidal coordinate system with variables
u, v, and phi and default parameter value 1.  Toroidal[u, v, phi, a]
represents the toroidal coordinate system with variables u, v, phi and
parameter a."
	(* Toroidal coordinates: These are obtained by rotating bipolar
	coordinates about an axis perpendicular to the line joining the
	two foci which are separated by a distance 2*a.  The third variable
	parameterizes this rotation. *)
	
Conical::usage = "Conical represents the conical coordinate system with \
default
variables Llambda, Mmu, and Nnu and default parameter values 1 and 2.
Conical[lambda, mu, nu] represents the conical coordinate system with \
variables
lambda, mu, and nu and default parameter values 1 and 2.  
Conical[lambda, mu, nu, a, b] represents the conical coordinate system with
variables lambda, mu, and nu and parameters a and b."
	(* Conical coordinates:  The coordinate surfaces are spheres
	with centers at the origin and radius Llambda (Llambda constant),
	cones with apexes at the origin and axes along the z-axis
	(Mmu constant), and cones with apexes at the origin and axes
	along the y-axis (Nnu constant). *)
	
ConfocalEllipsoidal::usage = "ConfocalEllipsoidal represents the confocal
ellipsoidal coordinate system with default variables Llambda, Mmu, and Nnu
and default parameter values 3, 2, and 1.  ConfocalEllipsoidal[lambda, mu, \
nu]  
represents the confocal ellipsoidal coordinate system with variables
lambda, mu, and nu and default parameter values 3, 2, and 1.
ConfocalEllipsoidal[lambda, mu, nu, a, b, c] represents the confocal
ellipsoidal coordinate system with variables lambda, mu, and nu and
parameters a, b, and c."
	(* Confocal ellipsoidal coordinates: The coordinate surfaces are
	ellipsoids (Llambda constant), hyperboloids of one sheet (Mmu
	constant), and hyperboloids of two sheets (Nnu constant). *)

ConfocalParaboloidal::usage = "ConfocalParaboloidal represents the confocal
paraboloidal coordinate system with default variables Llambda, Mmu, and
Nnu and default parameter values 2 and 1.  ConfocalParaboloidal[lambda, mu, \
nu] 
represents the confocal paraboloidal oordinate system with variables
lambda, mu, and nu and default parameter values 2 and 1.
ConfocalParaboloidal[lambda, mu, nu, a, b] represents the confocal
paraboloidal coordinate system with variables lambda, mu, and nu
and parameters a and b."
	(* Confocal paraboloidal coordinates: The coordinate surfaces are
	confocal families of elliptic paraboloids extending in the
	direction of the negative z-axis (Llambda constant), hyperbolic
	paraboloids (Mmu constant), and elliptic paraboloids extending in
	the direction of the positive z-axis (Nnu constant). *)


(* default coordinates *)

Xx::usage = Yy::usage =
"{Xx,Yy,Zz} are the default Cartesian coordinates."

Zz::usage =
"{Xx,Yy,Zz} are the default Cartesian coordinates; {Rr,Ttheta,Zz} are the
default Cylindrical coordinates; {Uu,Vv,Zz} are the default coordinates for \
the
ParabolicCylindrical, EllipticCylindrical, and Bipolar coordinate systems."

Rr::usage = Ttheta::usage =
"{Rr,Ttheta,Zz} are the default Cylindical coordinates; {Rr,Ttheta,Pphi} are \
the
default Spherical coordinates."

Pphi::usage =
"{Rr,Ttheta,Pphi} are the default Spherical coordinates; {Uu,Vv,Pphi} are the \

default coordinates for the Paraboloidal, Bispherical, and Toroidal \
coordinate 
systems; {Xxi,Eeta,Pphi} are the default coordinates for the Prolate \
Spheroidal 
and Oblate Spheroidal coordinate systems."

Uu::usage = Vv::usage =
"{Uu,Vv,Zz} are the default coordinates for the ParabolicCylindrical,
EllipticCylindrical, and Bipolar coordinate systems; {Uu,Vv,Pphi} are the 
default coordinates for the Paraboloidal,  Bispherical, and Toroidal \
coordinate
systems."

Xxi::usage = Eeta::usage =
"{Xxi,Eeta,Pphi} are the default coordinates for the ProlateSpheroidal and
OblateSpheroidal coordinate systems."

Llambda::usage = Mmu::usage = Nnu::usage =
"{Llambda,Mmu,Nnu} are the default coordinates for the Conical,
ConfocalEllipsoidal, and ConfocalParaboloidal coordinate systems."

Twirl::usage = "This performs the Transpose as per defintion (correctly). 
More precisely, Twirl[t, ord] performs the transpose of the list of indices \
ord for the tensor t.
Twirl[t, n] performs the transpose of the indice n with the last indice for \
the tensor t.
Twirl[t] performs the transpose of the 2 indices of the tensor t."

GDot::usage = "GDot[tf, tg, nf, ng] Forms the dot product over the indices nf \
and ng of two tensors
tf and tg."

GTr::usage = "GTr[t, n1, n2] performs the generalized trace (dot product) in 
the indices n1 and n2 of the tensor t "

IntGradParameters::usage = "IntegrationConstants is an option of IntegrateGrad[]  \
which defines the name of the integration constants "

Options[IntegrateGrad] = {IntegrationConstants -> C};

Begin["Calculus`Tensor2Analysis`Private`"]

Clear[$CoordinateSystem];  (* a person might read in this package twice! *)

$ScaQ[v_] := Not[ListQ[v]] || (VectorQ[v] && Length[v] == 1);
$VecQ[v_] := VectorQ[v] && Length[v] == 3;
$TenQ[v_] := MatrixQ[v] && Dimensions[v] == {3,3};

    (* these are the coordinate systems that are supported *)
$CoordSysList = {Cartesian, Cylindrical, Spherical, ParabolicCylindrical,
    Paraboloidal, EllipticCylindrical, ProlateSpheroidal, OblateSpheroidal,
    Bipolar, Bispherical, Toroidal, Conical, ConfocalEllipsoidal,
    ConfocalParaboloidal}

Coordinates::invalid = "`1` is not a valid coordinate system specification."

$BadCoordSys[coordsys_] := (Message[Coordinates::invalid, coordsys]; $Failed)

$ExpandCoordSys[coordsys_Symbol] :=
	If[MemberQ[$CoordSysList, coordsys],
		Apply[coordsys, Flatten[{Coordinates[coordsys],
				Parameters[coordsys]}]],
	    (* else *)
		$BadCoordSys[coordsys]];

$ExpandCoordSys[coordsys_[vp___]] :=
    Module[{vpars = {vp}, np, len},
	If[!MemberQ[$CoordSysList, coordsys],
		Return[$BadCoordSys[coordsys[vp]]]];
	len = Length[vpars];
	np = Length[Parameters[coordsys]];
	Which[
	    len == 0,
		Apply[coordsys,
		    Flatten[{Coordinates[coordsys], Parameters[coordsys]}]],
	    len == 3,
		If[Head /@ vpars =!= {Symbol, Symbol, Symbol},
		    $BadCoordSys[coordsys[vp]],
		    Apply[coordsys, Flatten[{vpars, Parameters[coordsys]}]]],
	    len == 3 + np,
		If[(Head /@ Take[vpars, 3]) =!= {Symbol, Symbol, Symbol},
		    $BadCoordSys[coordsys[vp]], Apply[coordsys, vpars]],
	    True,
		$BadCoordSys[coordsys[vp]]
	    ]
	]

$PointRule[cs_, pt_] :=
	{cs[[1]] -> pt[[1]], cs[[2]] -> pt[[2]], cs[[3]] -> pt[[3]]}

$DAbsSign = {Sign'[_] -> 0, Abs'[x_] -> Sign[x]};

  (* ===================== Dot Product ====================== *)

DotProduct[v1_?$VecQ, v2_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], cv1, cv2},
	cv1 . cv2 /; (cs =!= $Failed &&
			(cv1 = $CTToCart[v1, cs]) =!= $Failed &&
			(cv2 = $CTToCart[v2, cs]) =!= $Failed)];

  (* ===================== Cross Product ====================== *)

CrossProduct[v1_?$VecQ, v2_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], cv1, cv2, m},
	(m = Minors[{cv1, cv2}, 2][[1]];
	$CTFromCart[{m[[3]], -m[[2]], m[[1]]}, cs]) /; (cs =!= $Failed &&
				(cv1 = $CTToCart[v1, cs]) =!= $Failed &&
				(cv2 = $CTToCart[v2, cs]) =!= $Failed)];

  (* ================== Scalar Triple Product ====================== *)

ScalarTripleProduct[v1_?$VecQ, v2_?$VecQ, v3_?$VecQ,
					coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], cv1, cv2, cv3},
	Det[{cv1, cv2, cv3}] /; (cs =!= $Failed &&
				(cv1 = $CTToCart[v1, cs]) =!= $Failed &&
				(cv2 = $CTToCart[v2, cs]) =!= $Failed &&
				(cv3 = $CTToCart[v3, cs]) =!= $Failed)];

  (* ==================== Coordinates ========================== *)

Coordinates[coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	(List @@ Take[cs, 3]) /; cs =!= $Failed]

  (* ==================== Parameters ========================== *)

Parameters[coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	(List @@ Drop[cs, 3]) /; cs =!= $Failed]

  (* ==================== SetCoordinates ========================== *)

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

Coordinates[Cartesian] ^= {Xx, Yy, Zz}
Coordinates[Cylindrical] ^= {Rr, Ttheta, Zz}
Coordinates[Spherical] ^= {Rr, Ttheta, Pphi} 
Coordinates[ParabolicCylindrical] ^= {Uu, Vv, Zz}
Coordinates[Paraboloidal] ^= {Uu, Vv, Pphi} 
Coordinates[EllipticCylindrical] ^= {Uu, Vv, Zz}
Coordinates[ProlateSpheroidal] ^= {Xxi, Eeta, Pphi} 
Coordinates[OblateSpheroidal] ^= {Eeta, Xxi, Pphi} 
Coordinates[Bipolar] ^= {Uu, Vv, Zz}
Coordinates[Bispherical] ^= {Uu, Vv, Pphi}
Coordinates[Toroidal] ^= {Vv, Uu, Pphi} 
Coordinates[Conical] ^= {Llambda, Mmu, Nnu} 
Coordinates[ConfocalEllipsoidal] ^= {Llambda, Mmu, Nnu} 
Coordinates[ConfocalParaboloidal] ^= {Llambda, Mmu, Nnu}

Parameters[Cartesian] ^= {}
Parameters[Cylindrical] ^= {}
Parameters[Spherical] ^= {} 
Parameters[ParabolicCylindrical] ^= {}
Parameters[Paraboloidal] ^= {} 
Parameters[EllipticCylindrical] ^= {1}
Parameters[ProlateSpheroidal] ^= {1}
Parameters[OblateSpheroidal] ^= {1} 
Parameters[Bipolar] ^= {1}
Parameters[Bispherical] ^= {1}
Parameters[Toroidal] ^= {1} 
Parameters[Conical] ^= {1, 2} 
Parameters[ConfocalEllipsoidal] ^= {3, 2, 1} 
Parameters[ConfocalParaboloidal] ^= {2, 1}

CoordinateSystem := $CoordinateSystem

SetCoordinates::range = "The parameter `1` is not positive."

SetCoordinates::ranges = "The parameters `1` do not satisfy the range
requirements `2` for `3`."

SetCoordinates[coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], hcs, rr, par},
	cs /; ((cs =!= $Failed) &&
		(rr = ParameterRanges[cs];
		hcs = Head[cs];
		par = Apply[List, Drop[cs, 3]];
		If[rr =!= Null && (Evaluate[rr]& @@ par === False),
		    If[Length[par] == 1,
			Message[SetCoordinates::range, par[[1]]],
			Message[SetCoordinates::ranges, par, rr, hcs]
			];
		    False,
		  (* else *)
		    Unprotect[Evaluate[hcs]];
		    Coordinates[hcs] ^= Apply[List, Take[cs, 3]];
		    Parameters[hcs] ^= par;
		    Protect[Evaluate[hcs]];
		    $CoordinateSystem = hcs;
		    True
		    ]))
	];

  (* ================== Coordinate Ranges ====================== *)

CoordinateRanges[coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	Switch[Head[cs],
	    Cartesian,
		{ -Infinity < cs[[1]] < Infinity,
		  -Infinity < cs[[2]] < Infinity,
		  -Infinity < cs[[3]] < Infinity },
	    Cylindrical,
		{ 0 <= cs[[1]] < Infinity,
		  -Pi < cs[[2]] <= Pi,
		  -Infinity < cs[[3]] < Infinity },
	    Spherical,
		{ 0 <= cs[[1]] < Infinity,
		  0 <= cs[[2]] <= Pi,
		  -Pi < cs[[3]] <= Pi },
	    ParabolicCylindrical,
		{ -Infinity < cs[[1]] < Infinity,
		  0 <= cs[[2]] < Infinity,
		  -Infinity < cs[[3]] < Infinity },
	    Paraboloidal,
		{ 0 <= cs[[1]] < Infinity,
		  0 <= cs[[2]] < Infinity,
		  -Pi <= cs[[3]] <= Pi },
	    EllipticCylindrical,
		{ 0 <= cs[[1]] < Infinity,
		  -Pi < cs[[2]] <= Pi,
		  -Infinity < cs[[3]] < Infinity },
	    ProlateSpheroidal,
		{ 0 <= cs[[1]] < Infinity,
		  0 <= cs[[2]] <= Pi,
		  -Pi < cs[[3]] <= Pi },
	    OblateSpheroidal,
		{ -Pi/2 <= cs[[1]] <= Pi/2,
		  0 <= cs[[2]] < Infinity,
		  -Pi < cs[[3]] <= Pi },
	    Bipolar,
		{ -Pi < cs[[1]] <= Pi,
		  -Infinity < cs[[2]] < Infinity,
		  -Infinity < cs[[3]] < Infinity },
	    Bispherical,
		{ 0 <= cs[[1]] <= Pi,
		  -Infinity < cs[[2]] < Infinity,
		  -Pi < cs[[3]] <= Pi },
	    Toroidal,
		{ 0 <= cs[[1]] < Infinity,
		  -Pi < cs[[2]] <= Pi,
		  -Pi < cs[[3]] <= Pi },
	    Conical,
		{ -Infinity < cs[[1]] < Infinity,
		  cs[[4]]^2 < cs[[2]]^2 < cs[[5]]^2,
		  cs[[3]]^2 < cs[[4]]^2 },
	    ConfocalEllipsoidal,
		{ -Infinity < cs[[1]] < cs[[6]]^2,
		  cs[[6]]^2 < cs[[2]] < cs[[5]]^2,
		  cs[[5]]^2 < cs[[3]] < cs[[4]]^2 },
	    ConfocalParaboloidal,
		{ -Infinity < cs[[1]] < cs[[5]]^2,
		  cs[[5]]^2 < cs[[2]] < cs[[4]]^2,
		  cs[[4]]^2 < cs[[3]] < Infinity }
	] /; (cs =!= $Failed)
    ]

  (* ================== Parameter Ranges ====================== *)

ParameterRanges[coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	Switch[Head[cs],
	    Conical,
		0 < #1 < #2 < Infinity,
	    ConfocalEllipsoidal,
		0 < #3 < #2 < #1 < Infinity,
	    ConfocalParaboloidal,
		0 < #2 < #1 < Infinity,
	    _ (* all others *),
		If[Length[cs] == 4, 0 < #1 < Infinity, Null]
	] /; (cs =!= $Failed)
    ]

    (* ==================== ScaleFactors ===========================*)

ScaleFactors[arg_:$CoordinateSystem] :=
    Module[{pt, cs},
	ScaleFactors[pt, cs] /; (If[$VecQ[arg],
					cs = $CoordinateSystem;
					pt = arg,
					cs = $ExpandCoordSys[arg];
					If[cs =!= $Failed,
					    pt = List @@ Take[cs, 3]]];
				cs =!= $Failed)];

ScaleFactors[pt_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], tmp},
	(tmp = Switch[Head[cs],
	    Cartesian, {1, 1, 1},
	    Cylindrical, {1, cs[[1]], 1},
	    Spherical, {1, cs[[1]], cs[[1]] Sin[cs[[2]]]},
	    ParabolicCylindrical, 
		tmp = Sqrt[cs[[1]]^2 + cs[[2]]^2];
		{tmp, tmp, 1},
	    Paraboloidal,
		tmp = Sqrt[cs[[1]]^2 + cs[[2]]^2];
		{tmp, tmp, cs[[1]] cs[[2]]},
	    EllipticCylindrical,
		tmp = cs[[4]] Sqrt[Sinh[cs[[1]]]^2 + Sin[cs[[2]]]^2];
		{tmp, tmp, 1},
	    ProlateSpheroidal,
		tmp = cs[[4]] Sqrt[Sinh[cs[[1]]]^2+Sin[cs[[2]]]^2];
		{tmp, tmp, cs[[4]] Sinh[cs[[1]]] Sin[cs[[2]]]},
	    OblateSpheroidal,
		tmp = cs[[4]] Sqrt[Sinh[cs[[2]]]^2+Sin[cs[[1]]]^2];
		{tmp, tmp, cs[[4]] Cosh[cs[[2]]] Cos[cs[[1]]]},
	    Bipolar,
		tmp = cs[[4]]/(Cosh[cs[[2]]] - Cos[cs[[1]]]);
		{tmp, tmp, 1},
	    Bispherical,
		tmp = cs[[4]]/(Cosh[cs[[2]]] - Cos[cs[[1]]]);
		{tmp, tmp, Sin[cs[[1]]] tmp},
	    Toroidal,
		tmp = cs[[4]]/(Cosh[cs[[1]]] - Cos[cs[[2]]]);
		{tmp, tmp, tmp Sinh[cs[[1]]]},
	    Conical,
		tmp = Abs[cs[[1]]] Sqrt[cs[[2]]^2-cs[[3]]^2];
		{1, tmp/(Sqrt[cs[[2]]^2-cs[[4]]^2] Sqrt[cs[[5]]^2-cs[[2]]^2]),
		   tmp/(Sqrt[cs[[3]]^2-cs[[4]]^2] Sqrt[cs[[3]]^2-cs[[5]]^2])},
	    ConfocalEllipsoidal,
		tmp = 2 Sqrt[(cs[[4]]^2 - #)(cs[[5]]^2 - #)(cs[[6]]^2 - #)]&;
		{Sqrt[(cs[[2]]-cs[[1]])(cs[[3]]-cs[[1]])]/tmp[cs[[1]]],
		 Sqrt[(cs[[3]]-cs[[2]])(cs[[1]]-cs[[2]])]/tmp[cs[[2]]],
		 Sqrt[(cs[[1]]-cs[[3]])(cs[[2]]-cs[[3]])]/tmp[cs[[3]]]},
	    ConfocalParaboloidal,
		tmp = 2 Sqrt[(cs[[4]]^2 - #)(cs[[5]]^2 - #)]&;
		{Sqrt[(cs[[2]]-cs[[1]])(cs[[3]]-cs[[1]])]/tmp[cs[[1]]],
		 Sqrt[(cs[[3]]-cs[[2]])(cs[[1]]-cs[[2]])]/tmp[cs[[2]]],
		 Sqrt[(cs[[1]]-cs[[3]])(cs[[2]]-cs[[3]])]/tmp[cs[[3]]]}
	    ];
	tmp /. $PointRule[cs, pt]) /; (cs =!= $Failed)
    ]

  (* ================= CoordinatesToCartesian ===================*)

CoordinatesToCartesian[v_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], cv},
	cv /; (cs =!= $Failed && (cv = $CTToCart[v, cs]) =!= $Failed)
    ]

CoordinatesToCartesian::range = "The point `1` does not satisfy the
range requirements `2` for `3`."

$CTToCart[{a_, b_, c_}, cs_] :=
    Module[{tmp, d},
	If[MemberQ[{Conical,ConfocalEllipsoidal,ConfocalParaboloidal},Head[cs]],
	    tmp = CoordinateRanges[cs];
	    tmp = And @@ (tmp /. {cs[[1]] -> a, cs[[2]] -> b, cs[[3]] -> c});
	    If[tmp === False,
		Message[CoordinatesToCartesian::range,
			{a, b, c}, CoordinateRanges[cs], cs];
		Return[$Failed]
		]
	    ];
	Switch[Head[cs],
	    Cartesian, {a, b, c},
	    Cylindrical, {a Cos[b], a Sin[b], c},
	    Spherical, {a Sin[b] Cos[c], a Sin[b] Sin[c], a Cos[b]},
	    ParabolicCylindrical, {(a^2 - b^2)/2, a b, c},
	    Paraboloidal, {a b Cos[c], a b Sin[c], (a^2 - b^2)/2},
	    EllipticCylindrical,
		{cs[[4]] Cosh[a] Cos[b], cs[[4]] Sinh[a] Sin[b], c},
	    ProlateSpheroidal,
		tmp = cs[[4]] Sinh[a] Sin[b];
		{tmp Cos[c], tmp Sin[c], cs[[4]] Cosh[a] Cos[b]},
	    OblateSpheroidal,
		tmp = cs[[4]] Cosh[b] Cos[a];
		{tmp Cos[c], tmp Sin[c], cs[[4]] Sinh[b] Sin[a]},
	    Bipolar,
		tmp = cs[[4]]/(Cosh[b]-Cos[a]);
		{tmp Sinh[b], tmp Sin[a], c},
	    Bispherical,
		tmp = cs[[4]]/(Cosh[b]-Cos[a]);
		d = Sin[a];
		{tmp d Cos[c], tmp d Sin[c], tmp Sinh[b]},
	    Toroidal,
		tmp = cs[[4]]/(Cosh[a]-Cos[b]);
		d = Sinh[a];
		{tmp d Cos[c], tmp d Sin[c], tmp Sin[b]},
	    Conical,
		tmp = 1/(cs[[4]]^2 - cs[[5]]^2);
		{ a Abs[b c] / (cs[[4]] cs[[5]]),
		  Abs[a] Sign[b] Sqrt[(b^2 - cs[[4]]^2)(c^2 - cs[[4]]^2) tmp]/
			cs[[4]],
		  Abs[a] Sign[c] Sqrt[(b^2 - cs[[5]]^2)(c^2 - cs[[5]]^2)(-tmp)]/
			cs[[5]] },
	    ConfocalEllipsoidal,
		tmp = (#^2 - a)(#^2 - b)(#^2 - c)&;
		d = 1/((cs[[#1]]^2 - cs[[#2]]^2)(cs[[#1]]^2 - cs[[#3]]^2))&;
		d = { Sqrt[tmp[cs[[4]]] d[4, 5, 6]],
		      Sqrt[tmp[cs[[5]]] d[5, 4, 6]],
		      Sqrt[tmp[cs[[6]]] d[6, 4, 5]] },
	    ConfocalParaboloidal,
		tmp = (#^2 - a)(#^2 - b)(#^2 - c)&;
		{ Sqrt[tmp[cs[[4]]]/(cs[[5]]^2 - cs[[4]]^2)],
		  Sqrt[tmp[cs[[5]]]/(cs[[4]]^2 - cs[[5]]^2)],
		  (cs[[4]]^2 + cs[[5]]^2 - a - b - c)/2}
	]
    ]

  (* ================= CoordinatesFromCartesian ===================*)

$ArcTan[x_, y_] := If[TrueQ[x == y == 0], 0, ArcTan[x, y]];

CoordinatesFromCartesian[v_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	$CTFromCart[v, cs] /; (cs =!= $Failed)
    ]

$CTFromCart[{x_, y_, z_}, cs_] :=
    Module[{tmp, a, b, c, d, e},
	Switch[Head[cs],
	    Cartesian, {x, y, z},
	    Cylindrical, {Sqrt[x^2 + y^2], $ArcTan[x, y], z},
	    Spherical,
		a = Sqrt[x^2 + y^2 + z^2];
		{a, If[TrueQ[a == 0], 0, ArcCos[z/a]], $ArcTan[x, y]},
	    ParabolicCylindrical,
		d = Sqrt[x^2 + y^2];
		If[TrueQ[d == 0],
		    a = b = d,
		    If[TrueQ[x < 0],
			a = Sqrt[d + x] Sign[y]; b = y/a,
			b = Sqrt[d - x]; a = y/b]];
		{a, b, z},
	    Paraboloidal,
		d = Sqrt[x^2 + y^2 + z^2];
		If[TrueQ[d == 0],
		    {d, d, d},
		    e = Sqrt[x^2 + y^2];
		    If[TrueQ[z >= 0],
			a = Sqrt[d + z]; b = e/a,
			b = Sqrt[d - z]; a = e/b];
		    {a, b, $ArcTan[x, y]}],
	    EllipticCylindrical,
		e = ArcCosh[(x + I y)/cs[[4]]];
		{Re[e], Im[e], z},
	    ProlateSpheroidal,
		e = ArcCosh[(z + I Sqrt[x^2 + y^2])/cs[[4]]];
		{Re[e], Im[e], $ArcTan[x, y]},
	    OblateSpheroidal,
		e = ArcCosh[(Sqrt[x^2 + y^2] + I z)/cs[[4]]];
		{Im[e], Re[e], $ArcTan[x, y]},
	    Bipolar,
		e = 2 ArcCoth[x + I y];
		{-Im[e], Re[e], z},
	    Bispherical,
		e = 2 ArcCoth[(z + I Sqrt[x^2 + y^2])/cs[[4]]];
		{-Im[e], Re[e], $ArcTan[x, y]},
	    Toroidal,
		e = 2 ArcCoth[(Sqrt[x^2 + y^2] + I z)/cs[[4]]];
		{Re[e], -Im[e], $ArcTan[x, y]},
	    Conical,
		a = Sqrt[x^2 + y^2 + z^2];
		If[TrueQ[a == 0],
		    {a, 0, 0},
		  (* else *)
		    c = b - cs[[4]]^2;
		    d = b - cs[[5]]^2;
		    e = b /. Solve[x^2 c d + y^2 b d + z^2 b c == 0, b];
		    Prepend[Sqrt[Reverse[Sort[e]]],a] {Sign[x],Sign[y],Sign[z]}
		],
	    ConfocalEllipsoidal,
		a = d - cs[[4]]^2;
		b = d - cs[[5]]^2;
		c = d - cs[[6]]^2;
		Sort[d /. Solve[x^2 b c + y^2 a c + z^2 a b + a b c == 0, d]],
	    ConfocalParaboloidal,
		a = d - cs[[4]]^2;
		b = d - cs[[5]]^2;
		Sort[d /. Solve[x^2 b + y^2 a - 2 z a b - d a b == 0, d]]
	]


    ]

(* === Here are definitions auxiliary operations for vectors and tensors === \
*)

Twirl[t_?ListQ, ord_?ListQ] :=
  Module[{rank = TensorRank[t], d = Dimensions[t], ind, perm},
    (ind = Map[ToExpression["i" <> ToString[#]] &, Range[rank]];
        perm = Map[{ToExpression["i" <> ToString[#]], d[[#]]} &, ord];
        Table[Part[t, Evaluate[Sequence @@ ind]], 
          Evaluate[Sequence @@ perm]] ) /;
      ( Signature[ord] =!= 0 && Length[ord] === rank)
    ]

Twirl[t_?ListQ, ord_?IntegerQ] := 
    Twirl[t, Prepend[Drop[Range[TensorRank[t]], {ord}], ord]]

Twirl[t_?ListQ] := Twirl[t, 2]

GDot[f_?ListQ,g_?ListQ,nf_?IntegerQ,ng_?IntegerQ] := \
Inner[Times,f,Twirl[g,ng],Plus,nf];

GTr[f_?ListQ,n1_?IntegerQ,n2_?IntegerQ] := Tr[Twirl[Twirl[f,n1],n2],Plus,2];

(* === Here are all of the rules for differential operations on vectors === \
*)

ArcLengthFactor[t_Symbol, coordsys_:$CoordinateSystem] :=
    Module[{pt, cs = $ExpandCoordSys[coordsys]},
	ArcLengthFactor[pt, t, cs] /;
	    (If[cs =!= $Failed, pt = List @@ Take[cs, 3]]; cs =!= $Failed)]

ArcLengthFactor[pt_?$VecQ, t_Symbol, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	Sqrt[(ScaleFactors[pt, cs]^2) . (Dt[pt, t]^2)] /; (cs =!= $Failed)]

JacobianDeterminant[arg_:$CoordinateSystem] :=
    Module[{pt, cs},
	JacobianDeterminant[pt, cs] /; (If[$VecQ[arg],
					cs = $CoordinateSystem;
					pt = arg,
					cs = $ExpandCoordSys[arg];
					If[cs =!= $Failed,
					    pt = List @@ Take[cs, 3]]];
					cs =!= $Failed)];

JacobianDeterminant[pt_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	(Times @@ ScaleFactors[pt, cs]) /; (cs =!= $Failed)]

JacobianMatrix[arg_:$CoordinateSystem] :=
    Module[{pt, cs},
	JacobianMatrix[pt, cs] /; (If[$VecQ[arg],
					cs = $CoordinateSystem;
					pt = arg,
					cs = $ExpandCoordSys[arg];
					If[cs =!= $Failed,
					    pt = List @@ Take[cs, 3]]];
					cs =!= $Failed)];

JacobianMatrix[pt_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], uvw},
	(uvw = List @@ Take[cs, 3];
	Outer[D, $CTToCart[uvw, cs], uvw] /. $DAbsSign /.
		$PointRule[cs, pt]) /; (cs =!= $Failed)]

Hessian3Tensor[arg_:$CoordinateSystem] :=
    Module[{pt, cs},
	Hessian3Tensor[pt, cs] /; (If[$VecQ[arg],
					cs = $CoordinateSystem;
					pt = arg,
					cs = $ExpandCoordSys[arg];
					If[cs =!= $Failed,
					    pt = List @@ Take[cs, 3]]];
					cs =!= $Failed)];

Hessian3Tensor[pt_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], sf, q, dq, uvw},
	(uvw = List @@ Take[cs, 3];
         sf = ScaleFactors[cs];
         q = Inner[ Divide, JacobianMatrix[uvw] , sf, List]; 
         dq = Inner[ Divide, Outer[ D, Transpose[ q], uvw ], sf, List ] /. $DAbsSign /. $PointRule[cs, pt]; 
         Transpose[ Simplify[ Inner[Times, dq, q, Plus, 2]]] ) /; (cs =!= $Failed)]

(* change of coordinates for fields *)

FieldFromCartesian[f_?$ScaQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	f /; (cs =!= $Failed)] 

FieldFromCartesian[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
         uvw = List @@ Take[cs, 3];
         sf = ScaleFactors[cs];
         q = Inner[ Divide, JacobianMatrix[uvw] , sf, List];
         Transpose[q] . f /; (cs =!= $Failed)] 

FieldFromCartesian[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
         uvw = List @@ Take[cs, 3];
         sf = ScaleFactors[cs];
         q = Inner[ Divide, JacobianMatrix[uvw] , sf, List];
         Transpose[q] . f . q /; (cs =!= $Failed)] 

FieldToCartesian[f_?$ScaQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	f /; (cs =!= $Failed)] 

FieldToCartesian[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
         uvw = List @@ Take[cs, 3];
         sf = ScaleFactors[cs];
         q = Inner[ Divide, JacobianMatrix[uvw] , sf, List];
         q . f /; (cs =!= $Failed)] 

FieldToCartesian[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
         uvw = List @@ Take[cs, 3];
         sf = ScaleFactors[cs];
         q = Inner[ Divide, JacobianMatrix[uvw] , sf, List];
         q . f . Transpose[q] /; (cs =!= $Failed)] 

(* Nabla is an operator which is equivalent to PreNable and in cartesian 
coordinate applies to the following indices: 
Nabla_a f_b = [f_b,a] e_a x e_b *)

PreNabla[f_?$ScaQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]}, Print[cs]; Print[ ScaleFactors[cs] ];
	Outer[D, {f}, List @@ Take[cs, 3]][[1]]/ScaleFactors[cs] /;
		(cs =!= $Failed)] 

PostNabla[f_?$ScaQ, coordsys_:$CoordinateSystem] :=
    PreNabla[f, coordsys]

Nabla[f_?$ScaQ, coordsys_:$CoordinateSystem] :=
    PreNabla[f, coordsys]

PreNabla[f_?$VecQ, coordsys_:$CoordinateSystem] :=  
    Module[{cs = $ExpandCoordSys[coordsys],sf,ht},
        (sf=ScaleFactors[cs]; ht=Hessian3Tensor[cs];
	Transpose[ Inner[ Divide, Outer[D, f, List @@ Take[cs, 3]], sf,List ]  ] ) \
/;
		(cs =!= $Failed)] 

PostNabla[f_?$VecQ, coordsys_:$CoordinateSystem] :=  
    Transpose[PreNabla[f, coordsys]]

Nabla[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    PreNabla[f, coordsys]

PreNabla[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys],sf,ht},
       (sf=ScaleFactors[cs]; ht=Hessian3Tensor[cs];
	Twirl[ Inner[ Divide, Outer[D, f, List @@ Take[cs, 3]], sf,List ] ,3 ]  ) /;
		(cs =!= $Failed)] 

PostNabla[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    Twirl[PreNabla[f, coordsys],{2, 3, 1}]

Nabla[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    PreNabla[f, coordsys]

(* Grad of Scalar is exactly the same as Nabla *)
PreGrad[f_?$ScaQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	Outer[D, {f}, List @@ Take[cs, 3]][[1]]/ScaleFactors[cs] /;
		(cs =!= $Failed)] 

PostGrad[f_?$ScaQ, coordsys_:$CoordinateSystem] :=
    PreGrad[f, coordsys]

Grad[f_?$ScaQ, coordsys_:$CoordinateSystem] :=
    PostGrad[f, coordsys]

PreGrad[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys],sf,ht},
        (sf=ScaleFactors[cs]; ht=Hessian3Tensor[cs];
	Nabla[f,coordsys]   +   GDot[ ht, f, 2, 1] ) /;
		(cs =!= $Failed)] 

PostGrad[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    Transpose[PreGrad[f, coordsys]]

Grad[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    PostGrad[f, coordsys]

PreGrad[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys],sf,ht},
       (sf=ScaleFactors[cs]; ht=Hessian3Tensor[cs];
	Nabla[f,coordsys] 
        +  GDot [ht, f, 2, 1]  
        +  Twirl[ GDot [ht, f, 2, 2], {1,3,2} ]  ) /;
		(cs =!= $Failed)] 

PostGrad[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    Twirl[PreGrad[f, coordsys],{2, 3, 1}]

Grad[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    PostGrad[f, coordsys]

PreDiv[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    GTr[PreGrad[f,coordsys],1,2]

PostDiv[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    PreDiv[f,coordsys]

Div[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    PreDiv[f,coordsys]

PreDiv[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    GTr[PreGrad[f,coordsys], 1,2 ] 

PostDiv[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    GTr[PostGrad[f,coordsys], 2,3 ] 

Div[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    PostDiv[f, coordsys]

Curl[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], sf,ht,epst,df,br},
	( sf=ScaleFactors[cs]; ht=Hessian3Tensor[cs];
          epst= Table[ Signature[{i, j, k}], {i, 3}, {j, 3}, {k, 3}];
          df=Nabla[ f, coordsys];
          br=df - GDot[ht,f,3,1];
        GTr[ GDot[ epst, br,1,1] ,1,3] 
                ) /; (cs =!= $Failed)]

Curl[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], sf,ht,epst,df,br},
	( sf=ScaleFactors[cs]; ht=Hessian3Tensor[cs];
          epst= Table[ Signature[{i, j, k}], {i, 3}, {j, 3}, {k, 3}];
          df=Nabla[ f, coordsys];
          br=df + GDot[ht,f,2,1] + Twirl[ GDot[ht,f,2,2], {1,3,2} ];
        GTr[ GDot[epst,br,1,1], 1,3 ] 
                ) /; (cs =!= $Failed)]

Laplacian[f_, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	Div[Grad[f, coordsys], coordsys] /; (cs =!= $Failed)]

Biharmonic[f_, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	Laplacian[Laplacian[f, coordsys], coordsys] /; (cs =!= $Failed)]

Inc[f_?$TenQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys]},
	Simplify[Transpose[Curl[Transpose[Curl[f, coordsys]]]]] /; (cs =!= $Failed)]


IntegrateGrad[nf_?$VecQ, coordsys_:$CoordinateSystem] := 
  Module[{cs = $ExpandCoordSys[coordsys], sf, ht}, 
  (
    
    thef[v_] := f[v] ; 
    C2frule[v_] =  C[1] ->  thef[v] ; 

    var = cs[[1]];
    func = thef[var];

    eq = Thread[ Grad[func][[1]] == nf [[1]]];
    sol = DSolve[ eq, thef[var], var , GeneratedParameters -> C][[1]];

    var = cs[[2]];
    func = func /. sol /. C2frule[var] ; 
    eq = Thread[ Grad[func][[2]] == nf [[2]]];

    sol = DSolve[ eq, thef[var], var , GeneratedParameters -> C][[1]];

    var = cs[[3]];
    func = func /. sol /. C2frule[var];
    eq = Thread[ Grad[func][[3]] == nf [[3]]];

    sol = DSolve[ eq, thef[var], var , GeneratedParameters -> C][[1]];

    func = func /. sol ; func 


  )
      
      /; (cs =!= $Failed  || Curl[nf] =!= {0,0,0})]

IntegrateGrad[nf_?$TenQ, coordsys_:$CoordinateSystem] := 
  Module[{cs = $ExpandCoordSys[coordsys]}, 
  ( 
    
    thef[v_] := {f1[v], f2[v], f3[v]} ; 
    C2frule[v_] = Thread[ Table[ C[i], {i, 3}] ->  thef[v] ]; 

    var = cs[[1]];
    func = thef[var];

    eq = Thread[ Grad[func][[All, 1]] == nf [[All, 1]]];
    sol = DSolve[ eq, thef[var], var , GeneratedParameters -> C][[1]];

    var = cs[[2]];
    func = func /. sol /. C2frule[var] ; 
    eq = Thread[ Grad[func][[All, 2]] == nf [[All, 2]]];

    sol = DSolve[ eq, thef[var], var , GeneratedParameters -> C][[1]];

    var = cs[[3]];
    func = func /. sol /. C2frule[var];
    eq = Thread[ Grad[func][[All, 3]] == nf [[All, 3]]];

    sol = DSolve[ eq, thef[var], var , GeneratedParameters -> C][[1]];

    func = func /. sol ; func 

  )
      /; (cs =!= $Failed   ||  Curl[Transpose[nf]] =!= {{0,0,0},{0,0,0},{0,0,0}} )]





(* no newcommand after this line, to escape the default coordinate system *)

SetCoordinates[Cartesian[ ]];  (* default coordinate system *)


End[ ];		(* "Calculus`Tensor2Analysis`Private`" *)

Attributes[DotProduct] = {ReadProtected};
Attributes[CrossProduct] = {ReadProtected};
Attributes[ScalarTripleProduct] = {ReadProtected};
Attributes[SetCoordinates] = {ReadProtected};
Attributes[CoordinateSystem] = {ReadProtected};
Attributes[Coordinates] = {ReadProtected};
Attributes[Parameters] = {ReadProtected};
Attributes[CoordinateRanges] = {ReadProtected};
Attributes[CoordinatesToCartesian] = {ReadProtected};
Attributes[CoordinatesFromCartesian] = {ReadProtected};
Attributes[ScaleFactors] = {ReadProtected};
Attributes[ArcLengthFactor] = {ReadProtected};
Attributes[FieldToCartesian] = {ReadProtected};
Attributes[FieldFromCartesian] = {ReadProtected};
Attributes[JacobianDeterminant] = {ReadProtected};
Attributes[JacobianMatrix] = {ReadProtected};
Attributes[Hessian3Tensor] = {ReadProtected};
Attributes[Grad] = {ReadProtected};
Attributes[Div] = {ReadProtected};
Attributes[Curl] = {ReadProtected};
Attributes[Laplacian] = {ReadProtected};
Attributes[Biharmonic] = {ReadProtected};
Attributes[Inc] = {ReadProtected};
Attributes[IntegrateGrad] = {ReadProtected};

Protect[Cartesian, Cylindrical, Spherical, ParabolicCylindrical,
    Paraboloidal, EllipticCylindrical, ProlateSpheroidal, OblateSpheroidal,
    Bipolar, Bispherical, Toroidal, Conical, ConfocalEllipsoidal,
    ConfocalParaboloidal];

Protect[DotProduct, CrossProduct, ScalarTripleProduct, SetCoordinates,
    Coordinates, Parameters, CoordinateRanges, CoordinatesToCartesian,
    CoordinatesFromCartesian, FieldFromCartesian, FieldToCartesian, \
ScaleFactors, ArcLengthFactor,
    JacobianDeterminant, JacobianMatrix, Hessian3Tensor, Nabla, Grad, Div, \
Curl, Laplacian,
    Biharmonic, Inc, CoordinateSystem,Twirl,GDot,GTr,IntegrateGrad, $VecQ];



Protect[Eeta, Llambda, Mmu, Nnu, Pphi, Rr, Ttheta, Uu, Vv, Xx, Xxi, Yy, Zz]


EndPackage[ ]

