c  program recfunk_mwm
c  6/3/04 JJP  --- UPDATED 08/19/07 JJP
c  MIGRATION/MOVEOUT variant to compute cross-correlation in a moving time window
c  you tell the program the depth at which you want to focus the cross-correlation
c  and it computes the moving-window delay from an input model
c  using the epicentral distance to look up the horizontal slowness
c  multitaper program to invert for receiver functions
c  SINGLE STATION WITH SPECTRAL WEIGHTING
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 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_NNN.grid out[rt]_epi_NNN.grid 
c  where NNN corresponds to the target depth in kilometers
c
c  these files are overwritten everytime you run the program
c  so rename them to save them
c
c  xf77 -o /park/backup/bin/recfunk_mwm recfunk_mwm.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c  xf77 -o /Users/jjpark/bin/recfunk_mwm recfunk_mwm.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
      implicit real*4(a-h,o-z)
      real*8 ar,ai,el(12),aaa
      complex*8 zc,zero
      complex*8 afft,rf,rfs
      character*80 name,subtext,string,name_stack
      character*28 title(2)
      character*10 cmps(3)
      character comp(3)
      common/nnaamm/iwflag
      common/npiprol/anpi
      common/moovwin/nst(3)
      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),epic(799),s2n(799),ip(799),iflag(799)
      common/saveit2/crfts(4100,799,2)
      common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2)
c NOTE THAT MODEL PARAMETERS ARE REAL*4 IN THIS CODE, NOT REAL*8 AS IN ANIREC
c ONLY 30 LAYERS OVER HALFSPACE ARE ALLOWED
      common/model/z(31),dz(31),rho(31),vp(31),vs(31),vp2(31),
     x               vp4(31),vs2(31),vss(31)
      common/model2/xmu(31),xla(31),xmu2(31),xla2(31),xla4(31)
      common/modstk/tau(31),itau(31),psd(31),gam(31),xik(31)
      common/migstk/arm(4100,2,31),aim(4100,2,31),dam(4100,2,31)
      common/tt_p/d_p(200),d_pkp(200),s_p(200),s_pkp(200),n_p,n_pkp
      common/header/ahead(158)
      common/distrib/xi(72,36),yi(72,36),sumij(72,36)
      common/chisq/chi(160000),enchi(160000),bbaz(160000),
     x chim(200,2),enchim(200,2)
      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.
      zero=cmplx(0.,0.)
      fmax=5.
      print *,'velocity model for stacking (blank = stack_model)'
      read(5,101) name_stack
      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.30) then
        print *,'only 30 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)
      end do
      close(7)
      do i=2,nl
        dz(i)=z(i)-z(i-1)
      end do
      dz(1)=z(1)
c  print the organ-pipe mode count for 1Hz
c  the lowest layer (nl+1) is taken as evanescent region.
      sdelay=0.
      pdelay=0.
      do i=1,nl
        sdelay=sdelay+(dz(i)/vs(i))
        pdelay=pdelay+(dz(i)/vp(i))
        tau(i)=sdelay-pdelay
        print *,z(i),tau(i)
      end do
      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
      print *,'enter the depth at which you wish to focus the RFs (km)'
      read(5,*) itarget
      targm=itarget*1000.
      do i=nl,1,-1
        resz=targm-z(i)
	vdelay=tau(i)
        if(targm.gt.z(i)) go to 313
      end do
      i=0
      resz=targm
      vdelay=0.
  313 idelay=i+1
      vdelay=vdelay+resz/vs(idelay)-resz/vp(idelay)
      print *,'delay at vert-inc = ',vdelay,' sec, idelay=',idelay
  104 format(2f10.1)
      print *,'we will divide the spectrum later!'
      print *,'enter fmax (Hz)'
      read(5,*) fmax
c generate tapers for npad=16384 point records
      fmax0=fmax
      fmin=0
      anpi=2.5
      nwin=3
      npad=16384
      nnf=npad/2
      do i=1,1000
        crft(i,1)=0.
        crft(i,2)=0.
      end do
      irec=0
      call read_ttimes
 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)'
        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 *,'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) stop
        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
c  naturally, you are in deep poop if dt is not the same for all records
      npre=12./dt
      npost=20./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
      irec=irec+1
      iflag(irec)=1
      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
	  iflag(irec)=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
      if(irec.gt.799) then
        print *,'max records exceeded'
	stop
      endif
      baz=ahead(53)
      bazs(irec)=ahead(53)
      epic(irec)=ahead(54)
      if(imark.eq.0.and.epic(irec).gt.118.) iflag(irec)=2
      print *,baz,'= back azimuth'
c  ddf is cycles/day
      ddf=1./(dt*npad)
      nf1=1
      nf2=fmax/ddf+1
      nf=nf2-nf1+1
      if(nf.le.0) stop
      if(nf.gt.4100) then
        print *,'nf is too big', nf
	stop
      endif
      fr(1)=0.
      fr(2)=ddf
      fr(3)=ddf
c  get p-slowness -- iflag branches on the start time of the analysis window
      if(imark.eq.3) then    ! PP wave
        epic(irec)=epic(irec)/2.
      endif
      igcarc=epic(irec)
      fac=epic(irec)-igcarc
      print *,'igcarc,i,epic(i),fac,imark',igcarc,irec,epic(irec),fac,imark
      if(iflag(irec).eq.1) then		! P-phase
        if(igcarc.ge.180.or.igcarc.le.0) then
          print *,'WHOA! igcarc,gcarc = ',igcarc,epic(irec)
          pause
        endif
c  allow extreme P_diff phases (moveout is constant for P_diff)
        if(igcarc.gt.129) igcarc=129 
        slow=(1.-fac)*s_p(igcarc)+fac*s_p(igcarc+1)
      elseif(iflag(irec).eq.2) then	! PKP-phase
        if(igcarc.ge.180.or.igcarc.lt.90) then
          print *,'WHOA! igcarc,gcarc = ',igcarc,epic(irec)
          pause
        endif
        ig=igcarc-89
        slow=(1.-fac)*s_pkp(ig)+fac*s_pkp(ig+1)
      else
        print *,'WHOA! iflag(i),i = ',iflag(irec),irec
        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
c  if slowness of the P wave is too large to imply raybottom beneath target,
c  we skip the record
      if(slw*vp(idelay).ge.1.0) then
        irec=irec-1
	print *,'ray slowness too great',slw
	pause
	go to 10
      endif
      tshft=0.
      if(idelay.ge.1) then
        do k=1,idelay-1
          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
          tshft=tshft+dz(k)*(coss/vs(k)-cosp/vp(k))
        end do
      endif
      coss=sqrt(1.-vs(idelay)*vs(idelay)*slw*slw)
      cosp=sqrt(1.-vp(idelay)*vp(idelay)*slw*slw)
      tshft=tshft+resz*(coss/vs(idelay)-cosp/vp(idelay))
      s2n(irec)=tshft
      npts=postwin/dt
      nst(3)=pstart/dt
      nst(1)=(pstart+tshft)/dt
      nst(2)=(pstart+tshft)/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(j),j=1,3),npts
      print *,'dt,nscan,ddf',dt,nscan,ddf
c      compute spectrum of pre-event noise
      npts0=min0(nst(3)+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
      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
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,npts
          bb(i)=cs*a(nst(3)+i,3)+sn*a(nst(3)+i,1)
	  bb(i+npts)=cs*a(nst(1)+i,1)-sn*a(nst(1)+i,3)
        end do
        do i=1,npts
	  a(i+nst(3),3)=bb(i)
          a(i+nst(1),1)=bb(i+npts)
        end do
      endif
c  KLUGE to write out some data
c      if(epic(irec).gt.70.) then
c        open(9,file='data4plot.dat',form='formatted')
c	bmx1=0.
c	bmx2=0.
c	bmx3=0.
c	bmn1=0.
c	bmn2=0.
c	bmn3=0.
c	do i=1,nscan
c	  bmx1=amax1(bmx1,a(i,1))
c	  bmx2=amax1(bmx2,a(i,2))
c	  bmx3=amax1(bmx3,a(i,3))
c	  bmn1=amin1(bmn1,a(i,1))
c	  bmn2=amin1(bmn2,a(i,2))
c	  bmn3=amin1(bmn3,a(i,3))
c	end do
c	space2=1.1*(bmx1-bmn2)
c	space3=space2+1.1*(bmx2-bmn3)
c	do i=1,nscan
c	  tttt=i*dt
c	  ddd=a(i,1)
c	  write(9,*) tttt, ddd
c	end do
c	write(9,101) '>'
c	do i=1,nscan
c	  tttt=i*dt
c	  ddd=a(i,2)+space2
c	  write(9,*) tttt, ddd
c	end do
c	write(9,101) '>'
c	do i=1,nscan
c	  tttt=i*dt
c	  ddd=a(i,3)+space3
c	  write(9,*) tttt, ddd
c	end do
c	close(9)
c	pstart2=pstart+tshft
c	print *,space2,space3,nscan,pstart,tshft,pstart2
c        pause
c      endif
c  END KLUGE	
c  end alternate estimate of angle of incidence
      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(3)+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),kcmp)*tap(i,k)
            bb(i)=a(i+nst(kcmp),kcmp)*tap(i,k)
            ai(i)=0.d0
          end do
          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  compute a simple spectral  signal-to-noise measure 
c      s2n(irec)=0.
c      do n=1,nf
c        s2n(irec)=s2n(irec)+alog(sss(n,6)/sss(n,3))
c      end do
c      s2n(irec)=s2n(irec)/nf
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)
	  crfts(n,irec,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)
        end do
      end do
      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
c        bbmx=0.
c        do n=1,npad
c          bb(n)=dsqrt(ar(n)**2+ai(n)**2)
c          bbmx=amax1(bbmx,bb(n))
c        end do
c        do n=1,npad
c          if(bb(n).gt.1e-4) then
c            pdat(n)=con*(pdat(n)/bb(n)) 
c          else
c            pdat(n)=360.
c          endif
c          bb(n)=con*datan2(ai(n),ar(n))
c        end do
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
c        fac=(4./3.)*float(nnf)/(float(npad)*float(nf))
        fac=2.*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
        end do
        do i=1,npre
          bb(npre+1-i)=ar(npad+1-i)*fac
        end do
        do i=1,ntot 
          rft(i,irec,l)=baz+50.*bb(i)
        end do
      end do
      go to 10
  111 continue
      close(10)
c      call plotit(epic,s2n,dum,irec,'moveout chk','gcarc','delay (sec)',
c     x  2,0,0.1,3,1)
c  average coherence as function of frequency
c      do l=1,2
c        do i=1,nf
c          crft(i,l)=crft(i,l)/irec
c        end do
c      end do
c      call plotit_axes(0.,0.,0.,0.)
c      call plotit(fr,crft(1,1),dum,nf,'Aggregate Radial Coherence',
c     x    'Freq(Hz)','Coherence',1,0,0.0,0,21)
c      call plotit_axes(0.,0.,0.,0.)
c      call plotit(fr,crft(1,2),dum,nf,'Aggregate Transv Coherence',
c     x    'Freq(Hz)','Coherence',1,0,0.0,0,22)
      print *,irec
  484 continue
  116 format(a,i3,a)
      write(name,116) 'outr_baz_',itarget,'.grid'
      print *, itarget
      print *,name
      if(itarget.lt.100) name(10:10)='0'
      if(itarget.lt.10) name(11:11)='0'
      open(12,file=name,form='formatted')
      write(name,116) 'outt_baz_',itarget,'.grid'
      if(itarget.lt.100) name(10:10)='0'
      if(itarget.lt.10) name(11:11)='0'
      open(13,file=name,form='formatted')
c      open(14,file='out_baz.chisq',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
      print *,irec
      kaz=0
      naz=0
      naz0=1
      maz=0
c  ndf is a freq-spacing for stepping thru the spectral arrays at approx Rayleigh freq
      ndf=float(npad)/float(npts)+0.25
      print *,nf1,nf2
      bazmax=355.
      bazmin=0.
      bazinc=-5.
      print *,'back-azimuth range to plot:'
      print *,bazmax, bazmin, bazinc
      print *,'change baz-interval or baz-spacing? (1=yes)'
      read(5,*) ickbaz
      if(ickbaz.eq.1) then
        print *,'enter back-azimuth range:'
        print *,'bazmax, bazmin, bazinc'
        print *,'will make bazinc negative to step back thru baz values'
        read (5,*) bazmax, bazmin, bazinc
        if((bazmax-bazmin)*bazinc.gt.0) bazinc=-bazinc
      endif
      abazinc=abs(bazinc)
      do baz=bazmax, bazmin, bazinc
        do l=1,2
          do n=nf1,nf2
            rf(n,l)=zero
            drf(n,l)=0.
          end do    
          jrec=0    
          do i=1,irec
            test=abs(baz-bazs(i))
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
c	    if(s2n(i).ge.s2n_floor) then
              if(test.lt.abazinc.or.test.gt.360.-abazinc) then
                jrec=jrec+1
		ip(jrec)=i
                do n=nf1,nf2
                  sig=drfs(n,i,l)**2
                  rf(n,l)=rf(n,l)+rfs(n,i,l)/sig
                  drf(n,l)=drf(n,l)+1.0/sig
                end do  
              endif
c            endif
          end do
          print *,baz,jrec
          if(jrec.ge.ibskip) then
            kaz=kaz+1
c  compute chi-squared misfit, avoid zero-freq
            ndf1=max0(ndf,nf1)
c	    print *,ndf1,nf2,ndf
	    do n=ndf1,nf2,ndf
	      naz=naz+1
	      chisq=0.
              do j=1,jrec
                sig=drfs(n,ip(j),l)**2
	        chisq=chisq+cabs(rf(n,l)/drf(n,l)-rfs(n,ip(j),l))**2/sig
	      end do
  	      chi(naz)=chisq
	      bbaz(naz)=baz
	      nchi=2*jrec-2
	      enchi(naz)=nchi
c	      write(14,114) nchi,chisq,baz,l
c	      if(chisq.gt.1000.) then
c	        print *,jrec,chisq,baz,n
c		pause
c	      endif
	    end do
c  find median misfit across all frequencies
            nnaz=naz-naz0+1
	    nnazh=nnaz/2
	    do nn=naz0,naz
	      ij=0
	      do nnn=naz0,naz
	        if(chi(nnn).le.chi(nn)) ij=ij+1
	      end do
c	      if(ij.eq.nnazh.and.jrec.gt.1.and.l.eq.2) then
	      if(ij.eq.nnazh.and.jrec.gt.1) then
	        if(l.eq.1) maz=maz+1
	        chim(maz,l)=chi(nn)
		enchim(maz,l)=enchi(nn)
		go to 787
	      endif
	    end do
  787       continue	    
c	    print *,(chi(nn),nn=naz0,naz)
c	    print *,(enchi(nn),nn=naz0,naz)
c	    pause
	    naz0=naz+1
            do n=1,npad
              ar(n)=0.d0
              ai(n)=0.d0
            end do
            do n=nf1,nf2
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
              rf(n,l)=rf(n,l)/drf(n,l)
              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*sqrt(1./drf(n,l))
            end do
            write(string,109) baz
c            bbmx=0.
c            do n=1,nf
c              bb(n)=dsqrt(ar(n)**2+ai(n)**2)
c              bbmx=amax1(bbmx,bb(n))
c              pdat(n)=drf(n,l)/sq2
c            end do
c            do n=1,npad
c              if(bb(n).gt.1e-4) then
c                pdat(n)=con*(pdat(n)/bb(n)) 
c              else
c                pdat(n)=360.
c              endif
c              bb(n)=con*datan2(ai(n),ar(n))
c            end do
c            do n=1,nf
c              bb(n)=ar(n)
c            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))
c            fac=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
c  note that the transverse for plotit plotting is boosted by factor of 2
c but the transverse RF written to disk (for pswiggle) is not scaled
            do i=1,ntot 
              rft(i,kz,l)=baz+50.*bb(i)*l
              t3=-time(i)
              iunit=11+l
              write(iunit,1022) t3,baz,bb(i)
            end do
            write(iunit,101) '>'
          endif
        end do
      end do  
      kaz=kaz/2
      print *,kaz,' traces'   
      close(12)
      close(13)
c      close(14)
c      enmax=1
c      do i=1,maz
c        enmax=amax1(enmax,enchim(i,1))
c        enmax=amax1(enmax,enchim(i,2))
c      end do
c      bb(1)=1.
c      bb(2)=enmax
c      call plotit(bb,bb,dum,2,' ',' ',' ',2,0,0.0,0,0)
c      call plotit(enchim(1,1),chim(1,1),dum,maz,'chi-squared test',
c     x            'deg of freed','median chi-squared value',2,0,0.3,1,0)
c      call plotit(enchim(1,2),chim(1,2),dum,maz,'chi-squared test',
c     x           'deg of freed', 'median chi-squared value',2,0,0.3,2,1)
c      call plotit(enchi,chi,dum,naz,'chi-squared test','deg of freedom',
c     x            'chi-squared value',2,0,0.3,3,1)
c      call plotit(enchi,chi,dum,naz,'chi-squared test','deg of freedom',
c     x            'chi-squared value',2,1,0.3,3,1)
c      print *,'Try a different frequency interval? (0=no)'
c      read(5,*) ick
c      if(ick.ne.0) then
c        print *,fmin,fmax
c        print *,'enter new fmin & fmax (Hz)'
c	read(5,*) fmin,fmax
c	fmax=amin1(fmax,fmax0)
c	nf1=max0(int(fmin/ddf),1)
c	nf2=fmax/ddf+1
c	nf=nf2-nf1+1
c	print *,fmin,fmax,nf1,nf2,ddf,nf
c	go to 484
c      endif
c  return to default freq interval
      fmin=0.
      fmax=fmax0
      nf1=1
      nf2=fmax/ddf+1
      nf=nf2
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
  494 continue
      write(name,116) 'outr_epi_',itarget,'.grid'
      if(itarget.lt.100) name(10:10)='0'
      if(itarget.lt.10) name(11:11)='0'
      open(12,file=name,form='formatted')
      write(name,116) 'outt_epi_',itarget,'.grid'
      if(itarget.lt.100) name(10:10)='0'
      if(itarget.lt.10) name(11:11)='0'
      open(13,file=name,form='formatted')
c      open(14,file='out_epi.chisq',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
 2020   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) then
          print *,'HEY! There are ',nep,' traces!'
          go to 2020
	endif
      endif
c  end code for changing the epicentral bin width
      kaz=0
      naz=0
      naz0=1
      maz=0
      do epi=ep1,ep2,dep
        do l=1,2
          do n=nf1,nf2
            rf(n,l)=zero
            drf(n,l)=0.
          end do    
          jrec=0    
          do i=1,irec
            test=abs(epi-epic(i))
c	    if(s2n(i).ge.s2n_floor) then
              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
                if(ick1.eq.1) then
                  jrec=jrec+1
	 	  ip(jrec)=i
                  do n=nf1,nf2
                    sig=drfs(n,i,l)**2
                    rf(n,l)=rf(n,l)+rfs(n,i,l)/sig
                    drf(n,l)=drf(n,l)+1.0/sig
                  end do  
                endif
              endif
c            endif
          end do
          if(jrec.ge.ibskip) then
            print *,epi,jrec
            kaz=kaz+1
c  compute chi-squared misfit, avoid zero-freq
            ndf1=max0(ndf,nf1)
	    do n=ndf1,nf2,ndf
	      naz=naz+1
	      chisq=0.
              do j=1,jrec
                sig=drfs(n,ip(j),l)**2
	        chisq=chisq+cabs(rf(n,l)/drf(n,l)-rfs(n,ip(j),l))**2/sig
	      end do
  	      chi(naz)=chisq
	      bbaz(naz)=baz
	      nchi=2*jrec-2
	      enchi(naz)=nchi
c	      write(14,114) nchi,chisq,baz,l
c	      if(chisq.gt.1000.) then
c	        print *,jrec,chisq,baz,n
c		pause
c	      endif
            end do
c  find median misfit across all frequencies
            nnaz=naz-naz0+1
	    nnazh=nnaz/2
	    do nn=naz0,naz
	      ij=0
	      do nnn=naz0,naz
	        if(chi(nnn).le.chi(nn)) ij=ij+1
	      end do
c	      if(ij.eq.nnazh.and.jrec.gt.1.and.l.eq.2) then
	      if(ij.eq.nnazh.and.jrec.gt.1) then
	        if(l.eq.1) maz=maz+1
	        chim(maz,l)=chi(nn)
		enchim(maz,l)=enchi(nn)
		go to 797
	      endif
	    end do
  797       continue	    
c	    print *,(chi(nn),nn=naz0,naz)
c	    print *,(enchi(nn),nn=naz0,naz)
c	    pause
	    naz0=naz+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,npad
              ar(n)=0.d0
              ai(n)=0.d0
            end do
            do n=nf1,nf2
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
              rf(n,l)=rf(n,l)/drf(n,l)
              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*sqrt(1./drf(n,l))
            end do
c            write(string,1099) epi
c            bbmx=0.
c            do n=1,nf
c              bb(n)=dsqrt(ar(n)**2+ai(n)**2)
c              bbmx=amax1(bbmx,bb(n))
c              pdat(n)=drf(n,l)/sq2
c            end do
c            do n=1,npad
c              if(bb(n).gt.1e-4) then
c                pdat(n)=con*(pdat(n)/bb(n)) 
c              else
c                pdat(n)=360.
c              endif
c              bb(n)=con*datan2(ai(n),ar(n))
c            end do
c            do n=1,nf
c              bb(n)=ar(n)
c            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))
c            fac=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
              write(iunit,1022) t3,epi,bb(i)
            end do
            write(iunit,101) '>'
          endif
        end do
      end do  
 1022 format(f7.3,f6.1,f7.3)
      close(12)
      close(13)
c      close(14)
c      enmax=1
c      do i=1,maz
c        enmax=amax1(enmax,enchim(i,1))
c        enmax=amax1(enmax,enchim(i,2))
c      end do
c      bb(1)=1.
c      bb(2)=enmax
c      call plotit(bb,bb,dum,2,' ',' ',' ',2,0,0.0,0,0)
c      call plotit(enchim(1,1),chim(1,1),dum,maz,'chi-squared test',
c     x            'deg of freed','median chi-squared value',2,0,0.3,1,0)
c      call plotit(enchim(1,2),chim(1,2),dum,maz,'chi-squared test',
c     x           'deg of freed', 'median chi-squared value',2,0,0.3,2,1)
c      print *,'Try a different frequency interval? (0=no)'
c      read(5,*) ick
c      if(ick.ne.0) then
c        print *,fmin,fmax
c        print *,'enter new fmin & fmax (Hz)'
c	read(5,*) fmin,fmax
c	fmax=amin1(fmax,fmax0)
c	nf1=max0(int(fmin/ddf),1)
c	nf2=fmax/ddf+1
c	nf=nf2-nf1+1
c	print *,fmin,fmax,nf1,nf2,ddf,nf
c	go to 494
c      endif
      kaz=kaz/2
      print *,kaz,' traces'   
  109 format('Back Azimuth Centered on ',f4.0,' degrees')
  114 format(i5,g15.4,f6.0)
 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
      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,16)
      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
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
      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',
c     x  form='formatted')
      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',
c     x  form='formatted')
      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
c      call plotit_axes(0.,0.,0.,0.)
c      call plotit(d_p,s_p,dum,n_p,' ',' ',' ',2,0,0.0,0,0)
c      call plotit(d_pkp,s_pkp,dum,n_pkp,' ',
c     x 'Epicentral Distance (degrees)','Slowness (sec/km)',2,0,0.0,0,1)
      return
      end