subroutine line(x0,y0,z0,x1,y1,z1,x2,y2,z2,x,y,z,v1,v2,r)
c
c  Subroutine LINE determines the intersection (x,y,z) of two 
c  lines.  First line is defined by points (x1,y1,z1) and 
c  (x2,y2,z2).  Second line is perpendicular to the first and 
c  passes through point (x0,y0,z0).  Distance between (x,y,z) 
c  and (x0,y0,z0) is returned as r.  Computation is done by a 
c  transformation of coordinate systems.
c
      tx0=x0-x1
      ty0=y0-y1
      tz0=z0-z1
      tx2=x2-x1
      ty2=y2-y1
      tz2=z2-z1
      a=sqrt(tx2**2+ty2**2+tz2**2)
      call cross(tx2,ty2,tz2,tx0,ty0,tz0,cx,cy,cz,c)
      call cross(cx,cy,cz,tx2,ty2,tz2,dx,dy,dz,d)
      tt11=tx2/a
      tt12=ty2/a
      tt13=tz2/a
      tt21=dx/d
      tt22=dy/d
      tt23=dz/d
      tt31=cx/c
      tt32=cy/c
      tt33=cz/c
      u0=tt11*tx0+tt12*ty0+tt13*tz0
      r=tt21*tx0+tt22*ty0+tt23*tz0
      x=tt11*u0+x1
      y=tt12*u0+y1
      z=tt13*u0+z1
      v1=-u0
      v2=a-u0
      return
      end