program eigen
c        1         2         3         4         5         6        7          8
c2345678901234567890123456789012345678901234567890123456789012345678901234567890
c***********************************************************************
c* Calculates eigen-values and eigen-vectors for any symmetric second  *
c* order stress tensor, Sij, using the Jacobi method outlined in       *
c* Appplied Numerical Methods, B. Carnahan, H.A. Luther, J.O. Wikes,   *
c* John Wiley and Sons Inc., pp. 250-260, 1969.                        *
c***********************************************************************
c
      implicit real*8 (a-h,o-z)
      dimension a(20,20), eval(20), evect(20,20)
      dimension x(4), y(4), z(4), d(18)
      character  description(80)*1
c
c  open files for input and output
      open(1,file='eigen-vectors_1.out',status='unknown',err=8888)
      open(2,file='eigen-vectors_2.out',status='unknown',err=8888)
      open(3,file='eigen-vectors_3.out',status='unknown',err=8888)
      open(4,file='eigen-vectors_4.out',status='unknown',err=8888)
      open(5,file='sij.data',status='old',err=8888)
      open(6,file='eigen-calc-details.out',status='unknown',err=8888)
      open(7,file='eigen-stress.out',status='unknown',err=8888)
      open(8,file='eigen-vectors.out',status='unknown',err=8888)
      open(9,file='xyz_1.txt',status='unknown',err=8888)
      open(10,file='xyz_2.txt',status='unknown',err=8888)
      open(11,file='xyz_3.txt',status='unknown',err=8888)
      open(12,file='xyz_4.txt',status='unknown',err=8888)
      open(13,file='euler_angles.out',status='unknown',err=8888)
      open(14,file='translate_pnts.out',status='unknown',err=8888)
c
c  set tolerances for convergence in subroutine jacobi 
      eps1=1.00e-10
      eps2=1.00e-10
      eps3=1.00e-05
      itmax=50
      pi=acos(-1.)
c
c  read alpha-numeric description
      read(5,111)(description(i),i=1,80)
c
c  write descprition at the beginning of all output files
      write(1,222)(description(i),i=1,80)
      write(2,222)(description(i),i=1,80)
      write(3,222)(description(i),i=1,80)
      write(4,222)(description(i),i=1,80)
      write(6,222)(description(i),i=1,80)
      write(7,222)(description(i),i=1,80)
      write(8,222)(description(i),i=1,80)
      write(9,222)(description(i),i=1,80)
      write(10,222)(description(i),i=1,80)
      write(11,222)(description(i),i=1,80)
      write(12,222)(description(i),i=1,80)
      write(13,222)(description(i),i=1,80)
c
c  write titles for eigen-values and -vectors
      write(7,333)
      write(8,444)
c
c  loop thru nmat matrices
      nmat=4
      do 4 im=1,nmat
c
c  read and establish starting matrix a
      do 2 i=1,20
      do 2 j=1,20
      evect(i,j)=0.0
      eval(i)=0.0
    2 a(i,j)=0.0
c
c  all matrices are second order stress tensors, n=3
      n=3
      read(5,101)a(1,1),a(2,2),a(3,3),a(1,2),a(1,3),a(2,3),
     *x(im),y(im),z(im)
      a(2,1)=a(1,2)
      a(3,1)=a(1,3)
      a(3,2)=a(2,3)
      write(6,200) n,itmax,eps1,eps2,eps3
      do 3 i=1,n
    3 write(6,201) (a(i,j),j=1,n)
c
c  compute eigenvalues and eigenvectors
      call jacobi(n,a,eval,evect,itmax,eps1,eps2,eps3)
c
c  calculate angles from unprimed to primed axes using eigenvectors
      angx=180.0*acos(evect(1,1))/pi
      angy=180.0*acos(evect(2,2))/pi
      angz=180.0*acos(evect(3,3))/pi
c
c  calculate Euler angles (degrees) from eignevector (transformation matrix)
c
      aa=evect(3,1)/evect(3,3)
      eulerangy=180.0*atan(aa)/pi
      aa=-evect(3,2)
      eulerangx=180.0*asin(aa)/pi
      aa=evect(1,2)/evect(2,2)
      eulerangz=180.0*atan(aa)/pi
c
c  calculate Euler angles (radians)
c 
      eulernangyrad=eulerangy*pi/180.0
      eulernangxrad=eulerangx*pi/180.0
      eulernangzrad=eulerangz*pi/180.0
      write(13,102)eulernangxrad,eulernangyrad,eulernangzrad
c
c  write results  
      write(6,204)
      write(6,201)(eval(i),i=1,n)
      write(7,101)(eval(i),i=1,n)
c     write(6,205)
c     do 23 i=1,n
c  23 write(6,201)(a(i,j),j=1,n)
      write(6,206)
      do 24 i=1,n
   24 write(6,201)(evect(i,j),j=1,n)
      write(6,208)angx,angy,angz
      write(8,202)im,x(im),y(im),z(im)
      if(im.eq.1)write(1,202)im,x(im),y(im),z(im)
      if(im.eq.1)write(9,202)im,x(im),y(im),z(im)
      if(im.eq.2)write(2,202)im,x(im),y(im),z(im)
      if(im.eq.2)write(10,202)im,x(im),y(im),z(im)
      if(im.eq.3)write(3,202)im,x(im),y(im),z(im)
      if(im.eq.3)write(11,202)im,x(im),y(im),z(im)
      if(im.eq.4)write(4,202)im,x(im),y(im),z(im)
      if(im.eq.4)write(12,202)im,x(im),y(im),z(im)
      write(8,209)angx,angy,angz
      if(im.eq.1)write(1,209)angx,angy,angz
      if(im.eq.2)write(2,209)angx,angy,angz
      if(im.eq.3)write(3,209)angx,angy,angz
      if(im.eq.4)write(4,209)angx,angy,angz
      write(8,210)eulerangx,eulerangy,eulerangz
      if(im.eq.1)write(1,210)eulerangx,eulerangy,eulerangz
      if(im.eq.2)write(2,210)eulerangx,eulerangy,eulerangz
      if(im.eq.3)write(3,210)eulerangx,eulerangy,eulerangz
      if(im.eq.4)write(4,210)eulerangx,eulerangy,eulerangz
c
    4 continue
c
c  read and write bounding box [111] corner X-Y-Z position
c     read(5,300)xb,yb,zb
c     write(6,310)xb,yb,zb 
c
c  calculate scale factor
c
c   differences in x:
      d(1)=abs(x(2)-x(1))
      d(2)=abs(x(3)-x(1))
      d(3)=abs(x(4)-x(1))
      d(4)=abs(x(3)-x(2))
      d(5)=abs(x(4)-x(2))
      d(6)=abs(x(4)-x(3))
c   differences in y:
      d(7)=abs(y(2)-y(1))
      d(8)=abs(y(3)-y(1))
      d(9)=abs(y(4)-y(1))
      d(10)=abs(y(3)-y(2))
      d(11)=abs(y(4)-y(2))
      d(12)=abs(y(4)-y(3))
c   differences in z:
      d(13)=abs(z(2)-z(1))
      d(14)=abs(z(3)-z(1))
      d(15)=abs(z(4)-z(1))
      d(16)=abs(z(3)-z(2))
      d(17)=abs(z(4)-z(2))
      d(18)=abs(z(4)-z(3))
c
c   find maximum
      absmax=d(1)
      do 5 i=2,18
      if (d(i).gt.absmax) absmax=d(i)
    5 continue
c
c   find nonzero minimum
      absmin=absmax
      do 6 i=2,18
      if((d(i).gt.0.0).and.(d(i).lt.absmin)) absmin=d(i)
    6 continue
      scale=10.0/absmin
      write(6,555)absmax,absmin,scale
      do 7 i=1,4
      x(i)=scale*x(i)
      y(i)=scale*y(i)
      z(i)=scale*z(i)
      write(6,666)i,x(i),i,y(i),i,z(i)
      write(14,777)x(i),y(i),z(i)
    7 continue
c
c  formats for input and output statements
c        1         2         3         4         5         6        7
c23456789012345678901234567890123456789012345678901234567890123456789012
  111 format(80a1,/)
  222 format(80a1)
  333 format(35h Eigen-values (Principal stresses) )
  444 format(/,2x,46hEigen-vectors (Principal Stress Orientations) )
  555 format(/,2x,10hMaximum = ,e12.5,2x,10hMinimum = ,e12.5,2x,
     *15hScale Factor = ,e12.5)
  666 format(/,2x,2hx(,i1,2h)=,e12.5,2x,2hy(,i1,2h)=,e12.5,2x,2hz(,i1,2h
     *)=,e12.5)
  777 format(3(2x,e10.3))
  100 format ( )
  101 format (6(1x,e10.3),3(1x,e9.2))
  102 format (3(1x,e12.5))
  200 format ( 5x, 54h__________________________________________________
     *____,//, 
     * 4x, 55h determination of eigenvalues by jacobi's method, where,
     * //, 7x, 10h n      = ,i4,/
     * 7x, 10h itmax  = , i4/ 7x, 10h eps1   = , e12.3,/
     * 7x, 10h eps2   = , e12.3/ 7x, 10h eps3   = , e12.3,//,
     * 4x, 39h the starting matrix a(1,1)...a(n,n) is,/ )
  201 format (7x, 10e11.3)
  202 format (/,8x,35h-----------------------------------,
     */,8x,22h           LOCATION - ,i1,
     */,8x,35h-----------------------------------,
     */,8x,2hX=,e9.2,1x,2hY=,e9.2,1x,2hZ=,e9.2) 
  204 format (/,4x,36h eigenvalues eigen(1)...eigen(n) are,/ )
  205 format (/,4x,34h the final transformed matrix a is,/ )
  206 format (/,4x,66h eignevectors are in corresponding columns of the
     *following matrix,/ )
  208 format (/,4x,30h corresonding angles (degrees),/,7x,3f11.6,/)
  209 format (/,9x,30h angles (degrees) between axes,/, 
     *8x,33h (  X->x'      Y->y'      Z->z' ),/,
     *7x,3f11.6)
  210 format (/,13x,23h euler angles (degrees),/,
     *8x,33h (   x          y          z    ),/,
     *7x,3f11.6)
c 300 format(43x,3(1x,e9.2))
c 310 format(/,54h______________________________________________________
c    *,//,43hBounding Box Diagonal [111] Corner (X-Y-Z):3(1x,e9.2))
c
 8888 stop
      end
c
c#######################################################################
      subroutine jacobi(n,a,eigen,t,itmax,eps1,eps2,eps3)
c        1         2         3         4         5         6        7
c23456789012345678901234567890123456789012345678901234567890123456789012
c***********************************************************************
c* Appplied Numerical Methods, B. Carnahan, H.A. Luther, J.O. Wikes,   *
c* John Wiley and Sons Inc., pp. 250-260, 1969.                        *
c* Example 4.4: Jacobi's method for solving eigenvalues and            *
c*              eigenvectors from real symmetric matrices              *
c*                  (Modified by R.D. Kriz 9-27-03)                    *
c* Only the upper triangular part of the real symmetric starting       *
c* matrix a is read.   T is an n by n orthogonal matrix, being the     *
c* product of the sequence of transformation matrices used to ann-     *
c* ihilate successively the off-diagonal elements of A, and con-       *
c* sequently to reduce A iteratively to near-diagonal form.            *
c* Sigma1 and sigma2 are the values of sigma A(i,i)**2 before and      *
c* after one complete iteration.  Since annihilation of an element     *
c* may cause another already-annihilated element to assume a non-zero  *
c* value, the iterative process is repeated until the convergence test * 
c* is passed (1-sigma1/sigma2) less than eps3 or itmax complete        *
c* iterations have failed to produce convergence.   An off-diagonal    *
c* element smaller in magnitude than eps2 is ignored in the            *
c* annihilation process.   Eps1 is compared with (A(i,i)-A(j,j)) to    *
c* determine the method ofcomputing dsa and sna, the cosine and sine   *
c* of the rotation angle, respectively.  When A has been reduced to    *
c* diagonal form, the eigenvalues will be found in the diagonal        *
c* positions and are saved in eigen(1)... eigen(n).   The associated   *
c* eigenvectors are in corresponding columns of the final              *
c* transformation matrix T.  S is the sum of squares of all elements   *
c* in the original matrix A, and should theoretically equal the final  *
c* value of sigma2.                                                    *
c***********************************************************************
c
      implicit real*8 (a-h,o-z)
      dimension a(20,20),t(20,20),aik(20),eigen(20)
c
      nm1=n-1
c
c  set up initial matrix t, computer sigma1 and s
      sigma1=0.0
      offdsq=0.0
      do 5 i=1,n
      sigma1=sigma1+a(i,i)**2
      t(i,i)=1.0
      ip1=i+1
      if(i.ge.n)go to 6
      do 5 j=ip1,n
    5 offdsq = offdsq + a(i,j)**2
    6 s=2.0*offdsq+sigma1
c
c  begin jacobi iteration
      do 26 iter=1,itmax
      do 20 i=1,nm1
      ip1=i+1
      do 20 j=ip1,n
      q=dabs(a(i,i)-a(j,j))
c
c  computer sine and cosine of rotation angle
      if(q.le.eps1)go to 9
        if(dabs(a(i,j)).le.eps2)go to 20
        p=2.0*a(i,j)*q/(a(i,i)-a(j,j))
        spq=dsqrt(p*p+q*q)
        csa=dsqrt((1.0+q/spq)/2.0)
        sna=p/(2.0*csa*spq)
        go to 10
    9   csa=1.0/dsqrt(2.0D0)
        sna=csa
   10 continue
c
c  update colums i and j of t - equivalent to
c  multiplication of the annihilation matrix
      do 11 k=1,n
      holdki=t(k,i)
      t(k,i)=holdki*csa+t(k,j)*sna
   11 t(k,j)=holdki*sna-t(k,j)*csa
c
c  compute new elements of a in rows i and j
      do 16 k=1,n
      if(k.gt.j)go to 15
        aik(k)=a(i,k)
        a(i,k)=csa*aik(k)+sna*a(k,j)
        if(k.ne.j)go to 14
          a(j,k)=sna*aik(k)-csa*a(j,k)
   14   go to 16
   15   holdik=a(i,k)
        a(i,k)=csa*holdik+sna*a(j,k)
        a(j,k)=sna*holdik-csa*a(j,k)
   16 continue
c
c  compute new elements of matrix a in columns i and j
      aik(j)=sna*aik(i)-csa*aik(j)
c
c  when k is larger than i
      do 19 k=1,j
      if(k.le.i)go to 18
        a(k,j)=sna*aik(k)-csa*a(k,j)
        go to 19
   18 holdki=a(k,i)
      a(k,i)=csa*holdki+sna*a(k,j)
      a(k,j)=sna*holdki-csa*a(k,j)
   19 continue
   20 a(i,j)=0.0
c
c  find sigma2 for transformed a and test for convergence
      sigma2=0.0
      do 21 i=1,n
      eigen(i)=a(i,i)
   21 sigma2=sigma2+eigen(i)**2
      if(1.0-sigma1/sigma2.ge.eps3)go to 25
c
c  converged, write results 
        write(6,204) iter, s, sigma2
        return 
c
c  still converging, write results
   25 write(6,202) iter,sigma1,sigma2
      sigma1=sigma2
c
c  loop on variable iter ( do 26 iter=1,itmax )
   26 continue
c
c  iter exceeded itmax, no convergence write results and return
      write(6,203) iter,s,sigma1,sigma2
      do 27 i=1,n
   27 write(6,201) (a(i,j),j=1,n)
      write(6,207)
      do 28 i=1,n
   28 write(6,201) (t(i,j),j=1,n)
c
c  formats for output statements
  201 format (7x, 10fe11.3,/)
  202 format (6x, 10h iter   = , i5, 3x, 10h sigma1 = , e10.3,
     * 3x, 10h sigma2 = , e10.3)
  203 format (/, 4x, 21h no convergence, with// 7x, 10h iter   = , i5,
     * 5x, 10h s      = , e10.3, 5x, 10h sigma1 = , e10.3,
     * 5x, 10h sigma2 = , e10.3,//,5x, 24h the current a matrix is,/ )
  204 format (6x, 31h convergence has occured, where,/, 
     * 6x,10h iter   = , i5,8x,5h s = ,e10.3,3x,10h sigma2 = , e10.3,/
     * 5x, 20h____________________,/)
  207 format (/, 4x, 24h the current t matrix is,/ )
c
      return 
      end