      subroutine calc_cfr(npts,dt,ind,data,cfr)
c
c     subroutine for estimation of various characteristic functions
c
c     input:
c             npts   - number of points in time series
c             dt     - sampling rate in seconds
c             ind    - index of characteristic function of choice
c             data   - array containing 3 component time series
c     output:
c             cfr    - characteristic function
c
      parameter (maxp  = 65536)
c
      real*4 data(maxp,3)
      real*4 xn(maxp), yn(maxp), zn(maxp)
      real*4 ang(maxp,2), rn(maxp)
      real*4 cfr(*)
      real*4 k1, k2, k4, k8
c
      pi = 4.0*atan(1.0)
      cffc = 0.2
      k2 = 0.5/(dt*pi*cffc)
      k8 = 0.5
c
c     initialize values
c
      cfr(1) = 0.0
      sum1a = 0.0
      sum1b = 0.0
      sum4a = 0.0
      sum4b = 0.0
      savg = 0.0
      svar = 0.0
c
c Rex Allen
c
      if (ind .eq. 1) then
        do ii = 2, npts
          xx = data(ii,1)
          yy = data(ii,2)
          zz = data(ii,3)
          xdt = data(ii,1) - data(ii-1,1)
          ydt = data(ii,2) - data(ii-1,2)
          zdt = data(ii,3) - data(ii-1,3)
          sum1a = sum1a + abs(xx) + abs(yy) + abs(zz)
          sum1b = sum1b + abs(xdt) + abs(ydt) + abs(zdt)
          k1 = sum1a/sum1b
          cfr(ii) = xx*xx + yy*yy + zz*zz 
     .              + k1*(xdt*xdt + ydt*ydt + zdt*zdt)
        end do
c
c Bob Uhrhammer
c
      elseif (ind .eq. 2) then
        do ii = 2, npts
          xx = data(ii,1)
          yy = data(ii,2)
          zz = data(ii,3)
          xdt = data(ii,1) - data(ii-1,1)
          ydt = data(ii,2) - data(ii-1,2)
          zdt = data(ii,3) - data(ii-1,3)
          cfr(ii) = xx*xx + yy*yy + zz*zz
     .              + k2*(xdt*xdt + ydt*ydt + zdt*zdt)
        end do
c
c Rex Allen **2
c
      elseif (ind .eq. 3) then
        do ii = 2, npts
          xx = data(ii,1)
          yy = data(ii,2)
          zz = data(ii,3)
          xdt = data(ii,1) - data(ii-1,1)
          ydt = data(ii,2) - data(ii-1,2)
          zdt = data(ii,3) - data(ii-1,3)
          sum1a = sum1a + abs(xx) + abs(yy) + abs(zz)
          sum1b = sum1b + abs(xdt) + abs(ydt) + abs(zdt)
          k1 = sum1a/sum1b
          c1 = xx*xx + yy*yy + zz*zz
     .              + k1*(xdt*xdt + ydt*ydt + zdt*zdt)
          cfr(ii) = c1 * c1
        end do
c
c Baer envelop function
c
      elseif (ind .eq. 4) then
        do ii = 2, npts
          xx = data(ii,1)
          yy = data(ii,2)
          zz = data(ii,3)
          xdt = data(ii,1) - data(ii-1,1)
          ydt = data(ii,2) - data(ii-1,2)
          zdt = data(ii,3) - data(ii-1,3)
          sum4a = sum4a + xx*xx + yy*yy + zz*zz
          sum4b = sum4b + xdt*xdt + ydt*ydt + zdt*zdt
          k4 = sum4a/sum4b
          cfr(ii) = xx*xx + yy*yy + zz*zz
     .              + k4*(xdt*xdt + ydt*ydt + zdt*zdt)
        end do
c
c Baer preferred function
c -- must contain an error - to ugly to be true
c
      elseif (ind .eq. 5) then
        do ii = 2, npts 
          xx = data(ii,1) 
          yy = data(ii,2) 
          zz = data(ii,3) 
          xdt = data(ii,1) - data(ii-1,1) 
          ydt = data(ii,2) - data(ii-1,2) 
          zdt = data(ii,3) - data(ii-1,3) 
          sum4a = sum4a + xx*xx + yy*yy + zz*zz
          sum4b = sum4b + xdt*xdt + ydt*ydt + zdt*zdt   
          k4 = sum4a/sum4b
          c4 = xx*xx + yy*yy + zz*zz   
     .         + k4*(xdt*xdt + ydt*ydt + zdt*zdt) 
          sf = c4 * c4
          savg = savg + sf
          sfavg = savg/(3*ii)
          svar = svar + (sf - sfavg)**2
          sfvar = sqrt(svar)
          cfr(ii) = (sf - sfavg)/sfvar
        end do 
c
c Sarrin
c
      elseif (ind .eq. 6) then
        do ii = 2, npts  
          xx = data(ii,1)  
          yy = data(ii,2)  
          zz = data(ii,3)  
          cfr(ii) = abs(xx*yy*zz)
        end do
c
c  Uhrhammer angle function
c
      elseif (ind .eq. 7) then
        do ii = 1, npts
          xx = data(ii,1)  
          yy = data(ii,2)  
          zz = data(ii,3)  
          fmag = sqrt(xx*xx + yy*yy + zz*zz)
          rn(ii) = fmag
          xn(ii) = xx/fmag
          yn(ii) = yy/fmag
          zn(ii) = zz/fmag
        end do
        do ii = 2, npts - 1
          xnm = xn(ii-1)
          ynm = yn(ii-1)
          znm = zn(ii-1)
          xn0 = xn(ii)
          yn0 = yn(ii)
          zn0 = zn(ii)
          xnp = xn(ii+1)
          ynp = yn(ii+1)
          znp = zn(ii+1)
          dot_p = xnm*xnp + ynm*ynp + znm*znp
          ang(ii,1) = 0.5*asin(dot_p)
c
          xm=ynm*zn0-yn0*znm
          ym=znm*xn0-zn0*xnm
          zm=xnm*yn0-xn0*ynm
          xp=yn0*znp-ynp*zn0
          yp=zn0*xnp-znp*xn0
          zp=xn0*ynp-xnp*yn0
          cross_p = xm*xp + ym*yp + zm*zp
          ang(ii,2) = asin(cross_p)
        end do
        ang(1,1) = ang(2,1)
        ang(1,2) = ang(2,2)
        ang(npts,1) = ang(npts-1,1)
        ang(npts,2) = ang(npts-1,2)
c
        fx = 1.0/dt
        fx2 = fx*fx
        do ii = 2, npts - 1
          da1dt = ang(ii,1)*fx
          da2dt = ang(ii,2)*fx
          da1dt2 = (ang(ii+1,1)-2.0*ang(ii,1)+ang(ii-1,1))*fx2
          da2dt2 = (ang(ii+1,2)-2.0*ang(ii,2)+ang(ii-1,2))*fx2
          cfr(ii) = rn(ii)*(da1dt*da1dt+da2dt*da2dt
     .                + k2*(da1dt2*da1dt2+da2dt2*da2dt2))
        end do
        cfr(1) = cfr(2)
        cfr(npts) = cfr(npts-1)
c
c     Uhrhammer combined
c
      elseif (ind .eq. 8) then
        do ii = 1, npts
          xx = data(ii,1)
          yy = data(ii,2)
          zz = data(ii,3)
          fmag = sqrt(xx*xx + yy*yy + zz*zz)
          rn(ii) = fmag
          xn(ii) = xx/fmag
          yn(ii) = yy/fmag
          zn(ii) = zz/fmag
        end do
c
        do ii = 2, npts
          xnm = xn(ii-1)
          ynm = yn(ii-1)
          znm = zn(ii-1)
          xn0 = xn(ii)
          yn0 = yn(ii)
          zn0 = zn(ii)
          xnp = xn(ii+1)
          ynp = yn(ii+1)
          znp = zn(ii+1)
          dot_p = xnm*xnp + ynm*ynp + znm*znp
          ang(ii,1) = 0.5*asin(dot_p)
          xm=ynm*zn0-yn0*znm
          ym=znm*xn0-zn0*xnm
          zm=xnm*yn0-xn0*ynm
          xp=yn0*znp-ynp*zn0
          yp=zn0*xnp-znp*xn0
          zp=xn0*ynp-xnp*yn0
          cross_p = xm*xp + ym*yp + zm*zp
          ang(ii,2) = asin(cross_p)
        end do
        ang(1,1) = ang(2,1)
        ang(1,2) = ang(2,2)
        ang(npts,1) = ang(npts-1,1)
        ang(npts,2) = ang(npts-1,2)
c
        fx = 1.0/dt
        fx2 = fx*fx
        do ii = 2, npts - 1
          xx = data(ii,1)
          yy = data(ii,2)
          zz = data(ii,3)

          xdt = data(ii,1) - data(ii-1,1)
          ydt = data(ii,2) - data(ii-1,2)
          zdt = data(ii,3) - data(ii-1,3)
          cfr1 = xx*xx + yy*yy + zz*zz
     .              + k2*(xdt*xdt + ydt*ydt + zdt*zdt)

          da1dt = ang(ii,1)*fx
          da2dt = ang(ii,2)*fx
          da1dt2 = (ang(ii+1,1)-2.0*ang(ii,1)+ang(ii-1,1))*fx2
          da2dt2 = (ang(ii+1,2)-2.0*ang(ii,2)+ang(ii-1,2))*fx2
          cfr2 = rn(ii)*(da1dt*da1dt+da2dt*da2dt
     .                + k2*(da1dt2*da1dt2+da2dt2*da2dt2))

          cfr(ii) = k8*cfr1 + (1.0 - k8)*cfr2
        end do
        cfr(1) = cfr(2)
        cfr(npts) = cfr(npts-1)

      end if
c
c     all done
c
      return
      end
