c
      subroutine hires(ta,tai,el,nwin,nf,ares)
      real*8 el
      dimension ta(16400,1),tai(16400,1),el(1),ares(1)
      do j=1,nf
        ares(j)=0.
      end do
      do i=1,nwin
        a=1./(el(i)*nwin)
        do j=1,nf
          ares(j)=ares(j)+a*(ta(j,i)*ta(j,i)+tai(j,i)*tai(j,i))
        end do
      end do
      do j=1,nf
        ares(j)=sqrt(ares(j))
      end do
      return
      end
c
c
      subroutine adwait(ta,tai,dcf,el,nwin,nf,ares,degf,avar)
c
c  this version uses Thomson's algorithm for calculating 
c  the adaptive spectrum estimate
c
      real*8 avar,spw,as,das,tol,el,a1,bias,scale,ax,fn,fx
      dimension ta(16400,1),tai(16400,1),el(1),ares(1),degf(1)
      dimension spw(10),bias(10),dcf(16400,1)
c
c  set tolerance for iterative scheme exit
c
      tol=3.d-4
      jitter=0
      scale=avar
c
c  we scale the bias by the total variance of the frequency transform
c  from zero freq to the nyquist
c  in this application we scale the eigenspectra by the bias in order to avoid
c  possible floating point overflow
c
      do 200 i=1,nwin
  200 bias(i)=(1.d0-el(i))
      do 100 j=1,nf
        do 150 i=1,nwin
  150   spw(i)=(ta(j,i)*ta(j,i)+tai(j,i)*tai(j,i))/scale
c
c  first guess is the average of the two lowest-order eigenspectral estimates
c
        as=(spw(1)+spw(2))/2.d0
        do 300 k=1,20
c
c  find coefficients
c
          fn=0.d0
          fx=0.d0
          do 350 i=1,nwin
            a1=dsqrt(el(i))*as/(el(i)*as+bias(i))
            a1=a1*a1
            fn=fn+a1*spw(i)
            fx=fx+a1
  350     continue
          ax=fn/fx
          das=dabs(ax-as)
          if(das/as.lt.tol) go to 400
  300   as=ax
c
c  flag if iteration does not converge
c
      jitter=jitter+1
  400 continue
      ares(j)=as*scale
c
c  calculate degrees of freedom
c
      df=0.
      do 450 i=1,nwin
      dcf(j,i)=dsqrt(el(i))*as/(el(i)*as+bias(i))
  450 df=df+dcf(j,i)*dcf(j,i)
c
c  we normalize degrees of freedom by the weight of the first eigenspectrum
c  this way we never have fewer than two degrees of freedom
c
  100 degf(j)=df*2./(dcf(j,1)*dcf(j,1))
      return
      end
c
      subroutine regre(sr,si,nf,nwin,amu)
c
      real b
      real *8 junk
      dimension sr(16400,1),si(16400,1),amu(16400,1)
      common/tapsum/b(10),junk(10)
c
c     "b" is the DFT of Slepian eigentapers at zero frequency
c     "sr" and "si" are the eigenspectra
c     "amu" contains line frequency estimates and f-test parameter
c
      sum=0.
      do i=1,nwin
        sum=sum+b(i)*b(i)
      end do
      do i=1,nf
        amu(i,1)=0.
        amu(i,2)=0.
        do j=1,nwin
          amu(i,1)=amu(i,1)+sr(i,j)*b(j)
          amu(i,2)=amu(i,2)+si(i,j)*b(j)
        end do
        amu(i,1)=amu(i,1)/sum
        amu(i,2)=amu(i,2)/sum
        sum2=0.
        do j=1,nwin
          sumr=sr(i,j)-amu(i,1)*b(j)
          sumi=si(i,j)-amu(i,2)*b(j)
          sum2=sum2+sumr*sumr+sumi*sumi
        end do
        amu(i,3)=(nwin-1)*(amu(i,2)**2+amu(i,1)**2)*sum/sum2
      end do
      return
      end
c
c
      subroutine taper2(n,nwin,el)
c
c     generate slepian tapers
c     "ta" is a real*4 array
c
c     written by J. Park
c
      real*8 el,a,z,pi,ww,cs,ai,an,eps,rlu,rlb
      real*8 dfac,drat,gamma,bh,ell
      real dump
      common/npiprol/anpi
      common/tapsum/tapsum(10),ell(10)
      common/misc/npad
      common/work/ip(32800)        
      common/taperz/z(65536),dump(393728)
      common/stap2/ta(32800,16)
      dimension a(32800,8),el(10)
      data pi/3.14159265358979d0/
      equivalence (a(1,1),ta(1,1))
      an=dfloat(n)
      ww=dble(anpi)/an
      cs=dcos(2.d0*pi*ww)
c
c     note: 
c
c initialize matrix for eispack subroutine
c this matrix is not the bandwidth retention factor matrix A (for minimizing
c spectral leakage) described in various multitaper papers 
c --e.g. Thomson (1982) [proc ieee], Park et al (1987) [jgr]. 
c the bandwidth matrix x returns a cluster of eigenvectors 
c (which are the slepian tapers) with eigenvalues very close to unity.
c since the spacing of eigenvalues can be comparable to machine precision,
c numerical instability can occur for time-bandwidth product n .ge. 6
c (that is, 6pi-prolate tapers)
c also, the bandwidth retention matrix A is toeplitz, but full, and it is not
c feasible to calculate tapers explicitly for long time series (in an earlier
c code i interpolated an m-point taper from the eigenvectors of a 128x128
c bandwidth retention matrix). the following is cribbed from a 1978 paper
c by Slepian (Prolate spheroidal wave functions, Fourier analysis 
c and uncertainty, Bell System Tech Journal, v57, pp1371-1430, 1978)
c which gives a three-term recursion that is satisfied by the
c slepian tapers. the three-term recursion can be manipulated to show that the
c slepian tapers are the eigenvectors of a tridiagonal matrix a, with
c well-spaced eigenvalues that are (practically speaking) unrelated to the 
c eigenvalues of x. using this matrix a we obtain both numerical stability
c and speed (tridiagonal matrices can be decomposed rapidly enough to
c calculate m-point tapers on the fly).
c
      do i=0,n-1
        ai=dfloat(i)
        a(i+1,1)=-cs*((an-1.d0)/2.d0-ai)**2
        a(i+1,2)=-ai*(an-ai)/2.d0
c
c       next statement is eispack routine tridib, see its documentation
c
        a(i+1,3)=a(i+1,2)**2        
      end do
      eps=1.e-13
      m11=1
      call tridib(n,eps,a(1,1),a(1,2),a(1,3),rlb,rlu,m11,nwin,el,ip,
     x       ierr,a(1,4),a(1,5))
      call tinvit(n,n,a(1,1),a(1,2),a(1,3),nwin,el,ip,z,ierr,
     x            a(1,4),a(1,5),a(1,6),a(1,7),a(1,8)) 
c
c  note: 
c     
c  we calculate the eigenvalues of the dirichlet-kernel problem
c  i.e. the bandwidth retention factors
c  from slepian 1978 asymptotic formula, gotten from thomson 1982 eq 2.5
c  supplemented by the asymptotic formula for k near 2n from Slepian (1978) 
c  more precise values of these parameters, perhaps useful in adaptive
c  spectral estimation, can be calculated explicitly using the 
c  rayleigh-quotient formulas in Thomson (1982) and Park et al (1987)
c 
      dfac=an*pi*ww
      drat=8.d0*dfac
      dfac=4.d0*dsqrt(pi*dfac)*dexp(-2.d0*dfac)
      do k=1,nwin
        el(k)=1.d0-dfac
        dfac=dfac*drat/k  ! is this correct formula? yes,but fails as k -> 2n
      end do
      gamma=dlog(8.d0*an*dsin(2.d0*pi*ww))+0.5772156649d0
      do k=1,nwin
        bh=-2.d0*pi*(an*ww-dfloat(k-1)/2.d0-.25d0)/gamma
        ell(k)=1.d0/(1.d0+dexp(pi*bh))
      end do
      do i=1,nwin
        el(i)=dmax1(ell(i),el(i))
      end do 
c    
c  normalize the eigentapers to preserve power for a white process
c  i.e. they have rms value unity
c  "tapsum" is the average of the eigentaper, should be near zero for
c  antisymmetric tapers
c
      do k=1,nwin
        kk=(k-1)*n
        tapsum(k)=0.
        tapsq=0.
        do i=1,n
          aa=z(kk+i)
          ta(i,k)=aa
          tapsum(k)=tapsum(k)+aa
          tapsq=tapsq+aa*aa
        end do
        aa=sqrt(tapsq/n)
        tapsum(k)=tapsum(k)/aa
        do i=1,n
          ta(i,k)=ta(i,k)/aa
        end do
      end do
c
c  the real FFT will preserve amplitudes with zeropadding
c  for example, a(i)=1.,i=1,100 will transform at zero freq 
c  to b(f=0)=100 no matter how much zero padding is done
c  therefore we need not doctor the taper normalization,
c  but wait until the fft to force the desired units
c
      return
      end