***************************************************************
      subroutine polar(ndim,npts,dt,kwin,fw,fmin,fmax,nout)
***************************************************************
c  input
c   data in a(1...npts,k) k=1..3  where a is dimensioned a(nndim,3)
c...components vert,rad,trans corresponding to 1,2,3
c   dt is sample interval
c   kwin is desired # of windows (.le.8)
c   fw is time-bandwidth product of windows (usually 4)
c   fmin and fmax is freq range of calculation (in Hz)
c  output
c   all results in common/result/ans(1.....nout,1-20)
c   ans(..,1)=frequency (in Hz)
c   ans(..,2)=amplitude of component 1
c   ans(..,3)=amplitude of component 2
c   ans(..,4)=amplitude of component 3
c   ans(..,5)=smallest sing val
c   ans(..,6)=middle sing val
c   ans(..,7)=largest sing val
c   ans(..,8)=horiz azimuth
c   ans(..,9)=  and its error
c   ans(..,10)=horizontal phase
c   ans(..,11)= and its error
c   ans(..,12)=vertical azimuth
c   ans(..,13)= and its error
c   ans(..,14)=vertical phase
c   ans(..,15)= and its error
c...
c...*additional since july 1991, gabi
c...subroutine ang3d for calculation of position of major axis of ellipse
c...in 3-d space
c...ans(..,16)= tilt
c...ans(..,17)= horizontal projection of major axis from eigenvector
c...ans(..,18)= vertical angle of major axis        from eigenvector
c...ans(..,19)= ellipticity (1-b/a)
c...ans(..,20)= after ang3d:
c...            Paulssen et al. noise to signal
c...ans(..,20)= after eptis: 
c...            used by polart for wrapping ans(17)
c...*subroutine adwait changed so that calculations with one taper possible
c...*svd and calculation of angles extracted to subroutine decomp
c...
c  see comments in subroutine ang for definitions of angles
c  Note that common/result/ is used as work space and so must be dimensioned
c  2*kw*ipow points long where ipow is currently 8192 and kw is 10
c  nh=ipow/2
      parameter (nndim=5000,ipow=4096,icol=20,kw=10,nh=2048)
      real*8 tas(ipow,kw),el(ipow),w2(ipow),dfw,b(3,3,2),y(3)
      common/taperz/ta(nh,kw,3),ti(nh,kw,3)
      common/result/ans(ipow,icol)
      common/sdata/a(nndim,3)
      dimension w1(1),avar(3),dff(kw),dk(kw,3)
      equivalence(w1,w2),(tas,ans)
      data pi/3.14159265358979d0/
c*** compute fft parameters, number of tapers - currently no padding
      dfw=fw
      nwin=min0(kw,max0(kwin,1))
      nph=npts/2
      npad=min0(nph,ipow)
      call facdwn(npad)
      npad=2*npad
      df=1./(npad*dt)
      nf=npad/2+1.5
      nf2=nf+nf
c*** compute windows
      call staper(npts,dfw,nwin,tas,ipow,el,w2)
c*** loop over 3 components, window then fft and put in ta,ti
c*** avar is the power in w domain (with half contributions from
c*** zero  and nyquist frequencies (using fftl normalisation)
      np1=npts+1
      do 30 jser=1,3
      do 20 iwin=1,nwin
      do 5 i=1,npts
    5 w1(i)=a(i,jser)*tas(i,iwin)
      do 10 i=np1,npad
   10 w1(i)=0.
c...  fourier forward transform with exp(iwt)!
      call fftl(w1,npad,1,ierr)
      if(ierr.ne.0) print *,'problem in fft'
      j=0
      do 15 i=1,nf2,2
      j=j+1
      ta(j,iwin,jser)=w1(i)
   15 ti(j,iwin,jser)=w1(i+1)
   20 continue
      avar(jser)=0.d0
      do 25 i=1,npts
   25 avar(jser)=avar(jser)+a(i,jser)*a(i,jser)
   30 avar(jser)=2.*avar(jser)/npad
      nout=0

c*** perform polarisation analysis - loop over frequency

      do 100 l=1,nf
      frq=(l-1)*df
      if(frq.lt.fmin) goto 100
      if(frq.gt.fmax) return
c*** compute adaptive weights and effective degrees of freedom
      dfwin=0.d0
      do 35 j=1,3
   35 call adwait(l,ta(1,1,j),ti(1,1,j),dk(1,j),el,nwin,nf,avar(j),y(j)
     +                 ,nh)
      do 40 i=1,nwin
      dff(i)=amin1(dk(i,1),dk(i,2),dk(i,3))
      dff(i)=dff(i)*dff(i)
   40 dfwin=dfwin+dff(i)
      if(dfwin.eq.0.d0) goto 100
c*** fill spectral density matrix
      do 50 i=1,3
      do 50 j=i,3
      b(i,j,1)=0.d0
      b(i,j,2)=0.d0
c...  indices i and j are twisted for imaginary part (rows and columns twisted)?

      do 45 k=1,nwin
      b(i,j,1)=b(i,j,1)+(ta(l,k,i)*ta(l,k,j)+ti(l,k,i)*ti(l,k,j))*dff(k)
   45 b(i,j,2)=b(i,j,2)+(ta(l,k,i)*ti(l,k,j)-ti(l,k,i)*ta(l,k,j))*dff(k)
      b(j,i,1)=b(i,j,1)
   50 b(j,i,2)=-b(i,j,2)
      nout=nout+1
      call decomp(b,y,dfwin,nout,frq)
  100 continue
      return
      end

**********************************************************
      subroutine adwait(j,ta,ti,dk,el,nwin,nf,avar,amp,nh)
**********************************************************
c***  adaptive weights computed at frequency index j,avar should be total
c***  variance of process,amp is estimate of amplitude
      implicit real*8(a-h,o-z)
      parameter (kw=10)
      real*4 ta,ti,dk,avar
      dimension ta(nh,1),ti(nh,1),dk(1),spw(kw),el(1)
      data tol/3.d-4/
c...  clearing of array added, so that only 1 taper can be used, gabi
      do 30 i=1,kw
 30   spw(i)=0.
      do 5 i=1,nwin
    5 spw(i)=(ta(j,i)*ta(j,i)+ti(j,i)*ti(j,i))/avar
      as=0.5d0*(spw(1)+spw(2))
      ax=0.d0
      if(as.eq.0.d0) goto 20
      do 15 k=1,20
      fn=0.d0
      fx=0.d0
      do 10 i=1,nwin
      a1=dsqrt(el(i))*as/(el(i)*as+1.d0-el(i))
      a1=a1*a1
      fn=fn+a1*spw(i)
   10 fx=fx+a1
      ax=fn/fx
      if(dabs(ax-as)/as.lt.tol) go to 20
   15 as=ax
   20 sum1=0.d0
      sum2=0.d0
      do 25 i=1,nwin
      dk(i)=dsqrt(el(i))*ax/(el(i)*ax+1.d0-el(i))
      sum2=sum2+dk(i)*dk(i)
   25 sum1=sum1+dk(i)*dk(i)*spw(i)
      amp=0.d0
      if(sum2.eq.0.d0) return
      amp=dsqrt(avar*sum1/sum2)
      return
      end

****************************************
      subroutine decomp(b,y,dfwin,nout,frq)
****************************************
      parameter (ipow=4096,icol=20)
      real*8 b(3,3,2),z(3,3,2),y(3),
     + w(3),bcd(6),def(6),sig,haz,dhaz,hphs,dhphs,vaz,dvaz,vphs,dvphs
      common/result/ans(ipow,icol)
c*** use eispack routines to decompose b
c...  produce real symmetric tridiagonal matrix

      call htridi(3,3,b,b(1,1,2),w,bcd,bcd,def)
      do 60 i=1,3
      do 55 j=1,3
      z(j,i,2)=0.d0
   55 z(j,i,1)=0.d0
   60 z(i,i,1)=1.d0
c...  find eigenvalues and eigenvecotrs of symmetric tridiagonal matrix
      call tql2(3,3,w,bcd,z,ier)
      if(ier.ne.0) print *,'problem in tql2'
c...  backtransforming for hermitian matrix
c...  output vectors are stored in columns  of z
      call htribk(3,3,b,b(1,1,2),def,3,z,z(1,1,2))
c*** analyse eigenvectors for polarisation angles
      sig=(w(1)+w(2))*.75d0/(w(3)*dfwin)
      call ang(z,sig,haz,dhaz,hphs,dhphs,vaz,dvaz,vphs,dvphs)
      sum=0.
      do 65 i=1,3
      if(dabs(w(i)).eq.0.d0) goto 65
      w(i)=dsqrt(dabs(w(i)))
   65 sum=sum+w(i)

      ans(nout,1)=frq
      do 70 i=1,3
      ans(nout,i+1)=y(i)
   70 ans(nout,i+4)=w(i)/sum
      ans(nout,8)=haz
      ans(nout,9)=dhaz
      ans(nout,10)=hphs
      ans(nout,11)=dhphs
      ans(nout,12)=vaz
      ans(nout,13)=dvaz
      ans(nout,14)=vphs
      ans(nout,15)=dvphs
      call ang3d(z,nout)

      return
      end

*******************************************************************
      subroutine ang(z,sig,haz,dhaz,hphs,dhphs,vaz,dvaz,vphs,dvphs)
*******************************************************************
c  the third eigenvector is the principal polarization vector, it defines a plan
c  ar
c  ellipse in 3-space around which the particle motion rotates. in a perfect sei
c  smic
c  pulse, the ellipse will collapse to a line and motion with respect to each
c  coordinate will have relative phase lags of either zero or 180 degrees. we ca
c  n
c  divine such quantities as horizontal azimuth and angle of incidence. if there
c   is
c  an oddball phase lag twixt components we have elliptical particle motion, and
c   such
c  quantities become harder to quantify with certainty. the presence of noise, a
c  s
c  measured by the 1st and 2nd singular values, also lends uncertainty.
c   quantities computed are:
c    haz   = apparent horizontal azimuth
c    hphs  = phase lag twixt horizontal components
c    vaz   = apparent angle of incidence of vertical with the horiz
c    vphs  = phase lag twixt vertical and horizontal components
c  horizontal azimuth is measured with respect to the second component (assuming
c  a
c  clockwise rotation to the third component). this implies a right-handed coord
c  inate
c  system, which is not the typical seismic configuration. the quantity calculat
c  ed
c  here corresponds to an angle measured clockwise from the second component in
c  a
c  left-handed system. to let haz range over 360 degrees, we must specify the
c  orientation of the vertical motion w.r.t. horizontal. we use the sign of
c  re(z1*conjg(z2)) or re(z1*conjg(z3)), whichever is larger,to specify the prop
c  er
c  quadrant for haz,leaving vaz to vary between 0 and 90
c based on code by jeff park
      implicit real*8(a-h,o-z)
      dimension z(3,3,2),phi(3),zsz(3),dz(3,3,2),zsqr(3),zsqi(3)
      data pi/3.14159265358979d0/
      con=180.d0/pi
      p2=0.5d0*pi
      do 5 i=1,3
      zsz(i)=z(i,3,1)**2+z(i,3,2)**2
      phi(i)=datan2(z(i,3,2),z(i,3,1))
      zsqr(i)=z(i,3,1)**2-z(i,3,2)**2
    5 zsqi(i)=2.d0*z(i,3,1)*z(i,3,2)
c*** fill covariance matrix
      do 10 i=1,3
      do 10 j=1,3
      dz(j,i,1)=sig*(z(j,1,1)*z(i,1,1)+z(j,1,2)*z(i,1,2)+
     x               z(j,2,1)*z(i,2,1)+z(j,2,2)*z(i,2,2))
   10 dz(j,i,2)=sig*(z(j,1,2)*z(i,1,1)-z(j,1,1)*z(i,1,2)+
     x               z(j,2,2)*z(i,2,1)-z(j,2,1)*z(i,2,2))
c*** horiz azimuth
      zhr=zsqr(2)+zsqr(3)
      zhi=zsqi(2)+zsqi(3)
      zhsq=zhr*zhr+zhi*zhi
      zh=dsqrt(zhsq)
      thd=.5*datan2(zhi,zhr)
      r1=zsz(2)*dcos(thd-phi(2))**2+zsz(3)*dcos(thd-phi(3))**2
      if((zsz(2)+zsz(3)).gt.2.d0*r1) thd=thd+p2
      cs=dcos(thd)
      sn=dsin(thd)
      r1=z(3,3,1)*cs+z(3,3,2)*sn
      r2=z(2,3,1)*cs+z(2,3,2)*sn
      haz=con*datan2(r1,r2)
      if(haz.lt.0.d0) haz=haz+180.d0
      z1z2=z(1,3,1)*z(2,3,1)+z(1,3,2)*z(2,3,2)
      z1z3=z(1,3,1)*z(3,3,1)+z(1,3,2)*z(3,3,2)
      if(dabs(z1z2).gt.dabs(z1z3)) goto 15
      if(z1z3.ge.0.d0) haz=haz-180.d0
      go to 20
   15 if(z1z2.lt.0.d0.and.haz.gt.90.d0) haz=haz-180.d0
      if(z1z2.gt.0.d0.and.haz.le.90.d0) haz=haz-180.d0
   20 arg=dz(2,2,1)*zsz(3)+dz(3,3,1)*zsz(2)-2.d0*
     x  (dz(2,3,1)*(z(2,3,1)*z(3,3,1)+z(2,3,2)*z(3,3,2))
     x  +dz(2,3,2)*(z(2,3,2)*z(3,3,1)-z(2,3,1)*z(3,3,2)))
      dhaz=dmin1(con*dsqrt(dmax1(arg,0.d0)/(2.d0*zhsq)),360.d0)
c*** phase lag between horizontal components
      hphs=con*(phi(3)-phi(2))
      if(hphs.gt. 180.d0) hphs=hphs-360.d0
      if(hphs.lt.-180.d0) hphs=hphs+360.d0
      dhphs=dmin1(dhaz*dsqrt(zhsq/(zsz(2)*zsz(3))),360.d0)
c*** angle of incidence of vert with the horiz plane
      ztr=zhr+zsqr(1)
      zti=zhi+zsqi(1)
      ztsq=ztr*ztr+zti*zti
      thdv=.5*datan2(zti,ztr)
      r1=zsz(2)*dcos(thdv-phi(2))**2+zsz(3)*dcos(thdv-phi(3))**2
     x     +zsz(1)*dcos(thdv-phi(1))**2
      if(r1.lt.0.5d0) thdv=thdv+p2
      cs=dcos(thdv)
      sn=dsin(thdv)
      r3=z(3,3,1)*cs+z(3,3,2)*sn
      r2=z(2,3,1)*cs+z(2,3,2)*sn
      r1=z(1,3,1)*cs+z(1,3,2)*sn
      vaz=con*datan2(dsqrt(r2*r2+r3*r3),dabs(r1))
      arg=zhsq*dz(1,1,1)+zsz(1)*zsz(2)*dz(2,2,1)+zsz(1)*zsz(3)*dz(3,3,1)
     x         -2.d0*((zhr*z(1,3,1)+zhi*z(1,3,2))*(z(2,3,1)*dz(1,2,1)
     x      +z(2,3,2)*dz(1,2,2)+z(3,3,1)*dz(1,3,1)+z(3,3,2)*dz(1,3,2))
     x               -(zhi*z(1,3,1)-zhr*z(1,3,2))*(z(2,3,1)*dz(1,2,2)
     x      -z(2,3,2)*dz(1,2,1)+z(3,3,1)*dz(1,3,2)-z(3,3,2)*dz(1,3,1)))
     x  +2.d0*zsz(1)*(z(2,3,1)*(z(3,3,1)*dz(2,3,1)+z(3,3,2)*dz(2,3,2))
     x               -z(2,3,2)*(z(3,3,1)*dz(2,3,2)-z(3,3,2)*dz(2,3,1)))
      dvaz=dmin1(con*dsqrt(dmax1(arg,0.d0)/(2.d0*ztsq*zh)),360.d0)
c*** phase angle between vertical and horizontal components
      vphs=con*(thd-phi(1))
      if(vphs.gt. 180.d0) vphs=vphs-360.d0
      if(vphs.lt.-180.d0) vphs=vphs+360.d0
      if(vphs.lt. -90.d0) vphs=vphs+180.d0
      if(vphs.gt.  90.d0) vphs=vphs-180.d0
      dvphs=dmin1(dvaz*dsqrt(ztsq/(zsz(1)*zh)),360.d0)
      return
      end

**********************************************
      subroutine htridi(nm,n,ar,ai,d,e,e2,tau)
**********************************************
c  see tred1, num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
c  handbook for auto. comp., vol.ii-linear algebra, 212-226(1971).
c  reduces a complex hermitian matrix to a real symmetric tridiagonal matrix
c  using  unitary similarity transformations.
c     on input
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.
c        n is the order of the matrix.
c        ar and ai contain the real and imaginary parts,
c          respectively, of the complex hermitian input matrix.
c          only the lower triangle of the matrix need be supplied.
c     on output
c        ar and ai contain information about the unitary trans-
c          formations used in the reduction in their full lower
c          triangles.  their strict upper triangles and the
c          diagonal of ar are unaltered.
c        d contains the diagonal elements of the the tridiagonal matrix.
c        e contains the subdiagonal elements of the tridiagonal
c          matrix in its last n-1 positions.  e(1) is set to zero.
c        e2 contains the squares of the corresponding elements of e.
c        e2 may coincide with e if the squares are not needed.
c        tau contains further information about the transformations.
      implicit real*8(a-h,o-z)
      dimension ar(nm,n),ai(nm,n),d(n),e(n),e2(n),tau(2,n)
      tau(1,n) = 1.0d0
      tau(2,n) = 0.0d0
      do 100 i = 1, n
  100 d(i) = ar(i,i)
      do 300 ii = 1, n
         i = n + 1 - ii
         l = i - 1
         h = 0.0d0
         scale = 0.0d0
         if (l .lt. 1) go to 130
         do 120 k = 1, l
  120    scale = scale + dabs(ar(i,k)) + dabs(ai(i,k))
         if (scale .ne. 0.0d0) go to 140
         tau(1,l) = 1.0d0
         tau(2,l) = 0.0d0
  130    e(i) = 0.0d0
         e2(i) = 0.0d0
         go to 290
  140    do 150 k = 1, l
            ar(i,k) = ar(i,k) / scale
            ai(i,k) = ai(i,k) / scale
  150       h = h + ar(i,k) * ar(i,k) + ai(i,k) * ai(i,k)
         e2(i) = scale * scale * h
         g = dsqrt(h)
         e(i) = scale * g
         f = pythag(ar(i,l),ai(i,l))
c     .......... form next diagonal element of matrix t ..........
         if (f .eq. 0.0d0) go to 160
         tau(1,l) = (ai(i,l) * tau(2,i) - ar(i,l) * tau(1,i)) / f
         si = (ar(i,l) * tau(2,i) + ai(i,l) * tau(1,i)) / f
         h = h + f * g
         g = 1.0d0 + g / f
         ar(i,l) = g * ar(i,l)
         ai(i,l) = g * ai(i,l)
         if (l .eq. 1) go to 270
         go to 170
  160    tau(1,l) = -tau(1,i)
         si = tau(2,i)
         ar(i,l) = g
  170    f = 0.0d0
         do 240 j = 1, l
            g = 0.0d0
            gi= 0.0d0
c     .......... form element of a*u ..........
            do 180 k = 1, j
               g = g + ar(j,k) * ar(i,k) + ai(j,k) * ai(i,k)
  180          gi= gi- ar(j,k) * ai(i,k) + ai(j,k) * ar(i,k)
            jp1 = j + 1
            if (l .lt. jp1) go to 220
            do 200 k = jp1, l
               g = g + ar(k,j) * ar(i,k) - ai(k,j) * ai(i,k)
  200          gi= gi- ar(k,j) * ai(i,k) - ai(k,j) * ar(i,k)
c     .......... form element of p ..........
  220       e(j) = g / h
            tau(2,j) = gi / h
  240       f = f + e(j) * ar(i,j) - tau(2,j) * ai(i,j)
         hh = f / (h + h)
c     .......... form reduced a ..........
         do 260 j = 1, l
            f = ar(i,j)
            g = e(j) - hh * f
            e(j) = g
            fi = -ai(i,j)
            gi = tau(2,j) - hh * fi
            tau(2,j) = -gi
            do 260 k = 1, j
               ar(j,k) = ar(j,k) - f * e(k) - g * ar(i,k)
     x                           + fi * tau(2,k) + gi * ai(i,k)
  260          ai(j,k) = ai(j,k) - f * tau(2,k) - g * ai(i,k)
     x                           - fi * e(k) - gi * ar(i,k)
  270    do 280 k = 1, l
            ar(i,k) = scale * ar(i,k)
  280       ai(i,k) = scale * ai(i,k)
         tau(2,l) = -si
  290    hh = d(i)
         d(i) = ar(i,i)
         ar(i,i) = hh
  300    ai(i,i) = scale * dsqrt(h)
      return
      end

***********************************************
      subroutine htribk(nm,n,ar,ai,tau,m,zr,zi)
***********************************************
c  see trbak1, num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
c  handbook for auto. comp., vol.ii-linear algebra, 212-226(1971).
c  forms the eigenvectors of a complex hermitian matrix by back transforming  th
c  ose
c  of the corresponding real symmetric tridiagonal matrix determined by  htridi.

c     on input
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.
c        n is the order of the matrix.
c        ar and ai contain information about the unitary trans-
c          formations used in the reduction by  htridi  in their
c          full lower triangles except for the diagonal of ar.
c        tau contains further information about the transformations.
c        m is the number of eigenvectors to be back transformed.
c        zr contains the eigenvectors to be back transformed
c          in its first m columns.
c     on output
c        zr and zi contain the real and imaginary parts,
c          respectively, of the transformed eigenvectors
c          in their first m columns.
c     note that the last component of each returned vector
c     is real and that vector euclidean norms are preserved.
      implicit real*8(a-h,o-z)
      dimension ar(nm,n),ai(nm,n),tau(2,n),zr(nm,m),zi(nm,m)
      if (m .eq. 0) return
      do 50 k = 1, n
         do 50 j = 1, m
            zi(k,j) = -zr(k,j) * tau(2,k)
   50       zr(k,j) = zr(k,j) * tau(1,k)
      if (n .eq. 1) return
c     .......... recover and apply the householder matrices ..........
      do 140 i = 2, n
         l = i - 1
         h = ai(i,i)
         if (h .eq. 0.d0) go to 140
         do 130 j = 1, m
            s = 0.0d0
            si = 0.d0
            do 110 k = 1, l
               s = s + ar(i,k) * zr(k,j) - ai(i,k) * zi(k,j)
  110          si = si + ar(i,k) * zi(k,j) + ai(i,k) * zr(k,j)
            s = (s / h) / h
            si = (si / h) / h
            do 120 k = 1, l
               zr(k,j) = zr(k,j) - s * ar(i,k) - si * ai(i,k)
  120          zi(k,j) = zi(k,j) - si * ar(i,k) + s * ai(i,k)
  130    continue
  140 continue
      return
      end

**************************************
      subroutine tql2(nm,n,d,e,z,ierr)
**************************************
c  see tql2,num. math. 11, 293-306(1968) by bowdler et al
c  handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
c  finds the eigenvalues and eigenvectors of a symmetric tridiagonal matrix
c  by the ql method.
c     the eigenvectors of a full symmetric matrix can also
c     be found if  tred2  has been used to reduce this
c     full matrix to tridiagonal form.
c     on input
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.
c        n is the order of the matrix.
c        d contains the diagonal elements of the input matrix.
c        e contains the subdiagonal elements of the input matrix
c          in its last n-1 positions.  e(1) is arbitrary.
c        z contains the transformation matrix produced in the
c          reduction by  tred2, if performed.  if the eigenvectors
c          of the tridiagonal matrix are desired, z must contain
c          the identity matrix.
c      on output
c        d contains the eigenvalues in ascending order.  if an
c          error exit is made, the eigenvalues are correct but
c          unordered for indices 1,2,...,ierr-1.
c        e has been destroyed.
c        z contains orthonormal eigenvectors of the symmetric
c          tridiagonal (or full) matrix.  if an error exit is made,
c          z contains the eigenvectors associated with the stored
c          eigenvalues.
c        ierr is set to
c          zero       for normal return,
c          j          if the j-th eigenvalue has not been
c                     determined after 30 iterations.
      implicit real*8(a-h,o-z)
      dimension d(n),e(n),z(nm,n)
      ierr = 0
      if (n .eq. 1) go to 1001
      do 100 i = 2, n
  100 e(i-1) = e(i)
      f = 0.0d0
      b = 0.0d0
      e(n) = 0.0d0
      do 240 l = 1, n
         j = 0
         h = dabs(d(l)) + dabs(e(l))
         if (b .lt. h) b = h
c     .......... look for small sub-diagonal element ..........
         do 110 m = l, n
 110     if (b + dabs(e(m)) .eq. b) go to 120
c     .......... e(n) is always zero, so there is no exit
c                through the bottom of the loop ..........
  120    if (m .eq. l) go to 220
  130    if (j .eq. 30) go to 1000
         j = j + 1
c     .......... form shift ..........
         l1 = l + 1
         l2 = l1 + 1
         g = d(l)
         p = (d(l1) - g) / (2.0d0 * e(l))
         r = pythag(p,1.0d0)
         d(l) = e(l) / (p + dsign(r,p))
         d(l1) = e(l) * (p + dsign(r,p))
         dl1 = d(l1)
         h = g - d(l)
         if (l2 .gt. n) go to 145
         do 140 i = l2, n
  140    d(i) = d(i) - h
  145    f = f + h
c     .......... ql transformation ..........
         p = d(m)
         c = 1.0d0
         c2 = c
         el1 = e(l1)
         s = 0.0d0
         mml = m - l
c     .......... for i=m-1 step -1 until l do -- ..........
         do 200 ii = 1, mml
            c3 = c2
            c2 = c
            s2 = s
            i = m - ii
            g = c * e(i)
            h = c * p
            if (dabs(p) .lt. dabs(e(i))) go to 150
            c = e(i) / p
            r =dsqrt(c*c+1.0d0)
            e(i+1) = s * p * r
            s = c / r
            c = 1.0d0 / r
            go to 160
  150       c = p / e(i)
            r = dsqrt(c*c+1.0d0)
            e(i+1) = s * e(i) * r
            s = 1.0d0 / r
            c = c * s
  160       p = c * d(i) - s * g
            d(i+1) = h + s * (c * g + s * d(i))
c     .......... form vector ..........
            do 180 k = 1, n
               h = z(k,i+1)
               z(k,i+1) = s * z(k,i) + c * h
  180          z(k,i) = c * z(k,i) - s * h
  200    continue
         p = -s * s2 * c3 * el1 * e(l) / dl1
         e(l) = s * p
         d(l) = c * p
         if (b + dabs(e(l)) .gt. b) go to 130
  220    d(l) = d(l) + f
  240 continue
c     .......... order eigenvalues and eigenvectors ..........
      do 300 ii = 2, n
         i = ii - 1
         k = i
         p = d(i)
         do 260 j = ii, n
            if (d(j) .ge. p) go to 260
            k = j
            p = d(j)
  260    continue
         if (k .eq. i) go to 300
         d(k) = d(i)
         d(i) = p
         do 280 j = 1, n
            p = z(j,i)
            z(j,i) = z(j,k)
  280       z(j,k) = p
  300 continue
      go to 1001
c     .......... set error -- no convergence to an
c                eigenvalue after 30 iterations ..........
 1000 ierr = l
 1001 return
      end

      real*8 function pythag(a,b)
c  finds sqrt(a**2+b**2) without overflow or destructive underflow
      implicit real*8(a-h,o-z)
      p = dmax1(dabs(a),dabs(b))
      q = dmin1(dabs(a),dabs(b))
      if (q .eq. 0.0d0) go to 20
   10    r = (q/p)**2
         t = 4.0d0 + r
         if (t .eq. 4.0d0) go to 20
         s = r/t
         p = p + 2.0d0*p*s
         q = q*s
      go to 10
   20 pythag=p
      return
      end

*********************************************
      subroutine staper(nt,fw,nev,v,ndim,a,w)
*********************************************
c  Slepian - Thomson multi-taper procedure
c  Slepian, D.     1978  Bell Sys Tech J v57 n5 1371-1430
c  Thomson, D. J.  1982  Proc IEEE v70 n9 1055-1096
c    nt    the number of points in the series
c    fw    the time-bandwidth product (number of Rayleigh bins)
c    nev   the desired number of tapers
c    v     the eigenvectors (tapers) are returned in v(.,nev)
c    a and w are work arrays dimensioned at least nt long
c    a(1..nev) contains bandwidth retention factors on output
c  The tapers are the eigenvectors of the tridiagonal matrix sigma(i,j) [see Sle
c  pian(1978)
c  eq 14 and 25 ]hey are also the eigenvectors of the Toeplitz matrix eq. 18.
c  We solve
c  the tridiagonal system in tsturm for the tapers and use them in Slepians eq 1
c  8 to get
c  the bandwidth retention factors (i.e. the eigenvalues)
c  note that this only works for nt even ergo we ignore last point if nt is odd
c    tgm and jfg
      implicit real*8(a-h,o-z)
      dimension a(1),w(1),v(ndim,1)
      parameter (pi=3.14159265358979d0,r2=1.414213562373095d0)
      if(nt.lt.2) return
      nx=mod(nt,2)
      lh=(nt/2)+nx
      lp1=nt+1
      om=2.d0*pi*fw/nt
      com=dcos(om)
      hn=0.5d0*lp1
      do 10 i=1,lh
        a(i)=com*(i-hn)**2
   10   w(i)=0.5d0*float(i*(nt-i))
      if(nx.eq.0) then
        asav=a(lh)-w(lh)
        a(lh)=a(lh)+w(lh)
        rbd=1.d0/(a(lh)+w(lh-1))
      else
        asav=w(lh-1)
        rbd=1.d0/(w(lh)+w(lh-1))
        w(lh-1)=r2*w(lh-1)
      endif
      do 15 i=1,lh
        a(i+lh)=w(i)*rbd
        w(i)=a(i+lh)**2
   15   a(i)=a(i)*rbd
      neven=max0((nev+1)/2,1)
      nodd=nev-neven
c***    do the even tapers
      call tsturm(nt,lh,a,a(lh+1),w,neven,v,ndim,w(lh+1),0)
      do 20 i=1,neven
        k=2*i-1
        if(nx.eq.1) v(lh,k)=r2*v(lh,k)
          do 20 j=1,lh
   20     v(lp1-j,k)=v(j,k)
      if(nodd.le.0) goto 34
c***    do the odd tapers
      if(nx.eq.0) then
        a(lh)=asav*rbd
      else
        a(nt)=asav*rbd
        w(lh-1)=asav*asav
      endif
      call tsturm(nt,lh-nx,a,a(lh+1),w,nodd,v,ndim,w(lh+1),1)
      do 30 i=1,nodd
        k=2*i
        if(nx.eq.1) v(lh,k)=0.d0
          do 30 j=1,lh
   30     v(lp1-j,k)=-v(j,k)
   34 ntot=neven+nodd
c***    calculate bandwidth retention parameters
      dc=2.d0*com
      sm=0.d0
      s=dsin(om)
      w(1)=om/pi
      w(2)=s/pi
      do 35 j=3,nt
        sn=dc*s-sm
        sm=s
        s=sn
   35   w(j)=s/(pi*(j-1))
      do 55 m=1,ntot
        vmax=dabs(v(1,m))
        kmax=1
        do 40 kk=2,lh
          if(dabs(v(kk,m)).le.vmax) goto 40
          kmax=kk
          vmax=dabs(v(kk,m))
   40     continue
        a(m)=0.d0
        nlow=kmax-1
          do 45 j=1,nlow
   45     a(m)=a(m)+w(j+1)*v(nlow+1-j,m)
        nup=nt-nlow
          do 50 j=1,nup
   50     a(m)=a(m)+w(j)*v(nlow+j,m)
   55 a(m)=a(m)/v(kmax,m)
      return
      end

******************************************************
      subroutine tsturm(nt,n,a,b,w,nev,r,ndim,ev,ipar)
******************************************************
c    uses bisection and Sturm counting to isolate the eigenvalues of the sym-
c    metric tridiagonal matrix with main diagonal a(.) and sub/super diagonal b(
c   .).
c    Newtons method is used to refine the eigenvalue in subroutine root then
c    direct recursion is used to get the eigenvector as this is always stable.
c    Thomsen's normalisation is used with no attention to sign.
c    Note  ipar=0 for even tapers   =1 for odd tapers
      implicit real*8(a-h,o-z)
      parameter (eps=1.d-15,eps1=5.d-15)
      dimension a(1),b(1),ev(1),w(1),r(ndim,1)
      if(n.le.0.or.nev.le.0) return
      umeps=1.d0-eps
      do 5 i=1,nev
    5 ev(i)=-1.d0
      u=1.d0
      do 1000 ik=1,nev
      if(ik.gt.1) u=ev(ik-1)*umeps
      el=dmin1(ev(ik),u)
   10 elam=.5d0*(u+el)
      if(dabs(u-el).le.eps1) goto 35
      iag=0
      q=a(1)-elam
      if(q.ge.0.d0) iag=iag+1
      do 15 i=2,n
      if(q.eq.0.d0) x=dabs(b(i-1))/eps
      if(q.ne.0.d0) x=w(i-1)/q
      q=a(i)-elam-x
      if(q.ge.0.d0) iag=iag+1
      if(iag.gt.nev) goto 20
   15 continue
      if(iag.ge.ik) go to 20
      u=elam
      go to 10
   20 if(iag.eq.ik) go to 30
      m=ik+1
      do 25 i=m,iag
   25 ev(i)=elam
      el=elam
      go to 10
   30 el=elam
      call root(u,el,elam,a,b,w,n,ik)
   35 ev(ik)=elam
      jk=2*ik+ipar-1
      r(1,jk)=1.d0
      r(2,jk)=-(a(1)-ev(ik))/b(1)
      ddot=1.d0+r(2,jk)*r(2,jk)
      jm1=2
      do 45 j=3,n
      r(j,jk)=-((a(jm1)-ev(ik))*r(jm1,jk)+b(j-2)*r(j-2,jk))/b(jm1)
      ddot=ddot+r(j,jk)*r(j,jk)
   45 jm1=j
      rnorm=dsqrt(nt/(2.d0*ddot))
      do 50 j=1,n
   50 r(j,jk)=r(j,jk)*rnorm
 1000 continue
      return
      end

********************************************
      subroutine root(u,el,elam,a,bb,w,n,ik)
********************************************
c** newtons method to refine root
      implicit real*8(a-h,o-z)
      parameter (eps=1.d-15,eps1=5.d-15)
      dimension a(1),bb(1),w(1)
    5 elam=.5d0*(u+el)
   10 if(dabs(u-el).le.1.5d0*eps1) return
      an=a(1)-elam
      b=0.d0
      bn=-1.d0/an
      iag=0
      if(an.ge.0.d0) iag=iag+1
      do 20 i=2,n
      if(an.eq.0.d0) x=dabs(bb(i-1))/eps
      if(an.ne.0.d0) x=w(i-1)/an
      an=a(i)-elam-x
      if(an.eq.0.d0) an=eps
      bm=b
      b=bn
      bn=((a(i)-elam)*b-bm*x-1.d0)/an
      if(an.ge.0.d0) iag=iag+1
   20 continue
      if(iag.eq.ik) goto 25
      u=elam
      goto 30
   25 el=elam
   30 del=1.d0/bn
      if(dabs(del).le.eps1) del=dsign(eps1,del)
      elam=elam-del
      if(elam.ge.u.or.elam.le.el) goto 5
      goto 10
      end
