rVector.cpp

//#include <iostream.h>
#include <math.h>
#include "rPoint.h"
#include "rPointi.h"
#include "rPoint2D.h"
#include "rLatLon.h"
#include "rVector.h"
#include "rVectori.h"
#include "rLine.h"
#include "rMatrix.h"


/*
 * Class declaration^
 */
rVector::rVector () : rPoint () 
{}
rVector::rVector (rPoint2D pt) {
	set2d (pt);
}
rVector::rVector(rPoint p)
  : rPoint (p)
{ }
/*
rVector::rVector(rVector &p)
  : rPoint (p.x,p.y,p.z)
{ }
*/
rVector::rVector(rPoint A,rPoint B)
  : rPoint (B.x-A.x,B.y-A.y,B.z-A.z)
{ }
rVector::rVector (double i,double j,double k) 
  : rPoint (i,j,k)
{ }
//-------------------------------------
rVector &rVector::operator=(rVector p) {
	x=p.x; y=p.y; z=p.z; 
	return (*this);
}

rVector &rVector::operator=(rVectori p) {
	x=p.x; y=p.y; z=p.z; 
	return (*this);
}

rVector &rVector::operator=(rPoint p) { 
	x=p.x; y=p.y; z=p.z; 
	return (*this);
}

rVector &rVector::operator=(rPoint2D p)
{

	set2d (p);
	return *this;
}
void rVector::set2d (rPoint2D p)
{
	/* y is equivalent to lattitude and x is longitutde. */
	x = sin (PI/2.0-p.y) * cos (p.x);
	y = sin (PI/2.0-p.y) * sin (p.x);
	z = cos (PI/2.0-p.y);
}
rVector &rVector::operator=(rLatLon p) { 
	// not implemented yet, see old mathpack.
	/* y is equivalent to lattitude and x is longitutde. */
	x = sin (PI/2.0-p.lat) * cos (p.lon);
	y = sin (PI/2.0-p.lat) * sin (p.lon);
	z = cos (PI/2.0-p.lat);
	
	return (*this);
}

rVector &rVector::operator=(rPointi p) { 
	x=p.x; y=p.y; z=p.z; 
	return (*this);
}

//rVector::operator rPoint ()
//{
//  rPoint c;
//  c = *this;
//  return c;
//}

rVector::operator double *()
{
  return &x;
}

rVector::operator rLine ()
{
  rLine l(*this);
  //l = (rLine)*this;
  return l;
}


void rVector::normalize () {
	if (iszero()) return;  
	double len = sqrt (x*x+y*y+z*z); 
	if (len==0) return; 
	
	x = x/len; y = y/len; z = z/len; 
}
rVector rVector::normal () {
	rVector v (0,0,0); 
	if (iszero())return v; 
	double len = sqrt (x*x+y*y+z*z);  
	if (len==0) return v; 
	v.set(x/len,y/len,z/len); 
	return v;
}
void rVector::unit () { 
	if (iszero())return ; 
	double len = sqrt (x*x+y*y+z*z); 
	if (len==0) return; 
	x = x/len; y = y/len; z = z/len; 
}




void rVector::scale (rVector b)
{
  x = x*b.x;
  y = y*b.y;
  z = z*b.z;
}
int rVector::near (rVector &v,double n)
{
   return ((v.x-n) < x) && ((v.x+n) > x) &&
   ((v.y-n) < y) && ((v.y+n) > y) &&
   ((v.z-n) < z) && ((v.z+n) > z) ;
}

int rVector::nearxy (rVector &v,double n)
{
   return ((v.x-n) < x) && ((v.x+n) > x) &&
   ((v.y-n) < y) && ((v.y+n) > y);
}
int rVector::near (rVector &v,rVector &n)
{
   return ((v.x-n.x) < x) && ((v.x+n.x) > x) &&
   ((v.y-n.y) < y) && ((v.y+n.y) > y) &&
   ((v.z-n.z) < z) && ((v.z+n.z) > z) ;
}


/*
	Cross Product
 */
rVector rVector::operator*(rVector b)
{
  rVector c;

  c.x = y*b.z-z*b.y;
  c.y = z*b.x-x*b.z;
  c.z = x*b.y-y*b.x;
  return (c);
}
rVector &rVector::operator*=(rVector b)
{
  rVector c;

  c.x = y*b.z-z*b.y;
  c.y = z*b.x-x*b.z;
  c.z = x*b.y-y*b.x;
  *this = c;
  return (*this);
}

rVector rVector::operator+(rVector b)
{
  rVector c;

  c.x = x+b.x;
  c.y = y+b.y;
  c.z = z+b.z;
  return (c);
}
rVector rVector::operator-(rVector b)
{
  rVector c;

  c.x = x-b.x;
  c.y = y-b.y;
  c.z = z-b.z;
  return (c);
}

/*
	Scale
rVector rVector::operator*(double b)
{
  rVector c;

  c.x = c.x*b;
  c.y = c.y*b;
  c.z = c.z*b;
  return (c);
}
rVector rVector::operator/(double b)
{
  rVector c;

  c.x = c.x/b;
  c.y = c.y/b;
  c.z = c.z/b;
  return (c);
}
 */
rVector rVector::cross(rVector b)
{
  rVector c;

  c.x = y*b.z-z*b.y;
  c.y = z*b.x-x*b.z;
  c.z = x*b.y-y*b.x;
  return c;
}
/*
	Dot Product
 */
double rVector::operator|(rVector b)
{
  return (x*b.x+y*b.y+z*b.z);
}
double rVector::dot(rVector u)
{
  return (x*u.x+y*u.y+z*u.z);
}
/*
 	One rVector divided by another.
 */
rVector rVector::operator/(rVector b)
{
  if (b.x==0 || b.y == 0 || b.z == 0) return b;
  rVector c (x/b.x,y/b.y,z/b.z);
  return (c);
}
/*
        Add rVectors by scalars.
 */
rVector operator+(double a,rVector b)
{
  rVector c(a+b.x,a+b.y,a+b.z);
  return c;
}
rVector operator+(rVector b,double a)
{
  rVector c(a+b.x,a+b.y,a+b.z);
  return c;
}
/*
        Subtract rVectors by scalars.
 */
rVector operator-(double a,rVector b)
{
  rVector c(a-b.x,a-b.y,a-b.z);
  return c;
}
rVector operator-(rVector b,double a)
{
  rVector c(a-b.x,a-b.y,a-b.z);
  return c;
}
/*
        Multiply rVectors by scalars.
 */
rVector operator*(double a,rVector b)
{
  rVector c(a*b.x,a*b.y,a*b.z);
  return c;
}
rVector operator*(rVector b,double a)
{
  rVector c(a*b.x,a*b.y,a*b.z);
  return c;
}
/*
        Divide rVectors by scalars.
 */
rVector operator/(double a,rVector b)
{
  if (b.x == 0 || b.y == 0 || b.z == 0) {
    
    return b;
  }
  rVector c(a/b.x,a/b.y,a/b.z);
  return c;
}
rVector operator/(rVector b,double a)
{
  if (b.x == 0 || b.y == 0 || b.z == 0) return b;

  rVector c(b.x/a,b.y/a,b.z/a);
  return c;
}
/*
        Dot two rVectors.
 */

double operator|(rVector &a,rVector &b)
{
  return (a.x*b.x+a.y*b.y+a.z*b.z);
}

/*
 	One rVector divided by another.
 */

rVector operator/(rVector &a,rVector &b)
{
  if (b.x == 0 || b.y == 0 || b.z == 0) return b;

  rVector c (a.x/b.x,a.y/b.y,a.z/b.z);
  return (c);
}

/*
	Add two rVectors
 */
rVector operator+(rVector a,rVector b)
{
  rVector c(a.x+b.x,a.y+b.y,a.z+b.z);
  return c;
}
rPoint operator+(rPoint a,rVector b)
{
  rVector c(a.x+b.x,a.y+b.y,a.z+b.z);
  return c;
}

/*
        Subtract two rVectors.
 */
rVector operator-(rVector a,rVector b)
{
  rVector c(a.x-b.x,a.y-b.y,a.z-b.z);
  return c;
}
rPoint operator-(rPoint a,rVector b)
{
  rVector c(a.x-b.x,a.y-b.y,a.z-b.z);
  return c;
}




rVector& rVector::operator*=(double s)
{
  x *= s; y *= s; z *= s;
  return (*this);
}
rVector& rVector::operator/=(double s)
{
  if (s==0) return *this;
  x /= s; y /= s; z /= s;
  return (*this);
}
rVector& rVector::operator+=(double s)
{
  x += s; y += s; z += s;
  return (*this);
}
rVector& rVector::operator+=(rVector b)
{
  x += b.x; y += b.y; z += b.z;
  return (*this);
}
rVector& rVector::operator-=(rVector b)
{
  x -= b.x; y -= b.y; z -= b.z;
  return (*this);
}

rVector& rVector::operator-=(double s)
{
  x -= s; y -= s; z -= s;
  return (*this);
}
rVector& rVector::operator/=(rVector b)
{
  if (b.x == 0 || b.y == 0 || b.z == 0) return *this;
  x /= b.x; y /= b.y; z /= b.z;
  return (*this);
}
int rVector::is_acute_with (rVector &v)
{
  return (angle_with (v)==-1);
}
int rVector::is_obtuse_with (rVector &v)
{
  return (angle_with (v)==1);
}
/*
   this function determines whether the angle between two rVectors is acute,
   obtuse, or 180 degrees.  acuteness depends upon the direction of the cross 
   product in relation to the two rVectors.
 */
int rVector::angle_with (rVector v)
{
  rVector cp; cp  = *this*v;
  double thesin_plus = cp.magnitude ();
  cp = v* *this;
  double thesin_neg = cp.magnitude ();
  double thedot = this->dot (v);
  if (thedot == 1) return -1;
  if (thesin_plus == 0) return 0;
  if (thedot < 0) {
    if (thesin_plus < thesin_neg) return 1; // obtuse angle
    else if (thesin_plus > thesin_neg) return -1; // acute angle
    else return 0;// colinear
  }
  else if (thedot > 0) {
    if (thesin_plus > thesin_neg) return 1; // obtuse angle
    else if (thesin_plus < thesin_neg) return -1; // acute angle
    else return 0; //colinear
  }
  else {
    	return 90.0*PI/180.0;

    /* perturb one of the rVectors a little to make them not perpendicular, then
	do the test again.
     */
    //rVector tmp((v.x+x)/2.0,(v.y+y)/2.0,(v.z+z)/2.0);
   // return angle_with (tmp);
  }
}
int rVector::angle_with (rVector &v,rVector &cp,rVector &cn)
{
  double thesin_plus = cp.magnitude ();
  double thesin_neg = cn.magnitude ();
  double thedot = this->dot (v);
  if (thedot == 1) return -1;
  if (thesin_plus == 0) return 0;
  if (thedot < 0) {
    if (thesin_plus < thesin_neg) return 1; // obtuse angle
    else if (thesin_plus > thesin_neg) return -1; // acute angle
    else return 0;// colinear
  }
  else if (thedot > 0) {
    if (thesin_plus > thesin_neg) return 1; // obtuse angle
    else if (thesin_plus < thesin_neg) return -1; // acute angle
    else return 0; //colinear
  }
  else {
  	return 90.0*PI/180.0;
  	
    /* perturb one of the rVectors a little to make them not perpendicular, then
	do the test again.
     */
    //rVector tmp((v.x+x)/2.0,(v.y+y)/2.0,(v.z+z)/2.0);
   // return angle_with (tmp);
  }
}
int rVector::angle_with (rVector &v,double &thesin_plus,double &thesin_neg,double &thedot)
{
  if (thedot == 1) return -1;
  if (thesin_plus == 0) return 0;
  if (thedot < 0) {
    if (thesin_plus < thesin_neg) return 1; // obtuse angle
    else if (thesin_plus > thesin_neg) return -1; // acute angle
    else return 0;// colinear
  }
  else if (thedot > 0) {
    if (thesin_plus > thesin_neg) return 1; // obtuse angle
    else if (thesin_plus < thesin_neg) return -1; // acute angle
    else return 0; //colinear
  }
  else {
    	return 90.0*PI/180.0;

    /* perturb one of the rVectors a little to make them not perpendicular, then
	do the test again.
     */
    //rVector tmp((v.x+x)/2.0,(v.y+y)/2.0,(v.z+z)/2.0);
    //return angle_with (tmp);
  }
}