//******************************************************************************
//** SCATMECH: Polarized Light Scattering C++ Class Library
//** 
//** File: vector3d.h
//**
//** Thomas A. Germer
//** Optical Technology Division, National Institute of Standards and Technology
//** 100 Bureau Dr. Stop 8443; Gaithersburg, MD 20899-8443
//** Phone: (301) 975-2876; FAX: (301) 975-6991
//** Email: thomas.germer@nist.gov
//**
//** Version: 6.00 (February 2008)
//**
//******************************************************************************
#ifndef SCATMECH_VECTOR3D_H
#define SCATMECH_VECTOR3D_H

#include "scatmech.h"
#include <ios>

namespace SCATMECH {


    // 
    // template class Vector3D_Base is used by the class Vector3D
    //
    template <class TYPE>
    class Vector3D_Base {
        public:
            //
            // The components of the vector:
            //
            TYPE x,y,z;

            // Class constructors:
            Vector3D_Base<TYPE>() {}
            Vector3D_Base<TYPE>(const TYPE& a,const TYPE& b,const TYPE& c) 
                                {x=a;y=b;z=c;}      
        
            // Non-converting copy constructor:
            Vector3D_Base<TYPE>(const Vector3D_Base<TYPE>& a) 
                : x(a.x),y(a.y),z(a.z) {}     

            // Assignment operator:
            Vector3D_Base<TYPE>& operator=(const Vector3D_Base<TYPE>& a) 
                    { x=a.x; y=a.y; z=a.z; return *this; }   
        
            // Addition of vectors:
            Vector3D_Base<TYPE> operator+(const Vector3D_Base<TYPE>& a) const 
                    { return Vector3D_Base<TYPE>(x+a.x,y+a.y,z+a.z); }
        
            // Addition of vectors: (Added by TAG 25 NOV 03)
            const Vector3D_Base<TYPE>& operator+=(const Vector3D_Base<TYPE>& a)
                    { return *this = *this+a; }

            // Subtraction of vectors:
            Vector3D_Base<TYPE> operator-(const Vector3D_Base<TYPE>& a) const 
                    { return Vector3D_Base<TYPE>(x-a.x,y-a.y,z-a.z); }
        
            // Subtraction of vectors: (Added by TAG 25 NOV 03)
            const Vector3D_Base<TYPE>& operator-=(const Vector3D_Base<TYPE>& a) 
                    { return *this = *this-a; }

            // Unary minus sign:
            Vector3D_Base<TYPE> operator-() const 
                    { return Vector3D_Base<TYPE>(-x,-y,-z); }
        
            // Scalar product of two vectors: 
            TYPE operator*(const Vector3D_Base<TYPE>& a) const 
                    { return x*a.x+y*a.y+z*a.z; }
        
            // Product of vector with scalar:
            Vector3D_Base<TYPE> operator*(const TYPE& b) const 
                    { return Vector3D_Base<TYPE>(x*b,y*b,z*b); }
            friend Vector3D_Base<TYPE> operator*(const TYPE& b,
                                                 const Vector3D_Base<TYPE>& a)
                    {return Vector3D_Base<TYPE>(a.x*b,a.y*b,a.z*b); }
        
            // Division of a vector by a scalar:
            Vector3D_Base<TYPE> operator/(const TYPE& a) const 
                    { return Vector3D_Base<TYPE>(x/a,y/a,z/a); }
        
            // The cross product of two vectors:
            friend Vector3D_Base<TYPE> cross(const Vector3D_Base<TYPE>& a,
                                             const Vector3D_Base<TYPE>& b)
                    { return Vector3D_Base<TYPE>(a.y*b.z-a.z*b.y,
                                                 a.z*b.x-a.x*b.z,
                                                 a.x*b.y-a.y*b.x); }

            // Multiplication of vector by scalar: (Added by TAG 25 NOV 03)
            const Vector3D_Base<TYPE>& operator*=(const TYPE& a) 
                    { return *this = *this*a; }
        
            // Division of vector by scalar: (Added by TAG 25 NOV 03)
            const Vector3D_Base<TYPE>& operator/=(const TYPE& a) 
                    { return *this = *this/a; }

            // Comparison... (Added by TAG 25 NOV 03)
            bool operator==(const Vector3D_Base<TYPE>& b) const
                { return (x==b.x && y==b.y && z==b.z); }

            // Comparison... (Added by TAG 25 NOV 03)
            bool operator!=(const Vector3D_Base<TYPE>& b) const
                { return (x!=b.x || y!=b.y || z!=b.z); }

            // Send to output stream
            friend std::ostream& operator<<(std::ostream& os,const Vector3D_Base<TYPE>& v) 
                { return os << '(' << v.x << ',' << v.y << ',' << v.z << ')'; }

            friend std::istream& operator>>(std::istream& is,Vector3D_Base<TYPE>& v)
            {   TYPE vx,vy,vz;
                char c;
                is >> c;
                if (c=='(') {
                    is >> vx;
                    is >> c;
                    if (c==',') {
                        is >> vy;
                        is >> c;
                        if (c==',') {
                            is >> vz;
                            is >> c;
                            if (c==')') {
                                v.x=vx, v.y=vy, v.z=vz;
                                return is;
                            }
                        }
                    }
                }
                is.setstate(std::ios::failbit);
                return is;
            }

        protected:
        
            // A type-converting copy function:
            template <class TYPE2>
            void copy(const Vector3D_Base<TYPE2>& a) 
                { x=a.x; y=a.y; z=a.z; }

    };

    // 
    // Define templated Vector3D
    //
    template <class TYPE> class Vector3D: public Vector3D_Base<TYPE> {};

    //
    // Define Vector3D<double>
    // 
    template <>
    class Vector3D<double> : public Vector3D_Base<double> {
    public:
        // Constructors:
        Vector3D<double>() {};
        Vector3D<double>(const Vector3D<double>& a) 
            : Vector3D_Base<double>(a) {}
        Vector3D<double>(const Vector3D_Base<double>& a) 
            : Vector3D_Base<double>(a) {}
        Vector3D<double>(double ax,double ay,double az) 
            :  Vector3D_Base<double>(ax,ay,az) {}
    };

    // The norm of a vector:
    inline 
    double 
    Norm(const Vector3D_Base<double>& a) 
    { 
        return sqrt(sqr(a.x)+sqr(a.y)+sqr(a.z)); 
    }

    // A vector of unit length:
    inline
    Vector3D<double> 
    unit(const Vector3D_Base<double>& a) 
    { 
        return a/Norm(a); 

    }

    // A function returning a vector having specific polar coordinates:
    inline 
    Vector3D_Base<double> 
    polar(double r,double theta,double phi) 
    {
        return Vector3D_Base<double>(r*sin(theta)*cos(phi),
                              r*sin(theta)*sin(phi),
                              r*cos(theta));
    }


    //
    // Defining Vector3D<COMPLEX>
    //
    template <>
    class Vector3D<COMPLEX > : public Vector3D_Base<COMPLEX > {
    public:
        // Constructors:
        Vector3D<COMPLEX >() {};
        Vector3D<COMPLEX >(const Vector3D<COMPLEX >& a) 
            : Vector3D_Base<COMPLEX >(a)  {}
        Vector3D<COMPLEX >(const Vector3D_Base<COMPLEX >& a)  
            : Vector3D_Base<COMPLEX >(a)  {}
        Vector3D<COMPLEX >(const COMPLEX & ax,
                           const COMPLEX & ay,
                           const COMPLEX & az) 
                           :  Vector3D_Base<COMPLEX >(ax,ay,az) {}

        // Constructor which converts from a real vector to a complex vector:
        Vector3D<COMPLEX >(const Vector3D_Base<double>& a) { copy(a); }

    };

    inline
    Vector3D<COMPLEX>
    Conj(const Vector3D_Base<COMPLEX>& a)
    {
        using std::conj;
        return Vector3D<COMPLEX>(conj(a.x),conj(a.y),conj(a.z));
    }

    // The norm of a complex vector:
    inline
    double 
    Norm(const Vector3D_Base<COMPLEX >& a) 
    { 
        using std::norm;
        return sqrt(norm(a.x)+norm(a.y)+norm(a.z)); 
    }

    // A vector of unit length:
    inline
    Vector3D<COMPLEX > 
    unit(const Vector3D_Base<COMPLEX >& a) 
    { 
        return a/Norm(a); 
    }

    // A vector perpendicular to two vectors:
    Vector3D<double> perpto(const Vector3D<double>& a,
                            const Vector3D<double>& b);
    Vector3D<COMPLEX > perpto(const Vector3D<COMPLEX >& a, 
                              const Vector3D<COMPLEX >& b);               

    // Some useful typedefs for the SCATMECH library:
    typedef Vector3D<double> Vector;
    typedef Vector3D<COMPLEX > CVector;


} // namespace SCATMECH 

#endif // VECTOR3D_H