subroutine plane(x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3,x,y,z,r)
c
c  Subroutine plane computes the intersection (x,y,z) of a plane 
c  and a perpendicular line.  The plane is defined by three points 
c  (x1,y1,z1), (x2,y2,z2), and (x3,y3,z3).  The line passes through 
c  (x0,y0,z0).  Computation is done by a transformation and inverse 
c  transformation of coordinates systems.
c
      x2n=x2-x1
      y2n=y2-y1
      z2n=z2-z1
      x0n=x0-x1
      y0n=y0-y1
      z0n=z0-z1
      x3n=x3-x1
      y3n=y3-y1
      z3n=z3-z1
      call cross(x3n,y3n,z3n,x2n,y2n,z2n,cx,cy,cz,c)
      call cross(x2n,y2n,z2n,cx,cy,cz,dx,dy,dz,d)
      a=sqrt(x2n**2+y2n**2+z2n**2)
      t11=x2n/a
      t12=y2n/a
      t13=z2n/a
      t21=cx/c
      t22=cy/c
      t23=cz/c
      t31=dx/d
      t32=dy/d
      t33=dz/d
      tx0=t11*x0n+t12*y0n+t13*z0n
      tz0=t31*x0n+t32*y0n+t33*z0n
      r=t21*x0n+t22*y0n+t23*z0n
      x=t11*tx0+t31*tz0
      y=t12*tx0+t32*tz0
      z=t13*tx0+t33*tz0
      x=x+x1
      y=y+y1
      z=z+z1
      return
      end