c  program rfmig_boot
c  10/12/04 JJP -- adapted from rfmigrate  -- UPDATED 12/01/07 JJP
c  xf77 -o /Users/jjpark/bin/rfmig_boot rfmig_boot.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c  code computes frequency-domain stacks of receiver functions that follow a harmonic expansion in baz
c  for both radial and transverse RFs there are constant terms and sin/cos terms for 2- and 4-lobed
c  amplitude dependence.  The constant term should be zero for the transverse RF.  
c  The 2-lobed terms govern dipping interface effects and tilted symmetry-axis ansotropy. 
c  The 4-lobed term is anisotropy with a horizontal axis
c    The code regresses for the harmonic expansion, and bootstrap-resamples the data to estimate the 
c    uncertainty of the harmonic terms.
c  
c   output files are 	out[rt]_bexp.grid  -- harmonic-expansions of the RFs, in time domain
c     			out[rt]1_bexp.grid  -- harmonic-expansion RFs plus bootstrap uncertainty
c     			out[rt]2_bexp.grid  -- harmonic-expansion RFs minus bootstrap uncertainty
c     			out[rt]_bbaz.grid  -- harmonic-expansion RFs computed for ordered baz values
c
c  migrates MTC receiver functions in the frequency domain.
c  requires a stacking model in the anirec format
c  such a model may have anisotropy parameters, 
c  but migration code only uses the isotropic velocities. 
c
c  has kluge to cheat the pre-event noise for synthetic records  3/12/00 JJP
c  check to see if the kluge is commented out
c
c  this version of the RF code reads a file of data filenames 
c  you have two choices: either read the time intervals in the filename file
c  or read them in the sac header
c  the data must be binary SAC format
c  horizontals must be rotated to radial and transverse
c
c  for start times in the file:
c  the file is "in_recfunk" and has lines of the form:
c
c  1997.070.19.33.bh?  <-- code replaces ? with r,t,z
c  57 52 	       <-- start time of analysis window, window duration (sec)
c  1997.076.08.15.bh?
c  62 62 
c  ...
c  ...
c  ...
c  stop                <-- tells code that data is finished, 799 events max
c
c
c  for start times in the SAC header
c  reads seismic record start times from the sac header
c  will search the SAC header for specific markers of P phases
c  T1 - P, Pdiff    ahead(12)
c  T2 - PKP,PKIKP   ahead(13)
c  T3 - PP          ahead(14)
c  T1=T2=T3=0 ==> use original A-marker ahead(9)
c
c   code does NOT combine data with different sample rates
c  data files limited to 99K pnts. To increase, see common block /datastuff/
c
c   many intermediate quantities are plotted with PLOTIT as the code proceeds.
c  other intermediate quantities can be plotted by uncommenting calls to PLOTIT
c
c  The code computes a 
c  the code writes the BAZ-dependent RFs to files
c  in a format easily digested by GMT (traces are separated by ">" lines)
c
c  filenames: out[rt]_baz.grid oum[rt]_epi.grid oum[rt]_baz.grid
c
c  these files are overwritten everytime you run the program
c  so rename them to save them
c
      implicit real*4(a-h,o-z)
      implicit integer*4(i-n)
      real*8 ar,ai,el(12),aaa,bootgg,bootdd,bootmm
      complex*8 zc,zero,bootd,bootm,bootma,bootmv,cdev
      complex*8 afft,rf,rfs
      character*80 name,subtext,string,name_stack
      character*28 title(2)
      character*10 cmps(3)
      character*18 xlabel(2)
      character*14 xlabel2(2)
      character*12 xlabel3(2)
      character comp(3)
      common/nnaamm/iwflag
      common/npiprol/anpi
      common/stap2/tap(16400,16)
      common/taperz/ta(16400,8,3),tai(16400,8,3)
      common/stap0/ar(16400),ai(16400),bb(16400),time(16400)
      common/specs/afft(4100,8,3),rf(4100,2),sss(4100,6),crf(4100,2)
      common/saveit/rft(4100,799,2),crft(4100,2),rfs(4100,799,2),
     x   drfs(4100,799,2),bazs(799),gcarc(799),iflag(799),
     x                         epic(799),s2n(799),ip(799)
      common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2)
      common/boot/bootg(5,5,4100,2,21),bootd(5,4100,2,21),
     x            bootm(100,5,4100,2,21)
      common/boot1/bootma(5,4100,2,21),bootmv(5,5,4100,2,21)
      common/boot2/bootgg(5,5),bootdd(5),bootmm(5,2),bf(5),cdev(5)
c NOTE THAT MODEL PARAMETERS ARE REAL*4 IN THIS CODE, NOT REAL*8 AS IN ANIREC
c ONLY 20 LAYERS OVER HALFSPACE ARE ALLOWED
      common/model/z(21),dz(21),rho(21),vp(21),vs(21),vp2(21),
     x               vp4(21),vs2(21),vss(21)
      common/model2/xmu(21),xla(21),xmu2(21),xla2(21),xla4(21)
      common/modstk/tau(21),itau(21),psd(21),gam(21),xik(21)
      common/migstk/arm(4100,2,21),aim(4100,2,21),dam(4100,2,21)
      common/tt_p/d_p(200),d_pkp(200),s_p(200),s_pkp(200),n_p,n_pkp
      common/header/ahead(158)
      dimension iah(158),chead(158),mmonth(12),fr(3),tt(3),ttt(3),
     x  xx(2),yy(2)
      equivalence (ahead,iah),(ahead,chead)
      data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
      data comp/'r','t','z'/
      data cmps/'Radial    ','Transverse','Vertical  '/
      data title/'Radial Receiver Function    ',
     x 'Transverse Receiver Function'/
      data xlabel/'H\\sub{R}(f) phase','H\\sub{T}(f) phase'/
      data xlabel2/'|H\\sub{R}(f)|','|H\\sub{T}(f)|'/
      data xlabel3/'H\\sub{R}(t)','H\\sub{T}(t)'/
      pi=3.14159265358979
      con=180./pi
      pih=pi/2.
      zero=cmplx(0.,0.)
      fmax=5.
c  nboot is the number of random resamplings in the bootstrap
c  to get an answer quickly and ignore uncertainty, one chooses small nboot
c  max nboot is 100 (first dimension of bootm)
  321 print *,'enter number of bootstrap iterations for ',
     x 'RF uncertainty estimate (nboot.le.100)'
      print *,'nboot=0 skips bootstrap and compute simple regression'
      read(5,*) nboot
      if(nboot.gt.100.or.nboot.lt.0) go to 321
      print *,'enter velocity model for stacking'
      print *,'leading space = <stack_model>'
      read(5,101) name_stack
      print *,'enter sampling time e.g., 0.05 for 20 sps'
      read(5,*) dt
      if(name_stack(1:1).eq.' ') then
        open(7,file='stack_model',form='formatted')
      else
        open(7,file=name_stack,form='formatted')
      endif
      read(7,101) subtext
      print *,subtext
      read(7,*) nl
      if(nl.gt.20) then
        print *,'only ten layers in stack model are allowed!'
        stop
      endif
c  read in theta,phi in degrees - polar coords of fast axis
      nlp=nl+1
      nlm=nl-1
      do i=1,nlp
        read(7,*) theta,phi
c  read depth to ith interface, vp (m/sec), pk-to-pk cos(2th) relative P pert
c  pk-to-pk cos(4th) relative P pert, v_s,  pk-to-pk cos(2th) relative S pert
c  density (kg/m**3)
        read(7,*) z(i),vp(i),vp2(i),vp4(i),vs(i),vs2(i),rho(i)
c  recall that we interpret fractional values of b,c,e 
c  as peak-to-peak relative velocity perts.
c  therefore, e=0.02 is 2% pert to mu from slowest to fastest
c  DO NOT NORMALIZE THE MODEL
        xmu(i)=rho(i)*vs(i)**2
        xmu2(i)=vs2(i)*xmu(i)
        xla(i)=rho(i)*vp(i)**2
        xla2(i)=vp2(i)*xla(i)
        xla4(i)=vp4(i)*xla(i)
        vs(i)=vs(i)
        vp(i)=vp(i)
        rho(i)=rho(i)
        z(i)=z(i)
      end do
      close(7)
      do i=2,nl
        dz(i)=z(i)-z(i-1)
      end do
      dz(1)=z(1)
      sdelay=0.
      pdelay=0.
      do i=1,nl
        sdelay=sdelay+(dz(i)/vs(i))
        pdelay=pdelay+(dz(i)/vp(i))
        tau(i)=sdelay-pdelay
        itau(i)=tau(i)/dt
        print *,itau(i),tau(i)
      end do
      itau(nlp)=4100
      print *, 'organ-pipe mode count at 1 Hz in stack of layers: S & P'
      print *,'These are the S & P traveltimes from the basal interface'
      print 104,sdelay,pdelay
  104 format(2f10.1)
      print *,'enter fmax (Hz)'
      read(5,*) fmax
c  some default values
      anpi=2.5
      nwin=3
      npad=16384
      nnf=npad/2
      do i=1,4100
        crft(i,1)=0.
        crft(i,2)=0.
      end do
      irec=0
 1010 print *,'User-input time windows (0) or sac-header picks (1)'
      read(5,*) isegment
      if(isegment.eq.0) then
        print *,
     x'file of seismic records to read? (space-ret: in_recfunk)'
        read(5,101) name
        if(name(1:1).eq.' ') then
          open(10,file='in_recfunk',form='formatted')
        else
          open(10,file=name,form='formatted')
        endif
      elseif(isegment.eq.1) then
        print *,
     x'file of seismic records to read? (space-ret: in_recpick)'
        read(5,101) name
        if(name(1:1).eq.' ') then
          open(10,file='in_recpick',form='formatted')
        else
          open(10,file=name,form='formatted')
        endif
	print *,'duration of data windows for analysis'
	read(5,*) postwin
      else
        go to 1010
      endif
      print *,'rotate to LQT coordinates to isolate P and SV? (0=no)'
      read(5,*) lqt
   10 print *,'enter filename e.g. 1998.361.bhz'
      read(10,101) name
  101 format(a)
      if(name(1:4).eq.'stop') go to 111
c  we assume that the last character is either 'r','t',or 'z'
      do i=80,1,-1
        if(name(i:i).ne.' ') then
          namlen=i
          go to 98
        endif
      end do
      print *,name
      print *,'no filename?'
      stop
   98 continue
      do kcmp=1,3
        name(namlen:namlen)=comp(kcmp)
        print *,name
        call sacread(name,a(1,kcmp),ierr)
        dt=ahead(1) 
        tt(1)=0.
        tt(2)=dt
        tt(3)=dt
        nscan=iah(80)
        if(nscan.gt.99000) then
          print *,'CAREFUL! Data length is greater than 99000:', nscan
          pause
        endif
      end do
      irec=irec+1
      if(irec.gt.799) then
        print *,'max records exceeded'
	stop
      endif
      npre=10./dt
      npost=30./dt
      ntot=npre+npost
      ttt(1)=-npre*dt
      ttt(2)=dt
      ttt(3)=dt
      do i=1,ntot
        time(i)=-ttt(1)-(i-1)*ttt(2)
      end do
      baz=ahead(53)
      bazs(irec)=ahead(53)
      gcarc(irec)=ahead(54)
      epic(irec)=ahead(54)
c      do i=0,79,10
c        print *,(ahead(ii+i),ii=1,10)
c      end do
      print *,bazs(irec),'= back azim ',ahead(54),'= epicentral dist'
      if(ahead(53).lt.-0.1.or.ahead(53).gt.360.1) then
        print *,'HEY! THE BACK AZIMUTH OF THIS RECORD IS OUT OF RANGE!'
        stop
      elseif(ahead(54).lt.-0.1.or.ahead(54).gt.180.1) then
        print *,'HEY! THE GCARC OF THIS RECORD IS OUT OF RANGE!'
        stop
      endif
c  ddf is cycles/day
      ddf=1./(dt*npad)
      nf=fmax/ddf+1
      fr(1)=0.
      fr(2)=ddf
      fr(3)=ddf
C  pre and post-window durations are the same
c  program pads with zeroes if start time is too early
      if(isegment.eq.0) then
        print *,'P window start time, duration (sec)'
        read(10,*) pstart,postwin
	imark=0
      else
c   READING TIME INTERVAL FROM SAC HEADER
c  ahead(12) is SAC's T1 parameter, for arrival time.  Set by pressing T & 1 in plotpk
        if(ahead(12).gt.1) then
          tmark=ahead(12)
	  imark=1
c  ahead(13) is SAC's T2 parameter, for arrival time.  Set by pressing T & 2 in plotpk
        elseif(ahead(13).gt.1) then
          tmark=ahead(13)
	  imark=2
c  ahead(14) is SAC's T3 parameter, for arrival time.  Set by pressing T & 3 in plotpk
        elseif(ahead(14).gt.1) then
          tmark=ahead(14)
	  imark=3
c  ahead(9) is SAC's A parameter, for arrival time.  Set by pressing A in plotpk
        elseif(ahead(9).gt.1) then
          tmark=ahead(9)
	  imark=0
        else
          print *,' hey! no P-wave mark for file ',name
	  stop
        endif
        pstart=tmark-ahead(6)-3.
        print *,'sac params b and tmark, pstart',ahead(6),tmark,pstart
      endif
      npts=postwin/dt
      nst=pstart/dt
      call taper(npts,nwin,el)
      print *,'eigenvalues ',(el(k),k=1,nwin)
c  loop over components
      print *,'npad,nf,nst,npts',npad,nf,nst,npts
      print *,'dt,nscan,ddf',dt,nscan,ddf
c      pause
      npts0=min0(nst+1,npts)
      if(npts0.lt.npts) then
        call taper(npts0,nwin,el)
        print *,'eigenvalues ',(el(k),k=1,nwin)
      endif
      epipi=epic(irec)
c  preferred option -- use parameterized relation between GCARC and zrrot
c  compute approximate angle between Z and R components for P wave
c  if the phase is PP, its moveout is the same as a P wave at gcarc/2
      iflag(irec)=1
      if(imark.eq.3) then
	epipi=epipi/2.
	gcarc(irec)=epipi
      endif
      if(imark.eq.2) then
	zrrot=9.0
	iflag(irec)=2
      else
        if(epipi.lt.80.) then
          zrrot=42.+(20.-epipi)*(20./60.)
        elseif(epipi.lt.118.) then
          zrrot=16.+(6./900.)*(110.-epipi)**2
        elseif(imark.eq.0) then
          zrrot=9.0
	  iflag(irec)=2
	else
          zrrot=16.+(6./900.)*(110.-epipi)**2
	endif
      endif
c  here is the option to rotate the Radial and vertical Components to
c  a P and SV coordinate system    
      if(lqt.ne.0) then
        zrrot=zrrot/con
        cs=cos(zrrot)
        sn=sin(zrrot)
c  rotate from Z-R coords to L-Q, similar to P-SV
        do i=1,nst+npts
          plog=cs*a(i,3)+sn*a(i,1)
	  a(i,1)=cs*a(i,1)-sn*a(i,3)
   	  a(i,3)=plog
        end do
      endif
      do kcmp=1,3
c  first: compute spectrum of pre-event noise
        do i=1,nf
          sss(i,kcmp)=0.
        end do
        do k=1,nwin
          do i=1,npts0
c OK to reverse order of data here, as we only use mod-sq spectrum
            ar(i)=a(nst+1-i,kcmp)*tap(i,k)
            ai(i)=0.d0
          end do
          do i=npts0+1,npad
            ar(i)=0.d0
            ai(i)=0.d0
          end do   
c  we use the complex FFT routine fft2 -- uses real*8 input arrays
c      print *,'we use the complex FFT routine fft2'
          call fft2(ar,ai,npad)
          do i=1,nf
            sss(i,kcmp)=sss(i,kcmp)+(ar(i)**2+ai(i)**2)
          end do
        end do
c  if pre-event noise duration is truncated, boost the noise spectrum by
c  relative length factor
        if(npts0.lt.npts) then
          do i=1,nf
            sss(i,kcmp)=sss(i,kcmp)*float(npts)/float(npts0)
          end do
        endif
      end do
      if(npts0.lt.npts) then
        call taper(npts,nwin,el)
        print *,'eigenvalues ',(el(k),k=1,nwin)
      endif
      do kcmp=1,3
        do i=1,nf
          sss(i,kcmp+3)=0.
        end do
c  first: compute spectrum of P coda
        do k=1,nwin
          do i=1,npts
            ar(i)=a(i+nst,kcmp)*tap(i,k)
            bb(i)=a(i+nst,kcmp)*tap(i,k)
            ai(i)=0.d0
          end do
c          call plotit(tt,bb,dum,npts,'tapered data',
c     x     'time',' ',1,0,0.0,0,1)
          do i=npts+1,npad
            ar(i)=0.d0
            ai(i)=0.d0
          end do
c  we use the complex FFT routine fft2 -- uses real*8 input arrays
c      print *,'we use the complex FFT routine fft2'
          call fft2(ar,ai,npad)
          do i=1,nf
            afft(i,k,kcmp)=dcmplx(ar(i),ai(i))
            sss(i,kcmp+3)=sss(i,kcmp+3)+(ar(i)**2+ai(i)**2)
          end do
        end do
      end do
c  simple correlation coefficient with vertical (kcmp=3)
      do l=1,2
        do n=1,nf
          zc=zero
          do k=1,nwin
            zc=zc+conjg(afft(n,k,3))*afft(n,k,l)
          end do
          rf(n,l)=zc/(sss(n,6)+sss(n,3))
          crf(n,l)=cabs(zc)**2/(sss(n,6)*sss(n,3+l))
          drf(n,l)=sqrt((1.0-crf(n,l))/(crf(n,l)*(nwin-1)))
     x                                         *cabs(rf(n,l))
          crft(n,l)=crft(n,l)+crf(n,l)
c          print *,n,l,zc,aaa,rf(n,l)
          rfs(n,irec,l)=rf(n,l)
          drfs(n,irec,l)=drf(n,l)
c	  if(n.gt.100.and.n.lt.110) then
c	    print *,rfs(n,irec,l),drfs(n,irec,l)
c	  endif
        end do
      end do
c      pause
      sq2=sqrt(2.)
      do l=1,2
        do n=1,nf
c	      if(n.ge.nf/2) then
c                fac=cos(pi*(n-nf/2)/nf)**2
c	      else
c                fac=1.0
c	      endif
              fac=cos(pi*(n-1)/(2*nf))**2
          ar(n)=real(rf(n,l))*fac
          ar(npad+2-n)=real(rf(n,l))*fac
          ai(n)=-imag(rf(n,l))*fac
          ai(npad+2-n)=imag(rf(n,l))*fac
          pdat(n)=drf(n,l)*fac/sq2
          tdat(n)=fr(1)+(n-1)*fr(2)
        end do
        do n=nf+1,npad/2+2
          ar(n)=0.d0
          ar(npad+2-n)=0.d0
          ai(n)=0.d0
          ai(npad+2-n)=0.d0
        end do
        dtdat=(tdat(nf)-tdat(1))/30.
        tdata=tdat(1)-dtdat
        tdatb=tdat(nf)+dtdat
        bbmx=0.
        do n=1,npad
          bb(n)=dsqrt(ar(n)**2+ai(n)**2)
          bbmx=amax1(bbmx,bb(n))
        end do
        call plotit_axes(0.,0.,0.,0.)
        do n=1,npad
          if(bb(n).gt.1e-4) then
            pdat(n)=con*(pdat(n)/bb(n)) 
          else
            pdat(n)=360.
          endif
          bb(n)=con*datan2(ai(n),ar(n))
        end do
        call plotit_axes(0.,0.,0.,0.)
c  invert the Fourier transform
        call fft2(ar,ai,npad)
c  normalization factor:
c  divide by npad for fft2 routine normalization
c  mult by nnf/nf to compensate for losing high freqs
c  mult by 2 if cosine**2 taper is applied to spectral RF
        fac=2.*float(nnf)/(float(npad)*float(nf))
c        fac=(4./3.)*float(nnf)/(float(npad)*float(nf))
c        fac=float(nnf)/(float(npad)*float(nf))
        print *,'ddf,npad,nnf,nf,fac',ddf,npad,nnf,nf,fac
        do i=1,npost
          bb(npre+i)=ar(i)*fac
          time(npre+i)=(i-1)*dt
        end do
        do i=1,npre
          bb(npre+1-i)=ar(npad+1-i)*fac
          time(npre+1-i)=-i*dt
        end do
        ttt(1)=-npre*dt
        ttt(2)=dt
        ttt(3)=dt
        do i=1,ntot 
          rft(i,irec,l)=baz+50.*bb(i)
        end do
c        call plotit(ttt,bb,dum,400,'magnified version',
c     x    'time(sec)','H(t)',1,0,0.0,0,59+3*l)
        nmpt=40./dt
c        call plotit(ttt,bb,dum,ntot,title(l),
c     x    'time(sec)',xlabel3(l),1,0,0.0,0,60+3*l)
c        if(irec.eq.1) then
c        print *,irec,l
c	do n=101,109
c	  print *,rfs(n,irec,l),drfs(n,irec,l)
c	end do
c	pause
c	endif
      end do
 1076 format(3f10.4)
 1077 format(2f10.4)
 1022 format(f9.5,f6.1,f12.5)
  102 format(a)
      go to 10
  111 continue
      close(10)
c  plot the inferred angle between radial and vertical against epicentral dist
c      do i=1,irec
c        if(epic(i).lt.80.) then
c	  bb(i)=42.+(20.-epic(i))*(20./60.)
c	elseif(epic(i).lt.118.) then
c	  bb(i)=16.+(6./900.)*(110.-epic(i))**2
c	else
c	  bb(i)=9.0
c	endif
c      end do
c      call plotit(epic,bb,dum,irec,' ',' ',' ',2,0,0.0,0,0)
c      call plotit(epic,bbaz,dum,irec,'P-SV Rotation Check',
c     x 'Epicentral Distance','Inferred R/Z Angle',2,0,0.1,3,1)
c  average coherence as function of frequency
      do l=1,2
        do i=1,nf
          crft(i,l)=crft(i,l)/irec
        end do
      end do
      print *,irec
      call plotit_axes(tdata,tdatb,0.0,1.0)
      bb(1)=0.
      bb(2)=fmax
      bb(3)=1./float(nwin)
      bb(4)=bb(3)
      call plotit(bb,bb(3),dum,2,' ',' ',' ',2,0,0.0,0,0)
      call plotit(fr,crft(1,2),dum,nf,'Average P-Coda Coherence',
     x    'Freq(Hz)','Avg C\\sup{2}(f)',1,0,0.05,0,0)
      call plotit(fr,crft(1,1),dum,nf,'Average P-Coda Coherence',
     x    'Freq(Hz)','Avg C\\sup{2}(f)',1,0,0.0,0,1)
      print *,irec
      do i=1,nmpt
        time(i)=-ttt(1)-(i-1)*ttt(2)
      end do
      nnn=25./dt
      tmn=75*dt
      tmx=tmn-(50+nnn)*dt
C  KLUGE TO FOCUS THE PLOT
      tmx=tmx/2.
      tmn=tmn/2.
C  END KLUGE
      call plotit_axes(-10.,370.,tmx,tmn)
      irecd=1+irec/150  ! KLUGE TO AVOID TOO MANY POINTS
      iirec=((irec-1)/irecd)*irecd+1 
      do n=1,irec,irecd
        call plotit(rft(150,n,2),time(150),dum,nnn,
     x   'Transverse RF Section',
     x   'Back-Azimuth (degrees CW from N)',
     x  'time(sec)',2,0,0.0,0,21*(n/iirec))
      end do
      call plotit_axes(-10.,370.,tmx,tmn)
      do n=1,irec,irecd
        call plotit(rft(150,n,1),time(150),dum,nnn,
     x   'Radial RF Section',
     x   'Back-Azimuth (degrees CW from N)',
     x  'time(sec)',2,0,0.0,0,22*(n/iirec))
      end do
      call read_ttimes
      print *,'call read_ttimes'
c  here we begin the harmonic regression code.
c  the first section does a simple regression of the data RFs to 
c  obtain a "control" harmonic expansion, but does not write the result to disk.  
c  This control regression is not necessary 
c  if the bootstrap section is run, because that section computes a bootstrap average.
c  the section that follows line 737 is the bootstrap itself, 
c  which regresses for the harmonic expansion for nboot resamplings of the data set
c
      if(nboot.eq.0) then  ! skip over the regression check      
c  the output of this program is a least-square regression of RFs in the freq
c  domain, with complex-valued coefficients for the 
c  constant,cos(baz),sin(baz),cos(2*baz),sin(2*baz) variation
c
c  test inversion to verify the concept
      print *,irec,' seismic records'
      print *,nf,' frequencies'
      do kk=1,nlp
        do l=1,2
	  do n=1,nf
            do k=1,5
              do j=1,5
                bootg(j,k,n,l,kk)=0.
	      end do
     	      bootd(k,n,l,kk)=zero
	    end do
	  end do
	end do
      end do
      print *,' boot matrices zeroed'
      do i=1,irec
        print *,' processing record',i
c  get p-slowness -- iflag branches on the start time of the analysis window
        igcarc=gcarc(i)
        fac=gcarc(i)-igcarc
c        print *,'igcarc,i,gcarc(i),fac',igcarc,i,gcarc(i),fac
        if(iflag(i).eq.1) then		! P-phase
          if(igcarc.ge.130.or.igcarc.le.0) then
            print *,'WHOA! igcarc,gcarc = ',igcarc,gcarc(i)
c            pause
          endif
          slow=(1.-fac)*s_p(igcarc)+fac*s_p(igcarc+1)
        elseif(iflag(i).eq.2) then	! PKP-phase
          if(igcarc.ge.180.or.igcarc.lt.90) then
            print *,'WHOA! igcarc,gcarc = ',igcarc,gcarc(i)
c            pause
          endif
          ig=igcarc-89
          slow=(1.-fac)*s_pkp(ig)+fac*s_pkp(ig+1)
        else
          print *,'WHOA! iflag(i),i = ',iflag(i),i
c          pause
        endif
c  the stack-model velocities are in SI units, s_p is in sec/km
        slw=(slow/1000.)       
c  compute the slowness-dependent time delays at the interfaces of stack-model
        do k=1,nlp
          coss=sqrt(1.-vs(k)*vs(k)*slw*slw)
          cosp=sqrt(1.-vp(k)*vp(k)*slw*slw)
c  trapdoor to avoid incoming slownesses that are evanescent (e.g. head waves)
c  basically, one stacks with the last acceptable interval velocity
          if(vp(k)*slw.lt.1.) then
            gam(k)=(vp(k)*coss-vs(k)*cosp)/(vp(k)-vs(k))
          else
            gam(k)=gam(k-1)
          endif
        end do
        psd(1)=gam(1)*tau(1)
        xik(1)=0.
        do k=2,nlp
          psd(k)=psd(k-1)+gam(k)*(tau(k)-tau(k-1))
          xik(k)=psd(k-1)/gam(k)-tau(k-1)
        end do
c  note that we boost the constant term here to boost the tradeoff
c  between const and baz-dependent terms.
c  we bias the regression to a constant solution
c  must multiply the constant coefficient by identical factor	  
	bf(1)=3.0
        bf(2)=cos(bazs(i)/con)
        bf(3)=sin(bazs(i)/con)
	bf(4)=cos(2.*bazs(i)/con)
	bf(5)=sin(2.*bazs(i)/con)
c  spline-interpolate the RF in freq domain
        do l=1,2
          do n=1,nf
            sss(n,1)=real(rfs(n,i,l))
            sss(n,2)=imag(rfs(n,i,l))
            sss(n,3)=drfs(n,i,l)
	    sss(n,4)=0.
	    sss(n,5)=0.
	    sss(n,6)=0.
c	    if(n.gt.100.and.n.lt.110) then
c	      print *,(sss(n,j),j=1,3)
c	    endif
          end do
          call splneq(nf,sss(1,1),sss(1,4))
          call splneq(nf,sss(1,2),sss(1,5))
          call splneq(nf,sss(1,3),sss(1,6))
c	print *,i,(bf(k),k=1,5)
c	print *,(psd(kk),kk=1,nlp)
c	print *,(tau(kk),kk=1,nlp)
c	print *,(gam(kk),kk=1,nlp)
c	print *,(xik(kk),kk=1,nlp)
          do kk=1,nlp  ! loop over layers in the stacking model
            dfm=ddf/gam(kk)
            dxm=1.0/gam(kk)
            dcs=cos(2.*pi*dfm*xik(kk)) 
            dsn=sin(2.*pi*dfm*xik(kk)) 
c  we start at f=0., so that the migration phase factor = 1
c  csm and snm are real/imag parts of the migration phase factor
c  and the position in the spline is x=1  (fm is "migration frequency")
            csm=1.
            snm=0.
            xm=1.
	    do n=1,nf
              sig=evaleq(xm,nf,sss(1,3),sss(1,6),0,1.)
	      sig2=sig*sig
              dre=evaleq(xm,nf,sss(1,1),sss(1,4),0,1.)
              dim=evaleq(xm,nf,sss(1,2),sss(1,5),0,1.)
              dre=dre/(sig*gam(kk))
              dim=dim/(sig*gam(kk))
	      zc=cmplx(dre*csm+dim*snm,-dre*snm+dim*csm)
              csm=csm*dcs-snm*dsn
              snm=snm*dcs+csm*dsn
              xm=xm+dxm	      
	      do j=1,5
	        do k=1,5	     
	           bootg(k,j,n,l,kk)=bootg(k,j,n,l,kk)+bf(k)*bf(j)/sig2
   	        end do
		bootd(j,n,l,kk)=bootd(j,n,l,kk)+bf(j)*zc/sig
	      end do
c	      if(n.gt.100.and.n.lt.110) then
c	        print *,'data vector'
c	        print *,(bootd(j,n,l,kk),j=1,5)
c	        print *,'diagonal of gram matrix'
c	        print *,(bootg(1,j,n,l,kk),j=1,5)
c              endif
	    end do
c	    pause
	  end do
	end do
      end do
c  we regress the real and imag parts of the RF spectra separately
c  because the gram matrix is real-valued, the real and imag separate        	      
      do l=1,2
        do kk=1,nlp  ! loop over layers in the stacking model
	  do n=1,nf
	    do j=1,5
	      do k=1,5
	        bootgg(k,j)=bootg(k,j,n,l,kk)
	      end do
	      bootdd(j)=real(bootd(j,n,l,kk))
	    end do
c	    if(n.gt.100.and.n.lt.110) then
c	      print *,'data vector'
c	      print *,(bootd(j,n,l,kk),j=1,5)
c	      print *,'diagonal of gram matrix'
c	      print *,(bootg(j,j,n,l,kk),j=1,5)
c	    endif
	    call solve(5,bootgg,bootmm,bootdd)
	    do j=1,5
	      bootdd(j)=imag(bootd(j,n,l,kk))
	    end do
	    call solve(-5,bootgg,bootmm(1,2),bootdd)
	    do j=1,5
 	  bootma(j,n,l,kk)=cmplx(sngl(bootmm(j,1)),sngl(bootmm(j,2)))
	    end do
c  correct for bias factor in constant term	    
	    bootma(1,n,l,kk)=bootma(1,n,l,kk)*3.0
c	    if(n.gt.100.and.n.lt.110) then
c  	      print *,n,l,kk
c	      print *,(bootma(j,n,l,kk),j=1,5)
c	    endif
	  end do
c	  pause
	end do
      end do
	
cccccccccccccccccccccccccccccc   unmodified	
c      print *,'model parameters'
c      do j=1,5
c	do l=1,2
c	  do n=1,nf
c	    bb(n)=cabs(bootma(j,n,l,1))
c	  end do
c          call plotit_axes(0.,0.,0.,0.)
c	  call plotit(fr,bb,dum,nf,'regression term','freq(Hz)',
c     x       xlabel2(l),1,0,0.0,0,39+l*2)
c          do n=1,nf
c        bb(n)=con*atan2(imag(bootma(j,n,l,1)),real(bootma(j,n,l,1)))
c	  end do
c          call plotit_axes(0.,0.,0.,0.)
c	  call plotit(fr,bb,dum,nf,'regression term','freq(Hz)',
c     x       xlabel(l),1,0,0.1,2,40+l*2)
c        end do
c      end do
cccccccccccccccccccccccccccccc   unmodified	
c  the uncertainty on the RF expansion wont be computed or saved if nboot=0
      else      
  737 print *,'time for the bootstrap'
      print *,'nf,npad,nlp',nf,npad,nlp
c  the bootstrap is to obtain uncertainties, AND cross correlations
c  we regress nboot times -- resampling the records in the data set
c  uncertainties in the data set are accounted for in the regression
      do ib=1,nboot
        do kk=1,nlp
          do l=1,2
	    do n=1,nf
              do k=1,5
                do j=1,5
                  bootg(j,k,n,l,kk)=0.
	        end do
     	        bootd(k,n,l,kk)=zero
	      end do
	    end do
	  end do
        end do
	do i=1,irec
c          print *,' step',i
	  ii=1+rand(0)*irec*0.999999
c          print *,' processing record',ii
c  get p-slowness -- iflag branches on the start time of the analysis window
          igcarc=gcarc(ii)
          fac=gcarc(ii)-igcarc
c          print *,'igcarc,i,gcarc(i),fac',igcarc,i,gcarc(i),fac
          if(iflag(ii).eq.1) then		! P-phase
            if(igcarc.ge.130.or.igcarc.le.0) then
              print *,'WHOA! igcarc,gcarc = ',igcarc,gcarc(ii)
c            pause
            endif
            slow=(1.-fac)*s_p(igcarc)+fac*s_p(igcarc+1)
          elseif(iflag(ii).eq.2) then	! PKP-phase
            if(igcarc.ge.180.or.igcarc.lt.90) then
              print *,'WHOA! igcarc,gcarc = ',igcarc,gcarc(ii)
c              pause
            endif
            ig=igcarc-89
            slow=(1.-fac)*s_pkp(ig)+fac*s_pkp(ig+1)
          else
            print *,'WHOA! iflag(i),i = ',iflag(ii),ii
c          pause
          endif
c  the stack-model velocities are in SI units, s_p is in sec/km
          slw=(slow/1000.)       
c  compute the slowness-dependent time delays at the interfaces of stack-model
          do k=1,nlp
            coss=sqrt(1.-vs(k)*vs(k)*slw*slw)
            cosp=sqrt(1.-vp(k)*vp(k)*slw*slw)
c  trapdoor to avoid incoming slownesses that are evanescent (e.g. head waves)
c  basically, one stacks with the last acceptable interval velocity
            if(vp(k)*slw.lt.1.) then
              gam(k)=(vp(k)*coss-vs(k)*cosp)/(vp(k)-vs(k))
            else
              gam(k)=gam(k-1)
            endif
          end do
          psd(1)=gam(1)*tau(1)
          xik(1)=0.
          do k=2,nlp
            psd(k)=psd(k-1)+gam(k)*(tau(k)-tau(k-1))
            xik(k)=psd(k-1)/gam(k)-tau(k-1)
          end do
c  note that we boost the constant term here to boost the tradeoff
c  between const and baz-dependent terms.
c  we bias the regression to a constant solution
c  must multiply the constant coefficient by identical factor	  
	  bf(1)=3.0
          bf(2)=cos(bazs(ii)/con)
          bf(3)=sin(bazs(ii)/con)
	  bf(4)=cos(2.*bazs(ii)/con)
	  bf(5)=sin(2.*bazs(ii)/con)
          do l=1,2
c  spline-interpolate the RF in freq domain
            do n=1,nf
              sss(n,1)=real(rfs(n,ii,l))
              sss(n,2)=imag(rfs(n,ii,l))
              sss(n,3)=drfs(n,ii,l)
	      sss(n,4)=0.
	      sss(n,5)=0.
	      sss(n,6)=0.
c	      if(n.gt.100.and.n.lt.110) then
c	        print *,(sss(n,j),j=1,3)
c	      endif
            end do
c	    pause
            call splneq(nf,sss(1,1),sss(1,4))
            call splneq(nf,sss(1,2),sss(1,5))
            call splneq(nf,sss(1,3),sss(1,6))
            do kk=1,nlp  ! loop over layers in the stacking model
              dfm=ddf/gam(kk)
              dxm=1.0/gam(kk)
              dcs=cos(2.*pi*dfm*xik(kk)) 
              dsn=sin(2.*pi*dfm*xik(kk)) 
c  we start at f=0., so that the migration phase factor = 1
c  csm and snm are real/imag parts of the migration phase factor
c  and the position in the spline is x=1  (fm is "migration frequency")
              csm=1.
              snm=0.
              xm=1.
	      do n=1,nf
                sig=evaleq(xm,nf,sss(1,3),sss(1,6),0,1.)
	        sig2=sig*sig
                dre=evaleq(xm,nf,sss(1,1),sss(1,4),0,1.)
                dim=evaleq(xm,nf,sss(1,2),sss(1,5),0,1.)
                dre=dre/(sig*gam(kk))
                dim=dim/(sig*gam(kk))
	        zc=cmplx(dre*csm+dim*snm,-dre*snm+dim*csm)
                csm=csm*dcs-snm*dsn
                snm=snm*dcs+csm*dsn
                xm=xm+dxm	      
	        do j=1,5
	          do k=1,5	     
	        bootg(k,j,n,l,kk)=bootg(k,j,n,l,kk)+bf(k)*bf(j)/sig2
   	          end do
		  bootd(j,n,l,kk)=bootd(j,n,l,kk)+bf(j)*zc/sig
	        end do
	      end do
	    end do
	  end do
	end do
c  we regress the real and imag parts of the RF spectra separately
c  because the gram matrix is real-valued, the real and imag separate        	      
        do l=1,2
          do kk=1,nlp  ! loop over layers in the stacking model
	    do n=1,nf
	      do j=1,5
	        do k=1,5
	          bootgg(k,j)=bootg(k,j,n,l,kk)
	        end do
	        bootdd(j)=real(bootd(j,n,l,kk))
	      end do
	      call solve(5,bootgg,bootmm,bootdd)
	      do j=1,5
	        bootdd(j)=imag(bootd(j,n,l,kk))
	      end do
	      call solve(-5,bootgg,bootmm(1,2),bootdd)
	      do j=1,5
      bootm(ib,j,n,l,kk)=cmplx(sngl(bootmm(j,1)),sngl(bootmm(j,2)))
	      end do
c  correct for bias factor in constant term	    
	      bootm(ib,1,n,l,kk)=bootm(ib,1,n,l,kk)*3.0
	    end do
	  end do
	end do
	if(ib.eq.(ib/10)*10) print *,'iteration',ib
      end do

      print *,'compute bootstrap means and covariance matrices'
      do l=1,2
        do kk=1,nlp  ! loop over layers in the stacking model
	do n=1,nf
	  do j=1,5
	    bootma(j,n,l,kk)=zero
	    do k=1,5
	      bootmv(k,j,n,l,kk)=zero 
	    end do
            do ib=1,nboot      	    
	      bootma(j,n,l,kk)=bootma(j,n,l,kk)+bootm(ib,j,n,l,kk)  
	    end do
	    bootma(j,n,l,kk)=bootma(j,n,l,kk)/nboot  
	  end do
	  do ib=1,nboot
  	    do j=1,5
	      cdev(j)=bootm(ib,j,n,l,kk)-bootma(j,n,l,kk)
	    end do
c  conjugation conventions can be dicey -- check this before you use it
	    do j=1,5	  
	      do k=j,5
	  bootmv(k,j,n,l,kk)=bootmv(k,j,n,l,kk)+conjg(cdev(k))*cdev(j)
	      end do
	    end do
	  end do
c	  if(kk.eq.1) then
c	    print *,(bootmv(j,j,n,l,kk),j=1,5)
c	    pause
c          end if
	  do j=1,4
	    do k=j+1,5
	      bootmv(j,k,n,l,kk)=conjg(bootmv(k,j,n,l,kk))
	    end do
	  end do
	end do
	end do
      end do
c  note that here the two 737 and 747 branches converge.  
c  The average RF expansion is computed one way or the other, 
c  now we compute the predicted BAZ dependence of the RF and then the uncertainties
c  if the bootstrap was taken  --- the choice if nboot.gt.0 -- 
  747 endif      
      open(12,file='outr_bbaz.grid',form='formatted')
      open(13,file='outt_bbaz.grid',form='formatted')
      bb(1)=0.
      bb(2)=fmax
      bb(3)=1./3.
      bb(4)=bb(3)
c  there are ntot points available in the raw
      nss=7./dt
      npss=18./dt
      kaz=0
      gc1=0.
      gc2=180.
      do baz=5,360,10
c      do baz=150,265,5
        do l=1,2
	  bf(1)=1.
	  bf(2)=cos(baz/con)
	  bf(3)=sin(baz/con)
	  bf(4)=cos(2.*baz/con)
	  bf(5)=sin(2.*baz/con)
          kaz=kaz+1
	  do kk=1,nlp
            do n=1,nf
              rf(n,l)=zero
	      do j=1,5
	        rf(n,l)=rf(n,l)+bf(j)*bootma(j,n,l,kk)
              end do
            end do    
c          print *,baz
c  do the expectation of variance for the weighted sum yourself
c  if you doubt the formula here -- basically the weighted terms in rf-sum
c  have unit variance, so that variance of the total variance is jrec/drf**2
            do n=1,nf
c	      if(n.ge.nf/2) then
c                fac=cos(pi*(n-nf/2)/nf)**2
c	      else
c                fac=1.0
c	      endif
              fac=cos(pi*(n-1)/(2*nf))**2
              ar(n)=fac*real(rf(n,l))
              ai(n)=-fac*imag(rf(n,l))
              ar(npad+2-n)=fac*real(rf(n,l))
              ai(npad+2-n)=fac*imag(rf(n,l))
            end do
            do n=nf+1,npad-nf+1
              ar(n)=0.d0
              ai(n)=0.d0
            end do
            call fft2(ar,ai,npad)
c  normalization factor:
c  divide by npad for fft2 routine normalization
c  mult by nnf/nf to compensate for losing high freqs
c  mult by 2 if cosine**2 taper is applied to spectral RF
c            fac=(4./3.)*float(nnf)/(float(npad)*float(nf))
            fac=2.*float(nnf)/(float(npad)*float(nf))
            do i=1,npost
              bb(npre+i)=ar(i)*fac
            end do
            do i=1,npre
              bb(npre+1-i)=ar(npad+1-i)*fac
            end do
            kz=(kaz-l)/2+1
            do i=1,ntot 
              rft(i,kz+kk,l)=bb(i)
            end do
	  end do
          do i=1,ntot
            rft(i,kz,l)=0.
          end do         
          i1=1
          do kk=1,nlp
            i2=min(ntot,200+itau(kk))
            print *,i1,i2
            do i=i1,i2 
              rft(i,kz,l)=rft(i,kz+kk,l)
            end do
            i1=npre+1+itau(kk)
          end do
          iunit=11+l
          do i=1,ntot 
            t3=-time(i)
            write(iunit,1022) t3,baz,rft(i,kz,l)
c for PLOTIT display, not yet in this  code
	    rft(i,kz,l)=baz+50*rft(i,kz,l)*l 
          end do
          write(iunit,101) '>'
        end do
      end do  
      kaz=kaz/2
      print *,kaz,' traces'   
      close(12)
      close(13)
c  now write out the five traces for the truncated Fourier expansion in BAZ
      open(12,file='outr_bexp.grid',form='formatted')
      open(13,file='outt_bexp.grid',form='formatted')
      open(15,file='outr1_bexp.grid',form='formatted')
      open(16,file='outt1_bexp.grid',form='formatted')
      open(17,file='outr2_bexp.grid',form='formatted')
      open(18,file='outt2_bexp.grid',form='formatted')
      kaz=0
      do j=1,5
        space=10.*(5-j)
        do l=1,2
          kaz=kaz+1
          do kk=1,nlp
            do n=1,nf
	      rf(n,l)=bootma(j,n,l,kk)
            end do    
c          print *,baz
c  do the expectation of variance for the weighted sum yourself
c  if you doubt the formula here -- basically the weighted terms in rf-sum
c  have unit variance, so that variance of the total variance is jrec/drf**2
            do n=1,nf
c	      if(n.ge.nf/2) then
c                fac=cos(pi*(n-nf/2)/nf)**2
c	      else
c                fac=1.0
c	      endif
              fac=cos(pi*(n-1)/(2*nf))**2
              ar(n)=fac*real(rf(n,l))
              ai(n)=-fac*imag(rf(n,l))
              ar(npad+2-n)=fac*real(rf(n,l))
              ai(npad+2-n)=fac*imag(rf(n,l))
            end do
            do n=nf+1,npad-nf+1
              ar(n)=0.d0
              ai(n)=0.d0
            end do
            call fft2(ar,ai,npad)
c  normalization factor:
c  divide by npad for fft2 routine normalization
c  mult by nnf/nf to compensate for losing high freqs
c  mult by 2 if cosine**2 taper is applied to spectral RF
c            fac=(4./3.)*float(nnf)/(float(npad)*float(nf))
            fac=2.*float(nnf)/(float(npad)*float(nf))
            do i=1,npost
              bb(npre+i)=ar(i)*fac
            end do
            do i=1,npre
              bb(npre+1-i)=ar(npad+1-i)*fac
            end do
            kz=(kaz-l)/2+1
            iunit=11+l
            do i=1,ntot 
              rft(i,j+kk,l)=bb(i)
            end do
	  end do
          i1=1
          do kk=1,nlp
            i2=min(ntot,200+itau(kk))
            print *,i1,i2
            do i=i1,i2 
c              rft(i,j,l)=j-1+5.0*rft(i,j+kk,l)
              rft(i,j,l)=rft(i,j+kk,l)
              bb(i)=rft(i,j+kk,l)
            end do
            i1=npre+1+itau(kk)
          end do
          iunit=11+l
          do i=1,ntot 
            t3=-time(i)
            write(iunit,1022) t3,space,bb(i)
          end do
          write(iunit,101) '>'
c  now we compute bootstrap uncertainty in the time domain
c  this involves doing the inverse FFTs of the bootstrapped RFs
c  we have saved the bootstrap RFs in the frequency domain because there are fewer points in FFT
          if(nboot.gt.0) then  ! the code will likely crash if nboot=1, so dont try it!
	    do i=1,ntot
	      crft(i,l)=0.0
	    end do
	    do ib=1,nboot
              do kk=1,nlp
                do n=1,nf
	          rf(n,l)=bootm(ib,j,n,l,kk)
                end do    
c          print *,baz
c  do the expectation of variance for the weighted sum yourself
c  if you doubt the formula here -- basically the weighted terms in rf-sum
c  have unit variance, so that variance of the total variance is jrec/drf**2
                do n=1,nf
c	      if(n.ge.nf/2) then
c                fac=cos(pi*(n-nf/2)/nf)**2
c	      else
c                fac=1.0
c	      endif
              fac=cos(pi*(n-1)/(2*nf))**2
                  ar(n)=fac*real(rf(n,l))
                  ai(n)=-fac*imag(rf(n,l))
                  ar(npad+2-n)=fac*real(rf(n,l))
                  ai(npad+2-n)=fac*imag(rf(n,l))
                end do
                do n=nf+1,npad-nf+1
                  ar(n)=0.d0
                  ai(n)=0.d0
                end do
                call fft2(ar,ai,npad)
c  normalization factor:
c  divide by npad for fft2 routine normalization
c  mult by nnf/nf to compensate for losing high freqs
c  mult by 2 if cosine**2 taper is applied to spectral RF
c                fac=(4./3.)*float(nnf)/(float(npad)*float(nf))
                fac=2.*float(nnf)/(float(npad)*float(nf))
                do i=1,npost
                  bb(npre+i)=ar(i)*fac
                end do
                do i=1,npre
                  bb(npre+1-i)=ar(npad+1-i)*fac
                end do
                kz=(kaz-l)/2+1
                iunit=11+l
                do i=1,ntot 
                  rft(i,5+ib+kk,l)=bb(i)
                end do
	      end do
              i1=1
              do kk=1,nlp
                i2=min(ntot,200+itau(kk))
                print *,i1,i2
                do i=i1,i2 
                  rft(i,5+ib,l)=rft(i,5+ib+kk,l)
                end do
                i1=npre+1+itau(kk)
              end do
	      do i=1,ntot
	        crft(i,l)=crft(i,l)+(rft(i,5+ib,l)-rft(i,j,l))**2
	      end do
	      print *,'iteration',ib
            end do
c  note that the bootstrap standard deviation is the sqrt of the average summed variance of
c  the bootstrapped estimates about the bootstrap mean. Therefore the expectation of the 
c  bootstrap standard deviation does not scale up or down with nboot.
	    do i=1,ntot
	      bb(i)=sqrt(crft(i,l)/(nboot-1))
	      crft(i,1)=rft(i,j,l)+bb(i)
	      crft(i,2)=rft(i,j,l)-bb(i)
	    end do
            iunit=14+l
	    do i=1,ntot
              t3=-time(i)
              write(iunit,1022) t3,space,crft(i,1)
	    end do
            write(iunit,101) '>'
            iunit=16+l
	    do i=1,ntot
              t3=-time(i)
              write(iunit,1022) t3,space,crft(i,2)
            end do
            write(iunit,101) '>'
	  endif
        end do
      end do  
      close(12)
      close(13)
      close(15)
      close(16)
      close(17)
      close(18)
      kaz=kaz/2
      print *,kaz,' traces'   
      do iagain=1,2
      call plotit_axes(0.,0.,0.,0.)
      do n=1,kaz
        call plotit(rft(150,n,2),time(150),dum,nnn,
     x   'Composite Transverse RF Section',
     x   'Back-Azimuth (degrees CW from N)',
     x  'time(sec)',2,0,0.0,0,21*(n/kaz))
      end do
            
      call plotit_axes(0.,0.,0.,0.)
      do n=1,kaz
        call plotit(rft(150,n,1),time(150),dum,nnn,
     x   'Composite Radial RF Section',
     x   'Back-Azimuth (degrees CW from N)',
     x  'time(sec)',2,0,0.0,0,22*(n/kaz))
      end do
      end do

      stop
      end
c
      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
      common/nnaamm/iwflag
      common/npiprol/anpi
      common/tapsum/tapsum(20),ell(20)
      common/work/ip(16400)        
      common/taperzz/z(65536)  ! we use this common block for scratch space
      common/stap2/ta(16400,8)
      dimension a(16400,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      print *,'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      print *,ierr,rlb,rlu
      print *,'eigenvalues for the eigentapers'
c      print *,(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      print *,'eigenvalues for the eigentapers (small k)'
c      print *,(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
      do i=1,nwin
        el(i)=dmax1(ell(i),el(i))
      end do     
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      print *,'tapsum',(tapsum(i),i=1,nwin)
  101 format(80a)
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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c               routine
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine fft2(ar,ai,n)
c  fft routine with 2 real input arrays rather than complex - see ffttwo
c  for comments. 
c  fft2 subroutine mults by exp(i\omega t)
c  OUTPUT: f=0 in ar(1), f=f_N in ar(n/2+1)
c          ar(i)=ar(n+2-i), ai(i)=-ai(n+2-i)
c  fft2 is NOT a unitary transform, mults the series by sqrt(n)
c  the inverse FT can be effected by running fft2 on the conjugate of
c  the FFT-expansion, then taking the the conjugate of the output, 
c  and dividing thru by N. to wit:
c
c   assume Xr, Xi is in freq domain, xr, xi in the time domain
c
c   (Xr,Xi)=fft2(xr,xi,N)
c   (xr,-xi)=fft2(Xr,-Xi,N)/N
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      dimension ar(1),ai(1)
      mex=dlog(dble(float(n)))/.693147d0
      nv2=n/2
      nm1=n-1
      j=1
      do 7 i=1,nm1
      if(i .ge. j) go to 5
      tr=ar(j)
      ti=ai(j)
      ar(j)=ar(i)
      ai(j)=ai(i)
      ar(i)=tr
      ai(i)=ti
   5  k=nv2
   6  if(k .ge. j) go to 7
      j=j-k
      k=k/2
      go to 6
   7  j=j+k
      pi=3.14159265358979d0
      do 20 l=1,mex
      le=2**l
      le1=le/2
      ur=1.0
      ui=0.
      wr=dcos(pi/le1 )
      wi=dsin (pi/le1)
      do 20 j=1,le1
      do 10 i=j,n,le
      ip=i+le1
      tr=ar(ip)*ur - ai(ip)*ui
      ti=ai(ip)*ur + ar(ip)*ui
      ar(ip)=ar(i)-tr
      ai(ip)=ai(i) - ti
      ar(i)=ar(i)+tr
      ai(i)=ai(i)+ti
  10  continue
      utemp=ur
      ur=ur*wr - ui*wi
      ui=ui*wr + utemp*wi
  20  continue
      return
      end
      subroutine read_ttimes
      implicit real*4(a-h,o-z)
      implicit integer*4(i-n)
      common/tt_p/d_p(200),d_pkp(200),s_p(200),s_pkp(200),n_p,n_pkp
c  read in the P- and PKP-slowness values from ttimes
c  the code interpolates across the upper-mantle triplication for DELTA=19-27
c  and chooses P or PKP based on the start time of the data-analysis window
c  for DELTA.ge.105
      do i=1,200
        d_p(i)=0.
      end do
c  fill in the head-wave slowness for DELTA.lt.14
c  13.2 sec/degree is 8.4+eps km/sec
      do i=1,13
        d_p(i)=float(i)
        s_p(i)=13.2
      end do
c      open(7,file='/park/backup/Vary/Layers/ttimes_P.dat',
      open(7,file='/Users/jjpark/Research/Layers/ttimes_P.dat',
     x  form='formatted')
      do i=14,200
        read(7,*,end=110) d_p(i),s_p(i)
      end do
  110 close(7)
      print *,i,' = counter at end of ttimes file'
      do i=1,200
        if(d_p(i).eq.0.) go to 120
      end do
  120 n_p=i-1
c  interpolates the triplication
c  linearly interpolate slowness between Delta=17 and p=11.06 sec/degree
c                                    and Delta=28 and p=8.92 sec/degree
      do i=18,27
        d_p(i)=float(i)
        s_p(i)=11.06-(i-17)*(11.06-8.92)/11.
      end do
      do i=45,n_p
        d_p(i-18)=d_p(i)
        s_p(i-18)=s_p(i)
      end do
      n_p=n_p-18
      do i=1,200
        d_pkp(i)=0.
      end do
c      open(7,file='/park/backup/Vary/Layers/ttimes_PKP.dat',
      open(7,file='/Users/jjpark/Research/Layers/ttimes_PKP.dat',
     x  form='formatted')
      do i=1,200
        read(7,*,end=130) d_pkp(i),s_pkp(i)
      end do
  130 close(7)
      print *,i,' = counter at end of ttimes file'
      do i=1,200
        if(d_pkp(i).eq.0.) go to 140
      end do
  140 n_pkp=i-1
      factor=180./(3.14159265358979*6371.)
      do i=1,n_p
        s_p(i)=s_p(i)*factor
      end do
      do i=1,n_pkp
        s_pkp(i)=s_pkp(i)*factor
      end do
      print *,'n_p,n_pkp',n_p,n_pkp
      call plotit_axes(0.,0.,0.,0.)
      call plotit(d_p,s_p,dum,n_p,' ',' ',' ',2,0,0.0,0,0)
      call plotit(d_pkp,s_pkp,dum,n_pkp,' ',
     x 'Epicentral Distance (degrees)','Slowness (sec/km)',2,0,0.0,0,1)
      return
      end