// SpacePointVector.cpp

#include "SpacePointVector.h"

#include <cmath>

#ifndef DEFECT_NO_STDLIB_NAMESPACES
#ifndef DEFECT_CMATH_NOT_STD
using std::sqrt;
#endif
using std::ostream;
using std::endl;
#endif

// The rotation from Cartesian to cylindrical:
// |  cos(phi)  sin(phi)   0  |
// | -sin(phi)  cos(phi)   0  |
// |     0         0       1  |
//
// The rotation from cylindrical to Cartesian:
// |  cos(phi) -sin(phi)   0  |
// |  sin(phi)  cos(phi)   0  |
// |     0         0       1  |
//
// The rotation from Cartesian to spherical:
// |  sin(tht)*cos(phi)  sin(tht)*sin(phi)    cos(tht)  |
// |  cos(tht)*cos(phi)  cos(tht)*sin(phi)   -sin(tht)  |
// |          -sin(phi)           cos(phi)       0      |
//
// The rotation from spherical to Cartesian:
// |  sin(tht)*cos(phi)  cos(tht)*cos(phi)   -sin(phi)  |
// |  sin(tht)*sin(phi)  cos(tht)*sin(phi)    cos(phi)  |
// |  cos(tht)          -sin(tht)                0      |
//

//**********************************************************************
// SpacePointVector
//**********************************************************************

// Default constructor.
// Initial point is the origin with phi = theta = 0.
// Vector has zero length.

SpacePointVector::SpacePointVector( ) {
  _vx = _vy = _vz = 0.0;
}

//**********************************************************************

// Constructor from a space point.
// Vector has zero length.

SpacePointVector::SpacePointVector(const SpacePoint& spt)
: SpacePoint(spt) {
  _vx = _vy = _vz = 0.0;
}

//**********************************************************************

// Copy constructor.

SpacePointVector::SpacePointVector(const SpacePointVector& rhs)
: SpacePoint(rhs), _vx(rhs._vx),  _vy(rhs._vy),  _vz(rhs._vz) { }

//**********************************************************************

// Destructor.
SpacePointVector::~SpacePointVector() { }

//**********************************************************************

// Assignment operator.

SpacePointVector&
SpacePointVector::operator=(const SpacePointVector& rhs) {
  if ( this == &rhs ) return *this;
  ((SpacePoint&) *this) = rhs;
  _vx = rhs._vx;
  _vy = rhs._vy;
  _vz = rhs._vz;
  return *this;
}

//**********************************************************************

// Return the rxy component.

double SpacePointVector::v_rxy() const {
  return cos_phi()*_vx + sin_phi()*_vy;
}

//**********************************************************************

// Return the phi component.

double SpacePointVector::v_phi() const {
  return -sin_phi()*_vx + cos_phi()*_vy;
}

//**********************************************************************

// Return the rxyz component.

double SpacePointVector::v_rxyz() const {
  return sin_theta()*cos_phi()*_vx
       + sin_theta()*sin_phi()*_vy + cos_theta()*_vz;
}

//**********************************************************************

// Return the theta component.

double SpacePointVector::v_theta() const {
  return cos_theta()*cos_phi()*_vx
       + cos_theta()*sin_phi()*_vy - sin_theta()*_vz;
}

//**********************************************************************

// Return the magnitude.

double SpacePointVector::magnitude() const {
  return sqrt( _vx*_vx + _vy*_vy + _vz*_vz );
}

//**********************************************************************
// CartesianPointVector
//**********************************************************************

// Construct a Cartesian vector from coordinates and direction.

CartesianPointVector::
CartesianPointVector(double x, double y, double z,
                     double vx, double vy, double vz)
: SpacePointVector( CartesianPoint(x,y,z) ) {
  _vx = vx;
  _vy = vy;
  _vz = vz;
}

//**********************************************************************

// Construct a Cartesian vector from space point and direction.

CartesianPointVector::
CartesianPointVector(const SpacePoint& spt,
                                 double vx, double vy, double vz)
: SpacePointVector(spt) {
  _vx = vx;
  _vy = vy;
  _vz = vz;
}

//**********************************************************************
// CylindricalPointVector
//**********************************************************************

// Construct a cylindrical vector from coordinates and direction.

CylindricalPointVector::
CylindricalPointVector(double r, double phi, double z,
                       double vr, double vphi, double vz)
: SpacePointVector( CylindricalPoint(r,phi,z) ) {
  _vx = cos_phi()*vr - sin_phi()*vphi;
  _vy = sin_phi()*vr + cos_phi()*vphi;
  _vz = vz;
}

//**********************************************************************

// Construct a cylindrical vector from space point and direction.

CylindricalPointVector::
CylindricalPointVector(const SpacePoint& spt,
                                 double vr, double vphi, double vz)
: SpacePointVector(spt) {
  _vx = cos_phi()*vr - sin_phi()*vphi;
  _vy = sin_phi()*vr + cos_phi()*vphi;
  _vz = vz;
}

//**********************************************************************

// Construct a spherical vector from coordinates and direction.

SphericalPointVector::
SphericalPointVector(double r, double phi, double theta,
                               double vr, double vtheta, double vphi)
: SpacePointVector( SphericalPoint(r,phi,theta) ) {
  _vx = sin_theta()*cos_phi()*vr 
      + cos_theta()*cos_phi()*vtheta - sin_phi()*vphi;
  _vy = sin_theta()*sin_phi()*vr 
      + cos_theta()*sin_phi()*vtheta + cos_phi()*vphi;
  _vz = cos_theta()*vr - sin_theta()*vtheta;
}

//**********************************************************************

// Construct a spherical vector from space point and direction.

SphericalPointVector::
SphericalPointVector(const SpacePoint& spt,
                     double vr, double vtheta, double vphi)
: SpacePointVector(spt) {
  _vx = sin_theta()*cos_phi()*vr 
      + cos_theta()*cos_phi()*vtheta - sin_phi()*vphi;
  _vy = sin_theta()*sin_phi()*vr 
      + cos_theta()*sin_phi()*vtheta + cos_phi()*vphi;
  _vz = cos_theta()*vr - sin_theta()*vtheta;
}

//**********************************************************************
// Free functions.
//**********************************************************************

// Equality.

bool operator==(const SpacePointVector& lhs,
                const SpacePointVector& rhs) {
  return (const SpacePoint&) lhs == (const SpacePoint&) rhs &&
         SpacePoint::equal( lhs.v_x(), rhs.v_x() ) &&
         SpacePoint::equal( lhs.v_y(), rhs.v_y() ) &&
         SpacePoint::equal( lhs.v_z(), rhs.v_z() );
}

//**********************************************************************

// Inequality.

bool operator!=(const SpacePointVector& lhs,
                const SpacePointVector& rhs) {
  return ! (lhs == rhs);
}

//**********************************************************************

// Output stream operator.

ostream& operator<<(ostream & stream, const SpacePointVector& svec) {
  stream << (SpacePoint&) svec << endl;
  stream << "     V_x: " << svec.v_x() << endl;
  stream << "     V_y: " << svec.v_y() << endl;
  stream << "     V_z: " << svec.v_z() << endl;
  stream << "   V_rxy: " << svec.v_rxy() << endl;
  stream << "  V_rxyz: " << svec.v_rxyz() << endl;
  stream << "  V_dphi: " << svec.v_phi() << endl;
  stream << " V_theta: " << svec.v_theta() << endl;
  stream << " Magnitude: " << svec.magnitude();
  return stream;
}

//**********************************************************************