// *********************************************************************
//
// Vector: a C++ class of a vector 
//
// Last updated 25/05/2000
// Abdon Pijpelink, abdon@its.chem.uva.nl
//
// In the case of overloaded mathematical operators, no size-checks are
// performed to gain speed.
//
// Variables:
//   cv        array of Complex objects
//   size		length of the vector
//
// *********************************************************************

#include "DVector.h"
#include <assert.h>
#include <iostream>
#include <math.h>

    // ================================================================
	DVector::DVector(int sz) {
          init(0, sz);
        }
     // ================================================================
     DVector::DVector(const double *array, int sz) {
          init(array, sz);
     }

     // ================================================================
     DVector::DVector(const DVector &cV) {
          init(cV.cv, cV.size);
     }

     // ================================================================
	void DVector::clear() {
		for (int f1=0; f1<size; f1++) {
			cv[f1] = 0.;
		}
	}

     // ================================================================
        void DVector::unitVector() {
		double sum = 0.;
                for (int f1=0; f1<getSize(); f1++)
			sum += (cv[f1] * cv[f1]);
                sum = sqrt(sum);
                for (int f1=0; f1<getSize(); f1++) {
			if(sum != 0)   cv[f1] = cv[f1] / sum;
			else   cv[f1] = 0;
                }
		if(sum != 0) mag=1.0;
		else mag = 0.0;
        }
	
     // ================================================================
     void DVector::init(const double *array, int sz) {
          size = sz;
	  //cout<<"Create cv = new double[size]\n";
	  cv = new double[size];
		
          assert (cv != 0);
          if (array != 0) {
		mag = 0.0;
		for (int f1=0; f1<size; f1++) {
               	  cv[f1] = array[f1];
		  mag += (cv[f1] * cv[f1]);
          	}
		mag = sqrt(mag);
	  }
	  else {
		  for (int f1=0; f1<size; f1++) {
			cv[f1] = double(0.0);
		  }
		mag = 0.0;
	  }
     }

     // ================================================================
     DVector::~DVector() { 
       delete [] cv; 
       size =0;
       //cout<<"Delete [] cv\n";
     }      // destructor

     // ================================================================
     DVector& DVector::operator=( const DVector &cV) {
          if (this == &cV) return *this;

          delete [] cv;
          init(cV.cv, cV.size);
          return *this;
     }

     // ================================================================
     double& DVector::operator[](int index) {
          return cv[index];
     }

     // ================================================================
     int DVector::getSize() {
          return size;
     }

     // ================================================================
     double DVector::getMag() {
          return mag;
     }

     // ================================================================
     // Output stream
     // ================================================================
/*
     ostream& operator << (ostream& os, DVector& cV) {
          for(int f1=0; f1<cV.getSize(); f1++) {
                    os << cV[f1];
          }
               os << "\n";
          return os;
     }
*/
     // === Dot Product/Scalar Product ================================
	double DVector::operator * (DVector &vector) {
		double result = 0;
		
           /*	for(int f1=0; f1<=2; f1++) {
		   cout << cv[f1]<<"\n";
		   cout << vector[f1]<<"\n";
	        }
		
	    */	
		for(int f1=0; f1<getSize(); f1++) {
			result += cv[f1] * vector[f1];
		}
		//cout << result<<"\n";
		return result;
	}


   // === Cross Product =============================================
	DVector DVector::operator ^ (DVector &vector) {
	        DVector result(3);
	//cout<< cv[1]<<"*"<<vector[2]<<"="<<cv[1]*vector[2]<<"\n";
	//cout<< cv[2]<<"*"<<vector[1]<<"="<<cv[2]*vector[1]<<"\n";
	//cout<<cv[1]*vector[2] - cv[2]*vector[1]<<"\n";
                result[0]=cv[1]*vector[2] - cv[2]*vector[1];
                result[1]=cv[2]*vector[0] - cv[0]*vector[2];
                result[2]=cv[0]*vector[1] - cv[1]*vector[0];
	/*	for(int f1=0; f1<=2; f1++) {
		   cout << cv[f1]<<"\n";
		   cout << vector[f1]<<"\n";
		   cout << result[f1]<<"\n";
	        }
	*/	
	        return result;

        }
   //===============================================================	
   // Define operator %: vector A * B ->vector C
   // = Vector A times vector B to generate new vector C where C[i]=A[i]*B[i]
	DVector DVector::operator % (DVector &vector) {
	        DVector result(3);
                result[0]=cv[0]*vector[0];
                result[1]=cv[1]*vector[1];
                result[2]=cv[2]*vector[2];
	        return result;

        }

   //======Multiplication by a double constant
	DVector DVector::operator * (double dblnum) {
	        DVector result(3);
	//	cout <<"dblnum="<< dblnum<<"\n";
		for (int f=0; f<getSize();f++){
		  result[f]= dblnum *cv[f];
		}
		/*
		for(int f1=0; f1<=2; f1++) {
		   cout << "dblnum *cv[f] ="<<result[f1]<<"\n";
	        }
		*/
	        return result;
        }
         
   //======Devision by a double constant
	DVector DVector::operator / (double dblnum) {
	        DVector result(3);
	//	cout <<"dblnum="<< dblnum<<"\n";
		for (int f=0; f<getSize();f++){
		  result[f]= cv[f]/dblnum;
		}  
	/*	for(int f1=0; f1<=2; f1++) {
		   cout << result[f1]<<"\n";
	        }
	*/	
	        return result;
        }
   //======Addition with a vector
	DVector DVector::operator + (DVector &vector) {
	        DVector result(3);
		result[0]=0.;
		result[1]=0.;
		result[2]=0.;
		for (int f=0; f<getSize();f++){
		  result[f] = cv[f] + vector[f];
		} 
		/*
		for(int f=0; f<=2; f++) {
		   cout << "cv[f]+vector ="<<result[f]<<"\n";
	        }
		*/
	        return result;
        }
   //======Subtraction with a vector
	DVector DVector::operator - (DVector &vector) {
	        DVector result(3);
		for (int f=0; f<getSize();f++){
		  result[f] = cv[f] - vector.cv[f];
		}  
	        return result;
        }