v_vector.h

Go to the documentation of this file.

/*=========================================================================
  Original version -- Copyright (c) 2006 Sandia Corporation.
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

#ifndef VERDICT_VECTOR
#define VERDICT_VECTOR

#include "moab/verdict.h"
#include <cmath>
#include <cassert>

// computes the dot product of 3d vectors
inline double dot_product( double vec1[], double vec2[] )
{
  double answer =  vec1[0] * vec2[0] +
     vec1[1] * vec2[1] +
     vec1[2] * vec2[2];
  return answer;
}

// normalize a vector
inline void normalize( double vec[] )
{
  double x = sqrt( vec[0]*vec[0] +
             vec[1]*vec[1] +
             vec[2]*vec[2] );

  vec[0] /= x;
  vec[1] /= x;
  vec[2] /= x;

}

// computes the cross product
inline double * cross_product( double vec1[], double vec2[], double answer[] )
{
  answer[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
  answer[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
  answer[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
  return answer;
}

// computes the length of a vector
inline double length ( double vec[] )
{
  return sqrt ( vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] );
}

// computes the square length of a vector
inline double length_squared (double vec[] )
{
  return (vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] );
}

// computes the interior angle between 2 vectors in degrees
inline double interior_angle( double vec1[], double vec2[] )
{
  double cosAngle, angleRad;
  double length1 = length(vec1);
  double length2 = length(vec2);
  assert( (length1 > 0.0 && length2 > 0.0) );

  cosAngle = dot_product(vec1, vec2) / (length1 * length2);

  if ((cosAngle > 1.0) && (cosAngle < 1.0001))
  {
    cosAngle = 1.0;
  }
  else if (cosAngle < -1.0 && cosAngle > -1.0001)
  {
    cosAngle = -1.0;
  }
  else
  {
    assert(cosAngle < 1.0001 && cosAngle > -1.0001);
  }

  angleRad = acos(cosAngle);
  return( (angleRad * 180.) / VERDICT_PI );
}

#endif