c Zachary H. Levine    31 May 2001 -- 4 June 2001
c zlevine@nist.gov


c This code is in the public domain.  Do not remove comments.


c  This code is a product of the National Institute of Standards
c  and Technology.  As a courtesy, please acknowledge if used in
c  a publication or report and inform the author.  Thank you.


c  Test program for aneig.f  Expected answers then locally
c computed answers are given.


      implicit none


      integer n
      double precision pi
      parameter (n=10, pi=3.1415926535897932385d0)


      integer i, j
      double precision eigVal(2,n), eigVec(3,2,n), theta, phi, deg


      deg = 180.d0/pi


      i=0


      print 200
      print 100
      print 201
      print 202
      print 203
      print 204
      print 205
      print 206
      print 207
      print 208
      print 209
      print 210
 
      print 199
      print 100


c [001] direction
      i = i+1
      theta = 0.d0
      phi = 0.d0
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))


      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


c [100] direction
      i = i+1
      theta = pi/2
      phi = 0
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))
      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


c [010] direction
      i = i+1
      theta = pi/2
      phi = pi/2
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))
      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


c [111] direction
      i = i+1
      theta = acos(sqrt(1./3.))
      phi = pi/4
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))
      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


c   _
c [111] direction
      i = i+1
      theta = acos(sqrt(1./3.))
      phi = -pi/4
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))
      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


c [110] direction
      i = i+1
      theta = pi/2
      phi =  pi/4
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))
      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


c another direction
      i = i+1
      theta = 1.d0
      phi =  1.d0
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))
      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


c another direction
      i = i+1
      theta = -1.d0
      phi =  2.d0
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))
      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


c another direction
      i = i+1
      theta = 0.5
      phi =  0
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))
      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


c another direction
      i = i+1
      theta = 0.5
      phi =  pi/2
      call aneig (theta, phi, eigVal(1,i), eigVec(1,1,i))
      print 101, theta*deg, phi*deg, eigVal(1,i), (eigVec(j,1,i),j=1,3)
      print 102,                     eigVal(2,i), (eigVec(j,2,i),j=1,3)


      stop
  100 format('    Theta       Phi      EigVal      EigVec_x ',
     $       '   EigVec_y    EigVec_z ')
  101 format(2f10.5,f12.8,2x,3f12.8)
  102 format(20x,   f12.8,2x,3f12.8)
  199 format(//,'Locally computed output:')
  200 format('Expected output:')
  201 format('   0.00000   0.00000 -0.40000000',
     $       '    1.00000000  0.00000000  0.00000000',/,
     $       '                     -0.40000000',
     $       '    0.00000000  1.00000000  0.00000000')
  202 format('  90.00000   0.00000 -0.40000000',
     $       '    0.00000000  0.00000000 -1.00000000',/,
     $       '                     -0.40000000',
     $       '    0.00000000  1.00000000  0.00000000')
  203 format('  90.00000  90.00000 -0.40000000',
     $       '    0.00000000  0.00000000 -1.00000000',/,
     $       '                     -0.40000000 ',
     $       '  -1.00000000  0.00000000  0.00000000')
  204 format('  54.73561  45.00000  0.26666668 ',
     $       '   0.40824830  0.40824830 -0.81649657',/,
     $       '                      0.26666665 ',
     $       '  -0.70710678  0.70710678  0.00000000')
  205 format('  54.73561 -45.00000  0.26666668 ',
     $       '   0.40824830 -0.40824830 -0.81649657',/,
     $       '                      0.26666665 ',
     $       '   0.70710678  0.70710678  0.00000000')
  206 format('  90.00000  45.00000  0.60000000 ',
     $       '  -0.70710678  0.70710678  0.00000000',/,
     $       '                     -0.40000000 ',
     $       '   0.00000000  0.00000000 -1.00000000')
  207 format('  57.29578  57.29578  0.36967366 ',
     $       '  -0.29038716  0.69131851 -0.66162981',/,
     $       '                      0.07169012 ',
     $       '  -0.84200347  0.14391229  0.51992250')
  208 format(' -57.29578 114.59156  0.39937504 ',
     $       '   0.19118453  0.62306585  0.75844408',/,
     $       '                     -0.08539472 ',
     $       '   0.91696543  0.16229108 -0.36446673')
  209 format('  28.64789   0.00000  0.30807342 ',
     $       '   0.87758256  0.00000000 -0.47942554',/,
     $       '                     -0.40000000 ',
     $       '   0.00000000  1.00000000  0.00000000')
  210 format('  28.64789  90.00000  0.30807342 ',
     $       '   0.00000000  0.87758256 -0.47942554',/,
     $       '                     -0.40000000 ',
     $       '  -1.00000000  0.00000000  0.00000000')
      end