Fortran code for finding the distance from a point to a plane.

distance2plane.f77

      subroutine distance2plane(x, y, z, flat)
!
!     Return the distance of a point to a plane. The plane is created by
!     the first three coordinates. The point of interest is the fourth.  
!
          real x(4), y(4), z(4) ! the fourth coordinate is the coordinate 
                                ! which is off of the plane
          real distance         ! distance of central atom from plane
          real threshhold       ! threshhold for testing planarity
          real vect1(3)         ! vector created by point 3 and point 1
          real vect2(3)         ! vector created by point 2 and point 1
          real normvect(3)      ! cross product of vector 1 and vector 2
          real magnormvect      ! magnitude of normal vector
          real unitvect(3)      ! unit vector of normal vector
          logical flat          ! true if below threshhold
          
          threshhold = 0.20e0   ! below this threshhold is considered flat
!
!     Create the two vectors which will be used to find the normal to
!     the plane.
!
          vect1 = [x(3)-x(1), y(3)-y(1), z(3)-z(1)]
          vect2 = [x(2)-x(1), y(2)-y(1), z(2)-z(1)]
!
!     Compute their normal (cross product).
!
          normvect = [vect1(2)*vect2(3) - vect1(3)*vect2(2),
     &                vect1(3)*vect2(1) - vect1(1)*vect2(3),
     &                vect1(1)*vect2(2) - vect1(2)*vect2(1)]
!     
!     Compute the magnitude of the normal.
!
          magnormvect = sqrt(dot_product(normvect, normvect))
!
!     Turn it into a unit vector.
!
          unitvect = normvect / magnormvect
!
!     The distance from the plane to the point is the projection of the
!     vector created between the point of interest off of the plane
!     and any point on the plane onto the unit normal vector. 
!     The first point in the plane [x(1), y(1), z(1)] is used.
!
          distance = dot_product(unitvect,
     &                           [x(4)-x(1),y(4)-y(1),z(4)-z(1)])
!
!     Don't want negative distances as an artifact of the directionality
!     of the unit vector
!
          distance = abs(distance)

          if (distance.le.threshhold) then
              flat = .true.
          else
              flat = .false.
          end if

          return
          
      end subroutine distance2plane