text: 
  
  /* Ray-Triangle Intersection Test Routines          */
  /* Different optimizations of my and Ben Trumbore's */
  /* code from journals of graphics tools (JGT)       */
  /* http://www.acm.org/jgt/                          */
  /* by Tomas Moller, May 2000                        */
  /* P.S: Modified for the Witch Engine */
  
  #include <math.h>
  
  #define EPSILON 0.000001
  #define CROSS(dest,v1,v2) \
            dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
            dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
            dest[2]=v1[0]*v2[1]-v1[1]*v2[0];
  #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
  #define SUB(dest,v1,v2) \
            dest[0]=v1[0]-v2[0]; \
            dest[1]=v1[1]-v2[1]; \
            dest[2]=v1[2]-v2[2]; 
  
  /* the original jgt code */
  int intersect_triangle(double orig[3], double dir[3],
  		       double vert0[3], double vert1[3], double vert2[3],
  		       double *t, double *u, double *v)
  {
     double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
     double det,inv_det;
  
     /* find vectors for two edges sharing vert0 */
     SUB(edge1, vert1, vert0);
     SUB(edge2, vert2, vert0);
  
     /* begin calculating determinant - also used to calculate U parameter */
     CROSS(pvec, dir, edge2);
  
     /* if determinant is near zero, ray lies in plane of triangle */
     det = DOT(edge1, pvec);
  
     if (det > -EPSILON && det < EPSILON)
       return 0;
     inv_det = 1.0 / det;
  
     /* calculate distance from vert0 to ray origin */
     SUB(tvec, orig, vert0);
  
     /* calculate U parameter and test bounds */
     *u = DOT(tvec, pvec) * inv_det;
     if (*u < 0.0 || *u > 1.0)
       return 0;
  
     /* prepare to test V parameter */
     CROSS(qvec, tvec, edge1);
  
     /* calculate V parameter and test bounds */
     *v = DOT(dir, qvec) * inv_det;
     if (*v < 0.0 || *u + *v > 1.0)
       return 0;
  
     /* calculate t, ray intersects triangle */
     *t = DOT(edge2, qvec) * inv_det;
  
     return 1;
  }
  
  
  /* code rewritten to do tests on the sign of the determinant */
  /* the division is at the end in the code                    */
  int intersect_triangle1(double orig[3], double dir[3],
  			double vert0[3], double vert1[3], double vert2[3],
  			double *t, double *u, double *v)
  {
     double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
     double det,inv_det;
  
     /* find vectors for two edges sharing vert0 */
     SUB(edge1, vert1, vert0);
     SUB(edge2, vert2, vert0);
  
     /* begin calculating determinant - also used to calculate U parameter */
     CROSS(pvec, dir, edge2);
  
     /* if determinant is near zero, ray lies in plane of triangle */
     det = DOT(edge1, pvec);
  
     if (det > EPSILON)
     {
        /* calculate distance from vert0 to ray origin */
        SUB(tvec, orig, vert0);
        
        /* calculate U parameter and test bounds */
        *u = DOT(tvec, pvec);
        if (*u < 0.0 || *u > det)
  	 return 0;
        
        /* prepare to test V parameter */
        CROSS(qvec, tvec, edge1);
        
        /* calculate V parameter and test bounds */
        *v = DOT(dir, qvec);
        if (*v < 0.0 || *u + *v > det)
  	 return 0;
        
     }
     else if(det < -EPSILON)
     {
        /* calculate distance from vert0 to ray origin */
        SUB(tvec, orig, vert0);
        
        /* calculate U parameter and test bounds */
        *u = DOT(tvec, pvec);
  /*      printf("*u=%f\n",(float)*u); */
  /*      printf("det=%f\n",det); */
        if (*u > 0.0 || *u < det)
  	 return 0;
        
        /* prepare to test V parameter */
        CROSS(qvec, tvec, edge1);
        
        /* calculate V parameter and test bounds */
        *v = DOT(dir, qvec) ;
        if (*v > 0.0 || *u + *v < det)
  	 return 0;
     }
     else return 0;  /* ray is parallell to the plane of the triangle */
  
  
     inv_det = 1.0 / det;
  
     /* calculate t, ray intersects triangle */
     *t = DOT(edge2, qvec) * inv_det;
     (*u) *= inv_det;
     (*v) *= inv_det;
  
     return 1;
  }
  
  /* code rewritten to do tests on the sign of the determinant */
  /* the division is before the test of the sign of the det    */
  int intersect_triangle2(double orig[3], double dir[3],
  			double vert0[3], double vert1[3], double vert2[3],
  			double *t, double *u, double *v)
  {
     double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
     double det,inv_det;
  
     /* find vectors for two edges sharing vert0 */
     SUB(edge1, vert1, vert0);
     SUB(edge2, vert2, vert0);
  
     /* begin calculating determinant - also used to calculate U parameter */
     CROSS(pvec, dir, edge2);
  
     /* if determinant is near zero, ray lies in plane of triangle */
     det = DOT(edge1, pvec);
  
     /* calculate distance from vert0 to ray origin */
     SUB(tvec, orig, vert0);
     inv_det = 1.0 / det;
     
     if (det > EPSILON)
     {
        /* calculate U parameter and test bounds */
        *u = DOT(tvec, pvec);
        if (*u < 0.0 || *u > det)
  	 return 0;
        
        /* prepare to test V parameter */
        CROSS(qvec, tvec, edge1);
        
        /* calculate V parameter and test bounds */
        *v = DOT(dir, qvec);
        if (*v < 0.0 || *u + *v > det)
  	 return 0;
        
     }
     else if(det < -EPSILON)
     {
        /* calculate U parameter and test bounds */
        *u = DOT(tvec, pvec);
        if (*u > 0.0 || *u < det)
  	 return 0;
        
        /* prepare to test V parameter */
        CROSS(qvec, tvec, edge1);
        
        /* calculate V parameter and test bounds */
        *v = DOT(dir, qvec) ;
        if (*v > 0.0 || *u + *v < det)
  	 return 0;
     }
     else return 0;  /* ray is parallell to the plane of the triangle */
  
     /* calculate t, ray intersects triangle */
     *t = DOT(edge2, qvec) * inv_det;
     (*u) *= inv_det;
     (*v) *= inv_det;
  
     return 1;
  }
  
  /* code rewritten to do tests on the sign of the determinant */
  /* the division is before the test of the sign of the det    */
  /* and one CROSS has been moved out from the if-else if-else */
  int intersect_triangle3(double orig[3], double dir[3],
  			double vert0[3], double vert1[3], double vert2[3],
  			double *t, double *u, double *v)
  {
     double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
     double det,inv_det;
  
     /* find vectors for two edges sharing vert0 */
     SUB(edge1, vert1, vert0);
     SUB(edge2, vert2, vert0);
  
     /* begin calculating determinant - also used to calculate U parameter */
     CROSS(pvec, dir, edge2);
  
     /* if determinant is near zero, ray lies in plane of triangle */
     det = DOT(edge1, pvec);
  
     /* calculate distance from vert0 to ray origin */
     SUB(tvec, orig, vert0);
     inv_det = 1.0 / det;
     
     CROSS(qvec, tvec, edge1);
        
     if (det > EPSILON)
     {
        *u = DOT(tvec, pvec);
        if (*u < 0.0 || *u > det)
  	 return 0;
              
        /* calculate V parameter and test bounds */
        *v = DOT(dir, qvec);
        if (*v < 0.0 || *u + *v > det)
  	 return 0;
        
     }
     else if(det < -EPSILON)
     {
        /* calculate U parameter and test bounds */
        *u = DOT(tvec, pvec);
        if (*u > 0.0 || *u < det)
  	 return 0;
        
        /* calculate V parameter and test bounds */
        *v = DOT(dir, qvec) ;
        if (*v > 0.0 || *u + *v < det)
  	 return 0;
     }
     else return 0;  /* ray is parallell to the plane of the triangle */
  
     *t = DOT(edge2, qvec) * inv_det;
     (*u) *= inv_det;
     (*v) *= inv_det;
  
     return 1;
  }