bool OBRing::findCenterAndNormal(vector3 & center, vector3 &norm1, vector3 &norm2)
  {
    OBMol *mol= this->_parent;
    int j= 0;
    const int nA= this->_path.size();
    vector3 tmp;

    center.Set(0.0,0.0,0.0);
    norm1.Set(0.0,0.0,0.0);
    norm2.Set(0.0,0.0,0.0);
    for (j = 0; j != nA; ++j)
      {
        center += (mol->GetAtom(_path[j]))->GetVector();
      }
    center/= double(nA);

    for (j = 0; j != nA; ++j)
      {
        vector3 v1= (mol->GetAtom(_path[j]))->GetVector() - center;
        vector3 v2= (mol->GetAtom(_path[j+1==nA?0:j+1]))->GetVector() - center;
        tmp= cross(v1,v2);
        norm1+= tmp;
      }
    norm1/= double(nA);
    norm1.normalize();
    norm2= norm1;
    norm2 *= -1.0;
    return(true);
  }

//for only world
void UITextBox::setFace(vector3 right, vector3 down){
	right.normalize();
	down.normalize();

	direction = right;
	this->down = down;
}

 Camera(vector3 p, vector3 d, vector3 u) {
   position = p;
   direction = d.normalized();
   up = u.normalized();
   up = (u - direction.projected(u)).normalized();
   right = (direction ^ up).normalized();
 }

	matrix4x4 translation_matrix(vector3 const& v)
	{
		return matrix4x4 (1, 0, 0, 0, 
			0, 1, 0, 0, 
			0, 0, 1, 0, 
			v.x(), v.y(), v.z(), 1);
	}

const vector3 vector3::operator / (const vector3 &v2) const
{
	vector3 r;
	r.x() = x() / v2.x();
	r.y() = y() / v2.y();
	r.z() = z() / v2.z();
	return r;
}

float vector3::AngleBetweenVectors(const vector3& v0, const vector3& v1)
{
	float d = v0.Dot(v1);
	if ( d >= -1.f && d <= 1.f )
	{
		return acosf( d / ( v0.Length() * v1.Length() ) );
	}

	return 0.f;
}

vector3<float> get_rotation_axis(vector3<float> u, vector3<float> v) {
	u.normalize();
	v.normalize();
	// fix linear case
	if (u == v || u == -v) {
		v[0] += 0.1;
		v[1] += 0.0;
		v[2] += 0.1;
		v.normalize();
	}
	return u.cross(v);
}

vector3 vector3::Transform(const vector3& v, const matrix& m)
{
	vector3 result = vector3(	
						(v.x() * m(0,0)) + (v.x() * m(1,0)) + (v.x() * m(2,0))  + m(3,0),
						(v.y() * m(0,1)) + (v.y() * m(1,1)) + (v.y() * m(2,1))  + m(3,1),
						(v.z() * m(0,2)) + (v.z() * m(1,2)) + (v.z() * m(2,2) ) + m(3,2));
	return result;
	
}

 //! Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#calculateOrthogonalisationMatrix">blue-obelisk:calculateOrthogonalisationMatrix</a>
 void OBUnitCell::SetData(const vector3 v1, const vector3 v2, const vector3 v3)
 {
   matrix3x3 m (v1, v2, v3);
   _mOrtho.FillOrth(vectorAngle(v2,v3), // alpha
                    vectorAngle(v1,v3), // beta
                    vectorAngle(v1,v2), // gamma
                    v1.length(),        // a
                    v2.length(),        // b
                    v3.length());       // c
   _mOrient = m.transpose() * _mOrtho.inverse();
   _spaceGroup = NULL;
   _spaceGroupName = "";
   _lattice = OBUnitCell::Undefined;
 }

double get_rotation_angle(vector3<float> u, vector3<float> v) {
	u.normalize();
	v.normalize();
	double cosine_theta = u.dot(v);
	// domain of arccosine is [-1, 1]
	if (cosine_theta > 1) {
		cosine_theta = 1;
	}
	if (cosine_theta < -1) {
		cosine_theta = -1;
	}
	double angle = acos(cosine_theta);
	return angle;
}

///
//	FacePlane()
//   
//		Which of the six face-plane(s) is point P outside of? 
//
static
int FacePlane(const vector3& p)
{
	int outcode;

	outcode = 0;
	if (p.x() >  .5) outcode |= 0x01;
	if (p.x() < -.5) outcode |= 0x02;
	if (p.y() >  .5) outcode |= 0x04;
	if (p.y() < -.5) outcode |= 0x08;
	if (p.z() >  .5) outcode |= 0x10;
	if (p.z() < -.5) outcode |= 0x20;
   
	return(outcode);
}

//Return size of stabilizer group of a Cartesian displacement (given Cartesian symmetry rotations)
inline int nStabilizer(const vector3<>& n, const std::vector< matrix3<> >& symCart)
{	int nStab = 0;
	for(const matrix3<>& m: symCart)
		if((n - m * n).length_squared() < symmThresholdSq * n.length_squared())
			nStab++;
	return nStab;
}

void matrix::setrot(vector3 p1,vector3 p2)
{
    vector3 _ax=p1*p2;
    _ax.selfnormalize();
    double _th=p1.angle(p2);
    setrot(_th,_ax);
}

void moon::display(const vector3 &moon_pos, const vector3 &sun_pos, double max_view_dist) const
{
	vector3 moon_dir = moon_pos.normal();
	double moon_size = max_view_dist/20;
	float moon_azimuth = atan2(-moon_dir.y, moon_dir.x);
	float moon_elevation = asin(moon_dir.z);

	glsl_moon->use();
	glsl_moon->set_gl_texture(*map_diffuse, loc_diffcol, 0);
	glsl_moon->set_gl_texture(*map_normal, loc_nrml, 1);
	//	transform light into object space
	matrix4 roth = matrix4::rot_z(-RAD_TO_DEG(moon_azimuth));
	matrix4 rotv = matrix4::rot_y(-RAD_TO_DEG(moon_elevation));
	matrix4 model_mat = roth*rotv;
	vector3 l = model_mat.inverse() * sun_pos;
	vector3 nl = vector3(-l.y, l.z, -l.x).normal();	//	OpenGL coordinates
	glsl_moon->set_uniform(loc_lightdir, nl);

	//	render moon
	glPushMatrix();
	model_mat.multiply_gl();
	glTranslated(0.95*max_view_dist, 0, 0);

	primitives::textured_quad(vector3f( 0,  moon_size,  moon_size),
				  vector3f( 0, -moon_size,  moon_size),
				  vector3f( 0, -moon_size, -moon_size),
				  vector3f( 0,  moon_size, -moon_size),
				  *map_diffuse).render_plain();
	glPopMatrix();
}

double controller::ray_traverse(vector3 const & direction) const
{
    float coords[4];
    direction.getAs4Values(coords);
    vector3 dir(coords[0] * cos(-rotation) - coords[2] * sin(-rotation),
                coords[1], 
                coords[0] * sin(-rotation) + coords[2] * cos(-rotation));
    boost::lock_guard<boost::mutex> lock(mutex_);
    core::line3d<f32> ray;
    ray.start = position_ + vector3(0, 5, 0);
    ray.end = position_  + vector3(0, 5, 0) + dir * 1000.0f;

    core::vector3df intersection;
    core::triangle3df hitTriangle;
    
    scene::ISceneNode * selectedSceneNode =
      collision_manager_->getSceneNodeAndCollisionPointFromRay(
        ray,
        intersection, // This will be the position of the collision
        hitTriangle, // This will be the triangle hit in the collision
        0, //IDFlag_IsPickable, // This ensures that only nodes that we have
                // set up to be pickable are considered
        0); // Check 
   return (intersection - position_).getLength();
}

	quat rotation_quat(vector3 const& axe, real const& angle )
	{
		// create (sin(a/2)*axis, cos(a/2)) quaternion
		// which rotates the point a radians around "axis"
		quat res;
		vector4 u = vector4(axe.x(), axe.y(), axe.z(), 0); u.normalize();
		real sina2 = sin(angle/2);
		real cosa2 = cos(angle/2);

		res.qx() = sina2 * u.x();
		res.qy() = sina2 * u.y();
		res.qz() = sina2 * u.z();
		res.qw() = cosa2;

		return res;
	}

  /*! This method employs the static method matrix3x3::jacobi(...)
    to find the eigenvalues and eigenvectors of a symmetric
    matrix. On entry it is checked if the matrix really is
    symmetric: if isSymmetric() returns 'false', an OBError is
    thrown.
 
    \note The jacobi algorithm is should work great for all
    symmetric 3x3 matrices. If you need to find the eigenvectors
    of a non-symmetric matrix, you might want to resort to the
    sophisticated routines of LAPACK.
 
    @param eigenvals a reference to a vector3 where the
    eigenvalues will be stored. The eigenvalues are ordered so
    that eigenvals[0] <= eigenvals[1] <= eigenvals[2].
 
    @return an orthogonal matrix whose ith column is an
    eigenvector for the eigenvalue eigenvals[i]. Here 'orthogonal'
    means that all eigenvectors have length one and are mutually
    orthogonal. The ith eigenvector can thus be conveniently
    accessed by the GetColumn() method, as in the following
    example.
    \code
    // Calculate eigenvectors and -values
    vector3 eigenvals;
    matrix3x3 eigenmatrix = somematrix.findEigenvectorsIfSymmetric(eigenvals);
  
    // Print the 2nd eigenvector
    cout << eigenmatrix.GetColumn(1) << endl;
    \endcode
    With these conventions, a matrix is diagonalized in the following way:
    \code
    // Diagonalize the matrix
    matrix3x3 diagonalMatrix = eigenmatrix.inverse() * somematrix * eigenmatrix;
    \endcode
  
  */
  matrix3x3 matrix3x3::findEigenvectorsIfSymmetric(vector3 &eigenvals) const
#ifdef OB_OLD_MATH_CHECKS
  throw(OBError)
#endif
  {
    matrix3x3 result;

#ifdef OB_OLD_MATH_CHECKS
    if (!isSymmetric())
      {
        OBError er("matrix3x3::findEigenvectorsIfSymmetric(vector3 &eigenvals) const throw(OBError)",
                   "The method was called on a matrix that was not symmetric, i.e. where isSymetric() == false.",
                   "This is a runtime or a programming error in your application.");
        throw er;
      }
#endif

    double d[3];
    matrix3x3 copyOfThis = *this;

    jacobi(3, copyOfThis.ele[0], d, result.ele[0]);
    eigenvals.Set(d);

    return result;
  }

		BOHGE_FORCEINLINE void RotateAxis(const vector3<T>& axis, T r )
		{
			//T sa, ca;
			//Math::SinCos(r * T(0.5), sa, ca);
			T sc[2];
			Math::SinCos( r * T(0.5), sc );
			if ( Math::isEqual( axis.Length(), T(0.0) ) )
			{
				m_x = sc[0];
				m_y = sc[0];
				m_z = sc[0];
				m_w = sc[1];
			}
			else
			{
				vector3<T> temp = axis;
				temp.NormalizeSelf();
				temp *= sc[0];
				//*this = Quaternion<T>( temp * sc[0], sc[1] );
				m_x = temp.m_x;
				m_y = temp.m_y;
				m_z = temp.m_z;
				m_w = sc[1];
			}
		}

void uvaxes::uv2axes( vector3 & u,  vector3 &v, vector3 &e1, vector3 &e2, vector3 &e3 ){
  e2 = u^v;
  double au = u.mod();
  double av = v.mod();
  double ae2 = e2.mod();
  if (myequal(au,0.0) || myequal(av,0.0)){
    cerr << "Error: zero vector not usable for uvaxes";
    exit(1);
  }
  if (myequal(ae2,0.0)){
    cerr << "Error: colinear vectors not usable for uvaxes";
    exit(1);
  }
  e1 = u/au;
  e2 = e2 / ae2;
  e3 = e1^e2;
  return;
}

        static vector3 slerp(vector3 start, vector3 end, T percent) {
            T dot = start.dot(end);
            dot = vector3::clamp(dot, -1.0f, 1.0f);

            T theta = acos(dot) * percent;
            vector3 relative = end - start*dot;
            relative.normalize();
            return ((start * cos(theta)) + (relative*sin(theta)));
        }