
(* :Name: Calculus`VectorAnalysis` *)

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

(* :Author: Stephen Wolfram, Bruce K. Sawhill, Jerry B. Keiper *) 

(* :Summary:
	This package performs standard vector differential operations in 
	orthogonal coordinate systems.  Coordinate transformations are 
	also provided.
*)  

(* :Context: Calculus`VectorAnalysis` *)

(* :Package Version: 2.0 *)

(* :Copyright: Copyright 1990  Wolfram Research, Inc.
	Permission is hereby granted to modify and/or make copies of
	this file for any purpose other than direct profit, or as part
	of a commercial product, provided this copyright notice is left
	intact.  Sale, other than for the cost of media, is prohibited.

	Permission is hereby granted to reproduce part or all of
	this file, provided that the source is acknowledged.
*)

(* :History:
	Extensively revised by Jerry B. Keiper, November, 1990.
	Revised by Bruce K. Sawhill, August, 1990.
	Version 1.0 by Stephen Wolfram, 1988.
*)

(* :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:
	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.
*)

(* :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 x, y, z > 0.  The transformation from
	the confocal paraboloidal coordinate system always gives a point
	in the quadrant with x, y > 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.
*)

BeginPackage["Calculus`VectorAnalysis`"]

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 variable names.
SetCoordinates[coordsys[x, y, z]] sets the default coordinate system to be
the specified coordinate system with variable names x, y, and z.  Certain
coordinate systems have parameters associated with them, and these
parameters can be set by specifying the full description of the coordinate
system, as in SetCoordinates[Conical[lambda, mu, nu, a, b]]."

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

Coordinates::usage = "Coordinates[ ] gives a list of the default names of
the coordinate variables in the default coordinate system.
Coordinates[coordsys] gives a list of the default names of 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 the 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 names of 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 names of coordinate variables are used."

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 names of 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 names of coordinate variables are used."

Grad::usage = "Grad[f] gives the gradient of the scalar function f in the
default coordinate system. Grad[f, coordsys] gives the gradient of f in the
coordinate system coordsys."

Div::usage = "Div[f] gives the divergence of the vector-valued function f
in the default coordinate system. Div[f, coordsys] gives the divergence of f
in the coordinate system coordsys."

Curl::usage = "Curl[f] gives the curl of the vector-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 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 function f in the default coordinate system.
Biharmonic[f, coordsys] gives the biharmonic of f in the coordinate system
coordsys."

Cartesian::usage = "Cartesian represents the Cartesian coordinate system
with default variable names.  Cartesian[x, y, z] represents the Cartesian
coordinate system with variable names x, y, and z."

Cylindrical::usage = "Cylindrical represents the cylindrical coordinate
system with default variable names. Cylindrical[r, theta, z] represents the
cylindrical coordinate system with variable names r, theta, and z."

Spherical::usage = "Spherical represents the spherical coordinate system
with default variable names.  Spherical[r, theta, phi] represents the
spherical coordinate system with variable names r, theta, and phi."

ParabolicCylindrical::usage = "ParabolicCylindrical represents the parabolic
cylindrical coordinate system with default variable names.
ParabolicCylindrical[u, v, z] represents the parabolic cylindrical coordinate
system with variable names 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 variable names.  Paraboloidal[u, v, phi] represents
the paraboloidal coordinate system with variable names 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 parameter and default variable
names.  EllipticCylindrical[u, v, z] represents the elliptic cylindrical
coordinate system with variable names u, v, and z and default parameter.
EllipticCylindrical[u, v, z, a] represents the elliptic cylindrical
coordinate system  with variable names 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 parameter and default variable
names.  ProlateSpheroidal[xi, eta, phi] represents the prolate spheroidal
coordinate system with variable names xi, eta, and phi and default parameter.
ProlateSpheroidal[xi, eta, phi, a] represents the prolate spheroidal
coordinate system with variable names 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 2a.  The third variable parameterizes the rotation. *)

OblateSpheroidal::usage = "OblateSpheroidal represents the oblate spheroidal
coordinate system with default parameter and default variable names.
OblateSpheroidal[xi, eta, phi] represents the oblate spheroidal coordinate
system with variable names xi, eta, and phi and default parameter.
OblateSpheroidal[xi, eta, phi, a] represents the oblate spheroidal
coordinate system with variable names 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 2a.  The third variable parameterizes the
    	rotation. *)
         
Bipolar::usage = "Bipolar represents the bipolar coordinate system with 
default parameter and default variable names.  Bipolar[u, v, z] represents
the bipolar coordinate system with variable names u, v, and z and default
parameter.  Bipolar[u, v, z, a] represents the  bipolar coordinate system
with variable names u, v, and z and parameter a."
	(* Bipolar coordinates: These coordinates are built around two
	foci, separated by a distance 2a.  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 parameter and default variable names.
Bispherical[u, v, phi] represents the bispherical coordinate system with
variable names u, v, and phi and default parameter.  Bispherical[u, v, phi, a]
represents the bispherical coordinate system with variable names 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 z. *)

Toroidal::usage = "Toroidal represents the toroidal coordinate system
with default parameter and default variable names.  Toroidal[u, v, phi]
represents the toroidal coordinate system with variable names u, v, and phi
and default parameter.  Toroidal[u, v, phi, a] represents the toroidal
coordinate system with variable names u, v, and 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 2a.  The third variable
	parameterizes this rotation. *)
	
Conical::usage = "Conical represents the conical coordinate system with 
default parameters and default variable names.  Conical[lambda, mu, nu]
represents the conical coordinate system with variable names lambda, mu,
and nu and default parameters.  Conical[lambda, mu, nu, a, b] represents
the conical coordinate system with variable names lambda, mu, and nu and
parameters a and b."
	(* Conical coordinates:  The coordinate surfaces are spheres
	with centers at the origin and radius lambda (lambda constant),
	cones with apexes at the origin and axes along the z-axis
	(mu constant), and cones with apexes at the origin and axes
	along the y-axis (nu constant). *)
	
ConfocalEllipsoidal::usage = "ConfocalEllipsoidal represents the confocal
ellipsoidal coordinate system with default parameters and default variable
names.  ConfocalEllipsoidal[lambda, mu, nu] represents the confocal ellipsoidal
coordinate system with variable names lambda, mu, and nu and default
parameters.  ConfocalEllipsoidal[lambda, mu, nu, a, b, c] represents the
confocal ellipsoidal coordinate system with variable names lambda, mu, and
nu and parameters a, b, and c."
	(* Confocal ellipsoidal coordinates: The coordinate surfaces are
	ellipsoids (lambda constant), hyperboloids of one sheet (mu
	constant), and hyperboloids of two sheets (nu constant). *)

ConfocalParaboloidal::usage = "ConfocalParaboloidal represents the confocal
paraboloidal coordinate system with default parameters and default variable
names.  ConfocalParaboloidal[lambda, mu, nu] represents the confocal
paraboloidal coordinate system with variable names lambda, mu, and nu and
default parameters.  ConfocalParaboloidal[lambda, mu, nu, a, b] represents
the confocal paraboloidal coordinate system with variable names 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 (lambda constant), hyperbolic
	paraboloids (mu constant), and elliptic paraboloids extending in
	the direction of the positive z-axis (nu constant). *)


(* default coordinates *)

x::usage = y::usage =
"{x,y,z} are the default Cartesian coordinates."

z::usage =
"{x,y,z} are the default Cartesian coordinates; {r,theta,z} are the default
Cylindrical coordinates; {u,v,z} are the default coordinates for the
ParabolicCylindrical, EllipticCylindrical, and Bipolar coordinate systems."

r::usage = theta::usage =
"{r,theta,z} are the default Cylindical coordinates; {r,theta,phi} are the
default Spherical coordinates."

phi::usage =
"{r,theta,phi} are the default Spherical coordinates; {u,v,phi} are the default
coordinates for the Paraboloidal, Bispherical, and Toroidal coordinate systems;
{xi,eta,phi} are the default coordinates for the Prolate Spheroidal and
Oblate Spheroidal coordinate systems."

u::usage = v::usage =
"{u,v,z} are the default coordinates for the ParabolicCylindrical,
EllipticCylindrical, and Bipolar coordinate systems; {u,v,phi} are the default
coordinates for the Paraboloidal,  Bispherical, and Toroidal coordinate
systems."

xi::usage = eta::usage =
"{xi,eta,phi} are the default coordinates for the ProlateSpheroidal and
OblateSpheroidal coordinate systems."

lambda::usage = mu::usage = nu::usage =
"{lambda,mu,nu} are the default coordinates for the Conical,
ConfocalEllipsoidal, and ConfocalParaboloidal coordinate systems."


Begin["Calculus`VectorAnalysis`Private`"]

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

$VecQ[v_] := VectorQ[v] && Length[v] == 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] ^= {x, y, z}
Coordinates[Cylindrical] ^= {r, theta, z}
Coordinates[Spherical] ^= {r, theta, phi} 
Coordinates[ParabolicCylindrical] ^= {u, v, z}
Coordinates[Paraboloidal] ^= {u, v, phi} 
Coordinates[EllipticCylindrical] ^= {u, v, z}
Coordinates[ProlateSpheroidal] ^= {xi, eta, phi} 
Coordinates[OblateSpheroidal] ^= {eta, xi, phi} 
Coordinates[Bipolar] ^= {u, v, z}
Coordinates[Bispherical] ^= {u, v, phi}
Coordinates[Toroidal] ^= {v, u, phi} 
Coordinates[Conical] ^= {lambda, mu, nu} 
Coordinates[ConfocalEllipsoidal] ^= {lambda, mu, nu} 
Coordinates[ConfocalParaboloidal] ^= {lambda, mu, nu}

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 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)]

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

Div[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], sf, psf},
	(sf = ScaleFactors[cs];
	psf = Times @@ sf;
	(Inner[D, f psf/sf, List @@ Take[cs, 3]] /. $DAbsSign)/psf) /;
		(cs =!= $Failed)]

Curl[f_?$VecQ, coordsys_:$CoordinateSystem] :=
    Module[{cs = $ExpandCoordSys[coordsys], sf},
	(sf = ScaleFactors[cs];
	{(D[sf[[3]] f[[3]], cs[[2]]] - D[sf[[2]] f[[2]], cs[[3]]])/
		(sf[[2]] sf[[3]]),
	 (D[sf[[1]] f[[1]], cs[[3]]] - D[sf[[3]] f[[3]], cs[[1]]])/
		(sf[[1]] sf[[3]]),
	 (D[sf[[2]] f[[2]], cs[[1]]] - D[sf[[1]] f[[1]], cs[[2]]])/
		(sf[[1]] sf[[2]])} /. $DAbsSign) /; (cs =!= $Failed)]

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

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

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

End[ ];		(* "Calculus`VectorAnalysis`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[JacobianDeterminant] = {ReadProtected};
Attributes[JacobianMatrix] = {ReadProtected};
Attributes[Grad] = {ReadProtected};
Attributes[Div] = {ReadProtected};
Attributes[Curl] = {ReadProtected};
Attributes[Laplacian] = {ReadProtected};
Attributes[Biharmonic] = {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, ScaleFactors, ArcLengthFactor,
    JacobianDeterminant, JacobianMatrix, Grad, Div, Curl, Laplacian,
    Biharmonic, CoordinateSystem];

EndPackage[ ]
