--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

)abbrev package COORDSYS CoordinateSystems
++ Author: Jim Wen
++ Date Created: 12 March 1990
++ Date Last Updated: 19 June 1990, Clifton J. Williamson
++ Basic Operations: cartesian, polar, cylindrical, spherical, parabolic, elliptic, 
++ parabolicCylindrical, paraboloidal, ellipticCylindrical, prolateSpheroidal,
++ oblateSpheroidal, bipolar, bipolarCylindrical, toroidal, conical
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description: CoordinateSystems provides coordinate transformation functions 
++ for plotting.  Functions in this package return conversion functions 
++ which take points expressed in other coordinate systems and return points 
++ with the corresponding Cartesian coordinates.
 
CoordinateSystems(R): Exports == Implementation where

  R : Join(Field,TranscendentalFunctionCategory,RadicalCategory)
  Pt ==> Point R

  Exports ==> with
    cartesian : Pt -> Pt
      ++ cartesian(pt) returns the Cartesian coordinates of point pt.
    polar: Pt -> Pt
      ++ polar(pt) transforms pt from polar coordinates to Cartesian 
      ++ coordinates: the function produced will map the point \spad{(r,theta)}
      ++ to \spad{x = r * cos(theta)} , \spad{y = r * sin(theta)}.
    cylindrical: Pt -> Pt
      ++ cylindrical(pt) transforms pt from polar coordinates to Cartesian 
      ++ coordinates: the function produced will map the point \spad{(r,theta,z)}
      ++ to \spad{x = r * cos(theta)}, \spad{y = r * sin(theta)}, \spad{z}.
    spherical: Pt -> Pt
      ++ spherical(pt) transforms pt from spherical coordinates to Cartesian 
      ++ coordinates: the function produced will map the point \spad{(r,theta,phi)}
      ++ to \spad{x = r*sin(phi)*cos(theta)}, \spad{y = r*sin(phi)*sin(theta)},
      ++ \spad{z = r*cos(phi)}.
    parabolic: Pt -> Pt
      ++ parabolic(pt) transforms pt from parabolic coordinates to Cartesian 
      ++ coordinates: the function produced will map the point \spad{(u,v)} to
      ++ \spad{x = 1/2*(u**2 - v**2)}, \spad{y = u*v}.
    parabolicCylindrical: Pt -> Pt
      ++ parabolicCylindrical(pt) transforms pt from parabolic cylindrical 
      ++ coordinates to Cartesian coordinates: the function produced will 
      ++ map the point \spad{(u,v,z)} to \spad{x = 1/2*(u**2 - v**2)}, 
      ++ \spad{y = u*v}, \spad{z}.
    paraboloidal: Pt -> Pt
      ++ paraboloidal(pt) transforms pt from paraboloidal coordinates to 
      ++ Cartesian coordinates: the function produced will map the point 
      ++ \spad{(u,v,phi)} to \spad{x = u*v*cos(phi)}, \spad{y = u*v*sin(phi)},
      ++ \spad{z = 1/2 * (u**2 - v**2)}.
    elliptic: R -> (Pt -> Pt)
      ++ elliptic(a) transforms from elliptic coordinates to Cartesian 
      ++ coordinates: \spad{elliptic(a)} is a function which will map the 
      ++ point \spad{(u,v)} to \spad{x = a*cosh(u)*cos(v)}, \spad{y = a*sinh(u)*sin(v)}.
    ellipticCylindrical: R -> (Pt -> Pt)
      ++ ellipticCylindrical(a) transforms from elliptic cylindrical coordinates 
      ++ to Cartesian coordinates: \spad{ellipticCylindrical(a)} is a function
      ++ which will map the point \spad{(u,v,z)} to \spad{x = a*cosh(u)*cos(v)},
      ++ \spad{y = a*sinh(u)*sin(v)}, \spad{z}.
    prolateSpheroidal: R -> (Pt -> Pt)
      ++ prolateSpheroidal(a) transforms from prolate spheroidal coordinates to 
      ++ Cartesian coordinates: \spad{prolateSpheroidal(a)} is a function 
      ++ which will map the point \spad{(xi,eta,phi)} to 
      ++ \spad{x = a*sinh(xi)*sin(eta)*cos(phi)}, \spad{y = a*sinh(xi)*sin(eta)*sin(phi)}, 
      ++ \spad{z = a*cosh(xi)*cos(eta)}.
    oblateSpheroidal: R -> (Pt -> Pt)
      ++ oblateSpheroidal(a) transforms from oblate spheroidal coordinates to 
      ++ Cartesian coordinates: \spad{oblateSpheroidal(a)} is a function which
      ++ will map the point \spad{(xi,eta,phi)} to \spad{x = a*sinh(xi)*sin(eta)*cos(phi)},
      ++ \spad{y = a*sinh(xi)*sin(eta)*sin(phi)}, \spad{z = a*cosh(xi)*cos(eta)}.
    bipolar: R -> (Pt -> Pt)
      ++ bipolar(a) transforms from bipolar coordinates to Cartesian coordinates:
      ++ \spad{bipolar(a)} is a function which will map the point \spad{(u,v)} to
      ++ \spad{x = a*sinh(v)/(cosh(v)-cos(u))}, \spad{y = a*sin(u)/(cosh(v)-cos(u))}.
    bipolarCylindrical: R -> (Pt -> Pt)
      ++ bipolarCylindrical(a) transforms from bipolar cylindrical coordinates 
      ++ to Cartesian coordinates: \spad{bipolarCylindrical(a)} is a function which 
      ++ will map the point \spad{(u,v,z)} to \spad{x = a*sinh(v)/(cosh(v)-cos(u))},
      ++ \spad{y = a*sin(u)/(cosh(v)-cos(u))}, \spad{z}.
    toroidal: R -> (Pt -> Pt)
      ++ toroidal(a) transforms from toroidal coordinates to Cartesian 
      ++ coordinates: \spad{toroidal(a)} is a function which will map the point 
      ++ \spad{(u,v,phi)} to \spad{x = a*sinh(v)*cos(phi)/(cosh(v)-cos(u))},
      ++ \spad{y = a*sinh(v)*sin(phi)/(cosh(v)-cos(u))}, \spad{z = a*sin(u)/(cosh(v)-cos(u))}.
    conical: (R,R) -> (Pt -> Pt)
      ++ conical(a,b) transforms from conical coordinates to Cartesian coordinates:
      ++ \spad{conical(a,b)} is a function which will map the point \spad{(lambda,mu,nu)} to
      ++ \spad{x = lambda*mu*nu/(a*b)},
      ++ \spad{y = lambda/a*sqrt((mu**2-a**2)*(nu**2-a**2)/(a**2-b**2))},
      ++ \spad{z = lambda/b*sqrt((mu**2-b**2)*(nu**2-b**2)/(b**2-a**2))}.

  Implementation ==> add

    cartesian pt ==
      -- we just want to interpret the cartesian coordinates
      -- from the first N elements of the point - so the
      -- identity function will do
      pt

    polar pt0 ==
      pt := copy pt0
      r := elt(pt0,1); theta := elt(pt0,2)
      pt.1 := r * cos(theta); pt.2 := r * sin(theta)
      pt

    cylindrical pt0 == polar pt0 
    -- apply polar transformation to first 2 coordinates

    spherical pt0 ==
      pt := copy pt0
      r := elt(pt0,1); theta := elt(pt0,2); phi := elt(pt0,3)
      pt.1 := r * sin(phi) * cos(theta); pt.2 := r * sin(phi) * sin(theta)
      pt.3 := r * cos(phi)
      pt

    parabolic pt0 ==
      pt := copy pt0
      u := elt(pt0,1); v := elt(pt0,2)
      pt.1 := (u*u - v*v)/(2::R) ; pt.2 := u*v
      pt

    parabolicCylindrical pt0 == parabolic pt0
    -- apply parabolic transformation to first 2 coordinates

    paraboloidal pt0 ==
      pt := copy pt0
      u := elt(pt0,1); v := elt(pt0,2); phi := elt(pt0,3)
      pt.1 := u*v*cos(phi); pt.2 := u*v*sin(phi); pt.3 := (u*u - v*v)/(2::R)
      pt

    elliptic a ==
      pt := copy(#1)
      u := elt(#1,1); v := elt(#1,2)
      pt.1 := a*cosh(u)*cos(v); pt.2 := a*sinh(u)*sin(v)
      pt

    ellipticCylindrical a == elliptic a
    -- apply elliptic transformation to first 2 coordinates

    prolateSpheroidal a ==
      pt := copy(#1)
      xi := elt(#1,1); eta := elt(#1,2); phi := elt(#1,3)
      pt.1 := a*sinh(xi)*sin(eta)*cos(phi)
      pt.2 := a*sinh(xi)*sin(eta)*sin(phi)
      pt.3 := a*cosh(xi)*cos(eta)
      pt

    oblateSpheroidal a ==
      pt := copy(#1)
      xi := elt(#1,1); eta := elt(#1,2); phi := elt(#1,3)
      pt.1 := a*sinh(xi)*sin(eta)*cos(phi)
      pt.2 := a*cosh(xi)*cos(eta)*sin(phi)
      pt.3 := a*sinh(xi)*sin(eta)
      pt

    bipolar a ==
      pt := copy(#1)
      u := elt(#1,1); v := elt(#1,2)
      pt.1 := a*sinh(v)/(cosh(v)-cos(u))
      pt.2 := a*sin(u)/(cosh(v)-cos(u))
      pt

    bipolarCylindrical a == bipolar a
    -- apply bipolar transformation to first 2 coordinates

    toroidal a ==
      pt := copy(#1)
      u := elt(#1,1); v := elt(#1,2); phi := elt(#1,3)
      pt.1 := a*sinh(v)*cos(phi)/(cosh(v)-cos(u))
      pt.2 := a*sinh(v)*sin(phi)/(cosh(v)-cos(u))
      pt.3 := a*sin(u)/(cosh(v)-cos(u))
      pt

    conical(a,b) ==
      pt := copy(#1)
      lambda := elt(#1,1); mu := elt(#1,2); nu := elt(#1,3)
      pt.1 := lambda*mu*nu/(a*b)
      pt.2 := lambda/a*sqrt((mu**2-a**2)*(nu**2-a**2)/(a**2-b**2))
      pt.3 := lambda/b*sqrt((mu**2-b**2)*(nu**2-b**2)/(b**2-a**2))
      pt

