      subroutine dataft(nsect,jj,nf1,nfinc,nfsect,nfft,noff,
     $                  nfil,nar,deltat,dpss,wsave,rdata,fpar,cpar,
     $                  arf,fft,ffts)
c
c  birrp
c
c  computes fourier transforms of data, corrects for prewhitening
c  and instrument response, and returns selected frequencies
c
c  input variables
c    nsect-number of sections
c    jj-series index in ffts
c    nf1-starting frequency index
c    nfinc-frequency index increment
c    nfsect-number of frequencies per section starting with nf1
c    nfft-length of fft
c    noff-offset between sections
c    nfil-number of filter parameters in fpar if .gt.0 and
c         indicates filter parameters are contained in file cpar if .lt.0
c    nar-length of ar filter less leading 1
c    deltat-data sample interval
c    dpss-double precision vector containing slepian sequence data window
c    wsave-vector containing sinusoids for fft
c    rdata-vector containing data sequence
c    fpar-vector of filter parameters
c    cpar-character string of filter parameters
c    arf-vector of ar prewhitening filter coefficients
c    fft-complex work space of length nfft
c
c  output variables
c    ffts-complex array containing fourier transforms for each section
c         and component at selected frequencies
c    
      include 'parameters.h'
      parameter (narm=100,nfilm=47)
      parameter (nserm=noutm+ninpm+nrefm)
      character cpar*(*)
      double precision dpss,darf,pi,tpi,om,ct,st
      complex fft,ffts,ar,fil,filter,t0,t1
      dimension dpss(*),wsave(*),rdata(*),fpar(*),arf(*),
     $          fft(*),ffts(nsectm,nserm,*)
      dimension ffreq(nfilm),ffreq1(nfilm),rfpar(nfilm),
     $          afpar(nfilm),s(2*nfilm,6),rbrk(2*nfilm),
     $          rcoef(4,2*nfilm),abrk(2*nfilm),acoef(4,2*nfilm)
      dimension darf(narm+1)
      parameter (pi=3.141592653589793d0,tpi=2.d0*pi)
      write(*,*)cpar
      xnfft=nfft
      if(nar.ne.0)then
        darf(1)=1.d0
        do 1 ii=2,nar+1
    1     darf(ii)=arf(ii-1)
      endif
      if(nfil.lt.0)then
        call strlen(cpar,i1,i2)
        open(unit=10,file=cpar(i1:i2),status='old')
        read(10,*)scal
        do 5 i=1,nfilm
    5     read(10,*,end=6)ffreq(i),rfpar(i),afpar(i)
    6   nfi=i-1
        close(unit=10)
        do 7 i=1,nfi
    7     ffreq1(i)=ffreq(i)
        call ssort(ffreq,rfpar,nfi,2)

        call tautsp(ffreq,rfpar,nfi,5.5,s,rbrk,rcoef,
     $              lr,kr,iflag)

        write(*,*)"rbrk=",rbrk
        write(*,*)"rcoef=",rcoef
        call ssort(ffreq1,afpar,nfi,2)
        call tautsp(ffreq,afpar,nfi,5.5,s,abrk,acoef,
     $              la,ka,iflag)
      endif
      k=1
      nn=nsect/2
      do 30 n=1,2*nn,2
        do 10 i=1,nfft
          fft(i)=dpss(i)*cmplx(rdata(k),rdata(k+noff))
   10     k=k+1
        call cfftf(nfft,fft,wsave)
        kk=1
        do 20 i=nf1,nf1+(nfsect-1)*nfinc,nfinc
          if(nfil.gt.0)then
            fil=filter(float(i)/(deltat*xnfft),fpar,nfil+1)
          elseif(nfil.lt.0)then
            fr=float(i)/(deltat*xnfft)
            fil=cmplx(ppvalu(rbrk,rcoef,lr,kr,fr,0),
     $                ppvalu(abrk,acoef,la,ka,fr,0))
            if(nfil.eq.-2)fil=1./fil
            fil=scal*fil
          else
            fil=1.
          endif
          if(nar.ne.0)then
            om=tpi*dfloat(i)/dble(xnfft)
            call sft(darf,nar+1,om,ct,st)
            ar=cmplx(ct,st)
          else
            ar=1.
          endif
          t0=fft(i+1)
          t1=fft(nfft-i+1)
          ffts(n,jj,kk)=0.5*fil*(t0+conjg(t1))/ar
          ffts(n+1,jj,kk)=-cmplx(0.,0.5)*fil*(t0-
     $                     conjg(t1))/ar
          kk=kk+1
   20     continue
   30   k=k-nfft+2*noff
      if(mod(nsect,2).ne.0)then
        do 40 i=1,nfft
          fft(i)=dpss(i)*rdata(k)
   40     k=k+1
        call cfftf(nfft,fft,wsave)
        kk=1
        do 60 i=nf1,nf1+(nfsect-1)*nfinc,nfinc
          if(nfil.gt.0)then
            fil=filter(float(i)/(deltat*xnfft),fpar,nfil)
          elseif(nfil.lt.0)then
            fr=float(i)/(deltat*xnfft)
            fil=cmplx(ppvalu(rbrk,rcoef,lr,kr,fr,0),
     $                ppvalu(abrk,acoef,la,ka,fr,0))
          else
            fil=1.
          endif
          if(nar.ne.0)then
            om=tpi*dfloat(i)/dble(xnfft)
            call sft(darf,nar+1,om,ct,st)
            ar=cmplx(ct,st)
          else
            ar=1.
          endif
          ffts(nsect,jj,kk)=fil*fft(i+1)/ar
          kk=kk+1
   60     continue
      endif
      return
      end
      subroutine rotate(nsect,theta1,theta2,phi,fftx,ffty)
c
c  birrp
c
c  rotate ffts to new coordinate system
c
c  input variables:
c
c  nsect--number of sections
c  theta1--angle of first component pos cw wrt geomagnetic north
c  theta2--angle of second component pos cw wrt geomagnetic north
c  phi--pos cw rotation angle of geomagnetic coordinate system
c  fftx--complex array containing nsect section fft's for
c        x component (overwritten on return)
c  ffty--complex array containing nsect section fft's for
c        y component (overwritten on return)
c
c  output variables:
c
c  fftx--complex array containing nsect section fft's for x component
c  ffty--complex array containing nsect section fft's for y component
c
c  if original data are collected in orthogonal coordinate
c  system, set theta1=0, theta2=90
c
      complex fftx,ffty,t0
      dimension fftx(*),ffty(*)
      parameter (pi=3.141592653589793,rad=pi/180.)
      if(theta1.eq.0.)then
        spt1=sin(rad*phi)
        cpt1=cos(rad*phi)
      else
        spt1=sin(rad*(phi-theta1))
        cpt1=cos(rad*(phi-theta1))
      endif
      if(theta2.eq.90.)then
        spt2=-cos(rad*phi)
        cpt2=sin(rad*phi)
      else
        spt2=sin(rad*(phi-theta2))
        cpt2=cos(rad*(phi-theta2))
      endif
      st2t1=sin(rad*(theta2-theta1))
      do 10 i=1,nsect
        t0=(-fftx(i)*cpt2+ffty(i)*cpt1)/st2t1
        fftx(i)=(-fftx(i)*spt2+ffty(i)*spt1)/st2t1
 10     ffty(i)=t0
      return
      end
