/* 
   ORSA - Orbit Reconstruction, Simulation and Analysis
   Copyright (C) 2002-2004 Pasquale Tricarico
   
   This program is free software; you can redistribute it and/or
   modify it under the terms of the GNU General Public License
   as published by the Free Software Foundation; either version 2
   of the License, or (at your option) any later version.
   
   As a special exception, Pasquale Tricarico gives permission to
   link this program with Qt commercial edition, and distribute the
   resulting executable, without including the source code for the Qt
   commercial edition in the source distribution.
   
   This program 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 General Public License for more details.
   
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

#ifndef _ORSA_COORD_H_
#define _ORSA_COORD_H_

#include <vector>
#include <cmath>
#include <limits>

namespace orsa {
  
  class Vector {
  public:
    // constructors
    // inline Vector() { x = y = z = 0; }
    inline Vector() : x(0), y(0), z(0) { }
    // inline Vector(const Vector& v) { x = v.x; y = v.y; z = v.z; }
    inline Vector(const Vector & v) : x(v.x), y(v.y), z(v.z) { }
    // inline Vector(double _x, double _y, double _z) { x = _x; y = _y; z = _z; }
    inline Vector(double _x, double _y, double _z) : x(_x), y(_y), z(_z) { }
    
    // utils
    inline double GetX() const { return x; }
    inline double GetY() const { return y; }
    inline double GetZ() const { return z; }
    
    // operators
    inline Vector & operator += (const Vector & v) {
      x += v.x;
      y += v.y;
      z += v.z;
      return *this;
    }
   
    inline Vector & operator -= (const Vector & v) {
      x -= v.x;
      y -= v.y;
      z -= v.z;
      return *this;
    }
    
    inline Vector & operator *= (const double f) {
      x *= f;
      y *= f;
      z *= f;
      return *this;
    }
    
    inline Vector & operator /= (const double f) {
      x /= f;
      y /= f;
      z /= f;
      return * this;
    }
    
    // sign
    inline Vector operator + () const { return Vector(x,y,z); }   
    inline Vector operator - () const { return Vector(-x,-y,-z); }    
    
    // moved outside the class
    // inline Vector operator * (const double f) const { return Vector(x*f, y*f, z*f); }   
    // inline Vector operator / (const double f) const { return Vector(x/f, y/f, z/f); }
    
    // rotation
    Vector & rotate (const double, const double, const double);
    
    // set
    inline void Set(Vector v) {
      x = v.x;
      y = v.y;
      z = v.z;
    }
    
    inline void Set(double _x, double _y, double _z) {
      x = _x;
      y = _y;
      z = _z;
    }
    
    // metrics
    inline double Length() const {
      using std::sqrt;
      return sqrt( (x*x) +
                   (y*y) +
                   (z*z) );
    }
    
    inline double LengthSquared() const {
      return (x*x) +
             (y*y) +
             (z*z);
    }
    
    inline double ManhattanLength() const {
      using std::fabs;
      return (fabs(x)+fabs(y)+fabs(z));
    }
    
    inline bool IsZero() const {
      return ((x*x) + (y*y) + (z*z)) < (std::numeric_limits<double>::min() * 1.0e3);
    }
    
    // normalization
    inline Vector Normalized() const {
      double l = Length();
      if (l > (std::numeric_limits<double>::min() * 1.0e3))
	return Vector(x/l, y/l, z/l);
      else
        return Vector(0.0, 0.0, 0.0);
    }
    
    inline Vector & Normalize() {
      double l = Length();
      if (l > (std::numeric_limits<double>::min() * 1.0e3)) {
	x /= l;
	y /= l;
	z /= l;
      } else {
	z = 0.0;
	y = 0.0;
	z = 0.0;
      }
      return *this;
    }
    
  public:
    double x, y, z;
  };
  
  inline Vector operator * (const double f, const Vector & v) {
    return Vector(v.x*f, v.y*f, v.z*f);
  }
  
  inline Vector operator * (const Vector & v, const double f) {
    return Vector(v.x*f, v.y*f, v.z*f);
  }
  
  inline Vector operator / (const Vector & v, const double f) {
    return Vector(v.x/f, v.y/f, v.z/f); 
  }
  
  inline Vector operator + (const Vector& u, const Vector& v) {
    return Vector(u.x+v.x,
                  u.y+v.y,
                  u.z+v.z);
  }
  
  inline Vector operator - (const Vector& u, const Vector& v) {
    return Vector(u.x-v.x,
                  u.y-v.y,
                  u.z-v.z);
  }
  
  inline Vector ExternalProduct (const Vector& u, const Vector& v) {
    return Vector (u.y*v.z-u.z*v.y,
                   u.z*v.x-u.x*v.z,
                   u.x*v.y-u.y*v.x); 
  }
  
  inline Vector Cross (const Vector& u, const Vector& v) {
    return Vector (u.y*v.z-u.z*v.y,
                   u.z*v.x-u.x*v.z,
                   u.x*v.y-u.y*v.x); 
  }
  
  // scalar product
  inline double operator * (const Vector& u, const Vector& v) {
    return (u.x*v.x+
            u.y*v.y+
            u.z*v.z);
  }  
  
  inline bool operator == (const Vector &v1, const Vector &v2) {
    if (v1.x != v2.x) return false;
    if (v1.y != v2.y) return false;
    if (v1.z != v2.z) return false;
    return true;
  }
  
  inline bool operator != (const Vector &v1, const Vector &v2) {
    return !(v1 == v2);
  }
  
  // container for intepolation
  class VectorWithParameter : public Vector {
  public:
    double par;
  };
  
  void Interpolate(const std::vector < VectorWithParameter > v_in, const double x, Vector & v_out, Vector & err_v_out);
  
} // namespace orsa

#endif // _ORSA_COORD_H_
