      subroutine est_cfr(head,idc,dtr,dtc,d1,dw,iem,ic,
     &   file_pick,io_out,ierr)
c
c     subroutine for estimation of various characteristic functions
c
c     input:
c             head   - station/channel header
c             idc    - characteristic function index
c             dtr    - detection threshold
c             dtc    - detection time constant
c             d1     - data array
c             dw     - spectral array
c             iem    - crude component counter
c             io_out - unit number for writing
c     output:
c             iere   - error return code
c
      parameter (maxp  = 65536)
c
c     include files for header and time
c
      include '/usr/local/include/qlib2.inc'
      record /INT_TIME/ t0_time, t1_time
      record /EXT_TIME/ t2_time
      include 'station_hdr.inc'
c
      real*4 d1(*)
      real*4 data(maxp,3), a0(3)
      real*4 cfr(maxp)
      real*4 lta(maxp), sta(maxp)
      real*4 azm1(maxp), aoi1(maxp)
c
      real*8 tslew
c
      integer*4 ic(*)
c
      complex dw(*)
c
      character*1 ipo, ipt, ipp, istrc, cid
      character*3 chan
      character*4 stat
      character*6 ipd
      character*(*) file_pick
c
      common /filter/ np0, fl0, fh0
c
      data icall /0/
      data ipo /'i'/
      data ipt /'P'/
      data ipp /'Z'/
      data ipd /'0.0000'/
c
c     initialize error flag
c
      iere = 0
c
c     set constants for sta/lta determination
c
      c3 = 0.67
      c4 = 0.015
      c5 = dtr
      c6 = 0.05*dtr
      indc = idc
c
c     set variables for filtering
c
      pi = 4.0*atan(1.0)
      ind = 0
      ind2 = 7
      if (np0 .eq. 0) then
        ind2 = 0
      else
        npbb = np0
        fl = fl0
        fh = fh0
      end if
c
c     station parameters
c
      stat  = head.station_id
      chan  = head.channel_id
      npts  = head.num_pts
      dt  = 1.0/head.sample_rate
c
c     save this data
c
      do ii = 1, npts
        data(ii,iem) = d1(ii)
      end do
      a0(iem) = head.gain
c
c     save reference time if this is the first call
c
      if (icall .eq. 0) then
        t0_time = head.begtime
        icall = icall + 1
      end if
c
c     if we don't have three-components yet, return
c
      if (iem .ne. 3) then
        return
      else
        do ii = 1, 3
          if (ic(ii) .ne. 1) then
            iere = -iem
            return
          end if
        end do
      end if
c
c     open an output file
c
      inquire(file=file_pick,exist=lexist) 
      if (lexist) then 
        open(unit=io_out,file=file_pick,fileopt='eof') 
      else
        open(unit=io_out,file=file_pick)
        write(io_out,'(a13,i1,1x,2(f8.3,1x),i2,1x,3(f8.3,1x))')
     .   'Pickfile V.2 ', idc,dtr,dtc,npbb,fl,fh
        call f_int_to_ext(t0_time,t2_time)
        write(io_out,101) 'T ',t2_time.year, t2_time.doy,
     .     t2_time.hour, t2_time.minute, t2_time.second,  
     .     t2_time.usec
101     format(a,i4.4,1x,i3.3,1x,i2.2,":",i2.2,":",i2.2,".",i4.4)
      end if
c      
c     correct for instrument gain factor and filter
c
      gain = (a0(1) + a0(2) + a0(3))/3.0
      a1 = gain/a0(1)
      a2 = gain/a0(2)
      a3 = gain/a0(3)
c
      do ii = 1, npts
        data(ii,1) = data(ii,1)*a1
        data(ii,2) = data(ii,2)*a2
        data(ii,3) = data(ii,3)*a3
      end do
c
      call coinst2(npts,dt,ind,ind2,data(1,1),dw,npbb,fl,fh,iere)
      if (iere .ne. 0) then
        return
      end if
c
      call coinst2(npts,dt,ind,ind2,data(1,2),dw,npbb,fl,fh,iere)
      if (iere .ne. 0) then
        return
      end if
c
      call coinst2(npts,dt,ind,ind2,data(1,3),dw,npbb,fl,fh,iere)
      if (iere .ne. 0) then
        return
      end if
c
c     form appropriate characteristic functions
c
      sta(1) = 0.0
      lta(1) = 0.0
      call calc_cfr(npts,dt,indc,data,cfr)
c
c     now form the sta/lta functions for detection
c
      do ii = 2, npts
        sta(ii) = sta(ii-1) + c3*(cfr(ii) - sta(ii-1))
        lta(ii) = lta(ii-1) + c4*(cfr(ii) - lta(ii-1))
      end do
c
c     compute the azim and angle of incidence this set
c
      call calc_azm(npts,data,azm1,aoi1)
c
c     check for detections and write them out
c     --ignore the part of the time series which was tapered
c
      fps = npts * 5.0 * 0.005
      j1 = int(fps)
      j2 = npts - j1
c
      i1 = 0
      ftx = dt/dtc
      do ii = j1, j2
c
c       this helps prevent immediate redetections
c
c       if (i1 .ne. 0) then
c         dtx = (ii - i1)*ftx
c         alpha = 1.0 - exp(-dtx)
c       else
c         alpha = 1.0
c       end if
c
        if (i1 .ne. 0) then
          dtx = (ii - i1)*ftx
          if (dtx .lt. 1.0) then
            alpha = dtx*dtx*dtx
          else
            alpha = 1.0
          end if
        else
          alpha = 1.0
        end if
c
c       go back and try to find something early
c
        if (alpha*sta(ii) .gt. c5*lta(ii)) then
          k2 = ii
c
c         define an band for refining the estimate
c
          do kk = k2, i1 + 1, -1
            if (alpha*sta(kk) .lt. c6*lta(k2)) then
              k1 = kk
              go to 150
            end if
            k1 = k2 - 1
          end do
150       continue
c
          call find_onset(k1,k2,cfr,k0)
          call find_azm(k0,k2,npa,azm1,aoi1,azm,seazm,aaoim,seaaoim)
          cflevel = cfr(k0)
c
c         compute travel time based on current header time
c
          i1 = k0
          t1 = (k0 - 1)*dt
          it1s = int(t1)
          it1t = (t1 - it1s)*usecs_per_sec
          call f_add_time(head.begtime, it1s, it1t, t1_time)
          call f_int_to_ext(t1_time,t2_time)
c
c         compute time offset to reference time
c
          tslew = f_tdiff(t1_time,t0_time)
          t1 = tslew/usecs_per_sec
          k0 = t1/dt + 1          
c
c         what have we picked?
c
          if(aaoim .le. 45.0) then
            ipcnt=ipcnt+1
            ipt='P'//' '
            ipp='U'
            istrc='Z'
          elseif(aaoim .ge. 135.0) then
            ipcnt=ipcnt+1
            ipt='P'//' '
            ipp='D'
            istrc='Z'
          else
            iscnt=iscnt+1
            ipt='S'//' '
            if(azm .le. 45.0) then
              ipp='N'
              istrc='N'
            elseif(azm .le. 135.0) then
              ipp='E'
              istrc='E'
            elseif(azm .le. 225.0) then
              ipp='S'
              istrc='N'
            elseif(azm .le. 315.0) then
              ipp='W'
              istrc='E'
            else
              ipp='N'
              istrc='N'
            endif
          endif
          write(io_out,200) stat,chan(1:2),istrc,
     .      k0,ipo,ipt,t1,ipp,ipd,
     .      t2_time.year, t2_time.doy, t2_time.hour, 
     .      ':',t2_time.minute, ':', t2_time.second, '.', 
     .      t2_time.usec,
     .      npa,cflevel,azm,aaoim
200       format(1x,a4,1x,a2,a1,
     .      1x,i5,1x,a1,1x,a1,1x,f12.4,1x,a1,1x,a6,1x,
     .      i4.4,1x,i3.3,1x,i2.2,a1,i2.2,a1,i2.2,a1,i4.4,1x,
     .      i3,1x,e15.5,1x,2(f8.2,1x))
        end if
      end do
      close(io_out)
c
c     all done
c
      return
      end
