#pragma ident "$Id: Position.hpp 693 2007-07-19 14:56:11Z pben $"



/**
 * @file Position.hpp
 * class gpstk::Position encapsulates 3-D positions, including geographic positions,
 *    expressed as geodetic (with respect to an ellipsoidal geoid), geocentric or
 *    Earth-centered, Earth-fixed (cartesian) coordinates, as well as ordinary
 *    positions defined by spherical or cartesian coordinates. Position inherits
 *    from class Triple.
 */

#ifndef GPSTK_POSITION_HPP
#define GPSTK_POSITION_HPP

//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 2.1 of the License, or
//  any later version.
//
//  The GPSTk is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GPSTk; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//  
//  Copyright 2004, The University of Texas at Austin
//
//============================================================================

#include "Exception.hpp"
#include "StringUtils.hpp"
#include "DayTime.hpp"  // for FormatException
#include "Triple.hpp"
#include "GeoidModel.hpp"
#include "Xvt.hpp"

namespace gpstk
{
      /** @addtogroup geodeticgroup */
      //@{

      // The following forward declaration of Position and range are the only
      // way I can get range to be a member of namespace gpstk.
   class Position;
   double range(const Position& A, const Position& B) throw(GeometryException);
   
      /**
       * A position representation class for common 3D geographic position formats,
       * including geodetic (geodetic latitude, longitude, and height above the geoid)
       * geocentric (geocentric latitude, longitude, and radius from Earth's center),
       * cartesian (Earth-centered, Earth-fixed) and spherical (theta,phi,radius).
       *
       * Internally, the representation of Position consists of three coordinate
       * values (double), two doubles from a geoid model (see below, storing these
       * doubles is preferred over adding GeoidModel to calling arguments everywhere),
       * a flag of type 'enum CoordinateSystem' giving the coordinate system, and a
       * tolerance for use in comparing Positions. Class Position inherits from class
       * Triple, which is how the coordinate values are stored (Triple actually uses
       * std::valarray<double> of length 3). It is important to note that
       * Triple:: routines are properly used by Positions ONLY in the Cartesian
       * coordinate system.
       *
       * Only geodetic coordinates depend on a geoid, and then (for an ellipsoidal
       * geoid) only on the semi-major axis of the Earth and the square of its
       * eccentricity. Input of this geoid information (usually a pointer to a
       * GeoidModel) is required by functions involving constructors of, or
       * transformation to or from, Geodetic coordinates. However since a default
       * is supplied (WGS84), the user need never deal with geiods unless desired.
       * In fact, if the geodetic coordinate system is avoided, the Position class
       * can be interpreted simply as 3D vectors in any context, particularly since
       * the class inherits from Triple, which includes many vector manipulation
       * routines (although the Triple:: routines assume Cartesian coordinates).
       * Even the requirement that lengths (radius, height and the cartesian
       * coordinates) have units of meters is required only if geodetic coordinates
       * are used (because the semi-major axis in GeoidModel is in meters);
       * without using Geodetic one could apply the class using any units for
       * length as long as setTolerance() is called appropriately.
       *
       * Position relies on a series of fundamental routines to transform from
       * one coordinate system to another, these include, for example
       * void Position::convertGeodeticToCartesian(const Triple& llh, Triple& xyz,
       *    const double A, const double eccSq);
       * void Position::convertSphericalToCartesian(const Triple& tpr, Triple& xyz);
       * These functions use Triple in the calling arguments.
       *
       * Position will throw exceptions (gpstk::GeometryException) on bad input
       * (e.g. negative radius or latitude > 90 degrees); otherwise the class
       * attempts to handle all points, even the pole and the origin, consistently
       * and without throwing exceptions.
       * At or very near the poles, the transformation routines will set
       * latitude = +/-90 degrees, which is theta = 0 or 180, and (arbitrarily)
       * longitude = 0. At or very near the origin, the transformation routines
       * will set latitude = 0, which is theta = 90, and (arbitrarily) longitude = 0;
       * radius will be set to zero and geodetic height will be set to
       * -radius(Earth) (= -6378137.0 in WGS84). The tolerance used in testing
       * 'at or near the pole or origin' is radius < POSITION_TOLERANCE/5.
       * Note that this implies that a Position that is very near the origin may
       * be SET to the exact origin by the transformation routines, and that
       * thereby information about direction (e.g. latitude and longitude)
       * may be LOST. The user is warned to be very careful when working
       * near either the pole or the origin.
       *
       * Position includes setToString() and printf() functions similar to those
       * in gpstk::DayTime; this allows flexible and powerful I/O of Position to
       * strings and streams.
       *
       * @sa positiontest.cpp for examples.
       */
   class Position : public Triple
   {
   public:
      // ----------- Part  1: coordinate systems --------------------------------
      //
         /// The coordinate systems supported by Position
      enum CoordinateSystem
      {
         Unknown=0,  ///< unknown coordinate system
         Geodetic,   ///< geodetic latitude, longitude, and height above geoid
         Geocentric, ///< geocentric (regular spherical coordinates)
         Cartesian,  ///< cartesian (Earth-centered, Earth-fixed)
         Spherical   ///< spherical coordinates (theta,phi,radius)
      };

         /// return string giving name of coordinate system
      std::string getSystemName()
         throw();

      // ----------- Part  2: member functions: tolerance -----------------------
      //
         /// One millimeter tolerance.
      static const double ONE_MM_TOLERANCE;
         /// One centimeter tolerance.
      static const double ONE_CM_TOLERANCE;
         /// One micron tolerance.
      static const double ONE_UM_TOLERANCE;
      
         /// Default tolerance for time equality in days.
      static double POSITION_TOLERANCE;

         /// Changes the POSITION_TOLERANCE for all Position objects
      static double setPositionTolerance(const double tol)
         { POSITION_TOLERANCE = tol;  return POSITION_TOLERANCE; }

         /// Returns the current POSITION_TOLERANCE.
      static double getPositionTolerance()
         { return POSITION_TOLERANCE; }
   
         /**
          * Sets the tolerance for output and comparisons, for this object only.
          * See the constants in this file (e.g. ONE_MM_TOLERANCE)
          * for some easy to use tolerance values.
          * @param tol Tolerance in meters to be used by comparison operators.
          * @sa Position-Specific Definitions
          */
      Position& setTolerance(const double tol)
         throw();

      // ----------- Part  3: member functions: constructors --------------------
      //
         /**
          * Default constructor.
          * Initializes to zero, Unknown coordinates
          */
      Position()
         throw();

         /**
          * Explicit constructor. Coordinate system may be specified on input,
          * but defaults to Cartesian. Pointer to GeoidModel may be specified,
          * but default is NULL (in which case WGS84 values will be used).
          * @param a first coordinate [ X(m), or latitude (degrees N) ]
          * @param b second coordinate [ Y(m), or longitude (degrees E) ]
          * @param c third coordinate [ Z, height above ellipsoid or radius, in m ]
          * @param s coordinate system
          * @param geoid pointer to GeoidModel
          * @throw GeometryException on invalid input.
          */
      Position(const double& a,
               const double& b,
               const double& c,
               CoordinateSystem s = Cartesian,
               GeoidModel *geoid = NULL)
         throw(GeometryException);

         /**
          * Explicit constructor. Coordinate system may be specified on input,
          * but defaults to Cartesian. Pointer to GeoidModel may be specified,
          * but default is NULL (in which case WGS84 values will be used).
          * @param ABC double array[3] coordinate values
          * @param s CoordinateSystem
          * @param geoid pointer to GeoidModel
          * @throw GeometryException on invalid input.
          */
      Position(const double ABC[3],
               CoordinateSystem s = Cartesian,
               GeoidModel *geoid = NULL)
         throw(GeometryException);

         /**
          * Explicit constructor. Coordinate system may be specified on input,
          * but defaults to Cartesian. Pointer to GeoidModel may be specified,
          * but default is NULL (in which case WGS84 values will be used).
          * @param ABC coordinate values
          * @param s CoordinateSystem
          * @param geoid pointer to GeoidModel
          * @throw GeometryException on invalid input.
          */
      Position(const Triple& ABC,
               CoordinateSystem s = Cartesian,
               GeoidModel *geoid = NULL)
         throw(GeometryException);

         /**
          * Explicit constructor from Xvt. The coordinate system is Cartesian,
          * and the velocity and time information in the input is ignored.
          * @param xvt Input Xvt object, xvt.x contains the Cartesian coordinates
          */
      Position(const Xvt& xvt)
         throw();

         /// Destructor.
      ~Position()
         throw()
         {}

      // ----------- Part  4: member functions: arithmetic ----------------------
      //
         /** Subtract a Position from this Position. Perform the subtraction in
          * Cartesian coordinates, but return this Position to the system it
          * had originally.
          * @param right Position to subtract from this one.
          * @return new Position, in the original system.
          */
      Position& operator-=(const Position& right)
         throw();

         /** Add a Position to this Position. Perform the addition in
          * Cartesian coordinates, but return this Position to the system it
          * had originally.
          * @param right Position to add to this one.
          * @return new Position, in the original system.
          */
      Position& operator+=(const Position& right)
         throw();

         /**
          * Difference two Positions, returning result as a Position in Cartesian
          * coordinates, the only system in which a position difference makes sense.
          * @param right Position to subtract from this one.
          * @return difference as Position.
          */
      friend Position operator-(const Position& left,
                                      const Position& right)
         throw();

         /**
          * Add two Positions, returning result as a Position in Cartesian
          * coordinates, the only system in which a position sum makes sense.
          * @param right Position to add to this one.
          * @return The new Position.
          */
      friend Position operator+(const Position& left,
                                      const Position& right)
         throw();

         /** Multiply a Position by a double scalar on the left.
          * @param right Position to be multiplied by the scalar
          * @param scale the (double) scalar
          * @return The new Position.
          */
      friend Position operator*(const double& scale,
                                const Position& right)
         {
            Position tmp(right);
            tmp.theArray *= scale;
            return tmp;
         }

         /** Multiply a Position by a double scalar on the right.
          * @param left Position to be multiplied by the scalar
          * @param scale the (double) scalar
          * @return The new Position.
          */
      friend Position operator*(const Position& left,
                                const double& scale)
         {
            return operator*(scale, left);
         }

         /** Multiply a Position by an integer scalar on the left.
          * @param right Position to be multiplied by the scalar
          * @param scale the (int) scalar
          * @return The new Position.
          */
      friend Position operator*(const int& scale,
                                const Position& right)
         {
            return operator*(double(scale), right);
         }

         /** Multiply a Position by an integer scalar on the right.
          * @param left Position to be multiplied by the scalar
          * @param scale the (int) scalar
          * @return The new Position.
          */
      friend Position operator*(const Position& left,
                                const int& scale)
         {
            return operator*(double(scale), left);
         }

      // ----------- Part  5: member functions: comparisons ---------------------
      //
         /// Equality operator. Return true if range between this Position and
         /// the input Position is less than tolerance. Return false if geoid
         /// values differ.
         /// @param right Position to be compared to this Position
      bool operator==(const Position &right) const
         throw();

         /// Inequality operator. Return true if range between this Position and
         /// the input Position is greater than tolerance. Return true if geoid
         /// values differ.
         /// @param right Position to be compared to this Position
      bool operator!=(const Position &right) const
         throw();

      // ----------- Part  6: member functions: coordinate transformations ------
      //
         /**
          * Transform coordinate system. Does nothing if sys already matches the
          * current value of member CoordinateSystem 'system'.
          * @param sys CoordinateSystem into which this Position is transformed.
          */
      Position transformTo(CoordinateSystem sys)
         throw();
  
         /// Convert to geodetic coordinates (does nothing if
         /// system == Geodetic already).
      Position asGeodetic()
         throw()
      { transformTo(Geodetic); return *this; }

         /// Convert to another geoid, then to geodetic coordinates.
         /// @return a reference to this.
         /// @throw GeometryException if input is NULL.
      Position asGeodetic(GeoidModel *geoid)
         throw(GeometryException)
      {
         try { setGeoidModel(geoid); }
         catch(GeometryException& ge) { GPSTK_RETHROW(ge); }
         transformTo(Geodetic);
         return *this;
      }

         /// Convert to cartesian coordinates (does nothing if
         /// system == Cartesian already).
      Position asECEF()
         throw()
      { transformTo(Cartesian); return *this; }


      // ----------- Part  7: member functions: get -----------------------------
      // 
      // These routines retrieve coordinate values in all coordinate systems.
      //
         /// return X coordinate (meters)
      double X() const
         throw();

         /// return Y coordinate (meters)
      double Y() const
         throw();

         /// return Z coordinate (meters)
      double Z() const
         throw();

         /// return geodetic latitude (degrees North).
      double geodeticLatitude() const
         throw();

         /// return geocentric latitude (degrees North);
         /// equal to 90 degress - theta in regular spherical coordinates.
      double geocentricLatitude() const
         throw();

         /// return spherical coordinate theta in degrees
      double theta() const
         throw();

         /// return spherical coordinate phi in degrees
      double phi() const
         throw();

         /// return longitude (degrees East);
         /// equal to phi in regular spherical coordinates.
      double longitude() const
         throw();

         /// return distance from the center of Earth (meters),
         /// Same as radius in spherical coordinates.
      double radius() const
         throw();

         /// return height above ellipsoid (meters) (Geodetic).
      double height() const
         throw();

         /// return the coordinate system for this Position
      CoordinateSystem getCoordinateSystem() const
         throw() 
      { return system; };

         /// return geodetic latitude (deg N)
      double getGeodeticLatitude() const
         throw()
      { return geodeticLatitude(); }

         /// return geocentric latitude (deg N)
      double getGeocentricLatitude() const
         throw()
      { return geocentricLatitude(); }

         /// return longitude (deg E) (either geocentric or geodetic)
      double getLongitude() const
         throw()
      { return longitude(); }

         /// return height above ellipsoid (meters)
      double getAltitude() const
         throw()
      { return height(); }

         /// return height above ellipsoid (meters)
      double getHeight() const
         throw()
      { return height(); }

         /// return ECEF X coordinate (meters)
      double getX() const
         throw()
      { return X(); }

         /// return ECEF Y coordinate (meters)
      double getY() const
         throw()
      { return Y(); }

         /// return ECEF Z coordinate (meters)
      double getZ() const
         throw()
      { return Z(); }

         /// return spherical coordinate angle theta (deg) (90 - geocentric latitude)
      double getTheta() const
         throw()
      { return theta(); }

         /// return spherical coordinate angle phi (deg) (same as longitude)
      double getPhi() const
         throw()
      { return phi(); }

         /// return radius
      double getRadius() const
         throw()
      { return radius(); }

      // ----------- Part  8: member functions: set -----------------------------
      //
         /**
          * Set the geoid values for this Position given a geoid.
          * @param geoid Pointer to the GeoidModel.
          * @throw GeometryException if input is NULL.
          */
      void setGeoidModel(const GeoidModel *geoid)
         throw(GeometryException);

         /**
          * Set the Position given geodetic coordinates; system is set to Geodetic.
          * @param lat geodetic latitude in degrees North
          * @param lon geodetic longitude in degrees East
          * @param ht height above the ellipsoid in meters
          * @return a reference to this object.
          * @throw GeometryException on invalid input
          */
      Position& setGeodetic(const double lat,
                            const double lon,
                            const double ht,
                            const GeoidModel *geoid = NULL)
         throw(GeometryException);

         /**
          * Set the Position given geocentric coordinates; system is set to Geocentric
          * @param lat geocentric latitude in degrees North
          * @param lon geocentric longitude in degrees East
          * @param rad radius from the Earth's center in meters
          * @return a reference to this object.
          * @throw GeometryException on invalid input
          */
      Position& setGeocentric(const double lat,
                              const double lon,
                              const double rad)
         throw(GeometryException);

         /**
          * Set the Position given spherical coordinates; system is set to Spherical
          * @param theta angle from the Z-axis (degrees)
          * @param phi angle from the X-axis in the XY plane (degrees)
          * @param rad radius from the center in meters
          * @return a reference to this object.
          * @throw GeometryException on invalid input
          */
      Position& setSpherical(const double theta,
                             const double phi,
                             const double rad)
         throw(GeometryException);

         /**
          * Set the Position given ECEF coordinates; system is set to Cartesian.
          * @param X ECEF X coordinate in meters.
          * @param Y ECEF Y coordinate in meters.
          * @param Z ECEF Z coordinate in meters.
          * @return a reference to this object.
          */
      Position& setECEF(const double X,
                        const double Y,
                        const double Z)
         throw();

         /**
          * Set the Position given an array of ECEF coordinates;
          * system is set to Cartesian.
          * @param XYZ array[3] ECEF X,Y,Z coordinate in meters.
          * @return a reference to this object.
          */
      Position& setECEF(const double XYZ[3])
         throw()
      { return setECEF(XYZ[0],XYZ[1],XYZ[2]); }

         /**
          * Set the Position given ECEF coordinates; system is set to Cartesian.
          * @param XYZ ECEF X,Y,Z coordinates in meters.
          * @return a reference to this object.
          */
      Position& setECEF(const Triple& XYZ)
         throw()
      { return setECEF(XYZ[0],XYZ[1],XYZ[2]); }

      // ----------- Part 9: member functions: setToString, printf -------------
      //
         /**
          * setToString, similar to scanf, this function takes a string and a
          * format describing string in order to define Position
          * values.  The parameters it can take are listed below and
          * described above with the printf() function.
          *
          * The specification must be sufficient to define a Position.
          * The following table lists combinations that give valid
          * Positions. Anything more or other combinations will give
          * unknown (read as: "bad") results so don't try it.  Anything
          * less will throw an exception.
          *
          * @code
          *  %X %Y %Z  (cartesian or ECEF in kilometers)
          *  %x %y %z  (cartesian or ECEF in meters)
          *  %a %l %r  (geocentric lat,lon,radius, longitude E, radius in meters)
          *  %A %L %h  (geodetic lat,lon,height, longitude E, height in meters)
          *  %a %w %R  (geocentric lat,lon,radius, longitude W, radius in kilometers)
          *  %A %W %H  (geodetic lat,lon,height, longitude W, height in kilometers)
          *  %t %p %r  (spherical theta, phi, radius, degrees and meters)
          *  %T %P %R  (spherical theta, phi, radius, radians and kilometers)
          * @endcode
          *
          * So
          * @code
          * pos.setToString("123.4342,9328.1982,-128987.399", "%X,%Y,%Z");
          * @endcode
          *
          * works but 
          *
          * @code
          * pos.setToString("123.4342,9328.1982", "%X,%Y");
          * @endcode
          * doesn't work (incomplete specification because it doesn't specify
          * a Position).
          *
          * Whitespace is unimportant here; the function will handle it.
          * The caller must ensure that that the extra characters in
          * the format string (ie '.' ',') are in the same relative
          * location as they are in the actual string; see the example above.
          *
          * @param str string from which to get the Position coordinates
          * @param fmt format to use to parse \c str.
          * @throw GeometryException if \c fmt is an incomplete or invalid
          *    specification
          * @throw FormatException if unable to scan \c str.
          * @throw StringException if an error occurs manipulating the
          * \c str or \c fmt strings.
          * @return a reference to this object.
          */
      Position& setToString(const std::string& str,
                            const std::string& fmt)
         throw(GeometryException,
               DayTime::FormatException,
               StringUtils::StringException);


         // if you can see this, ignore the \'s below, as they are for
         // the nasty html-ifying of doxygen.  Browsers try to
         // interpret the % and they get all messed up.
         /**
          * Format this Position into a string.
          *
          * Generate and return a string containing formatted
          * Position coordinates, formatted by the specification \c fmt.
          *
          * \li \%x   X() (meters)
          * \li \%y   Y() (meters)
          * \li \%z   Z() (meters)
          * \li \%X   X()/1000 (kilometers)
          * \li \%Y   Y()/1000 (kilometers)
          * \li \%Z   Z()/1000 (kilometers)
          * \li \%A   geodeticLatitude() (degrees North)
          * \li \%a   geocentricLatitude() (degrees North)
          * \li \%L   longitude() (degrees East)
          * \li \%l   longitude() (degrees East)
          * \li \%w   longitude() (degrees West)
          * \li \%W   longitude() (degrees West)
          * \li \%t   theta() (degrees)
          * \li \%T   theta() (radians)
          * \li \%p   phi() (degrees)
          * \li \%P   phi() (radians)
          * \li \%r   radius() meters
          * \li \%R   radius()/1000 kilometers
          * \li \%h   height() meters
          * \li \%H   height()/1000 kilometers
          *
          * @param fmt format to use for this time.
          * @return a string containing this Position in the
          * representation specified by \c fmt.
          */
      std::string printf(const char *fmt) const
         throw(StringUtils::StringException);

         /// Format this time into a string.
         /// @see printf(const char*)
      std::string printf(const std::string& fmt) const
         throw(StringUtils::StringException) 
      { return printf(fmt.c_str()); }

         /// Returns the string that operator<<() would print.
      std::string asString() const
         throw(StringUtils::StringException);

      // ----------- Part 10: functions: fundamental conversions ---------------
      // 
         /** Fundamental conversion from spherical to cartesian coordinates.
          * @param trp (input): theta, phi (degrees), radius
          * @param xyz (output): X,Y,Z in units of radius
          * Algorithm references: standard geometry.
          */
      static void convertSphericalToCartesian(const Triple& tpr,
                                              Triple& xyz)
         throw();

         /** Fundamental routine to convert cartesian to spherical coordinates.
          * The zero vector is converted to (90,0,0).
          * @param xyz (input): X,Y,Z
          * @param trp (output): theta, phi (degrees), radius (units of input)
          * Algorithm references: standard geometry.
          */
      static void convertCartesianToSpherical(const Triple& xyz,
                                              Triple& tpr)
         throw();


         /** Fundamental routine to convert ECEF (cartesian) to geodetic coordinates,
          * (Geoid specified by semi-major axis and eccentricity squared).
          * The zero vector is converted to (90,0,-R(earth)).
          * @param xyz (input): X,Y,Z in meters
          * @param llh (output): geodetic lat(deg N), lon(deg E),
          *                             height above ellipsoid (meters)
          * @param A (input) Earth semi-major axis
          * @param eccSq (input) square of Earth eccentricity
          * Algorithm references: Leick, "GPS Satellite Surveying," 2nd edition.
          */
      static void convertCartesianToGeodetic(const Triple& xyz,
                                             Triple& llh,
                                             const double A,
                                             const double eccSq)
         throw();

         /** Fundamental routine to convert geodetic to ECEF (cartesian) coordinates,
          * (Geoid specified by semi-major axis and eccentricity squared).
          * @param llh (input): geodetic lat(deg N), lon(deg E),
          *                             height above ellipsoid (meters)
          * @param A (input) Earth semi-major axis
          * @param xyz (output): X,Y,Z in meters
          * @param eccSq (input) square of Earth eccentricity
          * Algorithm references: Leick, "GPS Satellite Surveying," 2nd edition.
          */
      static void convertGeodeticToCartesian(const Triple& llh,
                                             Triple& xyz,
                                             const double A,
                                             const double eccSq)
         throw();


         /** Fundamental routine to convert cartesian (ECEF) to geocentric
          * The zero vector is converted to (0,0,0).
          * @param xyz (input): X,Y,Z
          * @param llr (output): geocentric lat(deg N), lon(deg E),
          *                              radius (units of input)
          */
      static void convertCartesianToGeocentric(const Triple& xyz,
                                               Triple& llr)
         throw();

         /** Fundamental routine to convert geocentric to cartesian (ECEF)
          * @param llr (input): geocentric lat(deg N),lon(deg E),radius
          * @param xyz (output): X,Y,Z (units of radius)
          */
      static void convertGeocentricToCartesian(const Triple& llr,
                                               Triple& xyz)
         throw();


         /** Fundamental routine to convert geocentric to geodetic
          * @param llr (input): geocentric lat(deg N),lon(deg E),radius (meters)
          * @param geodeticllh (output): geodetic latitude (deg N),
          *            longitude (deg E), and height above ellipsoid (meters)
          * @param A (input) Earth semi-major axis
          * @param eccSq (input) square of Earth eccentricity
          */
      static void convertGeocentricToGeodetic(const Triple& llr,
                                              Triple& geodeticllh,
                                              const double A,
                                              const double eccSq)
         throw();

         /** Fundamental routine to convert geodetic to geocentric 
          * @param geodeticllh (input): geodetic latitude (deg N),
          *            longitude (deg E), and height above ellipsoid (meters)
          * @param llr (output): geocentric lat (deg N),lon (deg E),radius (meters)
          * @param A (input) Earth semi-major axis
          * @param eccSq (input) square of Earth eccentricity
          */
      static void convertGeodeticToGeocentric(const Triple& geodeticllh,
                                              Triple& llr,
                                              const double A,
                                              const double eccSq)
         throw();

      // ----------- Part 11: operator<< and other useful functions -------------
      //
         /**
         * Stream output for Position objects.
         * @param s stream to append formatted Position to.
         * @param t Position to append to stream \c s.
         * @return reference to \c s.
         */
      friend std::ostream& operator<<(std::ostream& s,
                                      const Position& p);

         /**
         * Compute the range in meters between two Positions.
         * Input Positions are not modified.
         * @param A,B Positions between which to find the range
         * @return the range (in meters)
         * @throw GeometryException if geoid values differ.
         *        or if transformTo(Cartesian) fails
         */
      friend double range(const Position& A,
                          const Position& B)
         throw(GeometryException);

         /**
         * Compute the radius of the ellipsoidal Earth, given the geodetic latitude.
         * @param geolat geodetic latitude in degrees
         * @return the Earth radius (in meters)
         */
      static double radiusEarth(const double geolat,
                                const double A,
                                const double eccSq)
         throw();

         /**
         * A member function that calls the non-member radiusEarth() for
         * this Position.
         * @return the Earth radius (in meters)
         */
      double radiusEarth() const
         throw()
      {
         Position p(*this);
         p.transformTo(Position::Geodetic);
         return Position::radiusEarth(p.theArray[0], p.AEarth, p.eccSquared);
      }

         /**
         * A member function that computes the elevation of the input
         * (Target) position as seen from this Position.
         * @param Target the Position which is observed to have the
         *        computed elevation, as seen from this Position.
         * @return the elevation in degrees
         */
      double elevation(const Position& Target) const
         throw(GeometryException);

         /**
         * A member function that computes the elevation of the input
         * (Target) position as seen from this Position, using a Geodetic
         * (ellipsoidal) system.
         * @param Target the Position which is observed to have the
         *        computed elevation, as seen from this Position.
         * @return the elevation in degrees
         */
      double elevationGeodetic(const Position& Target) const
         throw(GeometryException);

         /**
         * A member function that computes the azimuth of the input
         * (Target) position as seen from this Position.
         * @param Target the Position which is observed to have the
         *        computed azimuth, as seen from this Position.
         * @return the azimuth in degrees
         */
      double azimuth(const Position& Target) const
         throw(GeometryException);

         /**
         * A member function that computes the azimuth of the input
         * (Target) position as seen from this Position, using a Geodetic
         * (ellipsoidal) system.
         * @param Target the Position which is observed to have the
         *        computed azimuth, as seen from this Position.
         * @return the azimuth in degrees
         */
      double azimuthGeodetic(const Position& Target) const
         throw(GeometryException);

         /**
         * A member function that computes the position at which a signal, which
         * is received at this Position and there is observed at the (input)
         * azimuth and elevation angles, crosses a model ionosphere that is
         * taken to be a thin shell at constant (input) height.
         * This function will not transform this Position, and it will return
         * a Position in the same system; the algorithm itself is done in the
         * geocentric coordinate system.
         * @param elev elevation angle in degrees of the signal at reception
         * @param azim azimuth angle in degrees of the signal at reception
         * @param ionoht height of the ionosphere, in meters
         * @return Position IPP the position of the ionospheric pierce point,
         *     in the same coordinate system as *this; *this is not modified.
         */
      Position getIonosphericPiercePoint(const double elev,
                                         const double azim,
                                         const double ionoht) const
         throw();

      // ----------- Part 12: private functions and member data -----------------
      //
   private:

         /** Initialization function, used by the constructors.
          * @param a coordinate [ X(m), or latitude (degrees N) ]
          * @param b coordinate [ Y(m), or longitude (degrees E) ]
          * @param c coordinate [ Z, height above ellipsoid or radius, in m ]
          * @param s CoordinateSystem, defaults to Cartesian
          * @param geiod pointer to a GeoidModel, default NULL (WGS84)
          * @throw GeometryException on invalid input.
          */
      void initialize(const double a,
                     const double b,
                     const double c,
                     CoordinateSystem s = Cartesian,
                     GeoidModel *geoid = NULL)
         throw(GeometryException);

         /* Values of the coordinates, defined for each system as follows;
         *    Cartesian  : X,Y,Z in meters
         *    Geocentric : Latitude(degrees N), Longitude(degrees E),
         *                    Radius (meters)
         *    Geodetic   : Latitude(degrees N), Longitude(degrees E),
         *                    Height above ellipsoid (meters)
         *    Spherical  : theta (degrees) - angle from the z axis
         *                 phi (degrees) - angle in xy plane from x axis toward
         *                                     y axis (same as longitude)
         *                 radius (meters?) - distance from origin
         */
      // use std::valarray<double> theArray;  -- inherit from Triple

         /// semi-major axis of Earth (meters)
      double AEarth;

         /// square of geoid eccentricity
      double eccSquared;

         /// see #CoordinateSystem
      CoordinateSystem system;

         /// tolerance used in comparisons
      double tolerance;

   };   // end class Position

   //@}

}  // namespace gpstk

#endif   // GPSTK_POSITION_HPP
