#ifndef _utl_ReferenceEllipsoid_h_
#define _utl_ReferenceEllipsoid_h_

#include <string>
#include <map>
#include <boost/tuple/tuple.hpp>
#include <utl/CoordinateSystem.h>
#include <utl/AugerException.h>

namespace utl {

  class Point;
  class TransformationMatrix;


  /**
     \class ReferenceEllipsoid ReferenceEllipsoid.h utl/ReferenceEllipsoid.h

     \brief Reference ellipsoids for UTM transformations

     We only use global geodetic data, so the ellipsoids satisfy

     * centre coincides with centre of the earth

     * rotational ellipsoids (two of the three major axis are the same)

     * axis of rotational symmetry aligned with axis of rotation of the earth.

     Uses standard auger units for length and angles.

     \author Lukas Nellen
     \date 03 Aug 2003
     \ingroup geometry

     \todo Calculate tangent plane at a given point in a given
     coordinate system.
  */

  class ReferenceEllipsoid {

  public:
    /// ID's of known reference ellipsoid's
    enum EllipsoidID { eWGS84 = 1 };

    /// Coordinate triple
    typedef boost::tuple<double, double, double> Triple;

    /// Construct reference ellipsoid from equatorial radius, eccentricity, and coordinate system
    /**
       Construct reference ellipsoid from

       - equatorial radius

       - eccentricity

       - coordinate system (ECEF at the center of the ellipsoid

       This constructor should only be used for tests. Regular code should get
       ReferenceEllipsoids from the registry using a well-known id or, in
       special cases, string.
     */
    ReferenceEllipsoid(double theEquatorialRadius,
                       double theEccentricity2,
                       const CoordinateSystemPtr theECEF) :
      fEquatorialRadius(theEquatorialRadius),
      fEccentricity2(theEccentricity2),
      fECEF(theECEF)
    { }

    /// Get equatorial radius (\f$ a \f$)
    double GetEquatorialRadius() const { return fEquatorialRadius; }
    /// Get Polar radius (\f$ b \f$)
    double GetPolarRadius() const;
    /// Get eccentricity
    double GetEccentricity2() const { return fEccentricity2; }
    /// Get flattening
    double GetFlattening() const;
    /// Get the ECEF
    const CoordinateSystemPtr GetECEF() const { return fECEF; }

    /// compare for equality
    bool
    operator==(const ReferenceEllipsoid& r)
      const
    {
      return fEquatorialRadius == r.fEquatorialRadius &&
             fEccentricity2 == r.fEccentricity2 &&
             fECEF == r.fECEF;
    }

    /// compare for inequality
    bool operator!=(const ReferenceEllipsoid& r) const
    { return !operator==(r); }

    /// Convert Point to Lat/Long/Height
    Triple PointToLatitudeLongitudeHeight(const Point& thePoint) const;

    /// Convert Lat/Long/Height to Point
    Point LatitudeLongitudeHeightToPoint(double const latitude,
                                         double const longitude,
                                         double const height) const;

    /// Convert Lat/Long/Height to Point
    Point LatitudeLongitudeHeightToPoint(const Triple& theGeodeticCoordinates) const;

    /// Translation and rotation to go to local coordinate system
    boost::tuple<TransformationMatrix, TransformationMatrix>
    TransformECEFToLocalSystem(const Point& theOrigin) const;

    /// Get known ellipsoid by registered ID
    static const ReferenceEllipsoid& Get(const EllipsoidID theID);

    /// Get the auger standard ellipsoid: wgs84
    static const ReferenceEllipsoid& GetWGS84()
    { return Get(eWGS84); }

    /// Initialise the registry specifying the ECEF (instead of the Root CS)
    static void InitWithECEF(const CoordinateSystemPtr& theECEF);

    /// Get EllipsoidID from string
    static EllipsoidID GetEllipsoidIDFromString(const std::string& str);

    /**
       \class InvalidEllipsoidException ReferenceEllipsoid.h utl/ReferenceEllipsoid.h
       \brief Report request for non-existent ellipsoid
    */
    class InvalidEllipsoidException;

  private:
    // don't implement - doesn't make sense to have ellipsoid w/o paramters
    ReferenceEllipsoid();

    /// Radius at equator
    double fEquatorialRadius;
    /// first eccentricity squared \f$ e^2\f$
    double fEccentricity2;
    /// Coordinate system to use as ECEF (Earth Centred Earth Fixed)
    CoordinateSystemPtr fECEF;

    typedef std::map<EllipsoidID, const ReferenceEllipsoid> Registry;
    static Registry* gfRegistry;

    static void RegisterOneEllipsoid(const EllipsoidID theID,
                                     double theEquatorialRadius,
                                     double theEccentricity,
                                     const CoordinateSystemPtr& theECEF);

    double GetRn(const double phi) const;

  };


  class ReferenceEllipsoid::InvalidEllipsoidException : public AugerException {

  public:
    /// Construct Invalid Ellipsoid exception with message
    InvalidEllipsoidException(const std::string& message = std::string())
      : AugerException(message) { }

    /// Retrieve verbose exception name
    virtual std::string GetExceptionName() const
    { return "Invalid Reference Ellipsoid requested exception"; }

  };

}


#endif

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: