c  program recfunk_estack
c  5/11/01 JJP -->  REV2 5/16/01
c  updated 08/17/07 JJP
c  xf77 -o /park/backup/bin/recfunk_estack recfunk_estack.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c  xf77 -o /Users/jjpark/bin/recfunk_estack recfunk_estack.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c  will compute RFs using single eigenspectrum estimate per record
c  instead of cross-correlating eigenspectra in one event
c  --> cross-correlate across events
c  steps:
c  LOOP1
c  read the filenames to analyse, compute 1-pi prolate spectra 
c					and coherence estimates
c  save baz/epicenter info
c  LOOP2
c  choose criteria for estack
c  for each group of events: loop over freq:
c			assemble the spectral vectors, cross-correlate
c   				with coherence weighting
c  compute the H(f), display
c
c
c  BASED ON
c  program recfunk_noplot
c  5/25/99 JJP
c  multitaper program to invert for receiver functions
c  SINGLE STATION WITH SPECTRAL WEIGHTING
c  
c  "Multiple Taper Correlation" (MTC) Receiver functions are described in 
c
c  Park, J., and V. Levin, Receiver functions from multiple-taper spectral 
c  correlation estimates, Bull. Seismo. Soc Amer., v90, pp1507-1520,  2000.
c  
c  spectral estimates of pre-event noise 
c  weight the inversion in the frequency domain, thereby avoiding the 
c  usual problem with microseism noise in 0.1-0.5 Hz
c
c  the main advantage of MTC is not the noise damping, however -- the
c  cross-correlation of the horizontal and vertical components is used
c  in place of spectral division when estimating the RF.
c  this approach appears to be more stable and gives RF-uncertainties
c  in the frequency domain that are useful for stacking.
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  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 program computes RFs in the Freq domain for multiple records
c   then it computes composite RFs that are dependent on either 
c   Back-azimuth or epicentral distance.  The standard bin width is 10 degrees
c   with stacked RFs computed every 5 degrees.  You can change this in the code.
c  RF are not migrated in this process, so the user is prompted for parameters 
c  that narrow the focus of the stacks 
c  e.g. a BAZ-sector over which to compute an epicentral RF profile.
c
c  the PLOTIT graphs of RF sweeps can be improved upon, natch
c  the code writes the BAZ- and EPICEN-dependent RFs to files
c  in a format easily digested by GMT (traces are separated by ">" lines)
c
c  filenames: out[rt]_baz.grid out[rt]_epi.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)
      real*8 ar,ai,el(12),aaa
      complex*8 zc,zero
      complex*8 afft,rf,rff
      character*80 name,subtext,string
      character*28 title(2)
      character*10 cmps(3)
      character comp(3)
      common/nnaamm/iwflag
      common/npiprol/anpi
      common/stap2/tap(16400,16)
      common/taperz/ta(16400,2,3),tai(16400,2,3)
      common/stap0/ar(16400),ai(16400),bb(16400),time(99000)
      common/specs/afft(4100,799,3),cft(4100,799,3)
      common/specs2/rft(4100,100,2),crft(4100,2),rf(4100,2),rff(4100,2)
      common/saveit/bazs(799),epic(799),bbaz(799)
      common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2)
      common/header/ahead(158)
      common/distrib/xi(72,36),yi(72,36),sumij(72,36)
      dimension iah(158),chead(158),mmonth(12),fr(3),tt(3),ttt(3)
      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'/
      con=180./3.14159265358979
      pi=3.14159265358979
      pih=pi/2.
      sq2=sqrt(2.)
      zero=cmplx(0.,0.)
      fmax=5.
      print *,'enter fmax (Hz)'
      read(5,*) fmax
c generate tapers for npts-point records, zero-pad to npad=4096
      itt=16384
      npad=16384
      nnf=npad/2
      do i=1,4100
        crft(i,1)=0.
        crft(i,2)=0.
      end do
      irec=0
 1010 print *,'time intervals in file (0) or in SAC-header (1)'
      read(5,*) isegment
      if(isegment.eq.0) then
        print *,
     x'file of seismic records to read? (space-ret: in_recfunk_estack)'
        read(5,101) name
        if(name(1:1).eq.' ') then
          open(10,file='in_recfunk_estack',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_estack)'
        read(5,101) name
        if(name(1:1).eq.' ') then
          open(10,file='in_recpick_estack',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 *,'minimum number of events in a binsum (eg 1 or 2)'
      read(5,*) ibskip
      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)
	if(ierr.eq.1) then
	  print *,'problem reading file', name
	  pause
	endif
        dt=ahead(1) 
        tt(1)=0.
        tt(2)=dt
        tt(3)=dt
        nscan=iah(80)
        if(nscan.gt.99000) then
          print *,'Careful! data length greater than 99000:', nscan
          pause
        endif
      end do
      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 *,'sach params b and tmark, pstart',ahead(6),tmark,pstart
      endif
c  naturally, you are in deep poop if dt is not the same for all records
      npre=10./dt
      npost=85./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
c  here we set a trap for deep earthquakes - if too deep, skip the record
      if(ahead(39).gt.50000) go to 10
      irec=irec+1
      if(irec.gt.799) then
        print *,'too many records! 799 MAX'
	stop
      endif
      baz=ahead(53)
      bazs(irec)=ahead(53)
      epic(irec)=ahead(54)
      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
      if(nf.gt.4100) then
        print *,'nf is too big', nf
	stop
      endif
      fr(1)=0.
      fr(2)=ddf
      fr(3)=ddf
      npts=postwin/dt
      nst=pstart/dt
      anpi=1.
      nwin=1
      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
      if(lqt.ne.0) then
c  here is the option to rotate the Radial and vertical Components to
c  a P and SV coordinate system    
c  preferred option -- use parameterized relation between GCARC and zrrot
c  compute approximate angle between Z and R components for P wave
        epipi=epic(irec)
	if(imark.eq.3) epipi=epipi/2.
        if(imark.eq.2) then
	  zrrot=9.0
	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
	  else
            zrrot=16.+(6./900.)*(110.-epipi)**2
	  endif
        endif
        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 P coda
        do k=1,nwin
          do i=1,npts
            ar(i)=a(i+nst,kcmp)*tap(i,k)
c            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 n=1,nf
            afft(n,irec,kcmp)=dcmplx(ar(n),ai(n))
          end do
c          call plotit(fr,bb,dum,nf,'fft',
c     x     'frequency (Hz)','H(f)',1,0,0.0,0,1)
        end do
      end do
c  compute coherency estimates for regression weighting
c  to do this without fudge factors, recompute the spectra for 3 npi=2.5 tapers
      anpi=2.5
      nwin=3
      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  first part of calculation: cft accumulates hi-res spectra of 3 comps
c			       rf accumulates unnormalized correlation
c    2nd part of calculation: cft(*,*,[12]) has [RT] coherence
c			      cft(*,*,3) has total spectrum for weighting
      do j=1,3
        do n=1,nf
          cft(n,irec,j)=0.
        end do
      end do
      do j=1,2
        do n=1,nf
          rf(n,j)=zero
        end do
      end do
c  first part: compute spectrum of P coda
      do k=1,nwin
        do kcmp=1,3
          do i=1,npts
            ar(i)=a(i+nst,kcmp)*tap(i,k)
c            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)
          if(kcmp.le.2) then
            do n=1,nf
              rff(n,kcmp)=dcmplx(ar(n),ai(n))
              cft(n,irec,kcmp)=cft(n,irec,kcmp)+(ar(n)**2+ai(n)**2)
            end do
          elseif(kcmp.eq.3) then
            do n=1,nf
              zc=dcmplx(ar(n),ai(n))
              rf(n,1)=rf(n,1)+conjg(zc)*rff(n,1)
              rf(n,2)=rf(n,2)+conjg(zc)*rff(n,2)
              cft(n,irec,kcmp)=cft(n,irec,kcmp)+(ar(n)**2+ai(n)**2)
            end do
          endif
c          call plotit(fr,bb,dum,nf,'fft',
c     x     'frequency (Hz)','H(f)',1,0,0.0,0,1)
        end do
      end do
c  second part: assemble quantities for later regression weighting
      do n=1,nf
        cohrad=cabs(rf(n,1))**2/(cft(n,irec,1)*cft(n,irec,3))        
        cohtrn=cabs(rf(n,2))**2/(cft(n,irec,2)*cft(n,irec,3))        
        sum=0.
        do kcmp=1,3
          sum=sum+cft(n,irec,kcmp)
        end do
        cft(n,irec,1)=cohrad
        cft(n,irec,2)=cohtrn
        cft(n,irec,3)=sum
      end do
      go to 10
  111 continue
      close(10)
c  average coherence as function of frequency
      open(12,file='outr_baz.grid',form='formatted')
      open(13,file='outt_baz.grid',form='formatted')
c  plot distribution of back-azimuth and epicentral distance
      do j=1,36
        do i=1,72
          sumij(i,j)=0.0
          xi(i,j)=(i-1)*5.
          yi(i,j)=(j-1)*5.
        end do
      end do
      do baz=0,355,5
        ii=(baz+0.00000001)/5+1
        do i=1,irec
          test=abs(baz-bazs(i))
          if(test.lt.5..or.test.gt.355.) then
            jj=epic(i)/5+1
            sumij(ii,jj)=sumij(ii,jj)+1.0
            sumij(ii,jj+1)=sumij(ii,jj+1)+1.0
          endif
        end do
      end do
      nij=72*36
      call plotit(xi,yi,sumij,nij,'Data Distribution','Back Azimuth',
     x  'Epicentral Distance',5,0,0.0,3,1)
      do j=1,36
        do i=1,72
          if(sumij(i,j).lt.0.05) then
            sumij(i,j)=0.
          else
            sumij(i,j)=0.5+alog(sumij(i,j))
          endif
        end do
      end do
      call plotit(xi,yi,sumij,nij,'Data Distribution','Back Azimuth',
     x  'Epicentral Distance',5,0,0.0,3,1)
      print *,irec,' events processed'
 1009 print *,'enter gcarc interval for summing (GCMIN,GCMAX)'
      print *,'enter integers between 0 and 180'
      read(5,*) igcmin,igcmax
      if(igcmin.ge.igcmax) go to 1009
      print *,'bin halfwidth in degrees (e.g. 10)'
      read(5,*) dbaz
      kaz=0
      gc1=igcmin
      gc2=igcmax
      dbazh=dbaz/2.
      do baz=360.-dbazh,0.,-dbaz
c  this loop allows for greater overlap in averaging binsS
c      do baz=355.,0.,-10
        do l=1,2		! loop over R,T
          do n=1,nf  		! loop over freq
            zc=zero
            sumz=0.
            sumh=0.
            jrec=0    		! count records in the binstack
            do i=1,irec		! loop over seismic events
              test=abs(baz-bazs(i))
              if(test.lt.dbaz.or.test.gt.360.-dbaz) then
c  branch   to limit the Epicentral distance range for this sum
                if(epic(i).ge.gc1.and.epic(i).le.gc2) then
                  jrec=jrec+1
                dcoh=sqrt((1.0-cft(n,i,l))/(2.*cft(n,i,l)))*cft(n,i,3)
                  zc=zc+conjg(afft(n,i,3))*afft(n,i,l)/dcoh
                  sumz=sumz+cabs(afft(n,i,3))**2/dcoh
                  sumh=sumh+cabs(afft(n,i,l))**2/dcoh
                endif
              endif
            end do
            rf(n,l)=zc/sumz
            crft(n,l)=cabs(zc)**2/(sumz*sumh)
            if(jrec.gt.1) then
              drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*(jrec-1)))
     x                                         *cabs(rf(n,l))
            else  ! KLUGE TO AVOID ZERODIVIDE
              drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*0.5))
     x                                         *cabs(rf(n,l))
            endif
          end do
          if(jrec.ge.ibskip) then
            print *,baz,jrec
            kaz=kaz+1
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))
              drf(n,l)=fac*drf(n,l)
            end do
            do n=nf+1,npad-nf+1
              ar(n)=0.d0
              ai(n)=0.d0
            end do
            write(string,109) baz
            bbmx=0.
            do n=1,nf
              bb(n)=dsqrt(ar(n)**2+ai(n)**2)
              bbmx=amax1(bbmx,bb(n))
              pdat(n)=drf(n,l)/sq2
            end do
c            call plotit_axes(tdata,tdatb,-0.3*bbmx,1.03*bbmx)
c            call plotit(tdat,bb,pdat,nf,'Frequency Domain',
c     x   'frequency (Hz)','|H(f)|',3,0,0.05,0,58+l*3)
            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
c            call plotit_axes(tdata,tdatb,-180.,180.)
c            call plotit(tdat,bb,pdat,nf,'Frequency Domain',
c     x   'frequency (Hz)','H(f) phase',3,0,0.05,0,59+l*3)
            call plotit_axes(0.,0.,0.,0.)
            do n=1,nf
              bb(n)=ar(n)
            end do
            call plotit_axes(0.,0.,0.,0.)
            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))
            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,l)=baz+50.*bb(i)*l
              t3=-time(i)
              iunit=11+l
c              if(5*((i-1)/5).eq.i-1) write(iunit,1022) t3,baz,bb(i)
              write(iunit,1022) t3,baz,bb(i)
            end do
c            call plotit_axes(0.,0.,0.,0.)
c            call plotit(ttt,bb,dum,npost,string,
c     x    'time(sec)','H(t)',1,0,0.0,0,60+3*l)
          endif
        end do
        write(12,101) '>'
        write(13,101) '>'
        call plotit_axes(0.,0.,0.,0.)
        if(jrec.ge.ibskip) then
          do i=1,ntot
            bb(i)=(rft(i,kz,1)-baz)/50.
          end do
c          call plotit(ttt,bb,dum,npost,string,
c     x    'time(sec)','H(t)',1,0,0.0,0,0)
          do i=1,ntot
            bb(i)=(rft(i,kz,2)-baz)/100.
          end do
c          call plotit(ttt,bb,dum,npost,string,
c     x    'time(sec)','H(t)',1,0,0.1,0,1)
        endif
      end do  
      kaz=kaz/2
      print *,kaz,' traces'   
      nss=7./dt
      npss=nss+81./dt
      do iagain=1,2
      call plotit_axes(-10.,370.,0.,0.)
      do n=1,kaz
        call plotit(rft(nss,n,2),time(nss),dum,npss,
     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(-10.,370.,0.,0.)
      do n=1,kaz
        call plotit(rft(nss,n,1),time(nss),dum,npss,
     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
      call plotit_axes(-10.,370.,0.,0.)
      do n=1,kaz
        call plotit(rft(nss,n,2),time(nss),dum,npost,
     x   'Composite Transverse RF Section',
     x   'Back-Azimuth (degrees CW from N)',
     x  'time(sec)',2,0,0.0,0,(n/kaz))
      end do
            
      call plotit_axes(-10.,370.,0.,0.)
      do n=1,kaz
        call plotit(rft(nss,n,1),time(nss),dum,npost,
     x   'Composite Radial RF Section',
     x   'Back-Azimuth (degrees CW from N)',
     x  'time(sec)',2,0,0.0,0,(n/kaz))
      end do
      close(12)
      close(13)
cc  NOW PLOT TRACES THAT ARE BINNED WITH EPICENTRAL DISTANCE
      print *,'compute RFs binned with epicentral distance'
      print *,'enter back-azimuth limits ib1,ib2 (integers!)'
      print *,' ib1=ib2 -> 0,360 and 360-wraparound if ib1 > ib2 '
      read(5,*) ib1,ib2
      if(ib1.eq.ib2) then
        ick=1
        baz1=0.
        baz2=360.
      elseif(ib1.lt.ib2) then
        ick=1
        baz1=ib1
        baz2=ib2
      else
        ick=2
        baz1=ib1
        baz2=ib2
      endif
      kaz=0
      open(12,file='outr_epi.grid',form='formatted')
      open(13,file='outt_epi.grid',form='formatted')
c  code for changing the epicentral bin width
      ep1=0.
      ep2=175.
      dep=5.
      print *,'epicentral bins from',ep1,' to',ep2
      print *,'with halfbin width',dep
      print *,' do you want to change this spacing? (1=yes)'
      read(5,*) ickk
      if(ickk.eq.1) then
        print *,'enter delta range and spacing (ep1,ep2,dep)'
        read(5,*) ep1,ep2,dep
        nep=(ep2-ep1)/dep+1
        if(nep.gt.500.or.nep.lt.1) 
     x     print *,'HEY! There are ',nep,' traces!'
      endif
c  end code for changing the epicentral bin width
      do n=1,nf  		
        tdat(n)=fr(1)+(n-1)*fr(2)
      end do
      do epi=ep2,ep1,-dep
        do l=1,2		! loop over R,T
          do n=1,nf  		! loop over freq
            zc=zero
            sumz=0.
            sumh=0.
            jrec=0    		! count records in the binstack
            do i=1,irec		! loop over seismic events
              test=abs(epi-epic(i))
              if(test.lt.dep) then
c  branch to limit the back-azimuth range for this sum
                ick1=0
                if(ick.eq.1) then
                  if(bazs(i).ge.baz1.and.bazs(i).le.baz2) ick1=1
                else
                  if(bazs(i).ge.baz1.or.bazs(i).le.baz2) ick1=1
                endif
c  normalize the records to enhance highly coherent records
c  and normalize to equalize the signal size, so large events do not dominate
                if(ick1.eq.1) then
                  jrec=jrec+1
                  dcoh=sqrt((1.0-cft(n,i,l))/(2.*cft(n,i,l)))*cft(n,i,3)
                  zc=zc+conjg(afft(n,i,3))*afft(n,i,l)/dcoh
                  sumz=sumz+cabs(afft(n,i,3))**2/dcoh
                  sumh=sumh+cabs(afft(n,i,l))**2/dcoh
                endif
              endif
            end do
            rf(n,l)=zc/sumz
            crft(n,l)=cabs(zc)**2/(sumz*sumh)
            if(jrec.gt.1) then
              drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*(jrec-1)))
     x                                         *cabs(rf(n,l))
            else  ! KLUGE TO AVOID ZERODIVIDE
              drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*0.5))
     x                                         *cabs(rf(n,l))
            endif
          end do
          if(jrec.gt.ibskip) then
            print *,epi,jrec
            kaz=kaz+1
c  no weighted sum yet  -- just a crude sum with simple error estimate
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))
              drf(n,l)=fac*drf(n,l)
            end do
            do n=nf+1,npad-nf+1
              ar(n)=0.d0
              ai(n)=0.d0
            end do
            write(string,1099) epi
            bbmx=0.
            do n=1,nf
              bb(n)=dsqrt(ar(n)**2+ai(n)**2)
              bbmx=amax1(bbmx,bb(n))
              pdat(n)=drf(n,l)/sq2
            end do
            call plotit_axes(tdata,tdatb,-0.3*bbmx,1.03*bbmx)
c            call plotit(tdat,bb,pdat,nf,'Frequency Domain',
c     x                  'frequency (Hz)','|H(f)|',3,0,0.05,0,58+l*3)
            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(tdata,tdatb,-180.,180.)
c            call plotit(tdat,bb,pdat,nf,'Frequency Domain',
c     x   'frequency (Hz)','H(f) phase',3,0,0.05,0,59+l*3)
            call plotit_axes(0.,0.,0.,0.)
            do n=1,nf
              bb(n)=ar(n)
            end do
            call plotit_axes(0.,0.,0.,0.)
            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))
            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,l)=epi+50.*bb(i)*l
              t3=-time(i)
              iunit=11+l
c              if(5*((i-1)/5).eq.i-1) write(iunit,1022) t3,epi,bb(i)
              write(iunit,1022) t3,epi,bb(i)
            end do
            call plotit_axes(0.,0.,0.,0.)
c            call plotit(ttt,bb,dum,npost,string,
c     x    'time(sec)','H(t)',1,0,0.0,0,60+3*l)
          endif
        end do
        call plotit_axes(0.,0.,0.,0.)
        if(jrec.ge.ibskip) then
          do i=1,ntot
            bb(i)=(rft(i,kz,1)-epi)/50.
          end do
c          call plotit(ttt,bb,dum,npost,string,
c     x    'time(sec)','H(t)',1,0,0.0,0,0)
          do i=1,ntot
            bb(i)=(rft(i,kz,2)-epi)/100.
          end do
c          call plotit(ttt,bb,dum,npost,string,
c     x    'time(sec)','H(t)',1,0,0.1,0,1)
        endif
        write(12,101) '>'
        write(13,101) '>'
      end do  
 1022 format(f7.3,f6.1,f7.3)
      close(12)
      close(13)
      kaz=kaz/2
      print *,kaz,' traces' 
      nss=7./dt
      npss=nss+81./dt
      do iagain=1,2
      call plotit_axes(ep1-dep,ep2+dep,0.,0.)
      do n=1,kaz
        call plotit(rft(nss,n,2),time(nss),dum,npss,
     x   'Composite Transverse RF Section',
     x   'Epicentral Distance (degrees)',
     x  'time(sec)',2,0,0.0,0,21*(n/kaz))
      end do
            
      call plotit_axes(ep1-dep,ep2+dep,0.,0.)
      do n=1,kaz
        call plotit(rft(nss,n,1),time(nss),dum,npss,
     x   'Composite Radial RF Section',
     x   'Epicentral Distance (degrees)',
     x  'time(sec)',2,0,0.0,0,22*(n/kaz))
      end do
      end do
      call plotit_axes(ep1-dep,ep2+dep,0.,0.)
      do n=1,kaz
        call plotit(rft(nss,n,2),time(nss),dum,npost,
     x   'Composite Transverse RF Section',
     x   'Epicentral Distance (degrees)',
     x  'time(sec)',2,0,0.0,0,(n/kaz))
      end do
            
      call plotit_axes(ep1-dep,ep2+dep,0.,0.)
      do n=1,kaz
        call plotit(rft(nss,n,1),time(nss),dum,npost,
     x   'Composite Radial RF Section',
     x   'Epicentral Distance (degrees)',
     x  'time(sec)',2,0,0.0,0,(n/kaz))
      end do
  109 format('Back Azimuth Centered on ',f4.0,' degrees')
 1099 format('Epicentral Distance Centered on ',f4.0,' degrees')
      stop
      end




   
      subroutine taper(n,nwin,el)
c
c  to generate slepian tapers
c  ta is a real*4 array
c
c         j. park
c
      implicit real*4(a-h,o-z)
      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(270000)  ! we use this common block for scratch space
      common/stap2/ta(16400,16)
      dimension a(16400,8),el(10)
      equivalence (a(1,1),ta(1,1))
      data iwflag/0/,pi/3.14159265358979d0/
      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
c      print *,'eigenvalues for the eigentapers (k near 2n)'
c      print *,(ell(i),i=1,nwin)
      do i=1,nwin
        el(i)=dmax1(ell(i),el(i))
      end do     
c      print *,'composite asymptotics for taper eigenvalues'
c      print *,(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      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