subroutine kutl_intersect_2lines(x1, eta1, x2, eta2, xint)
*
* begin_doc
*
*  Returns intersection point in 3-d of two lines.
*
*  Input Parameters:
*  x1(3)          double precision array
*                 (x,y,z) of initial point defining line 1
*
*  eta1(3)        double precision array
*                 direction cosines of line 1 (do not have to be normalized)
*
*  x2(3)          double precision array
*                 (x,y,z) of initial point defining line 2
*
*  eta2(3)        double precision array
*                 direction cosines of line 2 (do not have to be normalized)
*
*  Output Parameters:
*  xint(3)        double precision array
*                 (x,y,z) intersection point. If cannot be determined, it is
*                 set to (0,0,0).
*
*  Other routines:
*
*  Notes:
*
*  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*  The lines are defined by
*
*                  x = x1 + eta1*r1
*                  x = x2 + eta2*r2
*
*    where r1, r2 are free parameters and eta1 and eta2 do not have to
*    be normalized to 1. The intersection point in 3-d is defined as
*    the midpoint between the lines at the distance of closest
*    approach. The intersection point is calculated from
*
*            xint(i) =(1/2) [x1(i) + x2(i)+ eta1(i)*r1 + eta2(i)*r2]
*
*                          delta*eta1 - (eta1*eta2)(delta*eta2)
*            where r1 =  --------------------------------------
*                                 1 - (eta1*eta2)**2
*
*                         -delta*eta2 + (eta1*eta2)(delta*eta1)
*                  r2 =  --------------------------------------
*                                 1 - (eta1*eta2)**2
*
*            delta(i) = x2(i) - x1(i)
*
*  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*  Author:   Paul R Avery      Created:  Sun Nov 16 18:54:29 EST 1997
*
*  Major revisions:
*
*
* end_doc