/*   []----------------------------------------[]
      |                Vector.cpp              |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       Vector routines...        |
      |                                        |
     []----------------------------------------[]
*/

// Copyright (C) 2005 M.F. Somers, Theoretical Chemistry Department, Leiden University
//
// This 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.
//
// http://www.gnu.org/licenses/lgpl.txt

#include "Global.h"
#include "Constants.h"
#include "Rand.h"
#include "Vector.h"

/* ------------------------------------------------------------------- */

double Angle( const Vector& A, const Vector& B )
{
 double cosangle, AbsA, AbsB;
 
 AbsA = Abs( A );
 AbsB = Abs( B );
 
 if( AbsA < CONST_EPSILON ) return( 0.0 );
 if( AbsB < CONST_EPSILON ) return( 0.0 );
 
 cosangle = Dot( A, B ) / ( AbsA * AbsB );

 if( fabs( cosangle ) > 1.0 - CONST_EPSILON )
  return( 0.0 );
 else
  return( acos( cosangle ) );
}

/* ------------------------------------------------------------------- */

Vector Cross( const Vector& A, const Vector& B )
{
 Vector C;
 
 C.X = A.Y*B.Z - A.Z*B.Y;
 C.Y = A.Z*B.X - A.X*B.Z;
 C.Z = A.X*B.Y - A.Y*B.X;
 
 return( C );
}

/* ------------------------------------------------------------------- */

Vector SphericalToCartesian( double r, double theta, double phi )
{
 Vector Cartesian;
 
 Cartesian.X = r * sin( theta ) * cos( phi );
 Cartesian.Y = r * sin( theta ) * sin( phi );
 Cartesian.Z = r * cos( theta );
 
 return( Cartesian );
}

/* ------------------------------------------------------------------- */

void   CartesianToSpherical( const Vector& V, double& r, double& theta, double& phi )
{
 r = Abs( V );
   
 if( r >= CONST_EPSILON )
 {
  phi   = acos( V.X / r );
  theta = acos( V.Z / r );
 } 
 else
  phi = theta = 0.0;
 
}

/* ------------------------------------------------------------------- */

void RotateVectorAroundAxis( Vector& V, Vector RotationAxis, double Angle )
{
 double   sRX = RotationAxis.X*RotationAxis.X;
 double   sRY = RotationAxis.Y*RotationAxis.Y;
 double   sRZ = RotationAxis.Z*RotationAxis.Z;
 double     N = (sRX+sRY+sRZ);

 if( N <= 0.0 ) return;

 double sqrtN = sqrt( N );
 double     D = Dot( V, RotationAxis );
 double     c = cos( Angle );
 double     s = sin( Angle );
 Vector     W; 
   
 W.X = RotationAxis.X*D + c*(V.X*(sRY+sRZ)-RotationAxis.X*(RotationAxis.Y*V.Y+RotationAxis.Z*V.Z)) +
       s*sqrtN*(V.Z*RotationAxis.Y-RotationAxis.Z*V.Y);
 
 W.Y = RotationAxis.Y*D + c*(V.Y*(sRX+sRZ)-RotationAxis.Y*(RotationAxis.X*V.X+RotationAxis.Z*V.Z)) +
       s*sqrtN*(V.X*RotationAxis.Z-RotationAxis.X*V.Z);        
       
 W.Z = RotationAxis.Z*D + c*(V.Z*(sRX+sRY)-RotationAxis.Z*(RotationAxis.X*V.X+RotationAxis.Y*V.Y)) +
       s*sqrtN*(V.Y*RotationAxis.X-RotationAxis.Y*V.X);              
       
 V = W / N;        
}

/* ------------------------------------------------------------------- */

inline double Remainder( double A, double B )
{
 double R = fmod( A, B );

 if( R >= 0.5 * B ) R = R - B;
 if( R <= -0.5 * B ) R = R + B;

 return( R );
}

/* ------------------------------------------------------------------- */

Vector MimimumSystemImageConventionVector( const Vector& V, const Vector& TheBox )
{
 Vector mV;
 
 mV.X = ( TheBox.X >= CONST_EPSILON ? Remainder( V.X, TheBox.X ) : 0.0 );
 mV.Y = ( TheBox.Y >= CONST_EPSILON ? Remainder( V.Y, TheBox.Y ) : 0.0 );
 mV.Z = ( TheBox.Z >= CONST_EPSILON ? Remainder( V.Z, TheBox.Z ) : 0.0 ); 

 return( mV );
}

/* ------------------------------------------------------------------- */

Vector ARandomNormalizedVector( void )
{
 Vector V;
 double AbsV;
 
 V.X = Rand();
 V.Y = Rand();
 V.Z = Rand();
 
 if( Rand() % 10000 > 5000 ) V.X = -V.X;
 if( Rand() % 10000 > 5000 ) V.Y = -V.Y;
 if( Rand() % 10000 > 5000 ) V.Z = -V.Z;
 
 AbsV = Abs( V ) ;
 if( AbsV >= CONST_EPSILON ) V /= AbsV;
 
 return( V );
} 

/* ------------------------------------------------------------------- */

ostream& operator<<( ostream& out, const Vector& A )
{
 return( out << "[ " << A.X << ", " << A.Y << ", " << A.Z << " ]" );
}

/* ------------------------------------------------------------------- */

istream& operator>>( istream& in, Vector& A )
{
 char LineRead[ 128 ];
 double X, Y, Z;
 
 in.getline( LineRead, sizeof( LineRead ), ']' );

 if( sscanf( LineRead, " [ %lf , %lf , %lf", &X, &Y, &Z ) == 3 )
  A.X = X, A.Y = Y, A.Z = Z;
   
 return( in ); 
}


/* ------------------------------------------------------------------- */