c$debug
C-----------------------------------------------------------------------------
c
c   QDIAG - Diagonalization of the inertial quadrupole tensor
c
C-----------------------------------------------------------------------------
c
c      All elements of the inertial tensor can be specified and the program
c      will calculate:
c
c      - principal tensor components and their error
c      - direction cosines and their errors
c      - rotation angles and their errors 
c
c      The results are identical to those obtained for the 2x2 block
c      using the tan(2 theta)=2 chiab/(chiaa-chibb) relation for the rotation
c      angle theta.za
c
C       ver. 12.II.2023                           ----- Zbigniew KISIEL -----
C                          __________________________________________________
C                         | Institute of Physics, Polish Academy of Sciences |
C                         | Al.Lotnikow 32/46, Warszawa, POLAND              |
C                         |                             kisiel@ifpan.edu.pl  |
C                         |     http://info.ifpan.edu.pl/~kisiel/prospe.htm  |
C_________________________/--------------------------------------------------
C
C
c      Default data file is QDIAG.INP
c
c      Output is to default output device (screen) - use pipeline QDIAG>OUT
c      etc. to get hardcopy
c
c      Use data file for template of input, only the first record will
c      be read in.  First twenty characters of data line are skipped, if there
c      are problems on reading the floating point calues then separate them
c      with commas.
c
c      The two control variables in the third line define:
c
c      ikappa = -1 for prolate prolate top, +1 for oblate top (this is only
c               used if idatyp=1
c      idatyp =  1 for the two determinable constants taken directly from
c               output of Pickett's fitting program (ie 1.5 chi.aa and
c               1/4(chibb-chicc) for prolate, or 1.5 chi.cc and
c               1/4(chibb-chiaa) for oblate)
c             =  2 for specification of diagonal terms as chiaa, chibb, chicc
c
c
C  Modification history:
C
c   5.07.97:  creation, partially from QTENS
c   1.01.98:  addition of SPFIT diagonal constants to output
c  12.02.23:  eliminating features deleted in Fortran2018 or earlier
c   
C-----------------------------------------------------------------------------
c
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*2 (I-N)
c
c      qin - quadrupole tensor in inertial axes
c    erqin - errors in elements of QIN
c      eig - the three nonzero components of the principal inertial tensor
c    ereig - errors in above
c   eigvec - eigenectors for inertial->principal quadrupole transformation
c            containing direction cosines of interaxial angles
c            (columnwise)
c   ereigv - errors in above
c      der - numerical derivatives of principal tensor components with
c            respect to all six determinable inertial tensor components
c   dereig - numerical derivatives of all direction cosines relative to
c            inertial tensor components
c      ang - rotation angles derived from direction cosines
c    erang - errors inthe above
C      map - i'th element defines whether i'th principal tensor eigenvalue
c            is x,y or z (1,2 or 3), this is to be provided on input
c   mapinv - elements 1,2,3 contain principal tensor eigenvalue numbers
c            for x,y and z resp.
c
      character cdummy*20,axis(3),paxis(3),coment*79,pidiag*10
      integer map(3),mapinv(3)
      real*8 qin(3,3),erqin(3,3),eig(3),eigvec(3,3),
     *   dere(3,6),dereig(3,3,6),ereig(3),ereigv(3,3),
     *   ang(3,3),erang(3,3)
      equivalence (qin(1,1),chiaa),(qin(2,2),chibb),(qin(3,3),chicc)
      equivalence (erqin(1,1),echiaa),(erqin(2,2),echibb),
     *  (erqin(3,3),echicc)
      equivalence (qin(1,2),chiab),(qin(1,3),chiac),(qin(2,3),chibc)
      equivalence (erqin(1,2),echiab),(erqin(1,3),echiac),
     *  (erqin(2,3),echibc)
c
      COMMON /HBLK/H(3,3)/RBLK/S(3,3)
      data axis/'a','b','c'/
      data paxis/'x','y','z'/
c
      conv=360.d0/(2.d0*3.14159265359d0)
c
      WRITE(*,3344)
3344  FORMAT('  ',76(1H_)/' |',T79,'|'/
     * ' |   Q D I A G  -  DIAGONALIZATION OF THE PRINCIPAL QUADRUPOLE ',
     * ' TENSOR',T79,'|'/
     * ' |',76(1H_),'|'/'  version 12.II.2023',T64,'Zbigniew KISIEL'/)
c
c
      open(1,file='qdiag.inp',status='old')
c
c...ikappa=-1 prolate, +1 oblate
c...idatyp= 1 SPFIT diagonal constants
c           2 Chiaa,Chibb,Chicc
c
        read(1,'(a)')coment
        write(*,'(1x/1x,a/)')coment
        read(1,'(a)')cdummy
        read(1,'(a,2i15)')cdummy,ikappa,idatyp
c
c...diagonals
c
      if(idatyp.eq.1)then
        read(1,2)cdummy,chi32,echi32
        read(1,2)cdummy,chi14,echi14
2       format(a,2f18.9)
        if(ikappa.eq.-1)then
          chiaa=chi32/1.5d0
          echiaa=echi32/1.5d0
          chbmc=chi14*4.d0
          echbmc=echi14*4.d0
          chicc=-(chiaa+chbmc)/2.d0
          echicc=dsqrt(echiaa**2+echbmc**2)/2.d0
          chibb=(chbmc-chiaa)/2.d0
          echibb=dsqrt(echiaa**2+echbmc**2)/2.d0
        else
          chicc=chi32/1.5d0
          echicc=echi32/1.5d0
          chbma=chi14*4.d0
          echbma=echi14*4.d0
          chiaa=-(chicc+chbma)/2.d0
          echiaa=dsqrt(echicc**2+echbma**2)/2.d0
          chibb=(chbma-chicc)/2.d0
          echibb=dsqrt(echicc**2+echbma**2)/2.d0
        endif
      endif
c
      if(idatyp.eq.2)then
        read(1,2)cdummy,chiaa,echiaa
        read(1,2)cdummy,chibb,echibb
        read(1,2)cdummy,chicc,echicc 
      endif
c
c...off-diagonals
c
      read(1,2)cdummy,chiab,echiab
      read(1,2)cdummy,chiac,echiac
      read(1,2)cdummy,chibc,echibc
c
c...xyz mapping
c
      read(1,'(a,3i15)')cdummy,map
c
      close(1)
c
      qin(2,1)=qin(1,2)
      qin(3,1)=qin(1,3)
      qin(3,2)=qin(2,3)
      erqin(2,1)=erqin(1,2)
      erqin(3,1)=erqin(1,3)
      erqin(3,2)=erqin(2,3)
c
      do 1 i=1,3
        do 101 j=1,3
          h(i,j)=qin(i,j)
101     continue
1     continue
c
c
      write(*,'('' Inertial tensor:''/1x,t14,''a'',t38,''b'',
     * t62,''c'')')
      do 6 i=1,3
        write(*,7)axis(i),(qin(i,j),erqin(i,j),j=1,3)
7       format(1x,a,3(f13.5,' +-',f8.5))
6     continue
c
      call jack(int2(3))
c
      write(*,'(1x)')     
      do 9 i=1,3
        eig(i)=h(i,i)
        do 10 j=1,3
          eigvec(j,i)=s(j,i)
10      continue
9     continue
c
c...Derivatives of eigenvalues wrt all quadrupole terms: numeric evaluation
c   by incrementing each quantity in turn with the increment DERINC
c
c   NDER=1,2,3,4,5,6 map to 1,1 1,2 1,3 2,2 2,3 3,3 locations in 3x3 arrays
c   ie. NDER=5 points to chibc etc.
c
      nder=0
      derinc=0.001d0
c
      do 4 i=1,3
        do 104 j=i,3
          nder=nder+1
c
          do 5 k=1,3
            do 105 l=k,3
              h(k,l)=qin(k,l)
              if(k.eq.i.and.j.eq.l)h(k,l)=h(k,l)+derinc
              if(k.ne.l)h(l,k)=h(k,l)
105         continue              
5         continue
c
          call jack(int2(3))
c
          do 11 k=1,3
            dere(k,nder)=(h(k,k)-eig(k))/derinc
            do 111 l=1,3
              dereig(l,k,nder)=(s(l,k)-eigvec(l,k))/derinc
111         continue
11        continue
c
104     continue
4     continue
c
c...errors in eigenvalues, using standard propagation
c
      do 14 i=1,3
        ereig(i)=0.d0
        do 114 j=1,3
          ereigv(i,j)=0.d0
114     continue          
14    continue
c
      do 15 l=1,3
        nder=0
        do 115 i=1,3
          do 215 j=i,3
            nder=nder+1
            ereig(l)=ereig(l)+dere(l,nder)**2*erqin(i,j)**2
215      continue
115     continue
15    continue
c
      do 16 i=1,3
        ereig(i)=dsqrt(ereig(i))
16    continue
c
c...errors in eigenvectors, using standard propagation
c
      do 20 jj=1,3
        do 120 ii=1,3
          nder=0
          do 220 i=1,3
            do 320 j=i,3
              nder=nder+1
              ereigv(ii,jj)=ereigv(ii,jj)+
     *          dereig(ii,jj,nder)**2*erqin(i,j)**2
320         continue
220       continue
120     continue
20    continue
c
c...errors in angles:  d( cos-1 (cos) )/d cos = -1/sqrt (1-cos**2)
c
      do 21 i=1,3
        do 121 j=1,3
          ereigv(i,j)=dsqrt(ereigv(i,j))        
          ang(i,j)=conv*dacos(eigvec(i,j))
          erang(i,j)= ereigv(i,j)**2 / (1.D0-eigvec(i,j)**2+1.d-14) 
          erang(i,j)=conv*dsqrt(erang(i,j))
121     continue
21    continue
c
c...output
c
      if(idatyp.eq.1)then
        write(*,'('' Principal tensor components: '',10x,
     *   ''SPFIT diagonal components:'')')
      else
        write(*,'('' Principal tensor components: '')')
      endif
      do 17 i=1,3
        if(idatyp.ne.1)then
          write(*,7)paxis(map(i)),eig(i),ereig(i)
        else
          if(i.eq.1)write(*,7)paxis(map(i)),eig(i),ereig(i)
          if(i.eq.2)then
            pidiag='3/2 chi.aa'
            if(ikappa.gt.0)pidiag='3/2 chi.cc'
            write(*,77)paxis(map(i)),eig(i),ereig(i),pidiag,chi32,echi32
          endif
          if(i.eq.3)then
            pidiag='1/4(bb-cc)'
            if(ikappa.gt.0)pidiag='1/4(bb-aa)'
            write(*,77)paxis(map(i)),eig(i),ereig(i),pidiag,chi14,echi14
          endif
        endif
        mapinv(map(i))=i
17    continue
77    format(1x,a,f13.5,' +-',f8.5,14x,a10,f13.5,' +-',f8.5)
c
      write(*,'(1x/'' Direction cosines:''/1x,t14,''a'',t38,''b'',
     * t62,''c'')')
      do 25 i=1,3
        write(*,8)paxis(map(i)),(eigvec(j,i),ereigv(j,i),j=1,3)
8       format(1x,a,3(f13.6,' +-',f8.6))
25    continue
c
      write(*,'(1x/'' Rotation angles:''/1x,t14,''a'',t38,''b'',
     * t62,''c'')')
      do 26 i=1,3
        write(*,22)paxis(map(i)),(ang(j,i),erang(j,i),j=1,3)
22      format(1x,a,3(f13.4,' +-',f8.4))
26    continue
c
c...eta =  (chixx-chiyy)/chizz
c
      chixx=eig(mapinv(1))
      chiyy=eig(mapinv(2))
      chizz=eig(mapinv(3))
      echixx=ereig(mapinv(1))
      echiyy=ereig(mapinv(2))
      echizz=ereig(mapinv(3))
c
      eta=  (chixx-chiyy)/chizz
      ereta=
     $ DSQRT(chixx**2*echizz**2-2.d0*chixx*chiyy*echizz**2
     $ +chiyy**2*echizz**2+chizz**2*(echixx**2+echiyy**2))
     $ /chizz**2
      write(*,'(1x)')
      if(ereta.le.0.0001d0)then
        write(*,37)'Eta = (chixx-chiyy)/chizz = ',eta,ereta
        stop
      endif
      if(ereta.le.0.00001d0)then
        write(*,38)'Eta = (chixx-chiyy)/chizz = ',eta,ereta
        stop
      endif       
      write(*,36)'Eta = (chixx-chiyy)/chizz = ',eta,ereta
36    format(1x,a,t30,f11.5,' +-',f8.5)
37    format(1x,a,t30,f12.6,' +-',f9.6)
38    format(1x,a,t30,f13.7,' +-',f10.7)
c
      stop
      end
C_____________________________________________________________________________
C
      SUBROUTINE JACK(N)
c
c   JACOBI type matrix diagonalization routine for real, symmetric matrices
c   (optimized version of AUTOValori from D.Lister's LINSHP program)
c
c     H - array to be diagonalized (contains the eigenvalues on the diagonal
c         on output
c     S - contains the eigenvectors on output (in order of eigenvalues)
c     N - dimension of the array to be diagonalized contained in the top
c         left hand corner of H
c
c
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*2 (I-N)
      PARAMETER (MAXDIM=3)
C
      COMMON /HBLK/H/RBLK/S
      REAL*8 H(MAXDIM,MAXDIM),S(MAXDIM,MAXDIM)
c
c ... initialization and setup of limits
c
      ss=0.7071067811865475d0
      n1=n-1
      sum=0.d0
      do 20 i=1,n1
        s(i,i)=1.d0
        k=i+1
        do 120 m=k,n
          sum=sum+h(i,m)**2
          s(i,m)=0.d0
          s(m,i)=0.d0
120     continue
20    continue
      s(n,n)=1.d0
      unu=dsqrt(2.d0*sum)
      ene=n
c
c ... convergence criterion
c
      unufi=1.d-12
c
c
c ... iteration sequence - determination of SIN and COS for rotation
c
54    unu=unu/ene
55    ind=0
      do 89 m=2,n
        l=m-1
        do 189 i=1,l
          if(h(i,m).eq.0.0d0)then
            goto 189
          else
            goto 60
          endif
60        q=dabs(h(i,m))
          if(q-unu.lt.0.d0)then
            goto 189
          else
            goto 62
          endif
62        ind=1
          if(h(i,i)-h(m,m).eq.0.d0)then
            goto 64
          else
            goto 68
          endif
64        si=dsign(ss,h(i,m))
          co=ss
          goto 77
68        x1=-h(i,m)*2.d0/(h(i,i)-h(m,m))
          x1a=dabs(x1)
c          
          if(1.d30-x1a.le.0.d0)then
            goto 66
          else
            goto 71
          endif
c
66        si=dsign(ss,x1)
          co=ss
          goto 77
71        if(1.d-30-x1a.lt.0.d0)then
            goto 72
          else
            goto 189
          endif
72        t=datan(x1)*0.5d0
          si=dsin(t)
          co=dcos(t)
c
c ... transformation loop
c
77        do 79 j=1,n
            hji=h(j,i)
            sji=s(j,i)
            h(j,i)=hji*co-h(j,m)*si
            h(j,m)=hji*si+h(j,m)*co
            s(j,i)=sji*co-s(j,m)*si
            s(j,m)=sji*si+s(j,m)*co
            if(j.eq.i.or.j.eq.m)goto 79
            h(i,j)=h(j,i)
            h(m,j)=h(j,m)
79        continue
c
          h(i,i)=co*h(i,i)-si*h(m,i)
          h(m,m)=co*h(m,m)+si*h(i,m)
          h(i,m)=0.d0
          h(m,i)=0.d0
189     continue
89    continue
      if(ind.gt.0.d0)then
        goto 55
      else
        if(ind.eq.0.d0)goto 91
        goto 92
      endif
c
91    if(unu-unufi.le.0.d0)then
        goto 92
      else
        goto 54
      endif
c
92    return
      end
C
C_____________________________________________________________________________
C_____________________________________________________________________________