c
         parameter (ntheta=360, nphi=360)
         parameter (nn=3, np=3, time=1.0)
         Dimension cc(3,3,3,3)  
         Dimension rr(3,nphi,nphi), rv(3,3,nphi,nphi)  
c
         do i=1,3
            do j=1,3
               do k=1,3
                  do l=1,3
                     cc(i,j,k,l) = 0.0
                  enddo
               enddo
            enddo
         enddo
         cc(1,1,1,1) = 5.063
         cc(2,2,2,2) = 5.063
         cc(3,3,3,3) = 4.000
         cc(1,1,2,2) = 1.938 
         cc(2,2,1,1) = 1.938
         cc(1,1,3,3) = 2.595 
         cc(3,3,1,1) = 2.595
         cc(2,2,3,3) = 2.595
         cc(3,3,2,2) = 2.595
c
         cc(2,3,2,3) = 1.103
         cc(2,3,3,2) = 1.103
         cc(3,2,3,2) = 1.103
         cc(3,2,2,3) = 1.103
c
         cc(1,3,1,3) = 1.103
         cc(1,3,3,1) = 1.103
         cc(3,1,3,1) = 1.103
         cc(3,1,1,3) = 1.103
c
         cc(1,2,1,2) = 1.562
         cc(1,2,2,1) = 1.562
         cc(2,1,2,1) = 1.562
         cc(2,1,1,2) = 1.562
c
         call wavefront(cc,time,rr,rv)
c
         stop
         end
c
c
c        "eigenv" is a subroutine from Press et al. (1986, pages 346 to 348). 
c        eigenvalues and eigenvectors of the matrix "g", which are located in
c        "dd" and "vv", respectively.
c
         subroutine wavefront(aa,time,rr,rv)
         parameter (nphi=360,nn=3,np=3)
         dimension g(3,3), dd(3), cn(3), vv(3,3), aa(3,3,3,3)  
         dimension rr(3,nphi,nphi), rv(3,3,nphi,nphi)
c
         pi2 = 8.*atan(1.)
         dpi = pi2/360.
         do ithe=1,nphi
            theta = (ithe-1)*dpi
            xxxs = sin(theta)
            cn(3) = cos(theta)
            do kphi=1,nphi
               phi   = (kphi-1)*dpi
               cn(1) = xxxs*cos(phi)
               cn(2) = xxxs*sin(phi)
               do i=1,3
                  do k=1,3
                     g(i,k) = 0.0
                     do j=1,3
                        do l=1,3
                           g(i,k) = g(i,k) + (aa(i,j,k,l)*cn(j)*cn(l))
                        enddo
                     enddo
                  enddo
               enddo
c
               call eigenv(g,nn,np,dd,vv,nrot)
c
               do m=1,3
                  rr(m,ithe,kphi) = sqrt(dd(m))*time
                  do n=1,3
                     rv(m,n,ithe,kphi) = vv(m,n) 
                  enddo 
               enddo 
c
            enddo 
         enddo 
c
         return
         end
c
c =========================================================================
c
         subroutine eigenv(aa,nn,np,dd,vv,nrot)
c
         parameter (nmax=100)
c
         dimension aa(np,np), dd(np), vv(np,np), bb(nmax), zz(nmax)
c
         do ip=1,nn
            do iq=1,nn
               vv(ip,iq) = 0.0
            enddo
            vv(ip,ip) = 1.0
         enddo
c
         do ip=1,nn
            bb(ip) = aa(ip,ip)
            dd(ip) = bb(ip)
            zz(ip) = 0.0
         enddo
c
         nrot = 0
         do i=1,50
            sm = 0.0
            do ip=1,nn-1
               do iq=ip+1,nn
                  sm = sm + abs(aa(ip,iq))
               enddo
            enddo
            if(sm.eq.0.) return
            if(i.lt.4) then
                 tresh = (0.2*sm)/(nn**2)
               else
                 tresh = 0.0
            endif 
            do ip=1,nn-1
               do iq=ip+1,nn
                  g = 100.*abs(aa(ip,iq))
                  if((i.gt.4).and.(abs(dd(ip))+g.eq.abs(dd(ip))).and.
     *               (abs(dd(iq))+g.eq.abs(dd(iq)))) then 
                        aa(ip,iq) = 0.0
                     else if(abs(aa(ip,iq)).gt.tresh) then
                        h = dd(iq)-dd(ip)
                        if(abs(h)+g.eq.abs(h)) then
                           t = aa(ip,iq)/h
                        else
                           theta = (0.5*h)/aa(ip,iq)
                           t = 1.0/(abs(theta)+sqrt(1.0+theta**2))
                           if(theta.lt.0.) t= -t
                        endif
                        c = 1./sqrt(1.0+t**2)
                        s = t*c
                        tau = s/(1.0+c)
                        h   = t*aa(ip,iq)
                        zz(ip) = zz(ip)-h
                        zz(iq) = zz(iq)+h
                        dd(ip) = dd(ip)-h
                        dd(iq) = dd(iq)+h
                        aa(ip,iq) = 0.0
                        do j=1,ip-1
                           g = aa(j,ip)
                           h = aa(j,iq)
                           aa(j,ip) = g - s*(h+g*tau)
                           aa(j,iq) = h + s*(g-h*tau)
                        enddo
                        do j=ip+1,iq-1
                           g = aa(ip,j)
                           h = aa(j,iq)
                           aa(ip,j) = g - s*(h+g*tau)
                           aa(j,iq) = h + s*(g-h*tau)
                        enddo
                        do j=iq+1,nn
                           g = aa(ip,j)
                           h = aa(iq,j)
                           aa(ip,j) = g - s*(h+g*tau)
                           aa(iq,j) = h + s*(g-h*tau)
                        enddo
                        do j=1,nn
                           g = vv(j,ip)
                           h = vv(j,iq)
                           vv(j,ip) = g - s*(h+g*tau)
                           vv(j,iq) = h + s*(g-h*tau)
                        enddo
                        nrot = nrot +1 
                     endif
                   enddo
                 enddo
                 do ip=1,nn
                    bb(ip) = bb(ip)+zz(ip)
                    dd(ip) = bb(ip)
                    zz(ip) = 0.0
                 enddo
               enddo
               Pause '50 iterations should never happen'
               return
               end