PROGRAM abc   

C     PROGRAM calculates principal A, B, C values             
C     from RAM A, B, C, Dab, Daci, Dbci values                       
C     f77 abc.f rs.f must be typed   
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(IDIM=3,NDIM=6)
      CHARACTER*4 TITLE(15)
      DIMENSION H(NDIM,NDIM),EI(NDIM,NDIM),E(NDIM),tmp(idim,idim)
      DIMENSION TOPCOS(IDIM),topangle(idim),test(idim,idim)
      dimension emhz(idim),fv1(idim),fv2(idim)
C
      do 10 nstruc=1,10
      READ(5,1001) (TITLE(it),it=1,15)
      if(title(1).eq.'-900') go to 999
      write(6,1002)
 1002 format(/)
      WRITE(6,1001)(TITLE(it),it=1,15)
      READ(5,1001) (TITLE(it),it=1,15)
 1001 FORMAT(15a4)
      write(6,1000)
 1000 format(/,' The RAM rotational constant matrix in abc order is',/)
      do 15 i=1,3
      READ(5,1010) (H(i,j),j=1,3)
 1010 format(3f15.11)
      write(6,1010) (H(i,j),j=1,3)
   15 continue
  
c     CALL DEVCSF(NDIM,H,IDIM,E,EI,IDIM)
c     call rs(nm,n,a,w,matz,z,fv1,fv2,ierr)
c      call rs(ndim,idim,h,e,1,ei,fv1,fv2,ierr)
       CALL TRED2(idim,ndim,H,E,WORK,EI)                              
         CALL TQL2(iDIM,ndim,E,WORK,EI,IERR)
c	CALL TRED2(NDIMTO,NDMX,H,EGVL,WORK,EGVC)                              
c         CALL TQL2(NDIMTO,NDMX,EGVL,WORK,EGVC,IERR)
cdebug  
c     write(6,9020) ierr
c9020 format(' ierr =',i5,'.  It should be 0.')
c     do 9000 i=1,3
c9000 write(6,1010) (H(i,j),j=1,3)
c end debug
c
c     change sign of a,b,c eigenvectors to get right-handed system
      if(ei(3,1).lt.0.d0) go to 21
      do 20 i=1,3
   20 ei(i,1)=-ei(i,1)
   21 if(ei(2,2).gt.0.d0) go to 23
      do 22 i=1,3
   22 ei(i,2)=-ei(i,2)
   23 if(ei(1,3).gt.0.d0) go to 25
      do 24 i=1,3
   24 ei(i,3)=-ei(i,3)
   25 continue
c
      write (6,1050)
 1050 format(/,' The column eigenvectors in C,B,A order are')
      do 30 i=1,3
   30 write(6,1010) (ei(i,j),j=1,3)
      WRITE(6,49)
 49   FORMAT(/,'EIGENVALUES = principal axis C,B,A in cm-1 and MHz')
      WRITE(6,50)(E(I), I=1,3)
 50   FORMAT(6F18.10)
      do 60 i=1,3
   60 emhz(i)=e(i)*29979.2458
      WRITE(6,50)(emhz(I), I=1,3)
C
c     check to see if eigenvectors are in columns
      do 70 i=1,3
      do 71 j=1,3
      tmp(i,j)=0.d0
      do 72 k=1,3
      tmp(i,j)=tmp(i,j)+H(i,k)*ei(k,j)
   72 continue
   71 continue
   70 continue
      do 80 i=1,3
      do 81 j=1,3
      test(i,j)=0.d0
      do 82 k=1,3
      test(i,j)=test(i,j)+ei(k,i)*tmp(k,j)
   82 continue
   81 continue
   80 continue
      write(6,1020)
 1020 format(/,' Utr*H*U below should agree with the eigenvalues above')
      do 90 i=1,3
      write(6,50) (test(i,j),j=1,3)
   90 continue

c     Begin computation of the top axis direction cosines in the PA system
c     using LAMPAM = normalization*BPAM**-1*RHOPAM = norm*BPAM**-1*Utr(.,1)
c     from notes EAA-I p. 62

      anorm=0.d0
      do 100 i=1,3
      topcos(i)=ei(1,i)/e(i)
      anorm=anorm+topcos(i)**2
  100 continue
      do 110 i=1,3
      topcos(i)=topcos(i)/dsqrt(anorm)
      topangle(i)=dacos(topcos(i))*180.d0/dacos(-1.d0)
  110 continue
      write(6,2000)
 2000 format(/,' The C,B,A direction cosines and angles of the top are')
      write(6,2010) (topcos(i),topangle(i),i=1,3)
 2010 format(' cos =',f12.8,'    angle =',f12.6,' degrees')
   10 continue
  999 STOP
      END