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)
      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 = 1
      if (t(2).gt.0.0) 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