      PROGRAM SAWave
c     USES GetSAW
c     :::
      call GetSAW
      stop
      END



      SUBROUTINE GetSAW
c     GIVES stdout
c     USES GetDat,ScanA
c     :::
      INTEGER angle(3,3),skip(3),till(3),ncall
      REAL rho,c(3,3,3,3),e(3,3,3),eps(3,3)
      COMMON ncall
      call GetDat(rho,c,e,eps,angle,skip,till)
      write(*,*) 'Calculations in progress.'
      ncall=0
      call ScanA(rho,c,e,eps,angle,skip,till)
      write(*,'(1x,''ncall='',i10)') ncall
      write(*,*) 'Calculations completed.'
      return
      END



      SUBROUTINE GetDat(rho,c,e,eps,angle,skip,till)
      INTEGER angle(3,3),skip(3),till(3)
      REAL rho,c(3,3,3,3),e(3,3,3),eps(3,3)
c     GETS stdin,files
c     GIVES rho,c,e,eps,angle,skip,till
c     USES SymMat,RetCE,CoeMat
c     Reads data from stdin and calculates the arguments.
c        The data (55 numbers): mass density rhoi (1 number), elastic
c     tensor ci(I,J) (21 numbers), piezoelectric tensor ei(i,J) (18
c     numbers), dielectric tensor epsi(i,j) (6 numbers), matrix of Euler
c     angles anglei(i,j) (9 numbers). Rows of the matrix anglei are
c     related to three Euler angles: lam,gam,the. Each row consists of
c     three values: first, last, and increment (these are three
c     parameters of the iteration loop in the subroutine ScanA).
c        In the above symbols, i,j=1,2,3 and I,J=1,...,6 , where I,J are
c     double indices (see the main comment of the subroutine RetCE).
c     Euler angles are given according to the convention (see the main
c     comment of the subroutine DefA).
c        The matrices ci,ei,epsi are represented in the form of a
c     symmetric matrix 9*9 (elasto-piezo-dielectric matrix according to
c     IEEE Standard on Piezoelectricity, see IEEE Trans. Ultrason.
c     Ferroel. Freq. Control, vol.43, 717 (1996)). Only elements on the
c     diagonal and below (45 numbers) are read (see the file saw-i.dgm).
c     All the matrices are read by rows.
c        The data are expressed in the following units: rhoi
c     [1e3*kg*m**(-3)]; ci [1e9*N*m**(-2)]; ei [C*m**(-2)]; epsi
c     [dimensionless]; anglei [degree].
c        The arguments are expressed in the following units: rho,c,e
c     [the same units as for rhoi,ci,ei]; eps [1e-9*F*m**(-1)]; angle
c     [the same unit as for anglei].
c        Field variables are expressed in adopted units as follows.
c     Independent variables: coordinates x,z [1e-3*m]; slownesses r,s
c     [1e-3*s*m**(-1)]; frequency omega [1e6*s**(-1)]. Dependent
c     variables: particle displacement vector u [1e-9*m]; electric
c     potential phi [V]; traction force vector T [1e3*N*m**(-2)];
c     electric displacement vector D [1e-6*C*m**(-2)].
c     :::
      INTEGER i,j,ios1,ios2
      REAL ci(6,6),cip(6,6),ei(3,6),epsi(3,3),epsip(3,3),EPS0
      PARAMETER (EPS0=8.85434e-3)
      read(*,*) rho,
     :   ((ci(i,j),j=1,i),i=1,6),
     :   ((ei(i,j),j=1,6),(epsi(i,j),j=1,i),i=1,3),
     :   ((angle(i,j),j=1,3),i=1,3)
      open(1,file='skip.ang',status='old',iostat=ios1)
      if (ios1.eq.0) then
       read(1,*) skip
      else
       skip(1)=angle(1,1)
       skip(2)=angle(2,1)
       skip(3)=angle(3,1)-angle(3,3)
      endif
      close(1)
      open(2,file='till.ang',status='old',iostat=ios2)
      if (ios2.eq.0) then
       read(2,*) till
      else
       till(1)=angle(1,2)
       till(2)=angle(2,2)
       till(3)=angle(3,2)
      endif
      close(2)
      call SymMat(ci,cip,6)
      call SymMat(epsi,epsip,3)
      call RetCE(cip,ei,c,e)
      call CoeMat(EPS0,epsip,eps,3)
      return
      END



      SUBROUTINE ScanA(rho,c,e,eps,angle,skip,till)
      INTEGER angle(3,3),skip(3),till(3)
      REAL rho,c(3,3,3,3),e(3,3,3),eps(3,3)
c     GETS rho,c,e,eps,angle,skip,till
c     GIVES stdout
c     USES Rotate,RootZ,WriteW,Offset
c     :::
      INTEGER a1,a2,a3,trial,fail
      REAL lam,gam,the,cp(3,3,3,3),ep(3,3,3),epsp(3,3),
     :   rq,rs,re,rz,epse,dinv,RAD,DEG
      LOGICAL done
      PARAMETER (RAD=3.14159265/180,DEG=180/3.14159265)
      do 13 a1=angle(1,1),angle(1,2),angle(1,3)
      do 12 a2=angle(2,1),angle(2,2),angle(2,3)
      do 11 a3=angle(3,1),angle(3,2),angle(3,3)
       if (a1.lt.skip(1)) goto 11
       if (a1.eq.skip(1)) then
        if (a2.lt.skip(2)) goto 11
        if (a2.eq.skip(2)) then
         if (a3.le.skip(3)) goto 11
        endif
       endif
       if (a1.gt.till(1)) return
       if (a1.eq.till(1)) then
        if (a2.gt.till(2)) return
        if (a2.eq.till(2)) then
         if (a3.gt.till(3)) return
        endif
       endif
       trial=0
       lam=a1*RAD
       gam=a2*RAD
       the=a3*RAD
  100  continue
       call Rotate(c,e,eps,lam,gam,the,cp,ep,epsp)
       call RootZ(rho,cp,ep,epsp,rq,rs,re,rz,epse,dinv,fail)
       if (fail.ne.0) goto 200
       call
     :    WriteW(rho,cp,ep,epsp,lam,gam,the,rq,rs,re,rz,epse,dinv,fail)
       if (fail.ne.0) fail=10000
  200  continue
       done=fail.le.0.or.trial.ge.7
       if (.not.done) then
        call Offset(a1,a2,a3,trial,lam,gam,the)
        goto 100
       endif
       if (fail.ne.0) then
        write(*,'(1x,3f8.2,2x,''**'',1x,i7,5x,''#0'')')
     :     lam*DEG,gam*DEG,the*DEG,fail
       endif
   11 continue
   12 continue
   13 continue
      return
      END



      SUBROUTINE WriteW
     :   (rho,cp,ep,epsp,lam,gam,the,rq,rs,re,rz,epse,dinv,fail)
      INTEGER fail
      REAL rho,cp(3,3,3,3),ep(3,3,3),epsp(3,3),lam,gam,the,
     :   rq,rs,re,rz,epse,dinv
c     GETS rho,cp,ep,epsp,lam,gam,the,rq,rs,re,rz,epse,dinv
c     GIVES stdout,fail
c     USES CalcH,FindVV,CalcPF
c     Performs additional calculations and writes results to stdin in
c     one line of the length of 157 characters (see the file saw-o.dgm).
c        The results are wave parameters (12 numbers): cutoff slowness
c     of bulk waves rq, slowness of secondary SAW rs, slowness of
c     open-circuit SAW re, slowness of short-circuit SAW rz; parameters
c     of short-circuit SAW: phase velocity vz=1/rz, beam steering angle
c     psi, piezoelectric coupling coefficient dinv, effective dielectric
c     permittivity epse, normal component of electric displacement D(3),
c     components of particle displacement u(1),u(2),u(3).
c        If rs=0 is output then the secondary SAW does not exist. If
c     re=0 is output then re has not been calculated. This occurs very
c     rarely, and only when dinv is small (rz-re<<0.01). If the error
c     message fail=-1 is output then SAW does not exist. Other error
c     messages mean that calculations have not been completed but SAW
c     may exist.
c        The results are expressed in the following units: lam,gam,the
c     [degree]; rq,rs,re,rz [1e-3*s*m**(-1)]; vz [1e3*m*s**(-1)]; psi
c     [degree]; dinv [dimensionless]; epse [dimensionless]; normalized
c     complex amplitudes D(3) [1e-6*C*m**(-2)/sp], and u(1),u(2),u(3)
c     [1e-9*m/sp], where sp=sqrt(1e-3*W*m**(-1)) (1e-3*W*m**(-1) is the
c     dimension of the Poynting vector integrated with respect to one
c     space variable).
c     :::
      INTEGER i
      REAL h(8,8),psi,vz,DEG,ua(4),up(4),ta(4),tp(4)
      COMPLEX q(8),f(8,8),utn(4),ttn(4)
      PARAMETER (DEG=180/3.14159265)
      call CalcH(rz,rho,cp,ep,epsp,h)
      call FindVV(h,q,f,fail)
      if (fail.ne.0) return
      call CalcPF(rz,cp,ep,epsp,q,f,psi,utn,ttn,fail)
      if (fail.ne.0) return
      do 10 i=1,4
       ua(i)=cabs(utn(i))
       if (real(utn(i)).eq.0.0.and.imag(utn(i)).eq.0.0) then
        up(i)=0.0
       else 
        up(i)=atan2(imag(utn(i)),real(utn(i)))
       endif
       ta(i)=cabs(ttn(i))
       if (real(ttn(i)).eq.0.0.and.imag(ttn(i)).eq.0.0) then
        tp(i)=0.0
       else 
        tp(i)=atan2(imag(ttn(i)),real(ttn(i)))
       endif
   10 continue
      vz=1/rz
      write(*,'(1x,3f8.2,2x,''*'',
     :   2x,4(f6.5,1x),f7.4,1x,f6.2,1x,f6.4,1x,f6.2,2x,''#1'',
     :   2x,f6.1,1x,f8.6,2x,3(1x,f6.1,1x,f7.5),2x,''#2'')')
     :   lam*DEG,gam*DEG,the*DEG,
     :   rq,rs,re,rz,vz,psi*DEG,dinv,epse,tp(4)*DEG,ta(4),
     :   (up(i)*DEG,ua(i),i=1,3)
      return
      END



      SUBROUTINE Balanc(nm,n,a,low,upp,scale)
      INTEGER nm,n,low,upp
      REAL a(nm,n),scale(n)
c     GETS nm,n,a
c     GIVES a,low,upp,scale
c     Balances a real general matrix in the 1-norm by exact similarity
c     transformations, and isolates eigenvalues whenever possible by
c     permutation similarity transformations.
c        The 1-norm of a matrix (or vector) is defined as the sum of the
c     absolute values of its elements. A matrix of rank n is said to be
c     balanced in the 1-norm if the 1-norms of the i-th row and i-th
c     column are equal for i=1,...,n.
c        Balancing reduces the 1-norm of the original matrix whenever
c     the 1-norms of some row and corresponding column are markedly
c     different, while at the same time leaving the eigenvalues
c     unchanged. Reducing the norm in this way can improve the accuracy
c     of the computed eigenvalues and/or eigenvectors.
c        Input: the array a (n must not be greater than nm) contains the
c     matrix of rank n to be balanced.
c        Output: the array a contains the transformed matrix (the
c     original matrix is overwritten).
c        Output: low,upp indicate the boundary indices for the balanced
c     matrix (low is not greater than upp). The (i,j)-element of the
c     balanced matrix is equal to zero if (1) i>j and (2) j<low or
c     i>upp. The eigenvalues are isolated in positions 1,...,low-1 and
c     upp+1,...,n.
c        Output: scale contains information about the similarity
c     transformations (determining the permutations used and the scaling
c     factors).
c        RADIX is the base (radix) of the number system used for
c     internal representation of numbers on the computer.
c     :::
      INTEGER i,j,m,p,q,iexch
      REAL c,f,g,r,s,swap,RADIX,SQRDX
      LOGICAL notcnv
      PARAMETER (RADIX=2.0,SQRDX=RADIX**2)
      p=1
      q=n
c     Search for rows isolating an eigenvalue and push them down.
  100 do 12 j=q,1,-1
      do 11 i=1,q
       if (i.eq.j) goto 11
       if (a(j,i).ne.0.0) goto 12
   11 continue
       m=q
       iexch=1
       goto 900
   12 continue
c     Search for columns isolating an eigenvalue and push them left.
  200 do 22 j=p,q
      do 21 i=p,q
       if (i.eq.j) goto 21
       if (a(i,j).ne.0.0) goto 22
   21 continue
       m=p
       iexch=2
       goto 900
   22 continue
c     Now balance the submatrix in rows p to q.
      do 23 i=p,q
       scale(i)=1.0
   23 continue
c     Iterative loop for norm reduction.
  300 notcnv=.false.
      do 38 i=p,q
       c=0.0
       r=0.0
      do 31 j=p,q
       if (j.eq.i) goto 31
       c=c+abs(a(j,i))
       r=r+abs(a(i,j))
   31 continue
       g=r/RADIX
       f=1.0
       s=c+r
   32  if (c.lt.g) then
        f=f*RADIX
        c=c*SQRDX
        goto 32
       endif
       g=r*RADIX
   34  if (c.ge.g) then
        f=f/RADIX
        c=c/SQRDX
        goto 34
       endif
c     Now balance.
       if ((c+r)/f.lt.0.95*s) then
        g=1.0/f
        scale(i)=scale(i)*f
        notcnv=.true.
       do 36 j=p,n
        a(i,j)=a(i,j)*g
   36  continue
       do 37 j=1,q
        a(j,i)=a(j,i)*f
   37  continue
       endif
   38 continue
      if (notcnv) goto 300
  500 low=p
      upp=q
      return
c     In-line procedure for row and column exchange.
  900 scale(m)=j
      if (j.ne.m) then
       do 90 i=1,q
        swap=a(i,j)
        a(i,j)=a(i,m)
        a(i,m)=swap
   90  continue
       do 95 i=p,n
        swap=a(j,i)
        a(j,i)=a(m,i)
        a(m,i)=swap
   95  continue
      endif
      if (iexch.eq.1) then
       if (q.eq.1) goto 500
       q=q-1
       goto 100
      else if (iexch.eq.2) then
       p=p+1
       goto 200
      endif
      END



      SUBROUTINE BalBak(nm,n,low,upp,scale,m,u)
      INTEGER nm,n,low,upp,m
      REAL scale(n),u(nm,m)
c     GETS nm,n,low,upp,scale,m,u
c     GIVES u
c     Transforms backward the right-hand eigenvectors of the balanced
c     matrix into the eigenvectors of the original matrix.
c        Input: low,upp indicate the boundary indices for the balanced
c     matrix (output from Balanc).
c        Input: scale contains information about the similarity
c     transformations (output from Balanc).
c        Input: m is equal to the number of eigenvectors to be
c     transformed.
c        Input: the array u (n must not be greater than nm), in the
c     first m columns, contains the real and imaginary parts of the
c     eigenvectors to be transformed.
c        Output: the array u, in the first m columns, contains the real
c     and imaginary parts of the transformed eigenvectors.
c     :::
      INTEGER i,j,k,l
      REAL s,swap
      if (upp.ne.low) then
       do 12 i=low,upp
        s=scale(i)
c     Left-hand eigenvectors are back transformed if the foregoing
c     statement is replaced by s=1.0/scale(i).
       do 11 j=1,m
        u(i,j)=u(i,j)*s
   11  continue
   12  continue
      endif
      do 22 l=1,n
       if (l.gt.upp) then
        i=l
       else if (l.lt.low) then
        i=low-l
       else
        goto 22
       endif
       s=scale(i)
       k=int(s)
       if (k.ne.i) then
       do 21 j=1,m
        swap=u(i,j)
        u(i,j)=u(k,j)
        u(k,j)=swap
   21  continue
       endif
   22 continue
      return
      END



      SUBROUTINE CalcH(r,rho,cp,ep,epsp,h)
      REAL r,rho,cp(3,3,3,3),ep(3,3,3),epsp(3,3),h(8,8)
c     GETS r,rho,cp,ep,epsp
c     GIVES h
c     USES DefU,DefV,DefW,DefX,DefY,CaMNPQ,DefH
c     :::
      REAL u(4,4),v(4,4),w(4,4),x(4,4),y(4,4),
     :   mt(4,4),nt(4,4),pt(4,4),qt(4,4)
      call DefU(cp,ep,epsp,u)
      call DefV(cp,ep,epsp,v)
      call DefW(rho,w)
      call DefX(cp,ep,epsp,x)
      call DefY(cp,ep,epsp,y)
      call CaMNPQ(u,v,w,x,y,r,mt,nt,pt,qt)
      call DefH(mt,nt,pt,qt,h)
      return
      END



      SUBROUTINE CalcPF(r,cp,ep,epsp,q,f,psi,utn,ttn,fail)
      INTEGER fail
      REAL r,cp(3,3,3,3),ep(3,3,3),epsp(3,3),psi
      COMPLEX q(8),f(8,8),utn(4),ttn(4)
c     GETS r,cp,ep,epsp,q,f
c     GIVES psi,utn,ttn,fail
c     USES Decomp,DefP
c     :::
      INTEGER a,k,l
      REAL OMEGA,TAU,MEPS,ps(3),psm,pss
      COMPLEX ru(4,4),lu(4,4),rl(4,4),ll(4,4),
     :   cl(4),pet(3,8,8),psl(3),utl(4)
      PARAMETER (OMEGA=1.0,TAU=-1.0,MEPS=0.596e-7)
      call Decomp(f,ru,lu,rl,ll,fail)
      if (fail.ne.0) return
      do 10 l=1,4
       cl(l)=-ll(l,4)*TAU
   10 continue
c     ...
      call DefP(r,cp,ep,epsp,q,f,pet)
      do 23 a=1,2
       psl(a)=0.0
      do 22 k=5,8
      do 21 l=5,8
       psl(a)=psl(a)-conjg(cl(k-4))*cl(l-4)*pet(a,k,l)/
     :    (cmplx(0.0,-r*OMEGA)*(conjg(q(k))-q(l)))
   21 continue
   22 continue
   23 continue
c     ...
      ps(1)=real(psl(1))
      ps(2)=real(psl(2))
      psm=sqrt(ps(1)**2+ps(2)**2)
      if (ps(1).eq.0.0.and.ps(2).eq.0.0) then
       psi=0.0
      else
       psi=atan2(ps(2),ps(1))
      endif
      pss=sqrt(psm)
c     ...
      do 32 k=1,4
       utl(k)=0.0
      do 31 l=1,4
       utl(k)=utl(k)-rl(k,l)*cl(l)
   31 continue
   32 continue
c     ...
      do 35 k=1,4
       utn(k)=utl(k)/(pss*cmplx(0.0,r*OMEGA))
       ttn(k)=0.0
   35 continue
       ttn(4)=TAU/pss
      return
      END



      SUBROUTINE CaMNPQ(u,v,w,x,y,r,mt,nt,pt,qt)
      REAL u(4,4),v(4,4),w(4,4),x(4,4),y(4,4),r,
     :   mt(4,4),nt(4,4),pt(4,4),qt(4,4)
c     GETS u,v,w,x,y,r
c     GIVES mt,nt,pt,qt
c     USES MatPro,MatSum,CoeMat,MatInv
c     :::
      REAL n(4,4),ym(4,4),OMEGA
      PARAMETER (OMEGA=1.0)
      call MatInv(v,n)
      call CoeMat(-1.0,n,nt,4)
      call MatPro(nt,u,mt)
      call MatPro(y,mt,ym)
      call MatSum(1/(r*OMEGA)**2,-1.0,-1.0,w,x,ym,pt)
      call MatPro(y,n,qt)
      return
      END



      SUBROUTINE CoeMat(c,a,ap,n)
      INTEGER n
      REAL c,a(n,n),ap(n,n)
c     GETS c,a,n
c     GIVES ap
c     :::
      INTEGER i,j
      do 12 i=1,n
      do 11 j=1,n
       ap(i,j)=c*a(i,j)
   11 continue
   12 continue
      return
      END



      SUBROUTINE Decomp(f,ru,lu,rl,ll,fail)
      INTEGER fail
      COMPLEX f(8,8),ru(4,4),lu(4,4),rl(4,4),ll(4,4)
c     GETS f
c     GIVES ru,lu,rl,ll,fail
c     USES MaInvC
c     :::
      INTEGER i,j
      COMPLEX su(4,4),sl(4,4)
      do 22 j=1,4
      do 21 i=1,4
       ru(i,j)=f(i,j)
       su(i,j)=f(i+4,j)
       rl(i,j)=f(i,j+4)
       sl(i,j)=f(i+4,j+4)
   21 continue
   22 continue
c     ...
      call MaInvC(su,lu,4,4,fail)
      if (fail.ne.0) return
      call MaInvC(sl,ll,4,4,fail)
      if (fail.ne.0) return
      return
      END



      SUBROUTINE DefA(lam,gam,the,a)
      REAL lam,gam,the,a(3,3)
c     GETS lam,gam,the
c     GIVES a
c     The matrix a is the transformation matrix for the Euler angles lam
c     (lambda), gam (gamma), the (theta) given in radians.
c        The transformation from the system of coordinates (X,Y,Z) to
c     the system of coordinates (x,y,z) consists of four consecutive
c     rotations: (0) by the angle of 180 degrees in the yz plane (i.e.
c     about the x axis; this is a new initial system: x=X,y=-Y,z=-Z),
c     (1) by the angle -lam in the xy plane (i.e about the z axis), (2)
c     by the angle gam in the yz plane (i.e. about the x axis), (3) by
c     the angle -the in the xy plane (i.e. about the z axis).
c        The names of Euler angles and the corresponding axes of
c     rotation agree with the convention (A.J.Slobodnik, in "Acoustic
c     Surface Waves", edited by A.A.Oliner, Springer-Verlag 1978, p.
c     269).
c        Alternatively, the transformation from the system (X,Y,Z) to
c     the system (x,y,z) may consist of the following four consecutive
c     rotations: (1) by the angle lam in the xy plane, (2) by the angle
c     gam in the yz plane, (3) by the angle the in the xy plane, (4) by
c     the angle of 180 degrees in the yz plane. The transformation
c     (1)+(2)+(3) gives the system of coordinates (x,y,z) which agrees
c     with the old convention (J.J.Campbell and W.R.Jones, IEEE Trans.
c     Son. Ultrason., vol.15, 209 (1968)). In that convention, the z
c     axis is perpendicular to the piezoelectric half-space and is
c     directed into it (in the new convention, the z axis is directed
c     oppositely).
c     :::
      a(1,1)=cos(the)*cos(lam)-sin(the)*cos(gam)*sin(lam)
      a(1,2)=cos(the)*sin(lam)+sin(the)*cos(gam)*cos(lam)
      a(1,3)=sin(the)*sin(gam)
      a(2,1)=sin(the)*cos(lam)+cos(the)*cos(gam)*sin(lam)
      a(2,2)=sin(the)*sin(lam)-cos(the)*cos(gam)*cos(lam)
      a(2,3)=-cos(the)*sin(gam)
      a(3,1)=-sin(gam)*sin(lam)
      a(3,2)=sin(gam)*cos(lam)
      a(3,3)=-cos(gam)
      return
      END



      SUBROUTINE DefH(mt,nt,pt,qt,h)
      REAL mt(4,4),nt(4,4),pt(4,4),qt(4,4),h(8,8)
c     GETS mt,nt,pt,qt
c     GIVES h
c     :::
      INTEGER i,j
      do 12 i=1,4
      do 11 j=1,4
       h(i,j)=mt(i,j)
       h(i,j+4)=nt(i,j)
       h(i+4,j)=pt(i,j)
       h(i+4,j+4)=qt(i,j)
   11 continue
   12 continue
      return
      END



      SUBROUTINE DefP(r,cp,ep,epsp,q,f,pet)
      REAL r,cp(3,3,3,3),ep(3,3,3),epsp(3,3)
      COMPLEX q(8),f(8,8),pet(3,8,8)
c     GETS r,cp,ep,epsp,q,f
c     GIVES pet
c     :::
      INTEGER a,i,j,k,l
      COMPLEX chi(3,4,4,8,8)
      do 15 a=1,3
      do 14 k=1,8
      do 13 l=1,8
      do 12 i=1,3
      do 11 j=1,3
       chi(a,i,j,k,l)=-cp(i,a,j,1)-q(l)*cp(i,a,j,3)
   11 continue
       chi(a,i,4,k,l)=-ep(1,i,a)-q(l)*ep(3,i,a)
     :    -ep(a,i,1)-conjg(q(k))*ep(a,i,3)
       chi(a,4,i,k,l)=0.0
   12 continue
       chi(a,4,4,k,l)=epsp(a,1)+conjg(q(k))*epsp(a,3)
   13 continue
   14 continue
   15 continue
      do 25 a=1,3
      do 24 k=1,8
      do 23 l=1,8
       pet(a,k,l)=0.0
      do 22 i=1,4
      do 21 j=1,4
       pet(a,k,l)=pet(a,k,l)-(0.5/r)*
     :     chi(a,i,j,k,l)*conjg(f(i,k))*f(j,l)
   21 continue
   22 continue
   23 continue
   24 continue
   25 continue
      return
      END



      SUBROUTINE DefU(cp,ep,epsp,u)
      REAL cp(3,3,3,3),ep(3,3,3),epsp(3,3),u(4,4)
c     GETS cp,ep,epsp
c     GIVES u
c     :::
      u(1,1)=cp(1,3,1,1)
      u(1,2)=cp(1,3,2,1)
      u(1,3)=cp(1,3,3,1)
      u(1,4)=ep(1,1,3)
      u(2,1)=cp(2,3,1,1)
      u(2,2)=cp(2,3,2,1)
      u(2,3)=cp(2,3,3,1)
      u(2,4)=ep(1,2,3)
      u(3,1)=cp(3,3,1,1)
      u(3,2)=cp(3,3,2,1)
      u(3,3)=cp(3,3,3,1)
      u(3,4)=ep(1,3,3)
      u(4,1)=ep(3,1,1)
      u(4,2)=ep(3,2,1)
      u(4,3)=ep(3,3,1)
      u(4,4)=-epsp(3,1)
      return
      END



      SUBROUTINE DefV(cp,ep,epsp,v)
      REAL cp(3,3,3,3),ep(3,3,3),epsp(3,3),v(4,4)
c     GETS cp,ep,epsp
c     GIVES v
c     :::
      v(1,1)=cp(1,3,1,3)
      v(1,2)=cp(1,3,2,3)
      v(1,3)=cp(1,3,3,3)
      v(1,4)=ep(3,1,3)
      v(2,1)=cp(2,3,1,3)
      v(2,2)=cp(2,3,2,3)
      v(2,3)=cp(2,3,3,3)
      v(2,4)=ep(3,2,3)
      v(3,1)=cp(3,3,1,3)
      v(3,2)=cp(3,3,2,3)
      v(3,3)=cp(3,3,3,3)
      v(3,4)=ep(3,3,3)
      v(4,1)=ep(3,1,3)
      v(4,2)=ep(3,2,3)
      v(4,3)=ep(3,3,3)
      v(4,4)=-epsp(3,3)
      return
      END



      SUBROUTINE DefW(rho,w)
      REAL rho,w(4,4)
c     GETS rho
c     GIVES w
c     :::
      INTEGER i,j
      REAL OMEGA
      PARAMETER (OMEGA=1.0)
      do 12 i=1,4
      do 11 j=1,4
       w(i,j)=0.0
   11 continue
   12 continue
      do 20 i=1,3
       w(i,i)=-(OMEGA**2)*rho
   20 continue
      return
      END



      SUBROUTINE DefX(cp,ep,epsp,x)
      REAL cp(3,3,3,3),ep(3,3,3),epsp(3,3),x(4,4)
c     GETS cp,ep,epsp
c     GIVES x
c     :::
      x(1,1)=-cp(1,1,1,1)
      x(1,2)=-cp(1,1,2,1)
      x(1,3)=-cp(1,1,3,1)
      x(1,4)=-ep(1,1,1)
      x(2,1)=-cp(2,1,1,1)
      x(2,2)=-cp(2,1,2,1)
      x(2,3)=-cp(2,1,3,1)
      x(2,4)=-ep(1,2,1)
      x(3,1)=-cp(3,1,1,1)
      x(3,2)=-cp(3,1,2,1)
      x(3,3)=-cp(3,1,3,1)
      x(3,4)=-ep(1,3,1)
      x(4,1)=-ep(1,1,1)
      x(4,2)=-ep(1,2,1)
      x(4,3)=-ep(1,3,1)
      x(4,4)=epsp(1,1)
      return
      END



      SUBROUTINE DefY(cp,ep,epsp,y)
      REAL cp(3,3,3,3),ep(3,3,3),epsp(3,3),y(4,4)
c     GETS cp,ep,epsp
c     GIVES y
c     :::
      y(1,1)=-cp(1,1,1,3)
      y(1,2)=-cp(1,1,2,3)
      y(1,3)=-cp(1,1,3,3)
      y(1,4)=-ep(3,1,1)
      y(2,1)=-cp(2,1,1,3)
      y(2,2)=-cp(2,1,2,3)
      y(2,3)=-cp(2,1,3,3)
      y(2,4)=-ep(3,2,1)
      y(3,1)=-cp(3,1,1,3)
      y(3,2)=-cp(3,1,2,3)
      y(3,3)=-cp(3,1,3,3)
      y(3,4)=-ep(3,3,1)
      y(4,1)=-ep(1,1,3)
      y(4,2)=-ep(1,2,3)
      y(4,3)=-ep(1,3,3)
      y(4,4)=epsp(1,3)
      return
      END



      SUBROUTINE DefZ(f,zet,fail)
      INTEGER fail
      COMPLEX f(8,8),zet
c     GETS f
c     GIVES zet,fail
c     USES Decomp
c     :::
      INTEGER i
      REAL EPS0
      COMPLEX ru(4,4),lu(4,4),rl(4,4),ll(4,4),zu,zl
      PARAMETER (EPS0=8.85434e-3)
      call Decomp(f,ru,lu,rl,ll,fail)
      if (fail.ne.0) return
      zu=0.0
      zl=0.0
      do 31 i=1,4
       zu=zu+ru(4,i)*lu(i,4)
       zl=zl+rl(4,i)*ll(i,4)
   31 continue
      zet=-zl
      return
      END



      FUNCTION Det(a,n)
      INTEGER n
      REAL Det,a(n,n)
c     GETS a,n
c     GIVES Det
c     Calculates the determinant Det of the matrix a for n=2,3,4.
c     :::
      INTEGER i(4),j,j1,j2,j3,j4,k,k2,k3,k4,l,l3,l4,swap
      REAL sj,sk,sl,sign
      Det=0.0
      sj=1.0
      j1=1
      j2=2
      j3=3
      j4=4
      do 13 j=1,n
       sign=sj
       i(1)=j1
       i(2)=j2
       i(3)=j3
       i(4)=j4
       swap=i(1)
       i(1)=i(j)
       i(j)=swap
       if (j.ne.1) sign=-sign
       if (n.eq.2) then
        Det=Det+sign*a(1,i(1))*a(2,i(2))
        goto 102
       endif
       sk=sign
       k2=i(2)
       k3=i(3)
       k4=i(4)
      do 12 k=2,n
       sign=sk
       i(2)=k2
       i(3)=k3
       i(4)=k4
       swap=i(2)
       i(2)=i(k)
       i(k)=swap
       if (k.ne.2) sign=-sign
       if (n.eq.3) then
        Det=Det+sign*a(1,i(1))*a(2,i(2))*a(3,i(3))
        goto 101
       endif
       sl=sign
       l3=i(3)
       l4=i(4)
      do 11 l=3,n
       sign=sl
       i(3)=l3
       i(4)=l4
       swap=i(3)
       i(3)=i(l)
       i(l)=swap
       if (l.ne.3) sign=-sign
       Det=Det+sign*a(1,i(1))*a(2,i(2))*a(3,i(3))*a(4,i(4))
   11 continue
  101 continue
   12 continue
  102 continue
   13 continue
      return
      END



      SUBROUTINE ElmHes(nm,n,low,upp,a,inter)
      INTEGER nm,n,low,upp,inter(n)
      REAL a(nm,n)
c     GETS nm,n,low,upp,a
c     GIVES a,inter
c     Reduces a real general matrix to upper Hessenberg form using
c     stabilized elementary similarity transformations. This reduced
c     form is used by other subroutines to find the eigenvalues and/or
c     eigenvectors of the original matrix.
c        Input: low,upp indicate the boundary indices for the balanced
c     matrix (output from Balanc; if Balanc is not used then low=1,
c     upp=n). ElmHes actually reduces the submatrix of rank upp-low+1 of
c     the balanced matrix (starting at the element (low,low) and
c     finishing at the element (upp,upp)).
c        Input: the array a (n must not be greater than nm) contains the
c     matrix of rank n to be reduced to the Hessenberg form.
c        Output: the array a contains the upper Hessenberg matrix as
c     well as (in the remaining triangle below) the multipliers used in
c     the reduction.
c        Output: inter identifies the rows and columns interchanged
c     during the reduction; the information is stored in the sub-array
c     inter(low:upp).
c     :::
      INTEGER i,j,m
      REAL x,y,swap
      do 20 m=low+1,upp-1
       x=0.0
       i=m
      do 10 j=m,upp
       if (abs(a(j,m-1)).gt.abs(x)) then
        x=a(j,m-1)
        i=j
       endif
   10 continue
       inter(m)=i
       if (i.ne.m) then
c     Interchange rows and columns of a.
        do 11 j=m-1,n
         swap=a(i,j)
         a(i,j)=a(m,j)
         a(m,j)=swap
   11   continue
        do 12 j=1,upp
         swap=a(j,i)
         a(j,i)=a(j,m)
         a(j,m)=swap
   12   continue
c     End interchange.
       endif
       if (x.ne.0.0) then
       do 15 i=m+1,upp
        y=a(i,m-1)
        if (y.ne.0.0) then
         y=y/x
         a(i,m-1)=y
        do 13 j=m,n
         a(i,j)=a(i,j)-y*a(m,j)
   13   continue
        do 14 j=1,upp
         a(j,m)=a(j,m)+y*a(j,i)
   14   continue
        endif
   15  continue
       endif
   20 continue
      return
      END



      SUBROUTINE ElmTra(nm,n,low,upp,a,inter,u)
      INTEGER nm,n,low,upp,inter(n)
      REAL a(nm,n),u(nm,n)
c     GETS nm,n,low,upp,a,inter
c     GIVES u
c     Accumulates the elementary similarity transformations used in the
c     reduction of a real general matrix to upper Hessenberg form by
c     ElmHes, and finds the resulting transformation matrix.
c        Input: low,upp indicate the boundary indices for the balanced
c     matrix (output from Balanc; if Balanc is not used then low=1,
c     upp=n).
c        Input: the array a (n must not be greater than nm) contains the
c     multipliers which were used in the reduction to the Hessenberg
c     form (in the lower triangle below the subdiagonal; the remaining
c     upper part of a may be arbitrary); this is output from ElmHes.
c        Input: inter identifies the rows and columns interchanged
c     during the reduction to the Hessenberg form (output from ElmHes).
c        Output: the array u (n must not be greater than nm) contains
c     the transformation matrix which performs the similarity reduction
c     from the original matrix to the Hessenberg matrix.
c     :::
      INTEGER i,j,k
c     Initialize u to identity matrix.
      do 12 i=1,n
      do 11 j=1,n
       u(i,j)=0.0
   11 continue
       u(i,i)=1.0
   12 continue
      do 25 i=upp-1,low+1,-1
      do 21 k=i+1,upp
       u(k,i)=a(k,i-1)
   21 continue
       j=inter(i)
       if (i.ne.j) then
       do 22 k=i,upp
        u(i,k)=u(j,k)
        u(j,k)=0.0
   22  continue
        u(j,i)=1.0
       endif
   25 continue
      return
      END



      SUBROUTINE FindVV(h,qn,fn,fail)
      INTEGER fail,NDIM
      PARAMETER (NDIM=8)
      REAL h(NDIM,NDIM)
      COMPLEX qn(NDIM),fn(NDIM,NDIM)
c     GETS h
c     GIVES qn,fn,fail
c     USES Balanc,ElmHes,ElmTra,HQR2,BalBak,UnPack,Renum
c     Finds the eigenvalues qn and the corresponding eigenvectors of the
c     matrix h. Eigenvectors of h are columns of the matrix fn.
c        The subroutines which are used here (except Renum) are taken
c     from EISPACK set of subroutines (B.T.Smith et al., Lect. Notes
c     Computer Science, vol.6 (1974)). Some changes have been made (the
c     original subroutines are written in Fortran 66). They consist
c     mainly in: (1) using the block statement "if" which eliminates a
c     great number of labels related to the statement "goto", (2) using
c     the statement "do" with a decreasing control variable, (3)
c     adjusting statements "do" to the standard of Fortran 77 (pre-test
c     instead of post-test).
c        A major change has been made in the subroutine HQR2, which
c     consists in inserting several lines of code (labelled with
c     five-digit numbers) with the statement "if". In these lines,
c     values of some variables are tested for underflow or overflow at
c     the processor level. This enables to stop calculations for a
c     particular crystal orientation without stopping the program.
c        Comments have been added which are based on the cited
c     publication and on a book by Wilkinson and Reinsch (J.H.Wilkinson
c     and C.Reinsch, "Linear Algebra", Springer-Verlag 1971). The book
c     contains algorithms and their implementations in the form of Algol
c     procedures (EISPACK subroutines are translations of these
c     procedures into Fortran 66).
c     :::
      INTEGER i,j,low,upp,inter(NDIM),ncall
      REAL ht(NDIM,NDIM),scale(NDIM),u(NDIM,NDIM),
     :   qr(NDIM),qi(NDIM),ur(NDIM,NDIM),ui(NDIM,NDIM)
      COMPLEX q(NDIM),f(NDIM,NDIM)
      COMMON ncall
c     ...
      ncall=ncall+1
c     ...
      do 12 i=1,NDIM
      do 11 j=1,NDIM
       ht(i,j)=h(i,j)
   11 continue
   12 continue
      call Balanc(NDIM,NDIM,ht,low,upp,scale)
      call ElmHes(NDIM,NDIM,low,upp,ht,inter)
      call ElmTra(NDIM,NDIM,low,upp,ht,inter,u)
      call HQR2(NDIM,NDIM,low,upp,ht,qr,qi,u,fail)
      if (fail.ne.0) return
      call BalBak(NDIM,NDIM,low,upp,scale,NDIM,u)
      call UnPack(NDIM,NDIM,qi,u,ur,ui)
      do 21 j=1,8
       q(j)=cmplx(qr(j),qi(j))
   21 continue
      do 32 i=1,8
      do 31 j=1,8
       f(i,j)=cmplx(ur(i,j),ui(i,j))
   31 continue
   32 continue
      call Renum(q,f,qn,fn,fail)
      return
      END



      SUBROUTINE GaussC(a,n,ns,b,m,ms,fail)
      INTEGER m,ms,n,ns,fail
      COMPLEX a(ns,ns),b(ns,ms)
c     GETS a,n,ns,b,m,ms
c     GIVES a,fail
c     :::
      INTEGER NMAX
      PARAMETER (NMAX=10)
      INTEGER i,j,k,l,p,icol,irow,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
      REAL big
      COMPLEX dum,pivinv
      fail=0
      do 11 j=1,n
       ipiv(j)=0
   11 continue
      do 22 i=1,n
       big=0.
      do 13 j=1,n
       if (ipiv(j).eq.1) goto 13
      do 12 k=1,n
       if (ipiv(k).eq.0) then
        if (cabs(a(j,k)).ge.big) then
         big=cabs(a(j,k))
         irow=j
         icol=k
        endif
       else if (ipiv(k).gt.1) then
        goto 999
       endif
   12 continue
   13 continue
       ipiv(icol)=ipiv(icol)+1
       if (irow.ne.icol) then
        do 14 l=1,n
         dum=a(irow,l)
         a(irow,l)=a(icol,l)
         a(icol,l)=dum
   14   continue
        do 15 l=1,m
         dum=b(irow,l)
         b(irow,l)=b(icol,l)
         b(icol,l)=dum
   15   continue
       endif
       indxr(i)=irow
       indxc(i)=icol
       if (a(icol,icol).eq.0.) goto 999 
       pivinv=1./a(icol,icol)
       a(icol,icol)=1.
      do 16 l=1,n
       a(icol,l)=a(icol,l)*pivinv
   16 continue
      do 17 l=1,m
       b(icol,l)=b(icol,l)*pivinv
   17 continue
      do 21 p=1,n
       if (p.eq.icol) goto 21
       dum=a(p,icol)
       a(p,icol)=0.
      do 18 l=1,n
       a(p,l)=a(p,l)-a(icol,l)*dum
   18 continue
      do 19 l=1,m
       b(p,l)=b(p,l)-b(icol,l)*dum
   19 continue
   21 continue
   22 continue
      do 24 l=n,1,-1
       if (indxr(l).eq.indxc(l)) goto 24
      do 23 k=1,n
       dum=a(k,indxr(l))
       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
   23 continue
   24 continue
      return
c     Set fail: Singular matrix occured.
  999 fail=1001
      return
      END



      SUBROUTINE HQR2(nm,n,low,upp,h,wr,wi,u,fail)
      INTEGER nm,n,low,upp,fail
      REAL h(nm,n),wr(n),wi(n),u(nm,n)
c     GETS nm,n,low,upp,h,u
c     GIVES wr,wi,fail
c     Computes the eigenvalues and eigenvectors of a real upper
c     Hessenberg matrix using the QR method. The eigenvectors of a real
c     general matrix can also be computed if ElmHes has been used to
c     reduce this general matrix to Hessenberg form and ElmTra has been
c     used to find the transformation matrix from the general to
c     Hessenberg form.
c        Input: low,upp indicate the boundary indices for the balanced
c     matrix (output from Balanc; if Balanc is not used then low=1,
c     upp=n).
c        Input: The array h (n must not be greater than nm), in its
c     relevant part, contains the upper Hessenberg matrix; this is
c     output from ElmHes. (Note: HQR2 destroys this matrix on output.)
c        Input: the array u contains the transformation matrix produced
c     by ElmTra (this matrix reduces the real general matrix to
c     Hessenberg form). If the original matrix is an upper Hessenberg
c     matrix then the array should contain the identity matrix of rank
c     n.
c        Output: the arrays wr,wi contain the real and imaginary parts,
c     respectively, of the eigenvalues of the original matrix. The
c     eigenvalues are unordered except that the eigenvalues which are a
c     complex conjugate pair appear consecutively, with the eigenvalue
c     having the positive imaginary part first.
c        Output: the array u contains the real and imaginary parts of
c     the eigenvectors of the original matrix (or balanced matrix if
c     Balanc is used). If the i-th eigenvalue is real, the i-th column
c     of u contains the corresponding real eigenvector. If the
c     eigenvalues i and i+1 are a complex conjugate pair, the columns i
c     and i+1 of u contain the real and imaginary part of the
c     eigenvector corresponding to the eigenvalue with the positive
c     imaginary part; the conjugate of this vector is the eigenvector of
c     the conjugate eigenvalue. The eigenvectors are not normalized.
c        Output: fail is set to zero if calculations are completed
c     successfully; if more than 30 iterations are required to determine
c     an eigenvalue, this subroutine terminates with fail set to the
c     index of the eigenvalue for which the failure occurs (the
c     eigenvalues in the arrays wr and wi should be correct for indices
c     fail+1,fail+2,...,fail+n but no eigenvectors are computed). If
c     underflow or overflow occurs then the subroutine terminates with
c     fail set to 1000.
c        MEPS is the relative machine precision defined as the smallest
c     positive number such that 1+MEPS>1 (when evaluated by the
c     computer). A (roughly) equivalent definition: MEPS=2**(-m) where m
c     is the amount of internal storage (in bits) allocated to the
c     mantissa of a number in the floating point representation.
c        MNUL is the machine null defined as the smallest positive
c     number such that MNUL*2>MNUL (when evaluated by the computer; the
c     inequality is not satisfied if underflow occurs on both sides).
c        MINF is the machine infinity defined as the greatest positive
c     number such that MINF/2<MINF (when evaluated by the computer; the
c     inequality is not satisfied if overflow occurs on both sides).
c     :::
      INTEGER i,j,k,l,m,na,nb,nc,iter
      REAL p,q,r,s,t,w,x,y,z,ra,sa,vr,vi,norm,d(2),
     :   MEPS,MNUL,MINF,nul,inf
      COMPLEX cdiv
      LOGICAL notlst
      PARAMETER (MEPS=0.596e-7,MNUL=0.118e-37,MINF=0.170e39)
      EQUIVALENCE (cdiv,d)
      nul=sqrt(MNUL)
      inf=sqrt(MINF)
      fail=0
c     Store roots isolated by Balanc.
      do 01 i=1,n
       if ((i.ge.low).and.(i.le.upp)) goto 01
       wr(i)=h(i,i)
       wi(i)=0.0
   01 continue
      nb=upp
      t=0.0
c     Search for next eigenvalues.
  100 if (nb.lt.low) goto 500
      iter=0
      na=nb-1
      nc=na-1
c     Look for single small sub-diagonal element.
  200 do 20 l=nb,low,-1
       if (l.eq.low) goto 210
       if (abs(h(l,l-1)).le.MEPS*(abs(h(l-1,l-1))+abs(h(l,l)))) goto 210
   20 continue
c     Form shift.
  210 x=h(nb,nb)
      if (l.eq.nb) goto 310
      y=h(na,na)
      w=h(nb,na)*h(na,nb)
      if (l.eq.na) goto 320
      if (iter.eq.30) goto 999
      if ((iter.eq.10).or.(iter.eq.20)) then
c     Form exceptional shift.
       t=t+x
       do 21 i=low,nb
        h(i,i)=h(i,i)-x
   21  continue
       s=abs(h(nb,na))+abs(h(na,nc))
       x=0.75*s
       y=x
       w=-0.4375*s**2
      endif
      iter=iter+1
c     Look for two consecutive small sub-diagonal elements.
      do 22 m=nc,l,-1
       z=h(m,m)
       r=x-z
       s=y-z
       p=(r*s-w)/h(m+1,m)+h(m,m+1)
       q=h(m+1,m+1)-z-r-s
       r=h(m+2,m+1)
       s=abs(p)+abs(q)+abs(r)
       p=p/s
       q=q/s
       r=r/s
       if (m.eq.l) goto 220
       if (abs(h(m,m-1))*(abs(q)+abs(r)).le.
     :    MEPS*abs(p)*(abs(h(m-1,m-1))+abs(z)+abs(h(m+1,m+1))))
     :    goto 220
   22 continue
  220 do 23 i=m+2,nb
       h(i,i-2)=0.0
       if (i.eq.m+2) goto 23
       h(i,i-3)=0.0
   23 continue
c     Double QR step involving rows l to nb and columns m to nb.
      do 29 k=m,na
       notlst=k.ne.na
       if (k.ne.m) then
        p=h(k,k-1)
        q=h(k+1,k-1)
        r=0.0
        if (notlst) r=h(k+2,k-1)
c     Prevent underflow.
10011   if (abs(p).le.nul) p=0.0
10012   if (abs(q).le.nul) q=0.0
10013   if (abs(r).le.nul) r=0.0
c     Prevent overflow.
10021   if (abs(p).ge.inf) goto 9999
10022   if (abs(q).ge.inf) goto 9999
10023   if (abs(r).ge.inf) goto 9999
        x=abs(p)+abs(q)+abs(r)
        if (x.eq.0.0) goto 29
        p=p/x
        q=q/x
        r=r/x
       endif
c     Prevent underflow.
20011  if ((abs(p).ne.0.0).and.(abs(p).le.nul)) goto 9999
20012  if ((abs(q).ne.0.0).and.(abs(q).le.nul)) goto 9999
20013  if ((abs(r).ne.0.0).and.(abs(r).le.nul)) goto 9999
       s=sign(sqrt(p**2+q**2+r**2),p)
       if (k.ne.m) then
        h(k,k-1)=-s*x
       else
        if (l.ne.m) h(k,k-1)=-h(k,k-1)
       endif
       p=p+s
       x=p/s
       y=q/s
       z=r/s
       q=q/p
       r=r/p
c     Row modification.
      do 24 j=k,n
       p=h(k,j)+q*h(k+1,j)
       if (notlst) then
        p=p+r*h(k+2,j)
        h(k+2,j)=h(k+2,j)-p*z
       endif
       h(k+1,j)=h(k+1,j)-p*y
       h(k,j)=h(k,j)-p*x
   24 continue
       j=min(nb,k+3)
c     Column modification.
      do 25 i=1,j
       p=x*h(i,k)+y*h(i,k+1)
       if (notlst) then
        p=p+z*h(i,k+2)
        h(i,k+2)=h(i,k+2)-p*r
       endif
       h(i,k+1)=h(i,k+1)-p*q
       h(i,k)=h(i,k)-p
   25 continue
c     Accumulate transformations.
      do 26 i=low,upp
       p=x*u(i,k)+y*u(i,k+1)
       if (notlst) then
        p=p+z*u(i,k+2)
        u(i,k+2)=u(i,k+2)-p*r
       endif
       u(i,k+1)=u(i,k+1)-p*q
       u(i,k)=u(i,k)-p
   26 continue
   29 continue
      goto 200
c     One root found.
  310 h(nb,nb)=x+t
      wr(nb)=h(nb,nb)
      wi(nb)=0.0
      nb=na
      goto 100
c     Two roots found.
  320 p=(y-x)/2.0
c     Prevent underflow.
30011 if ((abs(p).ne.0.0).and.(abs(p).le.nul)) goto 9999
30012 if (abs(w).le.MNUL) goto 9999
c     Prevent overflow.
30021 if (abs(p).ge.inf) goto 9999
30022 if (abs(w).ge.MINF) goto 9999
      q=p**2+w
      z=sqrt(abs(q))
      h(nb,nb)=x+t
      x=h(nb,nb)
      h(na,na)=y+t
      if (q.ge.0.0) then
c     Real pair.
       z=p+sign(z,p)
       wr(na)=x+z
       wr(nb)=wr(na)
       if (z.ne.0.0) wr(nb)=x-w/z
       wi(na)=0.0
       wi(nb)=0.0
       x=h(nb,na)
c     Prevent underflow.
40011  if ((abs(x).ne.0.0).and.(abs(x).le.nul)) goto 9999
40012  if ((abs(z).ne.0.0).and.(abs(z).le.nul)) goto 9999
c     Prevent overflow.
40021  if (abs(x).ge.inf) goto 9999
40022  if (abs(z).ge.inf) goto 9999
       r=sqrt(x**2+z**2)
       p=x/r
       q=z/r
c     Row modification.
       do 31 j=na,n
        z=h(na,j)
        h(na,j)=q*z+p*h(nb,j)
        h(nb,j)=q*h(nb,j)-p*z
   31  continue
c     Column modification.
       do 32 i=1,nb
        z=h(i,na)
        h(i,na)=q*z+p*h(i,nb)
        h(i,nb)=q*h(i,nb)-p*z
   32  continue
c     Accumulate transformations.
       do 33 i=low,upp
        z=u(i,na)
        u(i,na)=q*z+p*u(i,nb)
        u(i,nb)=q*u(i,nb)-p*z
   33  continue
      else
c     Complex pair.
       wr(na)=x+p
       wr(nb)=x+p
       wi(na)=z
       wi(nb)=-z
      endif
      nb=nc
      goto 100
c     All roots found. Backsubstitute to find vectors of upper
c     triangular form.
  500 norm=0.0
      k=1
      do 51 i=1,n
      do 50 j=k,n
       norm=norm+abs(h(i,j))
   50 continue
       k=i
   51 continue
      if (norm.eq.0.0) return
      do 59 nb=n,1,-1
       p=wr(nb)
       q=wi(nb)
       na=nb-1
       if (q.lt.0.0) goto 510
       if (q.gt.0.0) goto 59
c     Real vector.
       m=nb
       h(nb,nb)=1.0
      do 53 i=na,1,-1
       w=h(i,i)-p
       r=h(i,nb)
      do 52 j=m,na
       r=r+h(i,j)*h(j,nb)
   52 continue
       if (wi(i).lt.0.0) then
        z=w
        s=r
        goto 53
       endif
       m=i
       if (wi(i).eq.0.0) then
        t=w
        if (w.eq.0.0) t=MEPS*norm
        h(i,nb)=-r/t
        goto 53
       endif
c     Solve real equations.
       x=h(i,i+1)
       y=h(i+1,i)
       q=(wr(i)-p)**2+wi(i)**2
       t=(x*s-z*r)/q
       h(i,nb)=t
       if (abs(x).gt.abs(z)) then
        h(i+1,nb)=(-r-w*t)/x
        goto 53
       endif
       h(i+1,nb)=(-s-y*t)/z
   53 continue
c     End real vector.
       goto 59
c     Complex vector.
  510  m=na
c     Last vector component chosen imaginary so that eigenvector matrix
c     is triangular.
       if (abs(h(nb,na)).gt.abs(h(na,nb))) then
        h(na,na)=q/h(nb,na)
        h(na,nb)=-(h(nb,nb)-p)/h(nb,na)
       else
        cdiv=cmplx(0.0,-h(na,nb))/cmplx(h(na,na)-p,q)
        h(na,na)=d(1)
        h(na,nb)=d(2)
       endif
       h(nb,na)=0.0
       h(nb,nb)=1.0
       nc=na-1
      do 55 i=na-1,1,-1
       w=h(i,i)-p
       ra=0.0
       sa=h(i,nb)
      do 54 j=m,na
       ra=ra+h(i,j)*h(j,na)
       sa=sa+h(i,j)*h(j,nb)
   54 continue
       if (wi(i).lt.0.0) then
        z=w
        r=ra
        s=sa
        goto 55
       endif
       m=i
       if (wi(i).eq.0.0) then
        cdiv=cmplx(-ra,-sa)/cmplx(w,q)
        h(i,na)=d(1)
        h(i,nb)=d(2)
        goto 55
       endif
c     Solve complex equations.
       x=h(i,i+1)
       y=h(i+1,i)
       vr=(wr(i)-p)**2+wi(i)**2-q**2
       vi=(wr(i)-p)*2.0*q
       if ((vr.eq.0.0).and.(vi.eq.0.0))
     :    vr=MEPS*norm*(abs(w)+abs(q)+abs(x)+abs(y)+abs(z))
       cdiv=cmplx(x*r-z*ra+q*sa,x*s-z*sa-q*ra)/cmplx(vr,vi)
       h(i,na)=d(1)
       h(i,nb)=d(2)
       if (abs(x).gt.(abs(z)+abs(q))) then
        h(i+1,na)=(-ra-w*h(i,na)+q*h(i,nb))/x
        h(i+1,nb)=(-sa-w*h(i,nb)-q*h(i,na))/x
       else
        cdiv=cmplx(-r-y*h(i,na),-s-y*h(i,nb))/cmplx(z,q)
        h(i+1,na)=d(1)
        h(i+1,nb)=d(2)
       endif
   55 continue
c     End complex vector.
   59 continue
c     End back substitution.
c     Vectors of isolated roots.
      do 62 i=1,n
       if ((i.ge.low).and.(i.le.upp)) goto 62
      do 61 j=i,n
       u(i,j)=h(i,j)
   61 continue
   62 continue
c     Multiply by transformation matrix to give vectors of original full
c     matrix.
      do 73 j=n,low,-1
       m=min(j,upp)
      do 72 i=low,upp
       z=0.0
      do 71 k=low,m
       z=z+u(i,k)*h(k,j)
   71 continue
       u(i,j)=z
   72 continue
   73 continue
      return
c     Set fail: 30 iterations failed to isolate the current eigevalue.
  999 fail=nb
      return
c     Set fail: Underflow or overflow occured.
 9999 fail=1000
      return
      END



      SUBROUTINE Look(rho,cp,ep,epsp,rlow,rupp,rmid,Test,zinf,dzet,fail)
      INTEGER fail
      REAL rho,cp(3,3,3,3),ep(3,3,3),epsp(3,3),
     :   rlow,rupp,rmid,Test,zinf,dzet
      EXTERNAL Test
c     GETS rho,cp,ep,epsp,rlow,rupp,Test
c     GIVES rmid,zinf,dzet,fail
c     USES Seek
c     :::
      REAL rlo,rup,STEP
      PARAMETER (STEP=0.00001)
      rlo=rlow
      rup=rupp
   11 continue
       rlow=rlo
       rupp=rup
       call Seek(rho,cp,ep,epsp,rlow,rupp,rmid,Test,zinf,dzet,fail)
       if (fail.eq.0) goto 12
       rlo=rlo+STEP
       if (rlo.gt.rup) return
       goto 11
   12 continue
      return
      END



      SUBROUTINE MaInvC(a,ap,n,ns,fail)
      INTEGER n,ns,fail
      COMPLEX a(ns,ns),ap(ns,ns)
c     GETS a,n,ns,fail
c     GIVES ap,fail
c     USES GaussC
c     Calculates the inverse ap of the matrix a.
c        The subroutine GaussC is a complex version of the subroutine
c     GaussJ which has been taken from a collection of procedures by
c     Press and others (W.H.Press et al., "Numerical Recipes in Fortran.
c     The Art of Scientific Computing. Second Edition", Cambridge
c     University Press 1992).
c     :::
      INTEGER i,j,NMAX
      PARAMETER (NMAX=10)
      COMPLEX b(NMAX,1)
      do 12 i=1,n
      do 11 j=1,n
       ap(i,j)=a(i,j)
   11 continue
   12 continue
      do 20 i=1,n
       b(i,1)=0.0
   20 continue
      call GaussC(ap,n,ns,b,1,1,fail)
      return
      END



      SUBROUTINE MatInv(a,z)
      REAL a(4,4),z(4,4)
c     GETS a
c     GIVES z
c     USES Det
c     Calculates the inverse z of the matrix a.
c     :::
      INTEGER i,j,k,l
      REAL m(3,3),b(4,4),sign,Det,d
      LOGICAL row1,row2,col1,col2
      do 14 i=1,4
      do 13 j=1,4
      do 12 k=1,3
      do 11 l=1,3
       row1=k.lt.i
       row2=.not.row1
       col1=l.lt.j
       col2=.not.col1
       if (row1.and.col1) then
        m(k,l)=a(k,l)
       else if (row1.and.col2) then
        m(k,l)=a(k,l+1)
       else if (row2.and.col1) then
        m(k,l)=a(k+1,l)
       else if (row2.and.col2) then
        m(k,l)=a(k+1,l+1)
       endif
   11 continue
   12 continue
       sign=(-1.0)**(i+j)
       b(i,j)=sign*Det(m,3)
   13 continue
   14 continue
      d=Det(a,4)
      do 22 i=1,4
      do 21 j=1,4
       z(i,j)=b(j,i)/d
   21 continue
   22 continue
      return
      END



      SUBROUTINE MatPro(a1,a2,a)
      REAL a1(4,4),a2(4,4),a(4,4)
c     GETS a1,a2
c     GIVES a
c     Calculates the matrix product a=a1*a2.
c     :::
      INTEGER i,j,k
      do 13 i=1,4
      do 12 j=1,4
       a(i,j)=0.0
      do 11 k=1,4
       a(i,j)=a(i,j)+a1(i,k)*a2(k,j)
   11 continue
   12 continue
   13 continue
      return
      END



      SUBROUTINE MatSum(c1,c2,c3,a1,a2,a3,a)
      REAL c1,c2,c3,a1(4,4),a2(4,4),a3(4,4),a(4,4)
c     GETS c1,c2,c3,a1,a2,a3
c     GIVES a
c     :::
      INTEGER i,j
      do 12 i=1,4
      do 11 j=1,4
       a(i,j)=c1*a1(i,j)+c2*a2(i,j)+c3*a3(i,j)
   11 continue
   12 continue
      return
      END



      SUBROUTINE Offset(a1,a2,a3,trial,lam,gam,the)
      INTEGER a1,a2,a3,trial
      REAL lam,gam,the
c     GETS a1,a2,a3,trial
c     GIVES trial,lam,gam,the
c     :::
      REAL RAD,OFF
      PARAMETER (RAD=3.14159265/180,OFF=1e-3*RAD)
      trial=trial+1
      write(*,'(1x,''trial='',i1)') trial
      if (trial.eq.1) then
       lam=a1*RAD+OFF
       gam=a2*RAD+OFF
       the=a3*RAD+OFF
      else
     :if (trial.eq.2) then
       lam=a1*RAD
       gam=a2*RAD+OFF
       the=a3*RAD+OFF
      else
     :if (trial.eq.3) then
       lam=a1*RAD+OFF
       gam=a2*RAD
       the=a3*RAD+OFF
      else
     :if (trial.eq.4) then
       lam=a1*RAD+OFF
       gam=a2*RAD+OFF
       the=a3*RAD
      else
     :if (trial.eq.5) then
       lam=a1*RAD
       gam=a2*RAD
       the=a3*RAD+OFF
      else
     :if (trial.eq.6) then
       lam=a1*RAD
       gam=a2*RAD+OFF
       the=a3*RAD
      else
     :if (trial.eq.7) then
       lam=a1*RAD+OFF
       gam=a2*RAD
       the=a3*RAD
      endif
      return
      END



      SUBROUTINE Renum(q,f,qn,fn,fail)
      INTEGER fail
      COMPLEX q(8),f(8,8),qn(8),fn(8,8)
c     GETS q,f
c     GIVES qn,fn,fail
c     :::
      INTEGER i,j,k,l
      REAL pr(8),qi(8)
      LOGICAL pupp,qupp,plow,qlow
      fail=0
c     ...
      do 10 j=1,8
       pr(j)=-0.5*real(
     :    conjg(f(1,j))*f(5,j)+
     :    conjg(f(2,j))*f(6,j)+
     :    conjg(f(3,j))*f(7,j)+
     :    f(4,j)*conjg(f(8,j)))
       qi(j)=imag(q(j))
   10 continue
c     Note: (1/r)*pr(j) is the real part of the z component of the
c     Poynting vector for the j-th eigenwave at the plane z=0.
c     ...
      k=1
      l=5
      do 29 j=1,8
       pupp=(qi(j).eq.0.0).and.(pr(j).gt.0.0)
       qupp=qi(j).lt.0.0
       plow=(qi(j).eq.0.0).and.(pr(j).lt.0.0)
       qlow=qi(j).gt.0.0
       if (pupp.or.qupp) then
        qn(k)=q(j)
        do 21 i=1,8
         fn(i,k)=f(i,j)
   21   continue
        k=k+1
       else
     : if (plow.or.qlow) then
        qn(l)=q(j)
        do 22 i=1,8
         fn(i,l)=f(i,j)
   22   continue
        l=l+1
       else
        goto 999
       endif
   29 continue
      return
c     Set fail: Multiple eigenvalues occured.
  999 fail=-2
      return
      END



      SUBROUTINE RetCE(cip,ei,c,e)
      REAL cip(6,6),ei(3,6),c(3,3,3,3),e(3,3,3)
c     GETS cip,ei
c     GIVES c,e
c     Calculates c,e from cip,ei by applying the convention of double
c     indices: (1,1)<->1, (2,2)<->2, (3,3)<->3, (2,3)=(3,2)<->4,
c     (1,3)=(3,1)<->5, (1,2)=(2,1)<->6.
c     :::
      INTEGER i,j,k,l,y(3,3)
      y(1,1)=1
      y(2,2)=2
      y(3,3)=3
      y(2,3)=4
      y(3,2)=4
      y(1,3)=5
      y(3,1)=5
      y(1,2)=6
      y(2,1)=6
      do 14 i=1,3
      do 13 j=1,3
      do 12 k=1,3
       e(i,j,k)=ei(i,y(j,k))
      do 11 l=1,3
       c(i,j,k,l)=cip(y(i,j),y(k,l))
   11 continue
   12 continue
   13 continue
   14 continue
      return
      END



      SUBROUTINE RootZ(rho,cp,ep,epsp,rq,rs,re,rz,epse,dinv,fail)
      INTEGER fail
      REAL rho,cp(3,3,3,3),ep(3,3,3),epsp(3,3),rq,rs,re,rz,epse,dinv
c     GETS rho,cp,ep,epsp
c     GIVES rq,rs,re,rz,epse,dinv,fail
c     USES Seek,Look
c     :::
      REAL rlow,rupp,rmid,TestQ,TestE,TestZ,zinf,dzet,EPS0,
     :   MEPS,EPS,DELTA,RANGE
      PARAMETER (EPS0=8.85434e-3,
     :   MEPS=0.596e-7,EPS=15*MEPS,DELTA=0.01,RANGE=0.1)
      EXTERNAL TestQ,TestE,TestZ
      rlow=0.01
      rupp=1.0
      call Seek(rho,cp,ep,epsp,rlow,rupp,rmid,TestQ,zinf,dzet,fail)
      if (fail.ne.0) return
      rq=rmid
c     ...
      rlow=rq+EPS
      rupp=rq+DELTA
      call Seek(rho,cp,ep,epsp,rlow,rupp,rmid,TestZ,zinf,dzet,fail)
      rs=rmid
c     ...
      rlow=rq+DELTA
      rupp=rq+RANGE
      call Look(rho,cp,ep,epsp,rlow,rupp,rmid,TestZ,zinf,dzet,fail)
      if (fail.ne.0) return
      rz=rmid
      epse=1/(zinf*EPS0)
      dinv=zinf/(dzet*rz)
c     ...
      rlow=rz-DELTA
      rupp=rz
      call Look(rho,cp,ep,epsp,rlow,rupp,rmid,TestE,zinf,dzet,fail)
      fail=0
      re=rmid
c     ...
      return
      END



      SUBROUTINE Rotate(c,e,eps,lam,gam,the,cp,ep,epsp)
      REAL c(3,3,3,3),e(3,3,3),eps(3,3),lam,gam,the,
     :   cp(3,3,3,3),ep(3,3,3),epsp(3,3)
c     GETS c,e,eps,lam,gam,the
c     GIVES cp,ep,epsp
c     USES RotTn4,RotTn3,RotTn2
c     Transforms cp,ep,epsp into c,e,eps.
c     :::
      call RotTn4(c,lam,gam,the,cp)
      call RotTn3(e,lam,gam,the,ep)
      call RotTn2(eps,lam,gam,the,epsp)
      return
      END



      SUBROUTINE RotTn4(c,lam,gam,the,cp)
      REAL c(3,3,3,3),lam,gam,the,cp(3,3,3,3)
c     GETS c,lam,gam,the
c     GIVES cp
c     USES DefA
c     Transforms cp into c.
c     :::
      INTEGER i,j,k,l,ip,jp,kp,lp
      REAL a(3,3)
      call DefA(lam,gam,the,a)
      do 24 ip=1,3
      do 23 jp=1,3
      do 22 kp=1,3
      do 21 lp=1,3
       cp(ip,jp,kp,lp)=0.0
      do 14 i=1,3
      do 13 j=1,3
      do 12 k=1,3
      do 11 l=1,3
       cp(ip,jp,kp,lp)=cp(ip,jp,kp,lp)+
     :    a(ip,i)*a(jp,j)*a(kp,k)*a(lp,l)*c(i,j,k,l)
   11 continue
   12 continue
   13 continue
   14 continue
   21 continue
   22 continue
   23 continue
   24 continue
      return
      END



      SUBROUTINE RotTn3(e,lam,gam,the,ep)
      REAL e(3,3,3),lam,gam,the,ep(3,3,3)
c     GETS e,lam,gam,the
c     GIVES ep
c     USES DefA
c     Transforms ep into e.
c     :::
      INTEGER i,j,k,ip,jp,kp
      REAL a(3,3)
      call DefA(lam,gam,the,a)
      do 23 ip=1,3
      do 22 jp=1,3
      do 21 kp=1,3
       ep(ip,jp,kp)=0.0
      do 13 i=1,3
      do 12 j=1,3
      do 11 k=1,3
       ep(ip,jp,kp)=ep(ip,jp,kp)+
     :    a(ip,i)*a(jp,j)*a(kp,k)*e(i,j,k)
   11 continue
   12 continue
   13 continue
   21 continue
   22 continue
   23 continue
      return
      END



      SUBROUTINE RotTn2(eps,lam,gam,the,epsp)
      REAL eps(3,3),lam,gam,the,epsp(3,3)
c     GETS eps,lam,gam,the
c     GIVES epsp
c     USES DefA
c     Transforms epsp into eps.
c     :::
      INTEGER i,j,ip,jp
      REAL a(3,3)
      call DefA(lam,gam,the,a)
      do 22 ip=1,3
      do 21 jp=1,3
       epsp(ip,jp)=0.0
      do 12 i=1,3
      do 11 j=1,3
       epsp(ip,jp)=epsp(ip,jp)+a(ip,i)*a(jp,j)*eps(i,j)
   11 continue
   12 continue
   21 continue
   22 continue
      return
      END



      SUBROUTINE Seek(rho,cp,ep,epsp,rlow,rupp,rmid,Test,zinf,dzet,fail)
      INTEGER fail
      REAL rho,cp(3,3,3,3),ep(3,3,3),epsp(3,3),
     :   rlow,rupp,rmid,Test,zinf,dzet
      EXTERNAL Test
c     GETS rho,cp,ep,epsp,rlow,rupp,Test
c     GIVES rmid,zinf,dzet,fail
c     USES Test
c     :::
      INTEGER iter
      REAL rlast,rinf,zlow,zmid,zupp,MEPS,EPS
      PARAMETER (MEPS=0.596e-7,EPS=10*MEPS)
      rlast=rupp
c     ...
      zlow=Test(rlow,rho,cp,ep,epsp,fail)
      if (fail.ne.0) then
       rmid=0.99999
       return
      endif
      if (zlow.gt.0.0) then
       rmid=0.0
       goto 999
      endif
c     ...
      iter=0
   11 continue
       if (iter.gt.40) goto 9999
       rmid=(rlow+rupp)/2.0
       if (abs(rmid-rlast).le.EPS) goto 12
       zmid=Test(rmid,rho,cp,ep,epsp,fail)
       if (fail.ne.0) then
        rmid=9.99999
        return
       endif
       if (zmid.gt.0.0) then
        rupp=rmid
        zupp=zmid
       else
        rlow=rmid
        zlow=zmid
       endif
       rlast=rmid
       iter=iter+1
      goto 11
   12 continue
c     ...
      rinf=10.0
   21 continue
       if (rinf.gt.100.0) return
       zinf=Test(rinf,rho,cp,ep,epsp,fail)
       if (fail.eq.0) goto 22
       rinf=rinf*1.5
      goto 21
   22 continue
      dzet=(zupp-zlow)/(rupp-rlow)
      return
c     Set fail: Function Test has no zero in the given range.
  999 fail=-1
      return 
c     Set fail: More than 40 iterations attempted.
 9999 fail=500
      return
      END



      SUBROUTINE SymMat(a,s,n)
      INTEGER n
      REAL a(n,n),s(n,n)
c     GETS a,n
c     GIVES s
c     Calculates the symmetric matrix s such that its lower part is
c     equal to the lower part of the matrix a (i.e. elements on the
c     diagonal and below of the two matrices are equal).
c     :::
      INTEGER i,j
      do 12 i=1,n
      do 11 j=1,n
       if (i.ge.j) then
        s(i,j)=a(i,j)
        s(j,i)=a(i,j)
       endif
   11 continue
   12 continue
      return
      END



      FUNCTION TestE(r,rho,cp,ep,epsp,fail)
      INTEGER fail
      REAL TestE,r,rho,cp(3,3,3,3),ep(3,3,3),epsp(3,3)
c     GETS r,rho,cp,ep,epsp
c     GIVES TestE,fail
c     USES CalcH,FindVV,DefZ
c     :::
      REAL h(8,8),EPS0
      COMPLEX q(8),f(8,8),zet
      PARAMETER (EPS0=8.85434e-3)
      call CalcH(r,rho,cp,ep,epsp,h)
      call FindVV(h,q,f,fail)
      if (fail.ne.0) return
      call DefZ(f,zet,fail)
      if (fail.ne.0) return
      TestE=imag(zet)+1/EPS0
      return
      END



      FUNCTION TestQ(r,rho,cp,ep,epsp,fail)
      INTEGER fail
      REAL TestQ,r,rho,cp(3,3,3,3),ep(3,3,3),epsp(3,3)
c     GETS r,rho,cp,ep,epsp
c     GIVES TestQ,fail
c     USES CalcH,FindVV
c     :::
      INTEGER j
      REAL h(8,8)
      COMPLEX q(8),f(8,8)
      TestQ=-1.0
      call CalcH(r,rho,cp,ep,epsp,h)
      call FindVV(h,q,f,fail)
      if (fail.ne.0) return
      do 10 j=1,8
       if (imag(q(j)).eq.0.0) return
   10 continue
      TestQ=1.0
      return
      END



      FUNCTION TestZ(r,rho,cp,ep,epsp,fail)
      INTEGER fail
      REAL TestZ,r,rho,cp(3,3,3,3),ep(3,3,3),epsp(3,3)
c     GETS r,rho,cp,ep,epsp
c     GIVES TestZ,fail
c     USES CalcH,FindVV,DefZ
c     :::
      REAL h(8,8)
      COMPLEX q(8),f(8,8),zet
      call CalcH(r,rho,cp,ep,epsp,h)
      call FindVV(h,q,f,fail)
      if (fail.ne.0) return
      call DefZ(f,zet,fail)
      if (fail.ne.0) return
      TestZ=imag(zet)
      return
      END



      SUBROUTINE UnPack(nm,n,wi,u,ur,ui)
      INTEGER nm,n
      REAL wi(n),u(nm,n),ur(nm,n),ui(nm,n)
c     GETS nm,n,wi,u
c     GIVES ur,ui
c     Unpacks a matrix of complex eigenvectors into matrices of the real
c     and imaginary parts of the eigenvectors.
c        Input: the array wi contains the imaginary parts of the
c     eigenvalues of the original matrix (output from HQR2).
c        Input: the array u (n must not be greater than nm) contains the
c     eigenvectors to be unpacked (output from BalBak, or HQR2 if Balanc
c     is not used) corresponding to the eigenvalues with the imaginary
c     parts stored in the array wi.
c        Output: the arrays ur,ui contain the real and imaginary parts
c     of the eigenvectors.
c     :::
      INTEGER j,k
      do 20 k=1,n
       if (wi(k).eq.0.0) then
        do 11 j=1,n
         ur(j,k)=u(j,k)
         ui(j,k)=0.0
   11   continue
       else if (wi(k).gt.0.0) then
        do 12 j=1,n
         ur(j,k)=u(j,k)
         ui(j,k)=u(j,k+1)
   12   continue
       else if (wi(k).lt.0.0) then
        do 13 j=1,n
         ur(j,k)=u(j,k-1)
         ui(j,k)=-u(j,k)
   13   continue
       endif
   20 continue
      return
      END
