subroutine taper(n,nwin,el)
c
c  to generate slepian tapers
c  ta is a real*4 array
c
c         j. park
c
      real*8 el,a,z,pi,ww,cs,ai,an,eps,rlu,rlb
      real*8 dfac,drat,gamma,bh,ell,dumfil
      common/nnaamm/iwflag
      common/npiprol/anpi
      common/tapsum/tapsum(20),ell(20)
      common/work/ip(8194)        
      common/taperzz/z(65536)  ! we use this common block for scratch space
      common/stap2/ta(8192,16),dumfil(500)
      dimension a(8192,8),el(10)
      data iwflag/0/,pi/3.14159265358979d0/
      equivalence (a(1,1),ta(1,1))
      an=dfloat(n)
      ww=dble(anpi)/an
      cs=dcos(2.d0*pi*ww)
c initialize matrix for eispack subroutine
c      type *,'ww,cs,an',ww,cs,an
      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
        a(i+1,3)=a(i+1,2)**2        ! required by eispack routine
      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))
c      type *,ierr,rlb,rlu
c      call tnou('eigenvalues for the eigentapers')
c      type *,(el(i),i=1,nwin)
      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  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 eq 61
      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
c      call tnou('eigenvalues for the eigentapers (small k)')
c      type *,(el(i),i=1,nwin)
      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
c      call tnou('eigenvalues for the eigentapers (k near 2n)')
c      type *,(ell(i),i=1,nwin)
      do i=1,nwin
        el(i)=dmax1(ell(i),el(i))
      end do     
c      call tnou('composite asymptotics for taper eigenvalues')
c      type *,(el(i),i=1,nwin)  
c   normalize the eigentapers to preserve power for a white process
c   i.e. they have rms value unity
      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      type *,'tapsum',(tapsum(i),i=1,nwin)
c  refft preserves amplitudes with zeropadding
c  zum beispiel: a(i)=1.,i=1,100 will transform at zero freq to b(f=0.)=100
c  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
      iwflag=1
      return
      end
c
c
      subroutine envel(ff0,jsmoo)
c   10/25/89
c  made compartmental 7/19/90
c  made into a subroutine 11/20/90
      real ff0
      real*8 cs,sn,env0,amp,ell,dex,ex
      complex*16 g,d,env,dum,za,amp0,amp1,qrsave
      common/pltopt/npad,nst,npts,nmin,nmax,t(3),fmin,fmax,nf,
     x                                       dt,nwin,inc,tt(3)
      common/envel1/ell(20)
      common/taperz/ta(4100,16),tai(4100,16)
      common/staple/env(8200)
      common/stap2/tas(8192,16),ex(500)
      common/stap3/g(150000)
      common/stap4/cs(8200),sn(8200)
      common/stap5/d(100),za(100),qrsave(100),dum2(2000),dum(8200)
      dimension env0(2,1000),amp(2)
      equivalence (env0,env),(amp0,amp)
      tt(1)=0.
      tt(2)=dt
      pi=3.14159265358979
      rad=180./pi
      eps=1.e-4
      iboost=0
      fmin=t(1)
      fmax=t(1)+(nf-1)*t(2)
      qc=0.
      dex=dble(qc*pi/npts)
      call boost(iboost,npts,ex,dex)
c      per=1./ff0
      call fillup(ff0)
c  dot product with decaying envelope for fixed element
c  then integrate the kernels twice
      call premult(iboost,jsmoo)
      call zqrdc(g,npts,npts,nwin,qrsave,ip,dumm,0)
c  solve for the constant term
c  first we calculate double backtransform of data vector 
c  and the triangular matrix r
      call zbacktr(nwin,npts,g,d,dum,2)
c  mult by za to get a scalar     
      amp0=dcmplx(0.d0,0.d0)
      do k=1,nwin
        amp0=amp0+conjg(za(k))*dum(k)  ! note: conjugated
      end do
c  next we calculate double backtransform of za 
c  and the triangular matrix r
      call zbacktr(nwin,npts,g,za,dum,2)
c  mult by za to get a scalar     
      amp1=dcmplx(0.d0,0.d0)
      do k=1,nwin
        amp1=amp1+conjg(za(k))*dum(k)  ! note: conjugated
      end do
      amp0=amp0/amp1
      sum1=0.
      do k=1,nwin
        sum1=sum1+(cdabs(d(k)))**2
      end do
      do k=1,nwin
        d(k)=d(k)-za(k)*amp0
      end do
      sum2=0.
      do k=1,nwin
        sum2=sum2+(cdabs(d(k)))**2
      end do
      rat=sum2/sum1
c  solve for m-tilde
c  we backtransform once, then multiply by q
      do i=1,npts
        env(i)=dcmplx(0.,0.)
      end do
      call zbacktr(nwin,npts,g,d,env,1) 
      call zqrsl(g,npts,npts,nwin,qrsave,env,env,
     x dum,dum,dum,dum,10000,info)
      call postmult(jsmoo,npts,env,amp0,cs)
      npts2=npts*2
      end
c
      subroutine zbacktr(n,npts,g,d,dum,ick)
      complex*16 g(npts,1),dum(1),d(1)
c  g is assumed to be the qr decomposition of a matrix
c  we only use the upper triangle of g
c  ick.eq.1 => only perform backtransform with lower triangle
c  the lower triangular backtransform:
      do k=1,n
        dum(k)=d(k)
      end do
      do i=1,n
        i1=i-1
        if(i.gt.1) then
          do j=1,i1
            dum(i)=dum(i)-conjg(g(j,i))*dum(j)
c            dum(i)=dum(i)-(g(j,i))*dum(j)
          end do
        endif
        dum(i)=dum(i)/conjg(g(i,i))
c        dum(i)=dum(i)/(g(i,i))
      end do
      if(ick.eq.1) return
c  the upper triangular backtransform:
      do i=n,1,-1
        ip1=i+1
        if(i.lt.n) then
          do j=ip1,n
            dum(i)=dum(i)-g(i,j)*dum(j)
c            dum(i)=dum(i)-conjg(g(i,j))*dum(j)
          end do
        endif
        dum(i)=dum(i)/g(i,i)
c        dum(i)=dum(i)/conjg(g(i,i))
      end do
      return
      end
c
      subroutine boost(iboost,npts,ex,dex)
      implicit real*8 (a-h,o-z)
      dimension ex(1)
      if(iboost.eq.1) then
c  boost the envelope by multiplying tapers by decreasing exponential
        dex=dexp(-dex)
        ex(1)=1.d0
        do i=2,npts
          ex(i)=ex(i-1)*dex
        end do
      else  ! dont boost the envelope
        do i=1,npts
          ex(i)=1.d0
        end do        
      endif
      return
      end
c
      subroutine fillup(f0)
c  to fill data vector and kernel matrix
      real*8 cs,sn,ex
      real dumfil2
      complex*16 g,d,dumfil1,dumfil3
      common/pltopt/npad,nst,npts,nmin,nmax,t(3),fmin,fmax,nf,
     x                                dt,nwin,inc,tt(3)
      common/taperz/ta(4100,16),tai(4100,16)
      common/stap2/tas(8192,16),ex(500)
      common/stap3/g(150000)
      common/stap4/cs(8200),sn(8200)
      common/stap5/d(100),dumfil1(200),dumfil2(2000),dumfil3(8200)
      n1=(f0-fmin)/t(2)+1
      f1=(fmin+(n1-1)*t(2))
      df1=f0-f1
      if(df1*2.0.gt.t(2)) then
        n1=n1+1
        f1=f1+t(2)
        df1=f1-f0
      else
        df1=-df1    ! since formula uses f1-f0
      endif
c  ************* careful! make certain that freq and dt are reciprocal units!
c  ************* if dt is sec, df1 is mhz, need to scale by 1000.
c  ************* same if dt is kyr and df1 is cyc/myr
      call cossin(npts,cs,sn,1.d0,0.d0,dble(df1),dble(dt))
c      call cossin(npts,cs,sn,1.d0,0.d0,dble(df1),dble(dt/1000.))
c  assemble data vector
c  multiply by 2.0/inc to obtain proper normalization
c  2.0 from <cos(t)**2>=1/2
c  and inc from decimation of tapers in the data kernels 
      do k=1,nwin
        d(k)=dcmplx(ta(n1,k),tai(n1,k))*2.d0/inc
      end do
c  mult the tapers by sinusoid to get kernels
c  note that we use the complex conjugate
      do k=1,nwin
        ii=(k-1)*npts
        jj=(k-1)*npts
        j=0
        do i=1,npts
          j=j+1
          g(jj+j)=ex(i)*dcmplx(tas(i,k)*cs(j),-tas(i,k)*sn(j))        
        end do
      end do
      return
      end
c
      subroutine premult(iboost,jsmoo)
      real*8 cs,sn,cs0
      real*4 dumfil2
      complex*16 g,d,za,dumfil1,dumfil3
      common/pltopt/npad,nst,npts,nmin,nmax,t(3),fmin,fmax,nf,
     x                                dt,nwin,inc,tt(3)
      common/stap3/g(150000)
      common/stap4/cs(8200),sn(8200)
      common/stap5/d(100),za(100),dumfil1(100),dumfil2(2000),
     $       dumfil3(8200)
      pi=3.14159265358979
      eps=1.e-4
c  dot product with decaying envelope for fixed element
      if(iboost.eq.0) then
        cs0=qc*pi/(npts-1)
        cs0=dexp(-cs0)
        cs(1)=1.d0   
        do i=2,npts
          cs(i)=cs(i-1)*cs0
        end do   
      else   ! if we have boosted the kernels, the fixed element is constant
        do i=1,npts
          cs(i)=1.d0
        end do
      endif
      do k=1,nwin
        za(k)=dcmplx(0.,0.)
        kk=(k-1)*npts
        do i=1,npts
          za(k)=za(k)+cs(i)*g(kk+i)
        end do
        za(k)=conjg(za(k))  ! must conjugate to adhere to conventions
      end do
c  then integrate the kernels twice
      if(jsmoo.gt.0) then
        do ksmoo=1,jsmoo
          do k=1,nwin
            ii=(k-1)*npts
            nptsm=npts-1
            do i=nptsm,1,-1
              g(ii+i)=g(ii+i)+g(ii+i+1)
            end do
c  little fudge factor to ensure penalty of first derivative
            g(ii+1)=g(ii+1)/eps        
          end do
        end do
      endif
      return
      end
c
      subroutine postmult(jsmoo,npts,env,amp0,cs)
      real*8 cs(1)
      complex*16 env(1),amp0
      eps=1.e-4
      if(jsmoo.gt.0) then
c  post normalize to the envelope fluctuation
        do ksmoo=1,jsmoo
          env(1)=env(1)/eps
          do i=2,npts
            env(i)=env(i)+env(i-1)
          end do
        end do
      endif
c  add the decay envelope back in
      do i=1,npts
        env(i)=env(i)+amp0*cs(i)
      end do
      return
      end