// ------------------------------
//    vec3.h   3-vectors
// ------------------------------
#include "w.h"

#ifndef vec3_h
#define vec3_h

class vec
  { private:
	  double x,y,z;
	 public:
		vec(double xi = 0.0, double yi = 0.0, double zi = 0.0);
      vec(fptype1 fi, double xi, double yi, double zi); // f(v)
		friend void p(vec &u);
		friend void pl(vec &u);
		friend vec vsum(vec &u, vec &v);
		friend vec vdif(vec &u, vec &v);
		friend double vdot(vec &u, vec &v);
		friend double v3dot(vec &u, vec &v, vec &w);
		friend double vnorm(vec &u, vec &v);
		friend vec vcross(vec &u, vec &v);
      vec &operator=(const vec &);
		vec &operator=(const double &);  // v[i] = const
		double &operator[](int i);  // subscripts
		vec operator()(fptype1);    // v(f) = f(v)
		friend vec operator-(const vec &);
		friend vec operator+(const vec &, const vec &);
		friend vec operator-(const vec &, const vec &);
		friend vec operator*(const vec &, const vec &);

		};



// ------------------------------
//    vec3.cpp    3-vectors
// ------------------------------
	  vec::vec(double xi, double yi, double zi)
		  { x = xi; y = yi; z = zi; }

	  vec::vec(fptype1 fi, double xi, double yi, double zi)
		  { x = fi(xi); y = fi(yi); z = fi(zi); }

// -------------------------------------------------------------

		vec vsum(vec &u, vec &v)
			{ return vec(u.x + v.x, u.y + v.y, u.z + v.z); }

		vec vdif(vec &u, vec &v)
			{ return vec(u.x - v.x, u.y - v.y,u.z - v.z); }

		double vdot(vec &u, vec &v)
			{ return (u.x*v.x + u.y*v.y + u.z*v.z); }

		double vnorm(vec &u, vec &v)
			 { return sqrt(vdot(u,v)); }

		vec vcross(vec &u,  vec &v)
			{ return vec( u.y*v.z - u.z*v.y,
							  u.z*v.x - u.x*v.z,
							  u.x*v.y - u.y*v.x); }

		double v3dot(vec &u, vec &v, vec &w)
			{ vec t = vcross(v,w);
			  return vdot(u, t); }

		double &vec::operator[](int i)
			{ if (i == 1) return x;
			  else { if (i == 2) return y;
			         else return z; } }

		vec &vec::operator=(const vec &c)
			{ x = c.x; y = c.y; z = c.z;
			  return *this; }

		vec &vec::operator=(const double &a)
			{ x = y = z = a;
			  return *this; }

		vec vec::operator()(fptype1 fi)
			{
			  return (vec(fi(x),fi(y),fi(z))); }


		vec operator-(const vec &c)
			{ return vec(-c.x,-c.y,-c.z); }

  		vec operator+(const vec &c1, const vec &c2)
			{ return vec(c1.x+c2.x, c1.y+c2.y, c1.z+c2.z); }

		vec operator-(const vec &c1, const vec &c2)
			{ return vec(c1.x-c2.x, c1.y-c2.y, c1.z-c2.z); }

		vec operator*(const vec &c1, const vec &c2)
			{ return vec(c1.x*c2.x, c1.y*c2.y, c1.z*c2.z); }
 // --------------------------------------------------------
 // vector io

		void p(vec &u)
		  { cout << "(" << u.x << "," << u.y << "," << u.z << ")"; }

		void pl(vec &u) { p(u); nl(); }

  // ----------------------------------------------------------
 #endif