anchor.diff:  pnt(9) = m_enz(3,3)
   
   ! Eigenvalues are in the w array
  -pnt(10) = w(1)
  -pnt(11) = w(2)
  -pnt(12) = w(3)
  +pnt(10) = eigvals(1)
  +pnt(11) = eigvals(2)
  +pnt(12) = eigvals(3)
   
   #else
   
  -! Return zeros
  -write(stderr,*) 'mij2pnt: Hdef installed without LAPACK; returning pnt with zeros'
  -pnt = 0.0d0
  +! Calculate eigenvalues and eigenvectors with Numerical Recipes routine
  +call jacobi(m_enz,3,3,eigvals,eigvecs,nrot,ierr)
  +if (ierr.ne.0) then
  +    write(stderr,*) 'mij2pnt: error in jacobi; returning pnt with zeros'
  +    return
  +endif
  +
  +call eig_sort(eigvals,eigvecs,3,3)
  +
  +! Eigenvectors are in the mat array
  +pnt(1) = eigvecs(1,1)
  +pnt(2) = eigvecs(2,1)
  +pnt(3) = eigvecs(3,1)
  +pnt(4) = eigvecs(1,2)
  +pnt(5) = eigvecs(2,2)
  +pnt(6) = eigvecs(3,2)
  +pnt(7) = eigvecs(1,3)
  +pnt(8) = eigvecs(2,3)
  +pnt(9) = eigvecs(3,3)
  +
  +! Eigenvalues are in the vec array
  +pnt(10) = eigvals(1)
  +pnt(11) = eigvals(2)
  +pnt(12) = eigvals(3)
   
   #endif
   
  
positive.code: 
  use io, only: stderr
  
  implicit none
  
  ! Arguments
  double precision :: mrr, mtt, mpp, mrt, mrp, mtp, pnt(12)
  
  ! Local variables
  double precision :: m_enz(3,3)
  
  #ifdef USE_LAPACK
  
  integer :: info, lwork
  double precision :: w(3)
  double precision, allocatable :: work(:)
  
  #endif
  
  
negative.code: m_enz(1,2) = -mtp     ! mtp = -men
  m_enz(2,1) = -mtp
  
  
  ! Use LAPACK eigensystem solver if it is available. Otherwise, use Numerical Recipes routine.
  
  #ifdef USE_LAPACK
  
  ! Query the optimal workspace
  allocate(work(1))
  lwork = -1
  call dsyev('Vectors','Upper',3,m_enz,3,eigvals,work,lwork,info)
  lwork = int(work(1))
  deallocate(work)
  allocate(work(lwork))
  
  ! Solve eigenproblem
  call dsyev('Vectors','Upper',3,m_enz,3,eigvals,work,lwork,info)
  
  if (info.ne.0) then
      write(stderr,*) 'mij2pnt: error in dsyev; returning pnt with zeros'
  