c  program recfunk_svd
c  variant to compute receiver functions via SVD of all three components, 
c  rather than by pairwise crosscorrelation
c
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  original code 6/2/04; updated 08/19/07 JJP
c  multitaper program to compute receiver functions
c  SINGLE STATION WITH SPECTRAL WEIGHTING
c  
c  uses spectral estimates of pre-event noise 
c  to weight the inversion in the frequency domain, thereby avoiding the 
c  usual problem with microseism noise in 0.1-0.5 Hz
c IN PRACTICE, however, the regularization of what might be a spectral division is achieved
c  via multiple taper cross-correlation
c
c  Program has extra code to perform chi-squared tests, discussed in Park and Levin paper 
c
c  xf77 -o /park/backup/bin/recfunk_svd recfunk_svd.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c  xf77 -o /Users/jjpark/bin/recfunk_svd recfunk_svd.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,ssvd
      complex*8 zc,zero,uh1,vh1
      complex*8 afft,rf,rfs,asvd,usvd,vsvd
      complex*16 zsvd,zu,zv
      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,8,3),tai(16400,8,3)
      common/stap0/ar(16400),ai(16400),bb(16400),time(99000)
      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(2050,799,2),
     x   drfs(2050,799,2),bazs(799),epic(799),s2n(799),ip(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)
      common/chisq/chi(160000),enchi(160000),bbaz(160000),
     x chim(200,2),enchim(200,2)
      common/svdstuff/zsvd(3,3),zu(3,3),zv(3,3),ssvd(3),uh1(3),vh1(3),
     x                tnor(3)
      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 *,'This code reads start of interval from sac header'
      print *,'enter fmax (Hz)'
      read(5,*) fmax
c  CHANGES FOR READING TIME INTERVAL FROM SAC HEADER
      print *,'enter time duration for analysis (sec)'
      read(5,*) postwin
      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
C end changes
c some default values (often superceded)
      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
      open(10,file='in_recpick',form='formatted')
   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=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
      irec=irec+1
      if(irec.gt.799) then
        print *,'max records exceeded'
	stop
      endif
      baz=ahead(53)
      bazs(irec)=ahead(53)
      epic(irec)=ahead(54)
      print *,baz,'= back azimuth'
c  ddf is cycles/sec
      ddf=1./(dt*npad)
      nf1=1
      nf2=fmax/ddf+1
      nf=nf2-nf1+1
      if(nf.le.0) stop
      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 to provide a pre-event noise window of equal length
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
c  ahead(9) is SAC's A parameter, for arrival time.  Set by pressing A in plotpk
      pstart=tmark-ahead(6)-3.
      print *,'sach params b and a, pstart',ahead(6),ahead(9),pstart
      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      compute spectrum of pre-event noise
      npts0=min0(nst+1,npts)
      if(npts0.lt.npts) then
        call taper(npts0,nwin,el)
        print *,'eigenvalues ',(el(k),k=1,nwin)
      endif
      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 DELTA 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   compute spectrum of pre-event noise on each component
        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  compute K eigenspectra 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
          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  there is code (commented out) that allows the user to screen out events with low S2N
      s2n(irec)=0.
      do n=1,nf
        s2n(irec)=s2n(irec)+alog(sss(n,6)/sss(n,3))
      end do
      s2n(irec)=s2n(irec)/nf
c  SVD correlation of three components
c  at each frequency, we have nwin spectrum estimates on 3 differnt components
c  since nwin=3 typically, the correlation matric is 3x3 complex
c  form a matrix A with
c
c   A_ij = Y^{(j)}_i  -- the ith eignewindow (rows) for the jth component (columns)
c
c  in an ideal world the matrix A_ij would be a dyad u_i v_j,
c  with u_i a linear combination of eigenspectra representing the clean P wave
c  and  v_j the linear combination that maximizes variance.
c  in practice this means that we use the largest singular pair as a "clean filter"
c  for the P wave.
c
      do n=1,nf
        bb(1)=1.
	bb(2)=1.
	bb(3)=1.
	do l=1,3
	  tnor(l)=0.
          do k=1,nwin
	    zsvd(k,l)=afft(n,k,l)*bb(l)
	    tnor(l)=tnor(l)+cabs(afft(n,k,l))**2
	  end do
	end do
	call csvd(zsvd,3,3,3,3,0,3,3,ssvd,zu,zv)
c  isolate the first eigenvector  -- left eigenvectors must be conjugated
c        totvar=0.
        do k=1,3
          uh1(k)=conjg(zu(k,1))
	  vh1(k)=zv(k,1)*bb(k)
c	  totvar=totvar+ssvd(k)**2
	end do
c  there is another conjugation in the inverse FFT -- cancels this one
	do l=1,2    
	  rf(n,l)=conjg(vh1(l)/vh1(3))
c  this formula is by analogy to the ordinary coherence formula in recfunk
	  crf(n,l)=(cabs(vh1(l)*vh1(3))*ssvd(1)**2)**2/(tnor(l)*tnor(3))
          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)
          rfs(n,irec,l)=rf(n,l)
          drfs(n,irec,l)=drf(n,l)
        end do
      end do
      sq2=sqrt(2.)
      go to 10
  111 continue
      close(10)
c  plot the inferred angle between radial and vertical against epicentral dist
      do i=1,irec
        if(epic(i).lt.80.) then
	  bb(i)=42.+(20.-epic(i))*(20./60.)
	elseif(epic(i).lt.118.) then
	  bb(i)=16.+(6./900.)*(110.-epic(i))**2
	else
	  bb(i)=9.0
	endif
      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
      call plotit_axes(0.,0.,0.,0.)
      call plotit(fr,crft(1,1),dum,nf,'Aggregate SVD Coherence',
     x    'Freq(Hz)','Radial Coherence',1,0,0.0,0,21)
      call plotit_axes(0.,0.,0.,0.)
      call plotit(fr,crft(1,2),dum,nf,'Aggregate SVD Coherence',
     x    'Freq(Hz)','Transverse Coherence',1,0,0.0,0,22)
      print *,irec
  484 continue
      open(12,file='outr_baz.grid',form='formatted')
      open(13,file='outt_baz.grid',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
c      naz0=1
      maz=0
c  ndf is a freq-spacing for stepping thru the spectral arrays at approximately the 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(bazinc.gt.0) bazinc=-bazinc
      endif
      abazinc=abs(bazinc)
      do baz=bazmax, bazmin, bazinc
c      do baz=0,355,5
        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
          if(l.eq.2) print *,baz,jrec
          if(jrec.ge.ibskip) then
            kaz=kaz+1
c  THIS commented section of code will compute the chi-squared misfit of the frequency-domain RF estimate
c  an a posteriori check on the statistics can be performed
c  basic result is that the median chi-squared misfit is fairly well obtained,
c  but the histogram of chi-squared values is far broader than predicted by chi-sq distribution with 
c  N degrees of gfreedom.  A smaller DOF seems OK, suggesting that the misfits are actually correlated
c  The most likely explanation is that the misfits are correlated in ways that follow later-arriving 
c  Ps converted phases.  You can find the later-arriving phases with different stacks of the single-event RFs
cc  compute chi-squared misfit, avoid zero-freq
c            ndf1=max0(ndf,nf1)
cc	    print *,ndf1,nf2,ndf
c	    do n=ndf1,nf2,ndf
c	      naz=naz+1
c	      chisq=0.
c              do j=1,jrec
c                sig=drfs(n,ip(j),l)**2
c	        chisq=chisq+cabs(rf(n,l)/drf(n,l)-rfs(n,ip(j),l))**2/sig
c	      end do
c  	      chi(naz)=chisq
c	      bbaz(naz)=baz
c	      nchi=2*jrec-2
c	      enchi(naz)=nchi
c	      if(nchi.gt.0) write(14,114) nchi,chisq,baz,l
cc	      if(chisq.gt.1000.) then
cc	        print *,jrec,chisq,baz,n
cc		pause
cc	      endif
c	    end do
c  find median misfit across all frequencies
c            nnaz=naz-naz0+1
c	    nnazh=nnaz/2
c	    do nn=naz0,naz
c	      ij=0
c	      do nnn=naz0,naz
c	        if(chi(nnn).le.chi(nn)) ij=ij+1
c	      end do
cc	      if(ij.eq.nnazh.and.jrec.gt.1.and.l.eq.2) then
c	      if(ij.eq.nnazh.and.jrec.gt.1) then
c	        if(l.eq.1) maz=maz+1
c	        chim(maz,l)=chi(nn)
c		enchim(maz,l)=enchi(nn)
c		go to 787
c	      endif
c	    end do
c  787       continue	    
c	    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
            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  mult by 4/3 if cosine**2 taper is applied to second half of 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
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 amplified
            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      open(9,file='fig7_baz_sv.dat',form='formatted')
c      do i=1,maz
c        write(9,*) enchim(i,1),chim(i,1)
c      end do
c      close(9)
c      open(9,file='fig7_baz_sh.dat',form='formatted')
c      do i=1,maz
c        write(9,*) enchim(i,2),chim(i,2)
c      end do
c      close(9)
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      ick=0
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
      open(12,file='outr_epi.grid',form='formatted')
      open(13,file='outt_epi.grid',form='formatted')
c      open(14,file='out_epi.chisq',form='formatted')
c  code for changing the epicentral bin width
      ep1=0.
      ep2=120.
      dep=10.
      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
c      naz0=1
      maz=0
c      open(9,file='fig6.dat',form='formatted')
c      open(19,file='fig6avg.dat',form='formatted')
      print *,ep1,ep2,dep
      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
c	    write(11,*) epic(i)
            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
		    fff=fr(1)+(n-1)*fr(2)
                    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
c	  	    if(abs(epi-85).lt.0.1.and.l.eq.1) then
c		      write(9,*) fff,rfs(n,i,l),drfs(n,i,l)
c		    endif
                  end do  
                endif
              endif
c            endif
          end do
          if(jrec.ge.ibskip) then
            if(l.eq.2) print *,epi,jrec
            kaz=kaz+1
c  compute chi-squared misfit, avoid zero-freq
c            ndf1=max0(ndf,nf1)
c	    do n=ndf1,nf2,ndf
c	      naz=naz+1
c	      chisq=0.
c              do j=1,jrec
c                sig=drfs(n,ip(j),l)**2
c	        chisq=chisq+cabs(rf(n,l)/drf(n,l)-rfs(n,ip(j),l))**2/sig
c	      end do
c  	      chi(naz)=chisq
c	      bbaz(naz)=baz
c	      nchi=2*jrec-2
c	      enchi(naz)=nchi
c	      if(nchi.gt.0) write(14,114) nchi,chisq,epi,l
c            end do
c  find median misfit across all frequencies
c            nnaz=naz-naz0+1
c	    nnazh=nnaz/2
c	    do nn=naz0,naz
c	      ij=0
c	      do nnn=naz0,naz
c	        if(chi(nnn).le.chi(nn)) ij=ij+1
c	      end do
cc	      if(ij.eq.nnazh.and.jrec.gt.1.and.l.eq.2) then
c	      if(ij.eq.nnazh.and.jrec.gt.1) then
c	        if(l.eq.1) maz=maz+1
c	        chim(maz,l)=chi(nn)
c		enchim(maz,l)=enchi(nn)
c		go to 797
c	      endif
c	    end do
c  797       continue	    
c	    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)
c	      if(abs(epi-85).lt.0.1.and.l.eq.1) then
c		fff=fr(1)+(n-1)*fr(2)
c		write(19,*) fff,rf(n,l)
c              endif	      
              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
            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))
            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  
c      close(9)
c      close(19)
 1022 format(f7.3,f6.1,f7.3)
      close(12)
      close(13)
      close(14)
c      enmax=1
c      open(9,file='fig7_epi_sv.dat',form='formatted')
c      do i=1,maz
c        write(9,*) enchim(i,1),chim(i,1)
c      end do
c      close(9)
c      open(9,file='fig7_epi_sh.dat',form='formatted')
c      do i=1,maz
c        write(9,*) enchim(i,2),chim(i,2)
c      end do
c      close(9)
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
      ick=0
      if(ick.ne.0) then
        print *,fmin,fmax
        print *,'enter new fmin & fmax (Hz)'
	read(5,*) fmin,fmax
	fmax=amin1(fmax,fmax0)
	nf1=max0(int(fmin/ddf),1)
	nf2=fmax/ddf+1
	nf=nf2-nf1+1
	print *,fmin,fmax,nf1,nf2,ddf,nf
	go to 494
      endif
      kaz=kaz/2
      print *,kaz,' traces'   
  109 format('Back Azimuth Centered on ',f4.0,' degrees')
  114 format(i5,g15.4,f6.0,i5)
 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(262144)  ! 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 CSVD (A,MMAX,NMAX,M,N,IP,NU,NV,S,U,V)
C
C   Singular value decomposition of an  M by N  complex matrix  A,
C   where M .GT. N .  The singular values are stored in the vector
C   S. The first NU columns of the M by M unitary matrix U and the
C   first NV columns of the N by N unitary matrix V  that minimize
C   Det(A-USV*) are also computed.
C
C
C         P.A. Businger and G.H. Golub, "Singular Value Decomposition
C         of a Complex Matrix," Communications of the ACM, vol. 12,
C         pp. 564-565, October 1969.
C
C   This algorithm is reprinted by permission, Association for
C   Computing Machinery; copyright 1969.
C
      COMPLEX *16 A(MMAX,NMAX),U(MMAX,MMAX),V(NMAX,NMAX),Q,R
      REAL *8 S(NMAX),B(2000),C(2000),T(2000),zero,one,ETA,TOL,
     $        EPS
      ETA = 1.2d-7
      TOL = 2.4d-32
      zero = 0.d0
      one = 1.d0
      NP=N+IP
      N1=N+1
C   Householder reduction
      C(1)=zero
      K=1
10    K1=K+1
C   Elimination of A(I,K) , I=K+1,...,M
      Z=zero
      DO 20 I=K,M
20      Z=Z+dreal(A(I,K))**2+dIMAG(A(I,K))**2
      B(K)=zero
      IF (Z .LE. TOL)  GO TO 70
      Z=SQRT(Z)
      B(K)=Z
      W=CdABS(A(K,K))
      Q=dcmplx(one,zero)
      IF (W .NE. zero)  Q=A(K,K)/W
      A(K,K)=Q*(Z+W)
      IF (K .EQ. NP)  GO TO 70
      DO 50 J=K1,NP
        Q=dcmplx(zero,zero)
        DO 30 I=K,M
30        Q=Q+dCONJG(A(I,K))*A(I,J)
        Q=Q/(Z*(Z+W))
        DO 40 I=K,M
40        A(I,J)=A(I,J)-Q*A(I,K)
50      CONTINUE
C   Phase transformation
      Q=-dCONJG(A(K,K))/CdABS(A(K,K))
      DO 60 J=K1,NP
60      A(K,J)=Q*A(K,J)
C   Elimination of A(K,J) , J=K+2,...,N
70    IF (K .EQ. N)  GO TO 140
      Z=zero
      DO 80 J=K1,N
80      Z=Z+dreal(A(K,J))**2+dIMAG(A(K,J))**2
      C(K1)=zero
      IF (Z .LE. TOL)  GO TO 130
      Z=SQRT(Z)
      C(K1)=Z
      W=CdABS(A(K,K1))
      Q=dcmplx(one,zero)
      IF (W .NE. zero)  Q=A(K,K1)/W
      A(K,K1)=Q*(Z+W)
      DO 110 I=K1,M
        Q=dcmplx(zero,zero)
        DO 90 J=K1,N
90        Q=Q+dCONJG(A(K,J))*A(I,J)
        Q=Q/(Z*(Z+W))
        DO 100 J=K1,N
100       A(I,J)=A(I,J)-Q*A(K,J)
110     CONTINUE
C   Phase transformation
      Q=-dCONJG(A(K,K1))/CdABS(A(K,K1))
      DO 120 I=K1,M
120     A(I,K1)=A(I,K1)*Q
130   K=K1
      GO TO 10
C   Tolerance for negligible elements
140   EPS=zero
      DO 150 K=1,N
        S(K)=B(K)
        T(K)=C(K)
150   EPS=DMAX1(EPS,S(K)+T(K))
      EPS=EPS*ETA
C   Initialization of U and V
      IF (NU .EQ. 0)  GO TO 180
      DO 170 J=1,NU
        DO 160 I=1,M
160       U(I,J)=dcmplx(zero,zero)
170     U(J,J)=dcmplx(one,zero)
180   IF (NV .EQ. 0)  GO TO 210
      DO 200 J=1,NV
        DO 190 I=1,N
190       V(I,J)=dcmplx(zero,zero)
200     V(J,J)=dcmplx(one,zero)
C   QR diagonalization
210   DO 380 KK=1,N
        K=N1-KK
C   Test for split
220     DO 230 LL=1,K
          L=K+1-LL
          IF (ABS(T(L)) .LE. EPS)  GO TO 290
          IF (ABS(S(L-1)) .LE. EPS)  GO TO 240
230       CONTINUE
C   Cancellation of B(L)
240     CS=zero
        SN=one
        L1=L-1
        DO 280 I=L,K
          F=SN*T(I)
          T(I)=CS*T(I)
          IF (ABS(F) .LE. EPS)  GO TO 290
          H=S(I)
          W=SQRT(F*F+H*H)
          S(I)=W
          CS=H/W
          SN=-F/W
          IF (NU .EQ. 0)  GO TO 260
          DO 250 J=1,N
            X=dreal(U(J,L1))
            Y=dreal(U(J,I))
            U(J,L1)=dCMPLX(X*CS+Y*SN,0.)
250         U(J,I)=dCMPLX(Y*CS-X*SN,0.)
260       IF (NP .EQ. N)  GO TO 280
          DO 270 J=N1,NP
            Q=A(L1,J)
            R=A(I,J)
            A(L1,J)=Q*CS+R*SN
270         A(I,J)=R*CS-Q*SN
280       CONTINUE
C   Test for convergence
290     W=S(K)
        IF (L .EQ. K)  GO TO 360
C   Origin shift
        X=S(L)
        Y=S(K-1)
        G=T(K-1)
        H=T(K)
        F=((Y-W)*(Y+W)+(G-H)*(G+H))/(2.d0*H*Y)
        G=SQRT(F*F+one)
        IF (F .LT. zero)  G=-G
        F=((X-W)*(X+W)+(Y/(F+G)-H)*H)/X
C   QR step
        CS=one
        SN=one
        L1=L+1
        DO 350 I=L1,K
          G=T(I)
          Y=S(I)
          H=SN*G
          G=CS*G
          W=SQRT(H*H+F*F)
          T(I-1)=W
          CS=F/W
          SN=H/W
          F=X*CS+G*SN
          G=G*CS-X*SN
          H=Y*SN
          Y=Y*CS
          IF (NV .EQ. 0)  GO TO 310
          DO 300 J=1,N
            X=dreal(V(J,I-1))
            W=dreal(V(J,I))
            V(J,I-1)=dCMPLX(X*CS+W*SN,0.)
300         V(J,I)=dCMPLX(W*CS-X*SN,0.)
310       W=SQRT(H*H+F*F)
          S(I-1)=W
          CS=F/W
          SN=H/W
          F=CS*G+SN*Y
          X=CS*Y-SN*G
          IF (NU .EQ. 0)  GO TO 330
          DO 320 J=1,N
            Y=dreal(U(J,I-1))
            W=dreal(U(J,I))
            U(J,I-1)=dCMPLX(Y*CS+W*SN,0.)
320         U(J,I)=dCMPLX(W*CS-Y*SN,0.)
330       IF (N .EQ. NP)  GO TO 350
          DO 340 J=N1,NP
            Q=A(I-1,J)
            R=A(I,J)
            A(I-1,J)=Q*CS+R*SN
340         A(I,J)=R*CS-Q*SN
350       CONTINUE
        T(L)=zero
        T(K)=F
        S(K)=X
        GO TO 220
C   Convergence
360     IF (W .GE. zero)  GO TO 380
        S(K)=-W
        IF (NV .EQ. 0)  GO TO 380
        DO 370 J=1,N
370       V(J,K)=-V(J,K)
380     CONTINUE
C   Sort singular values
      DO 450 K=1,N
        G=-one
        J=K
        DO 390 I=K,N
          IF (S(I) .LE. G)  GO TO 390
          G=S(I)
          J=I
390       CONTINUE
        IF (J .EQ. K)  GO TO 450
        S(J)=S(K)
        S(K)=G
        IF (NV .EQ. 0)  GO TO 410
        DO 400 I=1,N
          Q=V(I,J)
          V(I,J)=V(I,K)
400       V(I,K)=Q
410     IF (NU .EQ. 0)  GO TO 430
        DO 420 I=1,N
          Q=U(I,J)
          U(I,J)=U(I,K)
420       U(I,K)=Q
430     IF (N .EQ. NP)  GO TO 450
        DO 440 I=N1,NP
          Q=A(J,I)
          A(J,I)=A(K,I)
440       A(K,I)=Q
450     CONTINUE
C   Back transformation
      IF (NU .EQ. 0)  GO TO 510
      DO 500 KK=1,N
        K=N1-KK
        IF (B(K) .EQ. zero)  GO TO 500
        Q=-A(K,K)/CdABS(A(K,K))
        DO 460 J=1,NU
460       U(K,J)=Q*U(K,J)
        DO 490 J=1,NU
          Q=dcmplx(zero,zero)
          DO 470 I=K,M
470         Q=Q+dCONJG(A(I,K))*U(I,J)
          Q=Q/(CdABS(A(K,K))*B(K))
          DO 480 I=K,M
480         U(I,J)=U(I,J)-Q*A(I,K)
490       CONTINUE
500     CONTINUE
510   IF (NV .EQ. 0)  GO TO 570
      IF (N .LT. 2)  GO TO 570
      DO 560 KK=2,N
        K=N1-KK
        K1=K+1
        IF (C(K1) .EQ. zero)  GO TO 560
        Q=-dCONJG(A(K,K1))/CdABS(A(K,K1))
        DO 520 J=1,NV
520       V(K1,J)=Q*V(K1,J)
        DO 550 J=1,NV
          Q=dcmplx(zero,zero)
          DO 530 I=K1,N
530         Q=Q+A(K,I)*V(I,J)
          Q=Q/(CdABS(A(K,K1))*C(K1))
          DO 540 I=K1,N
540         V(I,J)=V(I,J)-Q*dCONJG(A(K,I))
550       CONTINUE
560     CONTINUE
570   RETURN
      END