#ifndef  STARLAB_VECTOR_H
#  define  STARLAB_VECTOR_H
#define PR(x)  cerr << #x << " = " << x << " "
#define PRC(x) cerr << #x << " = " << x << ",  "
#define PRL(x) cerr << #x << " = " << x << "\n"

#ifndef real
# define real double
#endif
/*-----------------------------------------------------------------------------
 *  mvector  --  a class for 3-dimensional vectors
 *-----------------------------------------------------------------------------
 */
class mvector
    {
    private:
        real element[VLEN];

    public:
//	Default: initialize to zero.
        mvector(real c = 0)
	    {for (int i = 0; i<VLEN;i++)element[i] = c;}

//  []: the return type is declared as a reference (&), so that it can be used
//  on the left-hand side of an asignment, as well as on the right-hand side,
//  i.e.  v[1] = 3.14  and  x = v[2]  are both allowed and work as expected.

        real & operator [] (int i)       {return element[i];}

	inline void print() {for(int i = 0; i<VLEN;i++)
			     cout << element[i] << " ";
			     cout << "\n";}

//	Unary -

        const mvector operator - ()
	    {mvector v;for(int i = 0; i<VLEN;i++)v.element[i] = - element[i];
	          return v;}

        friend const mvector operator + (const mvector &, const mvector & );
        friend const mvector operator - (const mvector &, const mvector & );
        friend mvector operator * (real, const mvector & );
        friend real operator * (const mvector &, const mvector & ); 
        friend mvector operator * (const mvector &, real);
        friend mvector operator / (const mvector &, real);

//	Mvector +=, -=, *=, /=

        mvector& operator += (const mvector& b)
	    {for(int i = 0; i<VLEN;i++)element[i] += b.element[i];       
	     return *this;}

	mvector& operator -= (const mvector& b)
	    {for(int i = 0; i<VLEN;i++)element[i] -= b.element[i];
	     return *this;}

	mvector& operator *= (const real b)
	    {for(int i = 0; i<VLEN;i++)element[i] *= b;
	     return *this;}

	mvector& operator /= (const real b)
	    {for(int i = 0; i<VLEN;i++)element[i] /= b;
	     return *this;}
	
	//      Input / Output
	
        friend ostream & operator << (ostream & , const mvector & );
	
	friend istream & operator >> (istream & , mvector & );
    };

inline  ostream & operator << (ostream & s, const mvector & v)
{
    for(int i = 0; i<VLEN;i++)s << v.element[i] << "  " ;
    return s;
}
inline  istream & operator >> (istream & s, mvector & v)
{for(int i = 0; i<VLEN;i++)s >> v.element[i]; return s;}

inline  const mvector operator + (const mvector &v1, const mvector & v2)
{
    mvector v3;
    for(int i = 0; i<VLEN;i++)v3.element[i] = v1.element[i]+ v2.element[i];
    return v3;
}
inline  const mvector operator - (const mvector &v1, const mvector & v2)
{
    mvector v3;
    for(int i = 0; i<VLEN;i++)v3.element[i] = v1.element[i]- v2.element[i];
    return v3;
}

inline  mvector operator * (real b, const mvector & v)
{
    mvector v3;
    for(int i = 0; i<VLEN;i++)v3.element[i] = b* v.element[i];
    return v3;
}
inline  mvector operator * (const mvector & v, real b)
{
    mvector v3;
    for(int i = 0; i<VLEN;i++)v3.element[i] = b* v.element[i];
    return v3;
}

inline  real operator * (const mvector & v1, const mvector & v2)
{
    real x = 0;
    for(int i = 0; i<VLEN;i++) x+= v1.element[i]* v2.element[i];
    return x;
}


inline  mvector operator / (const mvector & v, real b)
{
    mvector v3;
    for(int i = 0; i<VLEN;i++)v3.element[i] =  v.element[i]/b;
    return v3;
}

typedef const mvector (*vfunc_ptr)(const mvector &); 
typedef const mvector (vfunc)(const mvector &); 

#endif

//=======================================================================//
//  +---------------+        _\|/_        +------------------------------\\
//  |  the end of:  |         /|\         |  inc/vector.h
//  +---------------+                     +------------------------------//
//========================= STARLAB =====================================\\