subroutine recon(a0,nwin0,npi,dt0,nscan,bndwdth,ff0,icon,rcon)
c
c     perform time-reconstruction of signal based on inversion
c     of multitaper envelope.
c
c     modified by M.E. Mann 12/96 based on original code of J. Park
c
c  program is limited to 16 windows
c  program is limited to 4100 frequency-ftest estimation pairs
c
      parameter (maxlen=50000,maxsignal=20)
      real a0(maxlen),dt0
      integer nwin0,nscan,npi
      real*8 ell,env0,cs,sn,ex
      real bndwdth,f0,ff0
      character*4 chead(158)
      complex zed
      common/npiprol/anpi
      common/header/iah(158)
      common/data/a(maxlen)
      common/nnaamm/iwflag
      common/work/b(8194)
      common/pltopt/npad,nst,mpts,nmin,nmax,t(3),fmin,fmax,nf,
     x                                ddt,nwin,inc,tt(3)
      common/taperz/ta(4100,16),tai(4100,16)
      common/staple/tad(4000)
      common/stap2/tas(8192,16),ex(500)
      common/stap4/cs(8200),sn(8200)
      common/envel1/ell(20)
      dimension dum(35),ahead(158),env0(2,1000)
      dimension reim(2),time(2)
      real *8 b1(8194),b2(8194)
      real env(3,500,2),aa(8200,3),acent(8200),rcon(maxlen)
      equivalence (zed,reim),(iy,dum),(iah,ahead),(iah,chead),
     x          (tad,env0)
      iwflag=0
      pi=3.14159265358979
c
      do i=1,maxlen
         a(i)=a0(i)
      end do
      time(1)=0.
      dt = dt0
      time(2)=dt/1.
      fzero=0.
      fny=0.5/dt
      npts=nscan
      nwin = nwin0
      anpi = float(npi)
c
c     set frequency range
c
      f1 = 0.0
      f2 = fny
c
      npad=2
      do while (npad.lt.nscan)
        npad=npad*2
      end do
      npad=8192
      ddf=1./(npad*dt)
      eps=.1*ddf
      nn1=(f1+eps)/ddf
      fmin=nn1*ddf
      nn2=(f2+eps)/ddf
      fmax=nn2*ddf
      nf=nn2-nn1+1
      nmin=nn1*2+1
      nmax=nn2*2+2
      t(1)=fmin
      t(2)=ddf
      t(3)=ddf
      ftfrq=ddf
  102 format(i6,g12.4,i6)
      n1=1
      n2=nscan
 2660 call taper1(npts,nwin,ell)
c
c     demean series
c 
      demn=0.
      do i=n1,n2
         demn=demn+a(i)
      end do
      demn=demn/(n2-n1+1)
      do i=n1,n2
         acent(i)=a(i)-demn
      end do
c
c  loop over # of windows
c  normalization anrm:  divide by npts as we wish spectrum per unit time
      anrm=npts
      icont=0
      do iwin=1,nwin
        do i=n1,n2
          b1(i-n1+1)=acent(i)*tas(i+1-n1,iwin)
          b2(i-n1+1)=0.0
        end do
        j=npts+1
        do i=j,npad
          b1(i)=0.
          b2(i)=0.
        end do
        call fft2(b1,b2,npad)
        j = 0
        do i=1,npad/2
          j = j+1
          b(i)=b1(i)
          b(i+1)=b2(i)
          ta(i,iwin)=b1(i)
          tai(i,iwin)=b2(i)
        end do
      end do
c
c     decimate the tapers
c
1660  inc=(npts-1)/500+1
      j=0
      do i=1,npts,inc
        j=j+1
        do k=1,nwin
          tas(j,k)=tas(i,k)
        end do
      end do
      mpts=j
      ddt=dt*inc
c
c     need to take into account bandwith retention factor
c     (leads to double counting of secular variance)
c
c     if specified carrier frequency is within one bandwidth
c     of f=0, then specify its carrier as zero to avoid
c     difficulty in estimating bandwidth retention factor
c     for 0<f<bndwdth
c
      f0=ff0
      if (f0.lt.bndwdth) f0=0.0
c
c     reconstruction
c
      anorm = 1.0
      if (f0.lt.eps) anorm=0.5
c
c
      jmax = 1
      if (icon.eq.0) jmax=3
      do i=1,jmax
        itype = icon-1
        if (icon.eq.0) itype=i-1
        call envel(f0,itype)
        do j=1,mpts
           env(i,j,1)=env0(1,j)
           env(i,j,2)=env0(2,j)
        end do
      end do
c
      call cossin(npts,cs,sn,1.d0,0.d0,dble(f0),dble(dt))
c
      finc=float(inc)
      do jj=1,jmax
        do i=1,mpts
           env0(1,i)=env(jj,i,1)
           env0(2,i)=env(jj,i,2)
        end do
c
c       linear interpolate the envelope from mpts (<500) to npts
c
        do i=1,mpts-1
          ii=(i-1)*inc
          do j=1,inc
            a1=(env0(1,i)*(inc+1-j)+env0(1,i+1)*(j-1))/finc
            a2=(env0(2,i)*(inc+1-j)+env0(2,i+1)*(j-1))/finc
            a(n1-1+nscan+ii+j)=a1*cs(ii+j)+a2*sn(ii+j)
            aa(n1-1+ii+j,jj)=anorm*a(n1-1+nscan+ii+j)
          end do
        end do
c
c       linear extrapolate envelope one point 
c     
        ii=(mpts-1)*inc
        e1=2.*env0(1,mpts)-env0(1,mpts-1)
        e2=2.*env0(2,mpts)-env0(2,mpts-1)
        do j=1,inc
          a1=(env0(1,mpts)*(inc+1-j)+e1*(j-1))/finc
          a2=(env0(2,mpts)*(inc+1-j)+e2*(j-1))/finc
          a(n1-1+nscan+ii+j)=a1*cs(ii+j)+a2*sn(ii+j)
          aa(n1-1+ii+j,jj)=anorm*a(n1-1+nscan+ii+j)
        end do
      end do
c
c     default constraint = pre-specified
c
      w1keep = 1.0
      w2keep = 0.0
      w3keep = 0.0
c
c     if indicated, determine weights on the 3 possible
c     reconstructions based on minimum-misfit criterion
c
      abest = 1.0e+36 
      step = 0.02
      if (icon.eq.0) then
        do w1=0,1,step
          do w2=0,1,step
            w3 = 1.0-w1-w2
            if ((w3.ge.0.0).and.(w3.le.1.0)) then
              amiss = 0.0
              do i=1,nscan
                signal = w1*aa(i,1)+w2*aa(i,2)+w3*aa(i,3)
                amiss = amiss + (signal-acent(i))**2
              end do
              if (amiss.lt.abest) then
                 abest=amiss
                 w1keep = w1
                 w2keep = w2
                 w3keep = w3
              endif
            endif
           end do
        end do
        do i=1,nscan
          rcon(i)=w1keep*aa(i,1)+w2keep*aa(i,2)+w3keep*aa(i,3)
        end do
      else
        do i=1,nscan
           rcon(i)=aa(i,1)
        end do
        do i=nscan+1,maxlen
           rcon(i)=0.0
        end do
      endif
      return
      end
c
c
      subroutine taper1(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,ex
      real*8 dfac,drat,gamma,bh,ell
      common/nnaamm/iwflag
      common/npiprol/anpi
      common/tapsum/tapsum(20),ell(20)
      common/work/ip(8194)        
      common/taperz/z(65600)  ! we use this common block for scratch space
      common/stap2/ta(8192,16),ex(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)
      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))
      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  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
      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  ! note that this doesn't hold 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
      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
  101 format(80a)
c
c     The real DFT preserves amplitudes with zeropadding
c     therefore we need not doctor the taper normalization,
c      but wait until the fft to force the desired units
c
      iwflag=1
      return
      end
c
      subroutine envel(f0,inorm)
c
      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(1000)
      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
      jsmoo = inorm
      iboost=0
      fmin=t(1)
      fmax=t(1)+(nf-1)*t(2)
      ff0 = f0
      qc=0.
      dex=dble(qc*pi/npts)
      call boost(iboost,npts,ex,dex)
      call fillup(f0)
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)
      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
c  subtract from eigenspectra
      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
      return
      end
c
      subroutine zbacktr(n,npts,g,d,dum,ick)
      complex*16 g(npts,1),dum(1),d(1)
c
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:
c
      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
c  to fill data vector and kernel matrix
c
      real*8 cs,sn,ex
      complex*16 g,d,junk1,junk3
      real *4 junk2
      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),junk1(200),junk2(2000),junk3(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
      complex*16 g,d,za,junk1,junk3
      real *4 junk2
      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),junk1(100),
     x    junk2(2000),junk3(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
c
      subroutine gfill(npts,nwin,ib,jb,df,ddt,g,ex,tas)
c  will fill ib,jb block of kernel matrix g with tapers
c  g has structure
c
c               |g11 g12|
c          g  = |g21 g22|
c
c  where each block is nptsxnwin, gij is ith carrier freq and jth fft freq
c  here ib is i jb is j
c
      implicit real*8 (a-h,o-z)
c      real*4 tas,tt(3),dum2
      real*4 tas,tt(3)
      real*8 junk1
      complex*16 g
      common/stap4/cs(500),sn(500),junk1(15400)
      dimension g(npts,2,nwin,2),ex(1),tas(8192,1)
c
      call cossin(npts,cs,sn,1.d0,0.d0,df,ddt)
c  mult the tapers by sinusoid to get kernels
c  note that we use the complex conjugate
      tt(1)=0.
      tt(2)=ddt
      do k=1,nwin
        do i=1,npts
          g(i,ib,k,jb)=ex(i)*dcmplx(tas(i,k)*cs(i),-tas(i,k)*sn(i))        
        end do
      end do
      return
      end