#ifndef _utl_BasicVector_h_
#define _utl_BasicVector_h_

/**
   \file
   Basic vector class template for geometry

   \author Lukas Nellen
   \version $Id$
   \date 23 May 2003
*/

#include <CLHEP/Geometry/Normal3D.h>
#include <CLHEP/Geometry/Point3D.h>
#include <CLHEP/Geometry/Vector3D.h>
#include <boost/tuple/tuple.hpp>
#include <utl/CoordinateSystem.h>
#include <utl/unconfig.h>
#include <utl/config.h>
#include <utl/GeometryException.h>
#include <utl/ConsecutiveEnumFactory.h>
#include <utl/ThrowPolicy.h>

static const char CVSId_utl_BasicVector[] =
  "$Id$";


namespace utl {

  /**
     \class BasicVector BasicVector.h utl/BasicVector.h

     \brief Basic vector class template for geometry

     This class wraps the concrete vector classes from CLHEP. Since most of
     the interface is the same for all types of vector-like objects.

     It acts as a base class for Vector, Point, and AxialVector.

     \author Lukas Nellen
     \date 23 May 2003
     \ingroup geometry
  */

  template<class HepVector>
  class BasicVector {

  public:
    typedef HepVector DataType;

    /// Base class for classes indicating coordinate types.
    class CoordinateType {
    public:
      virtual ~CoordinateType() { }
      virtual double x(double p1, double p2, double p3) const = 0;
      virtual double y(double p1, double p2, double p3) const = 0;
      virtual double z(double p1, double p2, double p3) const = 0;
      /// Coordinate type integral typification.
      /**
        Alternative typification via integral values for ease
        of streaming.
        Note that the values are explicitly defined with
        the integral value to be used in XML, binary files or 
        database.
      */   
      enum Kind { 
        eCartesian    = 0, 
        eCylindrical  = 1, 
        eSpherical    = 2
      };
      /// Coordinate type string typification.
      /**
        Alternative typification via strings 
        for ease of streaming and clarity.
       */  
      static const char * const KindTags[]; 
     /**
        \brief Convenience typedef for creation of Kind enumerators.
      */ 
      typedef utl::ConsecutiveEnumFactory<
          Kind, 
          eSpherical, 
          KindTags,
          utl::ThrowPolicy<Kind, CoordinateSystemException> >
        KindCreator;
      /// Factory method via enumeration.
      static const CoordinateType& Create(Kind k) {
        switch(k) {
          case eCartesian:
            return BasicVector<HepVector>::kCartesian;
          case eCylindrical:
            return BasicVector<HepVector>::kCylindrical;
          case eSpherical:
            return BasicVector<HepVector>::kSpherical;
        }
        // This shouldn't happen if the switch is kept in sync with the enum;
        // despite that is put here in order to avoid a compiler warning due
        // to the (supposedly) missing return at the end of the function 
        // (which is unreacheable under the former condition).
        std::ostringstream e;
        e << "No conversion from enum value " << k << " to utl::BasicVector<HepVector>::CoordinateType";
        throw CoordinateSystemException(e.str());
      }
    };

  protected:
    /// Class to have a type for indicating cartesian coordinate components
    /**
       Using a constant of this class allows to use virtual calls or
       overload resolution to select the correct constructor or transformation
       therein, when initialising the cartesian representation of a vector.
       Check the way the constructors of this class are defined.

       The concrete implementations of this class provide the functions
       \c x, \c y, and \c z to compute the transformation from a given
       coordinate representation to Cartesian coordinates.

       Only a single instance of this class exists as a public,
       static, constant member of this class.
    */

    class cCartesian : public CoordinateType {
    public:
      cCartesian() {};
      virtual double x(double x, double /*y*/, double /*z*/) const { return x; }
      virtual double y(double /*x*/, double y, double /*z*/) const { return y; }
      virtual double z(double /*x*/, double /*y*/, double z) const { return z; }
    };

    /// Class to have a type for indicating cylindrical coordinate components
    /**
       The arguments are in the order \c rho, \c phi, \c z.
     */
    class cCylindrical : public CoordinateType {
    public:
      cCylindrical() {};
      virtual double x(double rho, double phi, double z) const;
      virtual double y(double rho, double phi, double z) const;
      virtual double z(double rho, double phi, double z) const;
    };

    /// Class to have a type for indicating spherical coordinate components
    /**
       The arguments are in the order \c r, \c theta, \c phi.
     */
    class cSpherical : public CoordinateType {
    public:
      cSpherical() {};
      virtual double x(double r, double theta, double phi) const;
      virtual double y(double r, double theta, double phi) const;
      virtual double z(double r, double theta, double phi) const;
    };

  public:
    /// Coordinate triple for easy getting or setting of coordinates
    typedef typename boost::tuple<double, double, double> Triple;

  protected:
    /// Construct a point from cartesian coordinates
    /**
       The last argument of the constructor is used only for
       overloading the various constructors using three real
       arguments. The type of the last argument fixes the
       representation used (cartesian, cylindrical, or spherical).

       For construction from a cartesian representation, the last
       argument has a default value. That way, it can be omitted,
       making cartesian the default representation when constructing a
       vector.
    */
    BasicVector(const double x, const double y, const double z,
                const CoordinateSystemPtr& theCoordinateSystem)
      : fCoordinateSystem(theCoordinateSystem), fVector(x, y, z)
    { CoordinateSystemValid(theCoordinateSystem); }


    /// Construct a point from coordinates in arbitrary representation
    /**
       This constructor is used to construct points, vectors, and
       axial vectors from their representation in cartesian,
       cylindrical, or spherical coordinates.

       The actual transformation to the internal representation is
       implemented using a class implementing the CoordinateType
       interface.

       The pre-defined, static \c CoordinateType objects
       \c kCartesian, \c kSpherical, and \c kCylindrical
       can be used to construct an object in Cartesian, cylindrical, and
       spherical coordinates. The interpretation of the arguments is as
       follows:
       - Cartesian: \c x, \c y, \c z
       - Spherical: \c r (distance to origin), \c theta (zenith),
         \c phi (azimuth)
       - Cylindrical: \c rho (distance to z axis), \c phi (azimuth), \c z
    */
    BasicVector(double p1, double p2, double p3,
                const CoordinateSystemPtr& theCoordinateSystem,
                const CoordinateType& theType)
      : fCoordinateSystem(theCoordinateSystem),
        fVector(theType.x(p1, p2, p3),
                theType.y(p1, p2, p3),
                theType.z(p1, p2, p3))
    { CoordinateSystemValid(theCoordinateSystem); }

  public:
    bool operator==(const BasicVector& v) const
    { v.TransformTo(fCoordinateSystem); return this->fVector == v.fVector; }

    bool operator!=(const BasicVector& v) const
    { return !(*this == v); }

    double GetX(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.x(); }

    double GetY(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.y(); }

    double GetZ(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.z(); }

    /// Get coordinates as a tuple (x, y, z)
    Triple
    GetCoordinates(const CoordinateSystemPtr theCoordinateSystem)
      const
    {
      TransformTo(theCoordinateSystem);
      return boost::make_tuple(fVector.x(), fVector.y(), fVector.z());
    }

    /**
       \brief Get the coordinate system of the current internal representation

       \note This function is mainly for debugging purposes. In
       general, you \b can not \b rely on the coordinate system used
       for the internal representation having a particular value
       unless you set it first using \c TransformTo. In this case,
       though, you already have the information and you don't have to
       use this function.
    */
    CoordinateSystemPtr GetCoordinateSystem() const
    { return fCoordinateSystem; }

    // universally applicable operations

    // Assignment: rely on compiler to construct the correct operator

    /// Multiply by number
    BasicVector& operator*=(double a) { fVector *= a; return *this; }

    /// Divide by number
    BasicVector& operator/=(double a) { fVector /= a; return *this; }

    /// zenith (theta) angle in spherical coordinates
    double GetTheta(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.theta(); }
    /// cos of zenith (theta) angle
    double GetCosTheta(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.cosTheta(); }
    /// azimuth (phi) angle in spherical and cylindrical coordinates
    double GetPhi(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.phi(); }
    /// radius r in spherical coordinates coordinates (distance to origin)
    double GetR(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.mag(); }
    /// radius r^2 in spherical coordinates coordinates (distance to origin)^2
    double GetR2(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.mag2(); }
    /// radius r in cylindrical coordinates (distance to z axis)
    double GetRho(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.perp(); }
    /// radius r^2 in cylindrical coordinates (distance to z axis)^2
    double GetRho2(const CoordinateSystemPtr theCoordinateSystem) const
    { TransformTo(theCoordinateSystem); return fVector.perp2(); }
    /// normalize the vector
    void Normalize() { fVector = fVector.unit(); }

    /// Get spherical coordinates as a tuple (r, theta, phi)
    Triple
    GetSphericalCoordinates(const CoordinateSystemPtr theCoordinateSystem)
      const
    {
      TransformTo(theCoordinateSystem);
      return boost::make_tuple(fVector.r(), fVector.theta(), fVector.phi());
    }

    /// Get cylindrical coordinates as a tuple (rho, phi, z)
    Triple
    GetCylindricalCoordinates(const CoordinateSystemPtr theCoordinateSystem)
      const
    {
      TransformTo(theCoordinateSystem);
      return boost::make_tuple(fVector.perp(), fVector.phi(), fVector.z());
    }

    /// Check if two BasicVectors are nearly eaqual
    /**
      Two BasicVectors are nearly equal if either
      \f[
      |v_1-v_2| < \varepsilon
      \f]
      or
      \f[
      \frac{(v_1-v_2)^2}{(v_1+v_2)^2} < \varepsilon^2
      \f]
    */
    bool
    IsCloseTo(const BasicVector& v, const double eps = fgEpsilon)
      const
    {
      v.TransformTo(fCoordinateSystem);
      const double d = (fVector - v.fVector).mag2();
      if (d < eps*eps)
        return true;
      else
        return d / (fVector.mag2() + v.fVector.mag2()) < eps*eps;
    }

    /**
       \brief Transform vector to representation in new coordinate system

       This function is logically const, since it does not change the
       geometrical object, even though the internal representation changes.
    */
    void
    TransformTo(const CoordinateSystemPtr newCoordinateSystem)
      const
    {
      if (!SameCoordinateSystem(newCoordinateSystem) ) {
        fVector = TransformedVector(newCoordinateSystem);
        fCoordinateSystem = newCoordinateSystem;
      }
    }

    template<class V>
    bool SameCoordinateSystem(const BasicVector<V>& theVector) const
    { return fCoordinateSystem == theVector.fCoordinateSystem; }

    bool SameCoordinateSystem(const CoordinateSystemPtr theCoordinateSystem) const
    { return fCoordinateSystem == theCoordinateSystem; }

    static const cCartesian kCartesian;
    static const cCylindrical kCylindrical;
    static const cSpherical kSpherical;

  protected:
    /// Destructor - NOT virtual since class is NEVER used polymorphically
    ~BasicVector() { }

    BasicVector()
      : fCoordinateSystem(CoordinateSystem::GetRootCoordinateSystem()),
        fVector(0, 0, 0)
    { }

    /// Constructor from internal components for use by operators
    BasicVector(const DataType& theVector,
                const CoordinateSystemPtr& theCoordinateSystem)
      : fCoordinateSystem(theCoordinateSystem),
        fVector(theVector)
    { }

    DataType
    TransformedVector(const CoordinateSystemPtr newCoordinateSystem) const;

    mutable CoordinateSystemPtr fCoordinateSystem;
    mutable DataType fVector;

    static const double fgEpsilon;

  };


  typedef BasicVector<HepGeom::Normal3D<double> > AxialVectorBase;
  typedef BasicVector<HepGeom::Point3D<double> > PointBase;
  typedef BasicVector<HepGeom::Vector3D<double> > VectorBase;

} // utl


#endif // _utl_BasicVector_h_

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. -k"
// End: