c
c      QTRANS.f     Version, Feb 2019
c
c
c      Please cite as:  W. C. Bailey, Calculation of Nuclear
c      Quadrupole Coupling Constants in Gaseous State Molecules,
c      http://nqcc.wcbailey.net/index.html
c
c
c      Transforms qij and mu from Gaussian standard orientation axes to
c      inertial tensor a,b,c principal axes, to nqcc tensor x,y,z
c      principal axes.
c
c      INPUT file is qtrans.in, OUTPUT is qtrans.out
c     
c      INPUT: Line 1   Number of atoms, N
c             Lines 2,N+1   Isotopic mass; atomic x,y,z coordinates
c             Line N+2   mu_x, mu_y, mu_z, mu
c             Line N+3   qxx,qyy,qzz
c             Line N+4   qxy,qxz,qyz
c             Line N+5   eqQeff for conversion of qij to Xij
c
c      This program owes its existence to Dr. Zbigniew Kisiel for making
c      publicly available his version (JACK) of the jacobi matrix
c      diagonalization routine.  Visit http://info.ifpan.edu.pl/~kisiel/
c      rotlinks.htm.  See QDIAG.f.
c
c
c
       IMPLICIT REAL*8 (A-H,O-Z)
       IMPLICIT INTEGER (I-N)
c
       COMMON /HBLK/H(3,3)/RBLK/S(3,3)
       REAL*8 W(50),X(50),Y(50),Z(50),A(50),B(50),C(50)
       REAL*8 F(3,3),G(3,3),d(3),v(3,3),Aj(50),Bj(50),Cj(50)
       REAL*8 Q(3,3),aco(3,3),HO(3,3),SO(3,3)
       OPEN (11, FILE='qtrans.in')
       OPEN (12, FILE='qtrans.out')
       PI = (4.)*ATAN(1.)
    1  FORMAT(F12.6)
    2  FORMAT(3F12.6)
    3  FORMAT(4F12.6)
       WRITE (12,6)
    6  FORMAT (/' ------------------ INPUT ----------------'/)
       WRITE (12,10)
   10  FORMAT(' Isotopic mass and Gaussian x,y,z coordinates'/)
       READ(11,*) N
       DO 20 I=1,N
       READ (11,*) W(I),X(I),Y(I),Z(I)
   20  WRITE (12,3) W(I),X(I),Y(I),Z(I)
       READ (11,*) dipx, dipy, dipz, dip
       X(N+1)=dipx
       Y(N+1)=dipy
       Z(N+1)=dipz
       READ(11,*) qxx,qyy,qzz
       READ(11,*) qxy,qxz,qyz
       Q(1,1)=qxx
       Q(1,2)=qxy
       Q(1,3)=qxz
       Q(2,1)=qxy
       Q(2,2)=qyy
       Q(2,3)=qyz
       Q(3,1)=qxz
       Q(3,2)=qyz
       Q(3,3)=qzz
       WRITE(12,4)
    4  FORMAT(/' Electric Dipole Moments x, y, z, total'/)
       WRITE(12,3) X(N+1),Y(N+1),Z(N+1),dip
       WRITE(12,21)
   21  FORMAT(/' Gaussian Electric Field Gradient Tensor'/)
       WRITE(12,11)
   11  FORMAT('         x           y           z'/)
       WRITE(12,12) Q(1,1),Q(1,2),Q(1,3)
   12  FORMAT(' x', 3F12.6)
       WRITE(12,13) Q(2,1),Q(2,2),Q(2,3)
   13  FORMAT(' y', 3F12.6)
       WRITE(12,14) Q(3,1),Q(3,2),Q(3,3)
   14  FORMAT(' z', 3F12.6)
       READ(11,*) eQq
       WRITE(12,15) eQq
   15  FORMAT(/' Xij = eQeff/h*(-qij) where eQeff/h =', F12.6)
c
c      Compute center of mass
c
       XC = 0.
       YC = 0.
       ZC = 0.
       AM = 0.
       DO 30 I=1,N
       AM = AM + W(I)
       XC = XC + W(I)*X(I)
       YC = YC + W(I)*Y(I)
   30  ZC = ZC + W(I)*Z(I)
       XCM = XC/AM
       YCM = YC/AM
       ZCM = ZC/AM
       WRITE(12,7)
    7  FORMAT (//' ------------------ OUTPUT -------------------')
c      WRITE (12,35)
c  35  FORMAT(/'Coordinates of center of mass'/)
c      WRITE (12,2) XCM, YCM, ZCM
c
c      Translate to center of mass coordinates, center of mass at origin.
c
       DO 40 I=1,N
       X(I) = X(I) - XCM
       Y(I) = Y(I) - YCM
   40  Z(I) = Z(I) - ZCM
c
c    X,Y,Z are now center of mass atomic coordinates.
c
c      WRITE (12,45)
c  45  FORMAT (/'Atomic coordinates in center of mass system'/)
c      DO 50 I=1,N
c  50  WRITE (12,3) X(I),Y(I),Z(I)
c
c      Compute inertia tensor, F
c
       F(1,1)=0.
       F(1,2)=0.
       F(1,3)=0.
       F(2,2)=0.
       F(2,3)=0.
       F(3,3)=0.
       DO 55 I=1,N
       F(1,1) = F(1,1) + W(I)*(Y(I)*Y(I) + Z(I)*Z(I))
       F(1,2) = F(1,2) - W(I)*X(I)*Y(I)
       F(1,3) = F(1,3) - W(I)*X(I)*Z(I)
       F(2,2) = F(2,2) + W(I)*(X(I)*X(I) + Z(I)*Z(I))
       F(2,3) = F(2,3) - W(I)*Y(I)*Z(I)
   55  F(3,3) = F(3,3) + W(I)*(X(I)*X(I) + Y(I)*Y(I))
       F(2,1) = F(1,2)
       F(3,1) = F(1,3)
       F(3,2) = F(2,3)
c      WRITE (12,58)
c  58  FORMAT(/' Inertia tensor: Iij where i,j=x,y,z'/)
c      DO 60 I=1,3
c  60  WRITE (12,3) F(I,1), F(I,2), F(I,3)
c
       DO 56 I=1,3
       DO 56 J=1,3
   56  G(I,J)=F(I,J)
       DO 62 I=1,3
       DO 62 J=1,3
   62  H(I,J)=F(I,J)
c
       Call JACK(3)
       Call ORDER(H,S,HO,SO)
       DO 80 i=1,3
       H(i,i)=HO(i,i)
       DO 80 j=1,3
   80  S(i,j)=SO(i,j)
c
c      WRITE(12,57)
c  57  FORMAT(/' Eigenvalues: Iii where i=a,b,c (JACK)'/)
c      WRITE(12,3) H(1,1),H(2,2),H(3,3)
c      WRITE(12,59)
c  59  FORMAT(/' Eigenvectors, direction cosines (JACK)'/)
c      DO 61 I=1,3
c  61  WRITE(12,3) S(I,1), S(I,2), S(I,3)
c      DO 63 I=1,3
c      DO 63 J=1,3
c      aco(i,j) = acos(S(i,j))
c  63  aco(i,j) = (180./PI)*aco(i,j)
c      WRITE(12,64)
c  64  FORMAT(/' Rotation Angles'/)
c      DO 65 I=1,3
c  65  WRITE(12,3) aco(I,1), aco(I,2), aco(I,3)
c
c      Compute rotation angle and eigenvalues.
c
       qxx=F(1,1)
       qyy=F(2,2)
       qxy=F(1,2)
c  Rotation Angle
       t2t = 2.*qxy/(qxx-qyy)
       theta = 0.5*atan(t2t)
       theta = theta*(180./PI)
c
c  Compute for Cs symmetry rotation angle using dir cos
c
       thet = acos(S(1,1))
       thet = thet*(180./PI)
c  Eigenvalues
       ft = ((qxx+qyy)/2.)
       st = ((qxx-qyy)/2.)*sqrt(1.+t2t*t2t)
       q1 = ft+st
       q2 = ft-st
c  Write
c      write (12,33) qxx,qyy,qxy
c 33   format (///' Ixx,Iyy,Ixy = ', 3F12.6)
c      write (12,31) theta,thet
c 31   format(/'Cs or higher: Rotation Angle (x,u) =  ', 2F10.6)
c      write (12,32) q1,q2,F(3,3)
c 32   format (/'Cs or higher: Eigenvalues Iu,Iv,Iw = ', 3F12.6)
c
c  Calculate rotation constants
c
       conv=505379.0056
       AIa=conv/H(1,1)
       BIb=conv/H(2,2)
       CIc=conv/H(3,3)
       write (12,39)
  39   format(/' Rotational Constants (MHz)'/)
       write (12,71) AIa
  71   format('     A = ', F14.6)
       write (12,72) BIb
  72   format('     B = ', F14.6)
       write (12,73) CIc
  73   format('     C = ', F14.6)
  70   format(F15.6)
c
c  Rotate to a,b,c coordinates
c
       Call xyztoabc(N,X,Y,Z,S,Aj,Bj,Cj)
c
c      theta=(theta)*(PI/180.)
c      DO 27 I=1,N
c      A(I) = X(I)*cos(theta) + Y(I)*sin(theta)
c 27   B(I) = -X(I)*sin(theta) + Y(I)*cos(theta)
c      Write (12,28)
c 28   Format(//'Cs or higher: Atomic coordinates in u,v,w system'/)
c      Do 29 I=1,N
c 29   Write (12,3) A(I),B(I),Z(I)
c
       Call TQR(S,Q,eQq)
c
       CLOSE (11)
       CLOSE (12)
       END
c
c   -----------------------------------------------------------------
c
c   This subroutine is from Z.Kisiel's QDIAG.f program.
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
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER (NMAXDIM=3)
c
      COMMON /HBLK/H(3,3)/RBLK/S(3,3)
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 20 m=k,n
          sum=sum+h(i,m)**2
          s(i,m)=0.d0
          s(m,i)=0.d0
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 89 i=1,l
          if(h(i,m))60,89,60
60        q=dabs(h(i,m))
          if(q-unu)89,62,62
62        ind=1
          if(h(i,i)-h(m,m))68,64,68
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)
          if(1.d30-x1a)66,66,71
66        si=dsign(ss,x1)
          co=ss
          goto 77
71        if(1.d-30-x1a)72,89,89
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
89    continue
      if(ind)92,91,55
91    if(unu-unufi)92,92,54
c
92    return
      end
c
c  --------------------------------------------------------------
c
c    Coordinate transformation:  a(i) = T(i,j)*x(j)
c       where T is transpose of R (which is returned by JACK).
c
c    Performs transformation from Gaussian x,y,z standard
c    orientation coordinates to intertial a,b,c coordinates.
c    Called by MAIN.
c
        SUBROUTINE xyztoabc(N,XX,YY,ZZ,R,a,b,c)
c
        IMPLICIT REAL*8 (A-H,O-Z)
        IMPLICIT INTEGER (I-N)
c
        REAL*8 X(50,3),R(3,3),aa(50,3),XX(50),YY(50),ZZ(50)
        REAL*8 a(50),b(50),c(50),T(3,3)
        N=N+1
        DO 2 i=1,N
        X(i,1)=XX(i)
        X(i,2)=YY(i)
    2   X(i,3)=ZZ(i)
c       write(12,10)
c  10   format(/' Transformation Matrix, R (dir cos)'/)
c       do 3 j=1,3
c   3   write(12,7) R(j,1),R(j,2),R(j,3)
        do 12 i=1,3
        do 12 j=1,3
   12   T(i,j) = R(j,i)
c       write(12,14)
c  14   format(/' Transformation Matrix, T (R transpose)'/)
c       do 13 j=1,3
c  13   write(12,7) T(j,1),T(j,2),T(j,3)
        do 4 i=1,N
        do 4 j=1,3
    4   aa(i,j)=0.
        do 5 i=1,N
        do 5 j=1,3
        do 5 k=1,3
    5   aa(i,j)=aa(i,j)+T(j,k)*x(i,k)
        write(12,15)
   15   format(/' Atomic coordinates in a,b,c system'/)
        write(12,11)
   11   format('        a           b          c '/)
        do 8 i=1,N-1
    8   write(12,7) aa(i,1),aa(i,2),aa(i,3)
    7   format(4F12.6)
        write(12,20)
   20   format(/' Electric Dipole Moments a, b, c, total'/)
        totdip=sqrt(aa(N,1)**2 + aa(N,2)**2 + aa(N,3)**2)
       write(12,9) aa(N,1),aa(N,2),aa(N,3),totdip
   9   format(4F12.4)
        DO 6 i=1,N
        a(i)=aa(i,1)
        b(i)=aa(i,2)
    6   c(i)=aa(i,3)
        return
        END
c
c  ---------------------------------------------------------------
c
c     QI(i,j)=T(i,j)*Q(i,j)*R(i,j)
c
c     QI is efg tensor in a,b,c axes system.
c     Q is efg tensor in x,y,z; gaussian standard orientation.
c     R is transformation matrix from x,y,z to a,b,c.
c     T is transpose of R, which is returned by JACK.
c
c
c
c
       SUBROUTINE TQR(R,Q,eQq)
c
c
       IMPLICIT REAL*8 (A-H,O-Z)
       IMPLICIT INTEGER (I-N)
c
c
       COMMON /HBLK/H(3,3)/RBLK/S(3,3)
       REAL*8 T(3,3),Q(3,3),R(3,3),QI(3,3),A(3,3),QCC(3,3),
     & aco(3,3),qH(3,3)
c
       PI = (4.)*ATAN(1.)
c
       do 20 i=1,3
       do 20 j=1,3
   20  T(i,j)=R(j,i)
c      write(12,10)
c  10  format(/' T(3,3), transpose of R(3,3) '/)
c      do 1 i=1,3
c   1  write(12,*) T(i,1),T(i,2),T(i,3)
c      write(12,11)
c  11  format(/' Input: Q(3,3), Gaussian efg tensor '/)
c      do 2 i=1,3
c   2  write(12,*) Q(i,1),Q(i,2),Q(i,3)
       do 4 i=1,3
       do 4 j=1,3
    4  A(i,j)=0.
       do 5 i=1,3
       do 5 j=1,3
       do 5 k=1,3
    5  A(i,j)=A(i,j)+Q(i,k)*R(k,j)
       do 6 i=1,3
       do 6 j=1,3
    6  QI(i,j)=0.
       do 7 i=1,3
       do 7 j=1,3
       do 7 k=1,3
    7  QI(i,j)=QI(i,j)+T(i,k)*A(k,j)
       write(12,13)
   13  format(/' efg tensor in a,b,c system'/)
       WRITE(12,30)
   30  FORMAT('         a           b           c'/)
       WRITE(12,31) QI(1,1),QI(1,2),QI(1,3)
   31  FORMAT(' a', 3F12.6)
       WRITE(12,32) QI(2,1),QI(2,2),QI(2,3)
   32  FORMAT(' b', 3F12.6)
       WRITE(12,33) QI(3,1),QI(3,2),QI(3,3)
   33  FORMAT(' c', 3F12.6)
       eQq=-eQq
       do 9 i=1,3
       do 9 j=1,3
    9  QCC(i,j) = eQq*QI(i,j)
       write(12,40)
   40  format(/' nqcc tensor in a,b,c system'/)
       WRITE(12,41)
   41  FORMAT('         a           b           c'/)
       WRITE(12,82) QCC(1,1),QCC(1,2),QCC(1,3)
   82  FORMAT(' a', 3F12.6)
       WRITE(12,83) QCC(2,1),QCC(2,2),QCC(2,3)
   83  FORMAT(' b', 3F12.6)
       WRITE(12,84) QCC(3,1),QCC(3,2),QCC(3,3)
   84  FORMAT(' c', 3F12.6)
       WRITE(12,45)
   45  FORMAT(//'  --------------------------------------------'/)
    3  FORMAT(4F12.6)
       do 46 i=1,3
       do 46 j=1,3
   46  H(i,j)=QCC(i,j)
       Call JACK(3)
       Call ARRANGE(H,S,eQq)
       return
       end
c
c  --------------------------------------------------------------
c
c    This subroutine orders the eigenvalues of the intertia tensor
c    (and corresponding eigenvectors), H(1,1) < H(2,2) < H(3,3). 
c    It returns them in this order in HO and SO.  ORDER is
c    called by MAIN.
c
        SUBROUTINE ORDER(H,S,HO,SO)
c
c
        IMPLICIT REAL*8 (A-H,O-Z)
        IMPLICIT INTEGER (I-N)
c
        REAL*8 H(3,3),S(3,3),HO(3,3),SO(3,3)
c
c
       IF(H(1,1).LT.H(2,2).and.H(2,2).LT.H(3,3)) then
         HO(1,1)=H(1,1)
         HO(2,2)=H(2,2)
         HO(3,3)=H(3,3)
         DO 4 i=1,3
         DO 4 j=1,3
    4    SO(i,j)=S(i,j)
       else
       endif
       IF(H(2,2).LT.H(1,1).and.H(1,1).LT.H(3,3)) then
          HO(1,1)=H(2,2)
          HO(2,2)=H(1,1)
          HO(3,3)=H(3,3)
          DO 5 i=1,3
          SO(i,1)=S(i,2)
          SO(i,2)=S(i,1)
    5     SO(i,3)=S(i,3)
       else
       endif
       IF(H(1,1).LT.H(3,3).and.H(3,3).LT.H(2,2)) then
          HO(1,1)=H(1,1)
          HO(2,2)=H(3,3)
          HO(3,3)=H(2,2)
          DO 6 i=1,3
          SO(i,1)=S(i,1)
          SO(i,2)=S(i,3)
    6     SO(i,3)=S(i,2)
       else
       endif
       IF(H(3,3).LT.H(1,1).and.H(1,1).LT.H(2,2)) then
          HO(1,1)=H(3,3)
          HO(2,2)=H(1,1)
          HO(3,3)=H(2,2)
          DO 7 i=1,3
          SO(i,1)=S(i,3)
          SO(i,2)=S(i,1)
    7     SO(i,3)=S(i,2)
       else
       endif
       IF(H(2,2).LT.H(3,3).and.H(3,3).LT.H(1,1)) then
          HO(1,1)=H(2,2)
          HO(2,2)=H(3,3)
          HO(3,3)=H(1,1)
          DO 8 i=1,3
          SO(i,1)=S(i,2)
          SO(i,2)=S(i,3)
    8     SO(i,3)=S(i,1)
       else
       endif
       IF(H(3,3).LT.H(2,2).and.H(2,2).LT.H(1,1)) then
          HO(1,1)=H(3,3)
          HO(2,2)=H(2,2)
          HO(3,3)=H(1,1)
          DO 9 i=1,3
          SO(i,1)=S(i,3)
          SO(i,2)=S(i,2)
    9     SO(i,3)=S(i,1)
       else
       endif
       return
       end
c
c------------------------------------------------------------------------
c
c
c   SUBROUTINE ARRANGE(H,S,eQq)
c
c       This subroutine arrange the Xii
c       (and corresponding qii and eigenvectors) such that
c       |Xxx|<|Xyy|<|Xzz|, and writes the output.  ARRANGE is
c       called by TQR.
c
        SUBROUTINE ARRANGE(H,S,eQq)
c
        IMPLICIT REAL*8 (A-H,O-Z)
        IMPLICIT INTEGER (I-N)
c
        REAL*8 H(3,3),S(3,3),HO(3,3),SO(3,3),G(3,3),GO(3,3)
        REAL*8 q(3,3),aco(3,3)
c
       PI = (4.)*ATAN(1.)
c
       do 1 i=1,3
       do 1 j=1,3
    1  G(i,j)=0.
       do 2 i=1,3
    2  G(i,i)=H(i,i)
       If(G(1,1) .gt. 0.0 .and. G(2,2) .gt. 0.0) then
         continue
         If(G(1,1) .gt. G(2,2)) then
         HO(1,1)=G(2,2)
         HO(2,2)=G(1,1)
         HO(3,3)=G(3,3)
         DO 3 i=1,3
         SO(i,1)=S(i,2)
         SO(i,2)=S(i,1)
    3    SO(i,3)=S(i,3)
         else
         HO(1,1)=G(1,1)
         HO(2,2)=G(2,2)
         HO(3,3)=G(3,3)
         DO 4 i=1,3
         SO(i,1)=S(i,1)
         SO(i,2)=S(i,2)
    4    SO(i,3)=S(i,3)
         endif
       else
       endif
       If(G(1,1) .gt. 0.0 .and. G(3,3) .gt. 0.0) then
         continue
         If(G(1,1) .gt. G(3,3)) then
         HO(1,1)=G(3,3)
         HO(2,2)=G(1,1)
         HO(3,3)=G(2,2)
         DO 5 i=1,3
         SO(i,1)=S(i,3)
         SO(i,2)=S(i,1)
    5    SO(i,3)=S(i,2)
         else
         HO(1,1)=G(1,1)
         HO(2,2)=G(3,3)
         HO(3,3)=G(2,2)
         DO 6 i=1,3
         SO(i,1)=S(i,1)
         SO(i,2)=S(i,3)
    6    SO(i,3)=S(i,2)
         endif
       else
       endif
       If(G(2,2) .gt. 0.0 .and. G(3,3) .gt. 0.0) then
         continue
         If(G(2,2) .gt. G(3,3)) then
         HO(1,1)=G(3,3)
         HO(2,2)=G(2,2)
         HO(3,3)=G(1,1)
         DO 7 i=1,3
         SO(i,1)=S(i,3)
         SO(i,2)=S(i,2)
    7    SO(i,3)=S(i,1)
         else
         HO(1,1)=G(2,2)
         HO(2,2)=G(3,3)
         HO(3,3)=G(1,1)
         DO 8 i=1,3
         SO(i,1)=S(i,2)
         SO(i,2)=S(i,3)
    8    SO(i,3)=S(i,1)
         endif
       else
       endif
       If(G(1,1) .lt. 0.0 .and. G(2,2) .lt. 0.0) then
         continue
         If(G(1,1) .lt. G(2,2)) then
         HO(1,1)=G(2,2)
         HO(2,2)=G(1,1)
         HO(3,3)=G(3,3)
         DO 9 i=1,3
         SO(i,1)=S(i,2)
         SO(i,2)=S(i,1)
    9    SO(i,3)=S(i,3)
         else
         HO(1,1)=G(1,1)
         HO(2,2)=G(2,2)
         HO(3,3)=G(3,3)
         DO 10 i=1,3
         SO(i,1)=S(i,1)
         SO(i,2)=S(i,2)
   10    SO(i,3)=S(i,3)
         endif
       else
       endif
       If(G(1,1) .lt. 0.0 .and. G(3,3) .lt. 0.0) then
         continue
         If(G(1,1) .lt. G(3,3)) then
         HO(1,1)=G(3,3)
         HO(2,2)=G(1,1)
         HO(3,3)=G(2,2)
         DO 11 i=1,3
         SO(i,1)=S(i,3)
         SO(i,2)=S(i,1)
   11    SO(i,3)=S(i,2)
         else
         HO(1,1)=G(1,1)
         HO(2,2)=G(3,3)
         HO(3,3)=G(2,2)
         DO 12 i=1,3
         SO(i,1)=S(i,1)
         SO(i,2)=S(i,3)
   12    SO(i,3)=S(i,2)
         endif
       else
       endif
       If(G(2,2) .lt. 0.0 .and. G(3,3) .lt. 0.0) then
         continue
         If(G(2,2) .lt. G(3,3)) then
         HO(1,1)=G(3,3)
         HO(2,2)=G(2,2)
         HO(3,3)=G(1,1)
         DO 13 i=1,3
         SO(i,1)=S(i,3)
         SO(i,2)=S(i,2)
   13    SO(i,3)=S(i,1)
         else
         HO(1,1)=G(2,2)
         HO(2,2)=G(3,3)
         HO(3,3)=G(1,1)
         DO 14 i=1,3
         SO(i,1)=S(i,2)
         SO(i,2)=S(i,3)
   14    SO(i,3)=S(i,1)
         endif
       else
       endif
       do 15 i=1,3
   15  q(i,i)=HO(i,i)/eQq
       WRITE(12,21)
   21  FORMAT(/' Principal efg: qii where i=x,y,z'/)
       WRITE(12,22) q(1,1),q(2,2),q(3,3)
   22  FORMAT(' xx =',F12.6,'  yy =',F12.6,'  zz =',F12.6)
       WRITE(12,57)
   57  FORMAT(/' Principal nqcc: Xii where i=x,y,z'/)
       WRITE(12,85) HO(1,1),HO(2,2),HO(3,3)
   85  FORMAT(' xx =',F12.6,'  yy =',F12.6,'  zz =',F12.6)
       WRITE(12,59)
   59  FORMAT(/' Eigenvectors, direction cosines'/)
       WRITE(12,47)
   47  FORMAT('         x           y           z'/)
       WRITE(12,42) SO(1,1),SO(1,2),SO(1,3)
   42  FORMAT(' a', 3F12.6)
       WRITE(12,43) SO(2,1),SO(2,2),SO(2,3)
   43  FORMAT(' b', 3F12.6)
       WRITE(12,44) SO(3,1),SO(3,2),SO(3,3)
   44  FORMAT(' c', 3F12.6)
       DO 63 I=1,3
       DO 63 J=1,3
       aco(i,j) = acos(SO(i,j))
   63  aco(i,j) = (180./PI)*aco(i,j)
       WRITE(12,64)
   64  FORMAT(/' Rotation Angles (degrees)'/)
       WRITE(12,17) aco(1,1),aco(1,2),aco(1,3)
   17  FORMAT(' a', 3F12.4)
       WRITE(12,18) aco(2,1),aco(2,2),aco(2,3)
   18  FORMAT(' b', 3F12.4)
       WRITE(12,19) aco(3,1),aco(3,2),aco(3,3)
   19  FORMAT(' c', 3F12.4)
       ETA=(HO(1,1)-HO(2,2))/HO(3,3)
       WRITE(12,20) ETA
   20  FORMAT(/' (Xxx-Xyy)/Xzz = ', F12.6/)
       WRITE(12,45)
   45  FORMAT(/'  -------------------  fini  ------------------'/)
       return
       end
c
c  __________________________________________________________________
c









QTRANS.html



Last modified: 6 Feb 2019