subroutine TETRAQ(r,qg,qr)
c
c     To calculate the Glassmeier (QG) and Robert/Roux (QR) quality
c     factors for a tetrahedron.
c
c     Application: CLUSTER SCIENCE DATA SYSTEM
c                  (These are to be two auxiliary parameters)
c
c     Input:  R(3,4) = positions of the 4 points
c
c     Output: QG = Glassmeier factor
c             QR = Robert/Roux factor (their factor number 10)
c
c     Inputs and outputs are single precision, internal calculations
c       are double precision
c
c          True volume      True surface
c     QG = -----------  +   ------------  +   1.
c          Ideal vol        Ideal surf
c
c     where the ideals are volume and surface of a regular tetrahedron
c     of side length equal to the mean of the 6 sides.
c
c                   | True volume |(1/3)
c     QR = Factor * | ----------- |
c                   | Sphere vol  |
c
c     where the sphere is that circumscribing the tetrahedron (all 4
c     points on the surface), and the factor is such that QR=1
c     for a regular tetrahedron.
c
c*********************************************************************
c     Patrick W. Daly            daly@linmpi.mpg.de
c     Max-Planck-Institut        
c       fuer Aeronomie
c     D-37191 Katlenburg-Lindau
c     Germany
c                                  1994 June 7
c*********************************************************************
      implicit none
      real r(3,4),qg,qr
      double precision d(3,3),c(3,3),s1,s2,s3,s0,s,l1,l2,l3,vol
      double precision v(3),w,smean,vmean,lmean,rc
      double precision dot
      integer n,m,k
c
c     Find the differences
c
      do n=1,3
        w=dble(r(n,1))
        do m=1,3
          d(n,m)=dble(r(n,m+1))-w
        enddo      
      enddo
c
c     Find the average side length of all 6 sides
c
      l1=dsqrt(dot(d(1,1),d(1,1)))
      l2=dsqrt(dot(d(1,2),d(1,2)))
      l3=dsqrt(dot(d(1,3),d(1,3)))
      w= l1 + l2 + l3
      do n=1,3
        v(n)=d(n,2)-d(n,1)
      enddo
      w=w + dsqrt(dot(v,v))
      do n=1,3
        v(n)=d(n,3)-d(n,1)
      enddo
      w=w + dsqrt(dot(v,v))
      do n=1,3
        v(n)=d(n,3)-d(n,2)
      enddo
      w=w + dsqrt(dot(v,v))
      lmean=w/6.d0
c
c     Find the cross products
c
      do n=1,3
        m=mod(n,3) + 1
        k=mod(m,3) + 1
        call cross(d(1,m),d(1,k),c(1,n))
      enddo
c
c     Find the volume of the tetrahedron
c
      vol=dabs(dot(d(1,1),c(1,1)))/6.d0
c
c     Find the area of the 4 surfaces and their sum
c
      s1=0.5d0*dsqrt(dot(c(1,1),c(1,1)))      
      s2=0.5d0*dsqrt(dot(c(1,2),c(1,2)))      
      s3=0.5d0*dsqrt(dot(c(1,3),c(1,3)))      
      do n=1,3
        v(n)=c(n,1) + c(n,2) + c(n,3)
      enddo
      s0=0.5d0*dsqrt(dot(v,v))
      s=s0 + s1 + s2 + s3
c
c     Find the volume and total area for reg tetrahedron with mean side
c
      vmean=dsqrt(2.d0)*lmean*lmean*lmean/1.2d1
      smean=dsqrt(3.d0)*lmean*lmean
c
c     Calculate the Glassmeier factor
c
      w= vol/vmean + s/smean + 1.d0
      qg=sngl(w)
c
c     Find the center of the circumscribed circle
c
      w=1.d0/(1.2d1*vol)
      v(1)=w*(c(1,1)*l1*l1 + c(1,2)*l2*l2 + c(1,3)*l3*l3)
      v(2)=w*(c(2,1)*l1*l1 + c(2,2)*l2*l2 + c(2,3)*l3*l3)
      v(3)=w*(c(3,1)*l1*l1 + c(3,2)*l2*l2 + c(3,3)*l3*l3)
      rc=dsqrt(dot(v,v))
c
c     Calculate the Robert/Roux factor
c
      w=9.d0*dsqrt(3.d0)/8.d0
      w=(w*vol)**(1.d0/3.d0)/rc
      qr=sngl(w)
c
      return
      end

      function DOT(v,w)
c
c     To return the dot product of vectors V and W
c
      implicit none
      double precision dot,v(3),w(3)
      dot=v(1)*w(1) + v(2)*w(2) + v(3)*w(3)
      return
      end
      
      subroutine CROSS(v,w,x)
c
c     To return X as the cross product of vectors V and W
c
      implicit none
      double precision x(3),w(3),v(3)
      x(1)=v(2)*w(3) - v(3)*w(2)
      x(2)=v(3)*w(1) - v(1)*w(3)
      x(3)=v(1)*w(2) - v(2)*w(1)
      return
      end