	subroutine ezffti(n,wsave)
c       some minor modifications for ana interfaces, r. shine
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c                       version 2  february 1978
c
c          a package of fortran subprograms for the fast fourier
c           transform of periodic and other symmetric sequences
c
c                              by
c
c                       paul n swarztrauber
c
c       national center for atmospheric research  boulder,colorado 80307
c
c        which is sponsored by the national science foundation
c
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     this package consists of programs which perform fast fourier
c     transforms for both complex and real periodic sequences and
c     certain other symmetric sequences that are listed below.
c
c     1.   rffti     initialize  rfftf and rfftb
c     2.   rfftf     forward transform of a real periodic sequence
c     3.   rfftb     backward transform of a real coefficient array
c
c     4.   ezfftf    a simplified real periodic forward transform
c     5.   ezfftb    a simplified real periodic backward transform
c
c     6.   sinti     initialize sint
c     7.   sint      sine transform of a real odd sequence
c
c     8.   costi     initialize cost
c     9.   cost      cosine transform of a real even sequence
c
c     10.  sinqi     initialize sinqf and sinqb
c     11.  sinqf     forward sine transform with odd wave numbers
c     12.  sinqb     unnormalized inverse of sinqf
c
c     13.  cosqi     initialize cosqf and cosqb
c     14.  cosqf     forward cosine transform with odd wave numbers
c     15.  cosqb     unnormalized inverse of cosqf
c
c     16.  cffti     initialize cfftf and cfftb
c     17.  cfftf     forward transform of a complex periodic sequence
c     18.  cfftb     unnormalized inverse of cfftf
c
        implicit real*4 (a-h,o-z)
c---    8/19/82  this version modified to work only with even n
      dimension       wsave(*)
c
      if (n .le. 2) return
      ns2 = n/2
c      modn = n-ns2-ns2
c      if (modn .ne. 0) go to 101
      iw1 = n+1
      call rffti (n,wsave(iw1))
      return
c  101 iw1 = n+n+1
c      call cffti (n,wsave(iw1))
c      return
	end
c----------------------------------------------------------------------
      subroutine cffti (n,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine cffti initializes the array wsave which is used in
c     both cfftf and cfftb. the prime factorization of n together with
c     a tabulation of the trigonometric functions are computed and
c     stored in wsave.
c
c     input parameter
c
c     n       the length of the sequence to be transformed
c
c     output parameter
c
c     wsave   a work array which must be dimensioned at least 4*n+15
c             the same work array can be used for both cfftf and cfftb
c             as long as n remains unchanged. different wsave arrays
c             are required for different values of n. the contents of
c             wsave must not be changed between calls of cfftf or cfftb.
c
      dimension       wsave(1)
c
      if (n .eq. 1) return
      iw1 = n+n+1
      iw2 = iw1+n+n
      call cffti1 (n,wsave(iw1),wsave(iw2))
      return
      end
c----------------------------------------------------------------------
      subroutine cffti1 (n,wa,ifac)
        implicit real*4 (a-h,o-z)
      dimension       wa(1)      ,ifac(1)    ,ntryh(4)
      data ntryh(1),ntryh(2),ntryh(3),ntryh(4)/3,4,2,5/
      nl = n
      nf = 0
      j = 0
  101 j = j+1
      if (j-4) 102,102,103
  102 ntry = ntryh(j)
      go to 104
  103 ntry = ntry+2
  104 nq = nl/ntry
      nr = nl-ntry*nq
      if (nr) 101,105,101
  105 nf = nf+1
      ifac(nf+2) = ntry
      nl = nq
      if (ntry .ne. 2) go to 107
      if (nf .eq. 1) go to 107
      do 106 i=2,nf
         ib = nf-i+2
         ifac(ib+2) = ifac(ib+1)
  106 continue
      ifac(3) = 2
  107 if (nl .ne. 1) go to 104
      ifac(1) = n
      ifac(2) = nf
      tpi = 8.*atan(1.)
      arg1 = tpi/float(n)
      dc = cos(arg1)
      ds = sin(arg1)
      wa(1) = dc
      wa(2) = ds
      nt = n+n
      do 108 i=4,nt,2
         wa(i-1) = dc*wa(i-3)-ds*wa(i-2)
         wa(i) = ds*wa(i-3)+dc*wa(i-2)
  108 continue
      return
      end
c----------------------------------------------------------------------
      subroutine cfftb (n,c,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine cfftb computes the backward complex discrete fourier
c     transform (the fourier synthesis). equivalently , cfftb computes
c     a complex periodic sequence from its fourier coefficients.
c     the transform is defined below at output parameter c.
c
c     a call of cfftf followed by a call of cfftb will multiply the
c     sequence by n.
c
c     the array wsave which is used by subroutine cfftb must be
c     initialized by calling subroutine cffti(n,wsave).
c
c     input parameters
c
c
c     n      the length of the complex sequence c. the method is
c            more efficient when n is the product of small primes.
c
c     c      a complex array of length n which contains the sequence
c
c     wsave   a real work array which must be dimensioned at least 4n+15
c             in the program that calls cfftb. the wsave array must be
c             initialized by calling subroutine cffti(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c             the same wsave array can be used by cfftf and cfftb.
c
c     output parameters
c
c     c      for j=1,...,n
c
c                c(j)=the sum from k=1,...,n of
c
c                      c(k)*exp(i*j*k*2*pi/n)
c
c                            where i=sqrt(-1)
c
c     wsave   contains initialization calculations which must not be
c             destroyed between calls of subroutine cfftf or cfftb
      dimension       c(1)       ,wsave(1)
c
      if (n .eq. 1) return
      iw1 = n+n+1
      iw2 = iw1+n+n
      iw3 = iw2+1
      call cfftb1 (n,c,wsave,wsave(iw1),wsave(iw2))
      return
      end
c----------------------------------------------------------------------
      subroutine cfftb1 (n,c,ch,wa,ifac)
        implicit real*4 (a-h,o-z)
      dimension       ch(1)      ,c(1)       ,wa(1)      ,ifac(1)
      nf = ifac(2)
      l1 = 1
      do 105 k1=1,nf
         ip = ifac(k1+2)
         l2 = ip*l1
         ido = n/l2
         idot = ido+ido
         idl1 = idot*l1
         if (ip .ne. 4) go to 101
         ix2 = l1+l1
         ix3 = ix2+l1
         call passb4 (idot,l1,idl1,ix2,ix3,c,c,c,ch,ch,wa,wa,wa)
         go to 104
  101    if (ip .ne. 2) go to 102
         call passb2 (idot,l1,idl1,c,c,c,ch,ch,wa)
         go to 104
  102    if (ip .ne. 3) go to 103
         ix2 = l1+l1
         call passb3 (idot,l1,idl1,ix2,c,c,c,ch,ch,wa,wa)
         go to 104
  103    call passb (idot,ip,l1,idl1,c,c,c,ch,ch,wa)
  104    l1 = l2
  105 continue
      return
      end
c----------------------------------------------------------------------
      subroutine passb2 (ido,l1,idl1,cc,c1,c2,ch,ch2,wa1)
        implicit real*4 (a-h,o-z)
      dimension       cc(ido,2,l1)           ,c1(ido,l1,2)           ,
     1                c2(idl1,2) ,ch(ido,l1,2)           ,ch2(idl1,2),
     2                wa1(l1,1)
      idot = ido/2
      if (ido .lt. l1) go to 103
      do 102 k=1,l1
         do 101 i=1,ido
            ch(i,k,1) = cc(i,1,k)+cc(i,2,k)
            ch(i,k,2) = cc(i,1,k)-cc(i,2,k)
  101    continue
  102 continue
      go to 106
  103 do 105 i=1,ido
         do 104 k=1,l1
            ch(i,k,1) = cc(i,1,k)+cc(i,2,k)
            ch(i,k,2) = cc(i,1,k)-cc(i,2,k)
  104    continue
  105 continue
  106 do 107 ik=1,idl1
         c2(ik,1) = ch2(ik,1)
  107 continue
      do 108 k=1,l1
         c1(1,k,2) = ch(1,k,2)
         c1(2,k,2) = ch(2,k,2)
  108 continue
      if (ido .eq. 2) return
      if (idot .lt. l1) go to 111
      do 110 k=1,l1
         do 109 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
  109    continue
  110 continue
      return
  111 do 113 i=4,ido,2
         do 112 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
  112    continue
  113 continue
      return
      end
c----------------------------------------------------------------------
      subroutine passb3 (ido,l1,idl1,ix2,cc,c1,c2,ch,ch2,wa1,wa2)
        implicit real*4 (a-h,o-z)
      dimension       cc(ido,3,l1)           ,c1(ido,l1,3)           ,
     1                c2(idl1,3) ,ch(ido,l1,3)           ,ch2(idl1,3),
     2                wa1(l1,1)  ,wa2(ix2,1)
      data taur,taui /-.5,.866025403784439/
      idot = ido/2
      if (ido .lt. l1) go to 103
      do 102 k=1,l1
         do 101 i=1,ido
            ch(i,k,1) = cc(i,1,k)
            ch(i,k,2) = cc(i,2,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)-cc(i,3,k)
  101    continue
  102 continue
      go to 106
  103 do 105 i=1,ido
         do 104 k=1,l1
            ch(i,k,1) = cc(i,1,k)
            ch(i,k,2) = cc(i,2,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)-cc(i,3,k)
  104    continue
  105 continue
c
  106 do 107 ik=1,idl1
         c2(ik,1) = ch2(ik,1)+ch2(ik,2)
         c2(ik,2) = ch2(ik,1)+taur*ch2(ik,2)
         c2(ik,3) = taui*ch2(ik,3)
  107 continue
      do 108 ik=2,idl1,2
         ch2(ik-1,2) = c2(ik-1,2)-c2(ik,3)
         ch2(ik-1,3) = c2(ik-1,2)+c2(ik,3)
  108 continue
      do 109 ik=2,idl1,2
         ch2(ik,2) = c2(ik,2)+c2(ik-1,3)
         ch2(ik,3) = c2(ik,2)-c2(ik-1,3)
  109 continue
      do 111 j=2,3
         do 110 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  110    continue
  111 continue
      if (ido .eq. 2) return
      if (idot-1 .lt. l1) go to 114
      do 113 k=1,l1
         do 112 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)+
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)-
     1                    wa2(ix2,i-2)*ch(i,k,3)
  112    continue
  113 continue
      return
  114 do 116 i=4,ido,2
         do 115 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)+
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)-
     1                    wa2(ix2,i-2)*ch(i,k,3)
  115    continue
  116 continue
      return
      end
c----------------------------------------------------------------------
      subroutine passb4 (ido,l1,idl1,ix2,ix3,cc,c1,c2,ch,ch2,wa1,wa2,
     1                   wa3)
        implicit real*4 (a-h,o-z)
      dimension       cc(ido,4,l1)           ,c1(ido,l1,4)           ,
     1                c2(idl1,4) ,ch(ido,l1,4)           ,ch2(idl1,4),
     2                wa1(l1,1)  ,wa2(ix2,1) ,wa3(ix3,1)
      idot = ido/2
c
      if (ido .lt. l1) go to 106
      do 103 k=1,l1
         do 101 i=2,ido,2
            ch(i-1,k,4) = cc(i,4,k)-cc(i,2,k)
  101    continue
         do 102 i=2,ido,2
            ch(i,k,4) = cc(i-1,2,k)-cc(i-1,4,k)
  102    continue
  103 continue
      do 105 k=1,l1
         do 104 i=1,ido
            ch(i,k,2) = cc(i,1,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)+cc(i,4,k)
            ch(i,k,1) = cc(i,1,k)-cc(i,3,k)
  104    continue
  105 continue
      go to 111
  106 do 108 i=2,ido,2
         do 107 k=1,l1
            ch(i-1,k,4) = cc(i,4,k)-cc(i,2,k)
            ch(i,k,4) = cc(i-1,2,k)-cc(i-1,4,k)
  107    continue
  108 continue
      do 110 i=1,ido
         do 109 k=1,l1
            ch(i,k,2) = cc(i,1,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)+cc(i,4,k)
            ch(i,k,1) = cc(i,1,k)-cc(i,3,k)
  109    continue
  110 continue
  111 do 112 ik=1,idl1
         c2(ik,1) = ch2(ik,2)+ch2(ik,3)
  112 continue
      do 113 ik=1,idl1
         ch2(ik,3) = ch2(ik,2)-ch2(ik,3)
  113 continue
      do 114 ik=1,idl1
         ch2(ik,2) = ch2(ik,1)+ch2(ik,4)
  114 continue
      do 115 ik=1,idl1
         ch2(ik,4) = ch2(ik,1)-ch2(ik,4)
  115 continue
      do 117 j=2,4
         do 116 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  116    continue
  117 continue
      if (ido .eq. 2) return
      if (idot .lt. l1) go to 120
      do 119 k=1,l1
         do 118 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)+
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)-
     1                    wa2(ix2,i-2)*ch(i,k,3)
            c1(i,k,4) = wa3(ix3-1,i-2)*ch(i,k,4)+
     1                  wa3(ix3,i-2)*ch(i-1,k,4)
            c1(i-1,k,4) = wa3(ix3-1,i-2)*ch(i-1,k,4)-
     1                    wa3(ix3,i-2)*ch(i,k,4)
  118    continue
  119 continue
      return
  120 do 122 i=4,ido,2
         do 121 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)+
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)-
     1                    wa2(ix2,i-2)*ch(i,k,3)
            c1(i,k,4) = wa3(ix3-1,i-2)*ch(i,k,4)+
     1                  wa3(ix3,i-2)*ch(i-1,k,4)
            c1(i-1,k,4) = wa3(ix3-1,i-2)*ch(i-1,k,4)-
     1                    wa3(ix3,i-2)*ch(i,k,4)
  121    continue
  122 continue
      return
      end
c----------------------------------------------------------------------
      subroutine passb (ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa)
        implicit real*4 (a-h,o-z)
      dimension       ch(ido,l1,ip)          ,cc(ido,ip,l1)          ,
     1                c1(ido,l1,ip)          ,wa(1)      ,c2(idl1,ip),
     2                ch2(idl1,ip)
      idot = ido/2
      nt = ip*idl1
      ipp2 = ip+2
      ipph = (ip+1)/2
      l1t = l1+l1
c
      if (ido .lt. l1) go to 106
      do 103 j=2,ipph
         jc = ipp2-j
         do 102 k=1,l1
            do 101 i=1,ido
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
  101       continue
  102    continue
  103 continue
      do 105 k=1,l1
         do 104 i=1,ido
            ch(i,k,1) = cc(i,1,k)
  104    continue
  105 continue
      go to 112
  106 do 109 j=2,ipph
         jc = ipp2-j
         do 108 i=1,ido
            do 107 k=1,l1
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
  107       continue
  108    continue
  109 continue
      do 111 i=1,ido
         do 110 k=1,l1
            ch(i,k,1) = cc(i,1,k)
  110    continue
  111 continue
  112 do 113 ik=1,idl1
         c2(ik,1) = ch2(ik,1)
  113 continue
      idj = 0
      do 115 j=2,ipph
         jc = ipp2-j
         idj = idj+idl1
         do 114 ik=1,idl1
            c2(ik,j) = ch2(ik,1)+wa(idj-1)*ch2(ik,2)
            c2(ik,jc) = wa(idj)*ch2(ik,ip)
  114    continue
  115 continue
      do 117 j=2,ipph
         do 116 ik=1,idl1
            c2(ik,1) = c2(ik,1)+ch2(ik,j)
  116    continue
  117 continue
c
      idl = 0
      do 120 l=2,ipph
         lc = ipp2-l
         idl = idl+idl1
         idlj = idl
         do 119 j=3,ipph
            jc = ipp2-j
            idlj = idlj+idl
            if (idlj .gt. nt) idlj = idlj-nt
            war = wa(idlj-1)
            wai = wa(idlj)
            do 118 ik=1,idl1
               c2(ik,l) = c2(ik,l)+war*ch2(ik,j)
               c2(ik,lc) = c2(ik,lc)+wai*ch2(ik,jc)
  118       continue
  119    continue
  120 continue
c
      do 122 j=2,ipph
         jc = ipp2-j
         do 121 ik=2,idl1,2
            ch2(ik-1,j) = c2(ik-1,j)-c2(ik,jc)
            ch2(ik-1,jc) = c2(ik-1,j)+c2(ik,jc)
  121    continue
  122 continue
      do 124 j=2,ipph
         jc = ipp2-j
         do 123 ik=2,idl1,2
            ch2(ik,j) = c2(ik,j)+c2(ik-1,jc)
            ch2(ik,jc) = c2(ik,j)-c2(ik-1,jc)
  123    continue
  124 continue
c
      do 126 j=2,ip
         do 125 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  125    continue
  126 continue
      if (ido .eq. 2) return
      idj = 0
      if (idot .gt. l1) go to 130
      do 129 j=2,ip
         idj = idj+l1t
         idij = 0
         do 128 i=4,ido,2
            idij = idij+idj
            do 127 k=1,l1
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
  127       continue
  128    continue
  129 continue
      return
  130 do 134 j=2,ip
         idj = idj+l1t
         do 133 k=1,l1
            idij = 0
            do 131 i=4,ido,2
               idij = idij+idj
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
  131       continue
            idij = 0
            do 132 i=4,ido,2
               idij = idij+idj
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
  132       continue
  133    continue
  134 continue
      return
      end
c----------------------------------------------------------------------
      subroutine cfftf (n,c,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine cfftf computes the forward complex discrete fourier
c     transform (the fourier analysis). equivalently , cfftf computes
c     the fourier coefficients of a complex periodic sequence.
c     the transform is defined below at output parameter c.
c
c     the transform is not normalized. to obtain a normalized transform
c     the output must be divided by n. otherwise a call of cfftf
c     followed by a call of cfftb will multiply the sequence by n.
c
c     the array wsave which is used by subroutine cfftf must be
c     initialized by calling subroutine cffti(n,wsave).
c
c     input parameters
c
c
c     n      the length of the complex sequence c. the method is
c            more efficient when n is the product of small primes. n
c
c     c      a complex array of length n which contains the sequence
c
c     wsave   a real work array which must be dimensioned at least 4n+15
c             in the program that calls cfftf. the wsave array must be
c             initialized by calling subroutine cffti(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c             the same wsave array can be used by cfftf and cfftb.
c
c     output parameters
c
c     c      for j=1,...,n
c
c                c(j)=the sum from k=1,...,n of
c
c                      c(k)*exp(-i*j*k*2*pi/n)
c
c                            where i=sqrt(-1)
c
c     wsave   contains initialization calculations which must not be
c             destroyed between calls of subroutine cfftf or cfftb
c
      dimension       c(1)       ,wsave(1)
c
      if (n .eq. 1) return
      iw1 = n+n+1
      iw2 = iw1+n+n
      call cfftf1 (n,c,wsave,wsave(iw1),wsave(iw2))
      return
      end
c----------------------------------------------------------------------
      subroutine cfftf1 (n,c,ch,wa,ifac)
        implicit real*4 (a-h,o-z)
      dimension       ch(1)      ,c(1)       ,wa(1)      ,ifac(1)
      nf = ifac(2)
      l1 = 1
      do 105 k1=1,nf
         ip = ifac(k1+2)
         l2 = ip*l1
         ido = n/l2
         idot = ido+ido
         idl1 = idot*l1
         if (ip .ne. 4) go to 101
         ix2 = l1+l1
         ix3 = ix2+l1
         call passf4 (idot,l1,idl1,ix2,ix3,c,c,c,ch,ch,wa,wa,wa)
         go to 104
  101    if (ip .ne. 2) go to 102
         call passf2 (idot,l1,idl1,c,c,c,ch,ch,wa)
         go to 104
  102    if (ip .ne. 3) go to 103
         ix2 = l1+l1
         call passf3 (idot,l1,idl1,ix2,c,c,c,ch,ch,wa,wa)
         go to 104
  103    call passf (idot,ip,l1,idl1,c,c,c,ch,ch,wa)
  104    l1 = l2
  105 continue
      return
      end
c----------------------------------------------------------------------
      subroutine passf2 (ido,l1,idl1,cc,c1,c2,ch,ch2,wa1)
        implicit real*4 (a-h,o-z)
      dimension       cc(ido,2,l1)           ,c1(ido,l1,2)           ,
     1                c2(idl1,2) ,ch(ido,l1,2)           ,ch2(idl1,2),
     2                wa1(l1,1)
      idot = ido/2
      if (ido .lt. l1) go to 103
      do 102 k=1,l1
         do 101 i=1,ido
            ch(i,k,1) = cc(i,1,k)+cc(i,2,k)
            ch(i,k,2) = cc(i,1,k)-cc(i,2,k)
  101    continue
  102 continue
      go to 106
  103 do 105 i=1,ido
         do 104 k=1,l1
            ch(i,k,1) = cc(i,1,k)+cc(i,2,k)
            ch(i,k,2) = cc(i,1,k)-cc(i,2,k)
  104    continue
  105 continue
  106 do 107 ik=1,idl1
         c2(ik,1) = ch2(ik,1)
  107 continue
      do 108 k=1,l1
         c1(1,k,2) = ch(1,k,2)
         c1(2,k,2) = ch(2,k,2)
  108 continue
      if (ido .eq. 2) return
      if (idot .lt. l1) go to 111
      do 110 k=1,l1
         do 109 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
  109    continue
  110 continue
      return
  111 do 113 i=4,ido,2
         do 112 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
  112    continue
  113 continue
      return
      end
c----------------------------------------------------------------------
      subroutine passf3 (ido,l1,idl1,ix2,cc,c1,c2,ch,ch2,wa1,wa2)
        implicit real*4 (a-h,o-z)
      dimension       cc(ido,3,l1)           ,c1(ido,l1,3)           ,
     1                c2(idl1,3) ,ch(ido,l1,3)           ,ch2(idl1,3),
     2                wa1(l1,1)  ,wa2(ix2,1)
      data taur,taui /-.5,-.866025403784439/
      idot = ido/2
      if (ido .lt. l1) go to 103
      do 102 k=1,l1
         do 101 i=1,ido
            ch(i,k,1) = cc(i,1,k)
            ch(i,k,2) = cc(i,2,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)-cc(i,3,k)
  101    continue
  102 continue
      go to 106
  103 do 105 i=1,ido
         do 104 k=1,l1
            ch(i,k,1) = cc(i,1,k)
            ch(i,k,2) = cc(i,2,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)-cc(i,3,k)
  104    continue
  105 continue
c
  106 do 107 ik=1,idl1
         c2(ik,1) = ch2(ik,1)+ch2(ik,2)
         c2(ik,2) = ch2(ik,1)+taur*ch2(ik,2)
         c2(ik,3) = taui*ch2(ik,3)
  107 continue
      do 108 ik=2,idl1,2
         ch2(ik-1,2) = c2(ik-1,2)-c2(ik,3)
         ch2(ik-1,3) = c2(ik-1,2)+c2(ik,3)
  108 continue
      do 109 ik=2,idl1,2
         ch2(ik,2) = c2(ik,2)+c2(ik-1,3)
         ch2(ik,3) = c2(ik,2)-c2(ik-1,3)
  109 continue
      do 111 j=2,3
         do 110 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  110    continue
  111 continue
      if (ido .eq. 2) return
      if (idot-1 .lt. l1) go to 114
      do 113 k=1,l1
         do 112 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)-
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)+
     1                    wa2(ix2,i-2)*ch(i,k,3)
  112    continue
  113 continue
      return
  114 do 116 i=4,ido,2
         do 115 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)-
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)+
     1                    wa2(ix2,i-2)*ch(i,k,3)
  115    continue
  116 continue
      return
      end
c----------------------------------------------------------------------
      subroutine passf4 (ido,l1,idl1,ix2,ix3,cc,c1,c2,ch,ch2,wa1,wa2,
     1                   wa3)
        implicit real*4 (a-h,o-z)
      dimension       cc(ido,4,l1)           ,c1(ido,l1,4)           ,
     1                c2(idl1,4) ,ch(ido,l1,4)           ,ch2(idl1,4),
     2                wa1(l1,1)  ,wa2(ix2,1) ,wa3(ix3,1)
      idot = ido/2
c
      if (ido .lt. l1) go to 106
      do 103 k=1,l1
         do 101 i=2,ido,2
            ch(i-1,k,4) = cc(i,2,k)-cc(i,4,k)
  101    continue
         do 102 i=2,ido,2
            ch(i,k,4) = cc(i-1,4,k)-cc(i-1,2,k)
  102    continue
  103 continue
      do 105 k=1,l1
         do 104 i=1,ido
            ch(i,k,2) = cc(i,1,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)+cc(i,4,k)
            ch(i,k,1) = cc(i,1,k)-cc(i,3,k)
  104    continue
  105 continue
      go to 111
  106 do 108 i=2,ido,2
         do 107 k=1,l1
            ch(i,k,4) = cc(i-1,4,k)-cc(i-1,2,k)
            ch(i-1,k,4) = cc(i,2,k)-cc(i,4,k)
  107    continue
  108 continue
      do 110 i=1,ido
         do 109 k=1,l1
            ch(i,k,2) = cc(i,1,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)+cc(i,4,k)
            ch(i,k,1) = cc(i,1,k)-cc(i,3,k)
  109    continue
  110 continue
  111 do 112 ik=1,idl1
         c2(ik,1) = ch2(ik,2)+ch2(ik,3)
  112 continue
      do 113 ik=1,idl1
         ch2(ik,3) = ch2(ik,2)-ch2(ik,3)
  113 continue
      do 114 ik=1,idl1
         ch2(ik,2) = ch2(ik,1)+ch2(ik,4)
  114 continue
      do 115 ik=1,idl1
         ch2(ik,4) = ch2(ik,1)-ch2(ik,4)
  115 continue
      do 117 j=2,4
         do 116 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  116    continue
  117 continue
      if (ido .eq. 2) return
      if (idot .lt. l1) go to 120
      do 119 k=1,l1
         do 118 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)-
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)+
     1                    wa2(ix2,i-2)*ch(i,k,3)
            c1(i,k,4) = wa3(ix3-1,i-2)*ch(i,k,4)-
     1                  wa3(ix3,i-2)*ch(i-1,k,4)
            c1(i-1,k,4) = wa3(ix3-1,i-2)*ch(i-1,k,4)+
     1                    wa3(ix3,i-2)*ch(i,k,4)
  118    continue
  119 continue
      return
  120 do 122 i=4,ido,2
         do 121 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)-
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)+
     1                    wa2(ix2,i-2)*ch(i,k,3)
            c1(i,k,4) = wa3(ix3-1,i-2)*ch(i,k,4)-
     1                  wa3(ix3,i-2)*ch(i-1,k,4)
            c1(i-1,k,4) = wa3(ix3-1,i-2)*ch(i-1,k,4)+
     1                    wa3(ix3,i-2)*ch(i,k,4)
  121    continue
  122 continue
      return
      end
c----------------------------------------------------------------------
      subroutine passf (ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa)
        implicit real*4 (a-h,o-z)
      dimension       ch(ido,l1,ip)          ,cc(ido,ip,l1)          ,
     1                c1(ido,l1,ip)          ,wa(1)      ,c2(idl1,ip),
     2                ch2(idl1,ip)
      idot = ido/2
      nt = ip*idl1
      ipp2 = ip+2
      ipph = (ip+1)/2
      l1t = l1+l1
c
      if (ido .lt. l1) go to 106
      do 103 j=2,ipph
         jc = ipp2-j
         do 102 k=1,l1
            do 101 i=1,ido
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
  101       continue
  102    continue
  103 continue
      do 105 k=1,l1
         do 104 i=1,ido
            ch(i,k,1) = cc(i,1,k)
  104    continue
  105 continue
      go to 112
  106 do 109 j=2,ipph
         jc = ipp2-j
         do 108 i=1,ido
            do 107 k=1,l1
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
  107       continue
  108    continue
  109 continue
      do 111 i=1,ido
         do 110 k=1,l1
            ch(i,k,1) = cc(i,1,k)
  110    continue
  111 continue
  112 do 113 ik=1,idl1
         c2(ik,1) = ch2(ik,1)
  113 continue
      idj = 0
      do 115 j=2,ipph
         jc = ipp2-j
         idj = idj+idl1
         do 114 ik=1,idl1
            c2(ik,j) = ch2(ik,1)+wa(idj-1)*ch2(ik,2)
            c2(ik,jc) = -wa(idj)*ch2(ik,ip)
  114    continue
  115 continue
      do 117 j=2,ipph
         do 116 ik=1,idl1
            c2(ik,1) = c2(ik,1)+ch2(ik,j)
  116    continue
  117 continue
c
      idl = 0
      do 120 l=2,ipph
         lc = ipp2-l
         idl = idl+idl1
         idlj = idl
         do 119 j=3,ipph
            jc = ipp2-j
            idlj = idlj+idl
            if (idlj .gt. nt) idlj = idlj-nt
            war = wa(idlj-1)
            wai = wa(idlj)
            do 118 ik=1,idl1
               c2(ik,l) = c2(ik,l)+war*ch2(ik,j)
               c2(ik,lc) = c2(ik,lc)-wai*ch2(ik,jc)
  118       continue
  119    continue
  120 continue
c
      do 122 j=2,ipph
         jc = ipp2-j
         do 121 ik=2,idl1,2
            ch2(ik-1,j) = c2(ik-1,j)-c2(ik,jc)
            ch2(ik-1,jc) = c2(ik-1,j)+c2(ik,jc)
  121    continue
  122 continue
      do 124 j=2,ipph
         jc = ipp2-j
         do 123 ik=2,idl1,2
            ch2(ik,j) = c2(ik,j)+c2(ik-1,jc)
            ch2(ik,jc) = c2(ik,j)-c2(ik-1,jc)
  123    continue
  124 continue
c
      do 126 j=2,ip
         do 125 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  125    continue
  126 continue
      if (ido .eq. 2) return
      idj = 0
      if (idot .gt. l1) go to 130
      do 129 j=2,ip
         idj = idj+l1t
         idij = 0
         do 128 i=4,ido,2
            idij = idij+idj
            do 127 k=1,l1
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)+wa(idij)*ch(i,k,j)
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)-wa(idij)*ch(i-1,k,j)
  127       continue
  128    continue
  129 continue
      return
  130 do 134 j=2,ip
         idj = idj+l1t
         do 133 k=1,l1
            idij = 0
            do 131 i=4,ido,2
               idij = idij+idj
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)+wa(idij)*ch(i,k,j)
  131       continue
            idij = 0
            do 132 i=4,ido,2
               idij = idij+idj
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)-wa(idij)*ch(i-1,k,j)
  132       continue
  133    continue
  134 continue
      return
      end
c----------------------------------------------------------------------
      subroutine cosqi (n,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine cosqi initializes the array wsave which is used in
c     both cosqf and cosqb. the prime factorization of n together with
c     a tabulation of the trigonometric functions are computed and
c     stored in wsave.
c
c     input parameter
c
c     n       the length of the sequence to be transformed. n must be
c             even
c
c     output parameter
c
c     wsave   a work array which must be dimensioned at least 3.5*n+15.
c             the same work array can be used for both cosqf and cosqb
c             as long as n remains unchanged. different wsave arrays
c             are required for different values of n. the contents of
c             wsave must not be changed between calls of cosqf or cosqb.
c
      dimension       wsave(1)
c
      iw1 = n+1
      pih = 2.*atan(1.)
      dt = pih/float(n)
      dc = cos(dt)
      ds = sin(dt)
      wsave(1) = dc
      wsk = ds
      do 101 k=2,n
         wsave(k) = dc*wsave(k-1)-ds*wsk
         wsk = ds*wsave(k-1)+dc*wsk
  101 continue
      call rffti (n,wsave(iw1))
      return
      end
c------------------------------------------------------------------------
      subroutine cosqf (n,x,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine cosqf computes the fast fourier transform of quarter
c     wave data. that is , cosqf computes the coefficients in a cosine
c     series representation with only odd wave numbers. the transform
c     is defined below at output parameter x
c
c     cosqf is the unnormalized inverse of cosqb since a call of cosqf
c     followed by a call of cosqb will multiply the input sequence x
c     by 8*n.
c
c     the array wsave which is used by subroutine cosqf must be
c     initialized by calling subroutine cosqi(n,wsave).
c
c
c     input parameters
c
c     n       the length of the array x. n must be even and the method
c             is most efficient when n is a product of small primes.
c
c     x       an array which contains the sequence to be transformed
c
c     wsave   a work array which must be dimensioned at least 3.5n+15
c             in the program that calls cosqf. the wsave array must be
c             initialized by calling subroutine cosqi(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c
c     output parameters
c
c     x       for i=1,...,n
c
c                  x(i)=2*x(1) plus the sum from k=2 to k=n of
c
c                     4.*x(k)*cos((2*i-1)*(k-1)*pi/(2*n))
c
c                  a call of cosqf followed by a call of
c                  cosqb will multiply the sequence x by 8*n
c                  therefore cosqb is the unnormalized inverse
c                  of cosqf.
c
c     wsave   contains initialization calculations which must not
c             be destroyed between calls of cosqf or cosqb.
c
      dimension       x(1)       ,wsave(1)
      data tsq2/2.82842712474619/
c
C        type 2,n, x(1)
2       format(' cosqf, n =',i5, ' x(1) = ', g13.4)
      if (n .gt. 2) go to 101
      tx = x(1)+x(1)
      tsqx = tsq2*x(2)
      x(1) = tx+tsqx
      x(2) = tx-tsqx
      return
  101 iw1 = n+1
      iw2 = n/2+iw1
      call cosqf1 (n,x,wsave,wsave(iw1),wsave(iw2))
      return
      end
c----------------------------------------------------------------------
      subroutine cosqf1 (n,x,w,w1,xh)
        implicit real*4 (a-h,o-z)
      dimension       x(1)       ,xh(1)      ,w1(1)      ,w(1)
      ns2 = n/2
      nm1 = n-1
      np2 = n+2
      do 101 k=2,ns2
         kc = np2-k
         xh(k) = x(k)+x(kc)
         xh(kc) = x(k)-x(kc)
  101 continue
      xh(ns2+1) = x(ns2+1)+x(ns2+1)
      do 102 k=2,ns2
         kc = np2-k
         x(k) = w(k-1)*xh(kc)+w(kc-1)*xh(k)
         x(kc) = w(k-1)*xh(k)-w(kc-1)*xh(kc)
  102 continue
      x(ns2+1) = w(ns2)*xh(ns2+1)
      call rfftf (n,x,w1)
      xn = x(2)
      do 103 i=3,nm1,2
         x(i-1) = x(i)+x(i+1)
         x(i) = x(i)-x(i+1)
  103 continue
      x(n) = xn
      return
      end
c----------------------------------------------------------------------
      subroutine cosqb (n,x,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine cosqb computes the fast fourier transform of quarter
c     wave data. that is , cosqb computes a sequence from its
c     representation in terms of a cosine series with odd wave numbers.
c     the transform is defined below at output parameter x.
c
c     cosqb is the unnormalized inverse of cosqf since a call of cosqb
c     followed by a call of cosqf will multiply the input sequence x
c     by 8*n.
c
c     the array wsave which is used by subroutine cosqb must be
c     initialized by calling subroutine cosqi(n,wsave).
c
c
c     input parameters
c
c     n       the length of the array x. n must be even and the method
c             is most efficient when n is a product of small primes.
c
c     x       an array which contains the sequence to be transformed
c
c     wsave   a work array which must be dimensioned at least 3.5n+15
c             in the program that calls cosqb. the wsave array must be
c             initialized by calling subroutine cosqi(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c
c     output parameters
c
c     x       for i=1,...,n
c
c                  x(i)= the sum from k=1 to k=n of
c
c                    4*x(k)*cos((2*k-1)*(i-1)*pi/(2*n))
c
c                  a call of cosqb followed by a call of
c                  cosqf will multiply the sequence x by 8*n
c                  therefore cosqf is the unnormalized inverse
c                  of cosqb.
c
c     wsave   contains initialization calculations which must not
c             be destroyed between calls of cosqb or cosqf.
c
      dimension       x(1)       ,wsave(1)
      data tsq2/2.82842712474619/
c
      if (n .gt. 2) go to 101
      x1 = 4.*(x(1)+x(2))
      x(2) = tsq2*(x(1)-x(2))
      x(1) = x1
      return
  101 iw1 = n+1
      iw2 = n/2+iw1
      call cosqb1 (n,x,wsave,wsave(iw1),wsave(iw2))
      return
      end
c----------------------------------------------------------------------
      subroutine cosqb1 (n,x,w,w1,xh)
        implicit real*4 (a-h,o-z)
      dimension       x(1)       ,w(1)       ,w1(1)      ,xh(1)
      ns2 = n/2
      nm1 = n-1
      np2 = n+2
      xn = x(n)
      do 101 ii=3,nm1,2
         i = np2-ii
         x(i+1) = x(i-1)-x(i)
         x(i) = x(i)+x(i-1)
  101 continue
      x(1) = x(1)+x(1)
      x(2) = xn+xn
      call rfftb (n,x,w1)
      do 102 k=2,ns2
         kc = np2-k
         xh(k) = w(k-1)*x(kc)+w(kc-1)*x(k)
         xh(kc) = w(k-1)*x(k)-w(kc-1)*x(kc)
  102 continue
      x(ns2+1) = w(ns2)*(x(ns2+1)+x(ns2+1))
      do 103 k=2,ns2
         kc = np2-k
         x(k) = xh(k)+xh(kc)
         x(kc) = xh(k)-xh(kc)
  103 continue
      x(1) = x(1)+x(1)
      return
      end
c----------------------------------------------------------------------
      subroutine costi (n,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine costi initializes the array wsave which is used in
c     subroutine cost. the prime factorization of n together with
c     a tabulation of the trigonometric functions are computed and
c     stored in wsave.
c
c     input parameter
c
c     n       the length of the sequence to be transformed. n must be
c             odd and greater than 1.
c
c     output parameter
c
c     wsave   a work array which must be dimensioned at least 3*n+15.
c             different wsave arrays are required for different values
c             of n. the contents of wsave must not be changed between
c             calls of cost.
c
      dimension       wsave(1)
c
      if (n .le. 3) return
      nm1 = n-1
      ns2 = nm1/2
      ns2m = ns2-1
      iw1 = ns2+1
      pi = 4.*atan(1.)
      dt = pi/float(nm1)
      dcs = cos(dt)
      dss = sin(dt)
      wsave(1) = dss+dss
      wck = dcs+dcs
      if (ns2m .lt. 2) go to 102
      do 101 k=2,ns2m
         wsave(k) = dss*wck+dcs*wsave(k-1)
         wck = dcs*wck-dss*wsave(k-1)
  101 continue
  102 call rffti (nm1,wsave(iw1))
      return
      end
c----------------------------------------------------------------------
      subroutine cost (n,x,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine cost computes the discrete fourier cosine transform
c     of an even sequence x(i). the transform is defined below at output
c     parameter x.
c
c     cost is the unnormalized inverse of itself since a call of cost
c     followed by another call of cost will multiply the input sequence
c     x by 8*(n+1). the transform is defined below at output parameter x
c
c     the array wsave which is used by subroutine cost must be
c     initialized by calling subroutine costi(n,wsave).
c
c     input parameters
c
c     n       the length of the sequence x. n must be odd and greater
c             than 1. the method is most efficient when n-1 is a product
c             of small primes.
c
c     x       an array which contains the sequence to be transformed
c
c     wsave   a work array which must be dimensioned at least 3*n+15
c             in the program that calls cost. the wsave array must be
c             initialized by calling subroutine costi(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c
c     output parameters
c
c     x       for i=1,...,n
c
c                x(i)= 2*x(1)+2*(-1)**(i+1)*x(n)
c
c                  + the sum from k=2 to k=n-1
c
c                    4*x(k)*cos((k-1)*(i-1)*pi/(n-1))
c
c                  a call of cost followed by another call of
c                  cost will multiply the sequence x by 8(n+1)
c                  hence cost is the unnormalized inverse
c                  of itself.
c
c     wsave   contains initialization calculations which must not be
c             destroyed between calls of cost.
c
      dimension       x(1)       ,wsave(1)
c
      nm1 = n-1
      np1 = n+1
      ns2 = nm1/2
      iw1 = ns2+1
      if (n .gt. 3) go to 101
      x1px3 = x(1)+x(3)
      tx13 = x1px3+x1px3
      fx2 = 4.*x(2)
      x(2) = 2.*(x(1)-x(3))
      x(1) = tx13+fx2
      x(3) = tx13-fx2
      return
  101 c1 = x(1)-x(n)
      do 102 k=2,ns2
         kc = nm1-k
         ks = ns2-k
         c1 = c1+wsave(ks+1)*(x(k)-x(kc+2))
  102 continue
      c1 = c1+c1
      x(1) = x(1)+x(n)
      do 103 k=2,ns2
         kc = np1-k
         t1 = x(k)+x(kc)
         t2 = wsave(k-1)*(x(k)-x(kc))
         x(k) = t1-t2
         x(kc) = t1+t2
  103 continue
      x(ns2+1) = x(ns2+1)+x(ns2+1)
      call rfftf (nm1,x,wsave(iw1))
      cn = x(2)
      x(2) = c1
      do 104 i=4,nm1,2
         x(i) = x(i)+x(i-2)
  104 continue
      x(n) = cn
      return
      end
c----------------------------------------------------------------------
      subroutine rffti (n,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine rffti initializes the array wsave which is used in
c     both rfftf and rfftb. the prime factorization of n together with
c     a tabulation of the trigonometric functions are computed and
c     stored in wsave.
c
c     input parameter
c
c     n       the length of the sequence to be transformed. n must be
c             even
c
c     output parameter
c
c     wsave   a work array which must be dimensioned at least 2.5*n+15.
c             the same work array can be used for both rfftf and rfftb
c             as long as n remains unchanged. different wsave arrays
c             are required for different values of n. the contents of
c             wsave must not be changed between calls of rfftf or rfftb.
c
      dimension       wsave(1)
c
      ns2 = n/2
      nqm = (ns2-1)/2
      tpi = 8.*atan(1.)
      dt = tpi/float(n)
      dc = cos(dt)
      ds = sin(dt)
      wsave(1) = dc
      wsave(ns2-1) = ds
      if (nqm .lt. 2) go to 102
      do 101 k=2,nqm
         kc = ns2-k
         wsave(k) = dc*wsave(k-1)-ds*wsave(kc+1)
         wsave(kc) = ds*wsave(k-1)+dc*wsave(kc+1)
  101 continue
  102 iw1 = ns2+1
      call cffti (ns2,wsave(iw1))
      return
      end
c----------------------------------------------------------------------
      subroutine rfftf (n,r,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine rfftf computes the fourier coefficients of a real
c     perodic sequence (fourier analysis). the transform is defined
c     below at output parameter r.
c
c     input parameters
c
c     n       the length of the array r. n must be even and the method
c             is most efficient when n is a product of small primes.
c             n may change so long as different work arrays are provided
c
c     r       a real array of length n which contains the sequence
c             to be transformed
c
c     wsave   a work array which must be dimensioned at least 2.5*n+15
c             in the program that calls rfftf. the wsave array must be
c             initialized by calling subroutine rffti(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c             the same wsave array can be used by rfftf and rfftb.
c
c
c     output parameters
c
c     r       for k=2,...,n/2
c
c                  r(2*k-1)= the sum from i=1 to i=n of
c
c                       2.*r(i)*cos((k-1)*(i-1)*2*pi/n)
c
c                  r(2*k)= the sum from i=1 to i=n of
c
c                       2.*r(i)*sin((k-1)*(i-1)*2*pi/n)
c
c             also
c
c                  r(1)= the sum from i=1 to i=n of 2.*r(i)
c
c                  r(2)= the sum from i=1 to i=n of 2.*(-1)**(i+1)*r(i)
c
c      *****  note
c                  this transform is unnormalized since a call of rfftf
c                  followed by a call of rfftb will multiply the input
c                  sequence by 2*n.
c
c     wsave   contains results which must not be destroyed between
c             calls of rfftf or rfftb.
c
c
      dimension       r(2,1)     ,wsave(1)
c
c	type *,'n, r(1,1) =', n, r(1,1)
      if (n .gt. 2) go to 101
      r1 = 2.*(r(1,1)+r(2,1))
      r(2,1) = 2.*(r(1,1)-r(2,1))
      r(1,1) = r1
      return
  101 iw1 = n/2+1
      call rfftf1 (n,r,wsave(iw1),wsave)
      return
      end
c----------------------------------------------------------------------
      subroutine rfftf1 (n,x,xh,w)
        implicit real*4 (a-h,o-z)
      dimension       x(2,1)     ,xh(2,1)    ,w(1)
      ns2 = n/2
      ns2p2 = ns2+2
      nq = ns2/2
      ipar = ns2-nq-nq
      nqm = nq
      if (ipar .eq. 0) nqm = nqm-1
      nqp = nqm+1
      call cfftf (ns2,x,xh)
      if (nqp .lt. 2) go to 107
      do 101 k=2,nqp
         kc = ns2p2-k
         xh(1,kc) = x(2,k)+x(2,kc)
         xh(2,kc) = x(1,kc)-x(1,k)
  101 continue
      do 102 k=2,nqp
         kc = ns2p2-k
         xh(1,k) = x(1,k)+x(1,kc)
         xh(2,k) = x(2,k)-x(2,kc)
  102 continue
      do 103 k=2,nqp
         kc = ns2p2-k
         x(1,kc) = w(k-1)*xh(1,kc)+w(kc-1)*xh(2,kc)
         x(2,kc) = w(k-1)*xh(2,kc)-w(kc-1)*xh(1,kc)
  103 continue
      do 104 k=2,nqp
         kc = ns2p2-k
         xh(1,kc) = x(1,kc)
         xh(2,kc) = x(2,kc)
  104 continue
      do 105 k=2,nqp
         kc = ns2p2-k
         x(1,kc) = xh(1,k)-xh(1,kc)
         x(2,kc) = xh(2,k)-xh(2,kc)
  105 continue
      do 106 k=2,nqp
         kc = ns2p2-k
         x(1,k) = xh(1,k)+xh(1,kc)
         x(2,k) = -xh(2,k)-xh(2,kc)
  106 continue
      if (ipar .ne. 0) go to 108
  107 x(1,nqp+1) = x(1,nqp+1)+x(1,nqp+1)
      x(2,nqp+1) = x(2,nqp+1)+x(2,nqp+1)
  108 xhold1 = x(2,1)+x(1,1)
      xhold2 = x(1,1)-x(2,1)
      x(2,1) = xhold2+xhold2
      x(1,1) = xhold1+xhold1
      return
      end
c----------------------------------------------------------------------
      subroutine rfftb (n,r,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine rfftb computes the real perodic sequence from its
c     fourier coefficients (fourier synthesis). the transform is defined
c     below at output parameter r.
c
c     input parameters
c
c     n       the length of the array r. n must be even and the method
c             is most efficient when n is a product of small primes.
c             n may change so long as different work arrays are provided
c
c     r       a real array of length n which contains the sequence
c             to be transformed
c
c     wsave   a work array which must be dimensioned at least 2.5*n+15
c             in the program that calls rfftb. the wsave array must be
c             initialized by calling subroutine rffti(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c             the same wsave array can be used by rfftf and rfftb.
c
c
c     output parameters
c
c     r       for i=1,...,n
c
c                  r(i)=x(1)+(-1)**(i+1)*x(2)
c
c                       plus the sum from k=2 to k=n/2 of
c
c                         2*r(2k-1)*cos((k-1)*(i-1)*2*pi/n)
c
c                        +2*r(2k)*sin((k-1)*(i-1)*2*pi/n)
c
c      *****  note
c                  this transform is unnormalized since a call of rfftf
c                  followed by a call of rfftb will multiply the input
c                  sequence by 2*n.
c
c     wsave   contains results which must not be destroyed between
c             calls of rfftb or rfftf.
c
c
      dimension       r(2,1)     ,wsave(1)
c
      if (n .gt. 2) go to 101
      r1 = r(1,1)+r(2,1)
      r(2,1) = r(1,1)-r(2,1)
      r(1,1) = r1
      return
  101 iw1 = n/2+1
      call rfftb1 (n,r,wsave(iw1),wsave)
      return
      end
c----------------------------------------------------------------------
      subroutine rfftb1 (n,x,xh,w)
        implicit real*4 (a-h,o-z)
      dimension       x(2,1)     ,xh(2,1)    ,w(1)
      ns2 = n/2
      ns2p2 = ns2+2
      nq = ns2/2
      ipar = ns2-nq-nq
      nqm = nq
      if (ipar .eq. 0) nqm = nqm-1
      nqp = nqm+1
      xhold1 = x(1,1)-x(2,1)
      x(1,1) = x(2,1)+x(1,1)
      x(2,1) = xhold1
      if (ipar .ne. 0) go to 101
      x(1,nqp+1) = x(1,nqp+1)+x(1,nqp+1)
      x(2,nqp+1) = x(2,nqp+1)+x(2,nqp+1)
  101 if (nqp .lt. 2) go to 108
      do 102 k=2,nqp
         kc = ns2p2-k
         xh(1,k) = x(1,k)+x(1,kc)
         xh(2,k) = x(2,kc)-x(2,k)
  102 continue
      do 103 k=2,nqp
         kc = ns2p2-k
         xh(1,kc) = x(1,k)-x(1,kc)
         xh(2,kc) = -x(2,k)-x(2,kc)
  103 continue
      do 104 k=2,nqp
         kc = ns2p2-k
         x(1,kc) = xh(1,kc)
         x(2,kc) = xh(2,kc)
  104 continue
      do 105 k=2,nqp
         kc = ns2p2-k
         xh(1,kc) = w(k-1)*x(1,kc)-w(kc-1)*x(2,kc)
         xh(2,kc) = w(k-1)*x(2,kc)+w(kc-1)*x(1,kc)
  105 continue
      do 106 k=2,nqp
         kc = ns2p2-k
         x(1,k) = xh(1,k)-xh(2,kc)
         x(2,k) = xh(2,k)+xh(1,kc)
  106 continue
      do 107 k=2,nqp
         kc = ns2p2-k
         x(1,kc) = xh(1,k)+xh(2,kc)
         x(2,kc) = xh(1,kc)-xh(2,k)
  107 continue
  108 call cfftb (ns2,x,xh)
      return
      end
c----------------------------------------------------------------------
      subroutine sinqi (n,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine sinqi initializes the array wsave which is used in
c     both sinqf and sinqb. the prime factorization of n together with
c     a tabulation of the trigonometric functions are computed and
c     stored in wsave.
c
c     input parameter
c
c     n       the length of the sequence to be transformed. n must be
c             even
c
c     output parameter
c
c     wsave   a work array which must be dimensioned at least 3.5*n+15.
c             the same work array can be used for both sinqf and sinqb
c             as long as n remains unchanged. different wsave arrays
c             are required for different values of n. the contents of
c             wsave must not be changed between calls of sinqf or sinqb.
c
      dimension       wsave(1)
c
      call cosqi (n,wsave)
      return
      end
c----------------------------------------------------------------------
      subroutine sinqf (n,x,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine sinqf computes the fast fourier transform of quarter
c     wave data. that is , sinqf computes the coefficients in a sine
c     series representation with only odd wave numbers. the transform
c     is defined below at output parameter x.
c
c     sinqb is the unnormalized inverse of sinqf since a call of sinqf
c     followed by a call of sinqb will multiply the input sequence x
c     by 8*n.
c
c     the array wsave which is used by subroutine sinqf must be
c     initialized by calling subroutine sinqi(n,wsave).
c
c
c     input parameters
c
c     n       the length of the array x. n must be even and the method
c             is most efficient when n is a product of small primes.
c
c     x       an array which contains the sequence to be transformed
c
c     wsave   a work array which must be dimensioned at least 3.5n+15
c             in the program that calls sinqf. the wsave array must be
c             initialized by calling subroutine sinqi(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c
c     output parameters
c
c     x       for i=1,...,n
c
c                  x(i)=2*(-1)**(i+1)*x(n)
c
c                     + the sum from k=1 to k=n-1 of
c
c                     4*x(k)*sin((2i-1)*k*pi/(2*n))
c
c                  a call of sinqf followed by a call of
c                  sinqb will multiply the sequence x by 8*n
c                  therefore sinqb is the unnormalized inverse
c                  of sinqf.
c
c     wsave   contains initialization calculations which must not
c             be destroyed between calls of sinqf or sinqb.
c
      dimension       x(1)       ,wsave(1)
c
      ns2 = n/2
      do 101 k=1,ns2
         kc = n-k
         xhold = x(k)
         x(k) = x(kc+1)
         x(kc+1) = xhold
  101 continue
      call cosqf (n,x,wsave)
      do 102 k=2,n,2
         x(k) = -x(k)
  102 continue
      return
      end
c----------------------------------------------------------------------
      subroutine sinqb (n,x,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine sinqb computes the fast fourier transform of quarter
c     wave data. that is , sinqb computes a sequence from its
c     representation in terms of a sine series with odd wave numbers.
c     the transform is defined below at output parameter x.
c
c     sinqf is the unnormalized inverse of sinqb since a call of sinqb
c     followed by a call of sinqf will multiply the input sequence x
c     by 8*n.
c
c     the array wsave which is used by subroutine sinqb must be
c     initialized by calling subroutine sinqi(n,wsave).
c
c
c     input parameters
c
c     n       the length of the array x. n must be even and the method
c             is most efficient when n is a product of small primes.
c
c     x       an array which contains the sequence to be transformed
c
c     wsave   a work array which must be dimensioned at least 3.5n+15
c             in the program that calls sinqb. the wsave array must be
c             initialized by calling subroutine sinqi(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c
c     output parameters
c
c     x       for i=1,...,n
c
c                  x(i)= the sum from k=1 to k=n of
c
c                    4*x(k)*sin((2k-1)*i*pi/(2*n))
c
c                  a call of sinqb followed by a call of
c                  sinqf will multiply the sequence x by 8*n
c                  therefore sinqf is the unnormalized inverse
c                  of sinqb.
c
c     wsave   contains initialization calculations which must not
c             be destroyed between calls of sinqb or sinqf.
c
      dimension       x(1)       ,wsave(1)
c
      ns2 = n/2
      do 101 k=2,n,2
         x(k) = -x(k)
  101 continue
      call cosqb (n,x,wsave)
      do 102 k=1,ns2
         kc = n-k
         xhold = x(k)
         x(k) = x(kc+1)
         x(kc+1) = xhold
  102 continue
      return
      end
c----------------------------------------------------------------------
      subroutine sinti (n,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine sinti initializes the array wsave which is used in
c     subroutine sint. the prime factorization of n together with
c     a tabulation of the trigonometric functions are computed and
c     stored in wsave.
c
c     input parameter
c
c     n       the length of the sequence to be transformed. n must be
c             odd.
c
c     output parameter
c
c     wsave   a work array which must be dimensioned at least 3*n+18.
c             different wsave arrays are required for different values
c             of n. the contents of wsave must not be changed between
c             calls of sint.
c
      dimension       wsave(1)
c
      if (n .le. 1) return
      np1 = n+1
      ns2 = np1/2
      ns2m = ns2-1
      iw1 = ns2+1
      pi = 4.*atan(1.)
      dt = pi/float(np1)
      dcs = cos(dt)
      dss = sin(dt)
      wsave(1) = dss+dss
      wck = dcs+dcs
      if (ns2m .lt. 2) go to 102
      do 101 k=2,ns2m
         wsave(k) = dss*wck+dcs*wsave(k-1)
         wck = dcs*wck-dss*wsave(k-1)
  101 continue
  102 call rffti (np1,wsave(iw1))
      return
      end
c----------------------------------------------------------------------
      subroutine sint (n,x,wsave)
        implicit real*4 (a-h,o-z)
c     subroutine sint computes the discrete fourier sine transform
c     of an odd sequence x(i). the transform is defined below at
c     output parameter x.
c
c     sint is the unnormalized inverse of itself since a call of sint
c     followed by another call of sint will multiply the input sequence
c     x by 8*(n+1).
c
c     the array wsave which is used by subroutine sint must be
c     initialized by calling subroutine sinti(n,wsave).
c
c     input parameters
c
c     n       n is the length of x. n must be odd and the method is
c             most efficient when n+1 is a product of small primes
c
c     x       an array which contains the sequence to be transformed
c
c                   ************important*************
c
c                   x must be dimensioned at least n+1
c
c     wsave   a work array which must be dimensioned at least 3n+18
c             in the program that calls sint. the wsave array must be
c             initialized by calling subroutine sinti(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c
c     output parameters
c
c     x       for i=1,...,n
c
c                  x(i)= the sum from k=1 to k=n
c
c                      4*x(k)*sin(k*i*pi/(n+1))
c
c                  a call of sint followed by another call of
c                  sint will multiply the sequence x by 8(n+1)
c                  hence sint is the unnormalized inverse
c                  of itself.
c
c     wsave   contains initialization calculations which must not be
c             destroyed between calls of sint.
c
      dimension       x(1)       ,wsave(1)
c
      if (n .gt. 1) go to 101
      x(1) = 4.*x(1)
      return
  101 np1 = n+1
      ns2 = np1/2
      ns2m = ns2-1
      iw1 = ns2+1
      x1 = x(1)
      x(1) = 0.
      do 102 k=1,ns2m
         kc = np1-k
         t1 = x1-x(kc)
         t2 = wsave(k)*(x1+x(kc))
         x1 = x(k+1)
         x(k+1) = t1+t2
         x(kc+1) = t2-t1
  102 continue
      x(ns2+1) = 4.*x1
      call rfftf (np1,x,wsave(iw1))
      x(1) = .5*x(1)
      do 103 i=3,n,2
         x(i-1) = x(i+1)
         x(i) = x(i)+x(i-2)
  103 continue
      return
      end
c----------------------------------------------------------------------
        subroutine ezpower(n,r,p,wsave)
        implicit real*4 (a-h,o-z)
c       modified from ezfftf to just return power
      dimension       r(*),wsave(*),p(*)
c
        if (n.le.2) then
C        type 2,n
2       format(' a power spectrum for',i5,' points is silly !')
        return
        endif
  102 ns2 = n/2
c      modn = n-ns2-ns2
c      if (modn .ne. 0) go to 105
      iw1 = n+1
      do 103 i=1,n
         wsave(i) = r(i)
  103 continue
      call rfftf (n,wsave,wsave(iw1))
      cf = 1./float(n)
        cf2=cf*cf
c      azero = .5*cf*wsave(1)
        p(1)=.25*cf2*wsave(1)*wsave(1)
c      a(ns2) = .5*cf*wsave(2)
c      b(ns2) = 0.
        p(ns2+1)=cf2*.25*wsave(2)*wsave(2)
      ns2m = ns2-1
      do 104 i=1,ns2m
c         a(i) = cf*wsave(2*i+1)
c         b(i) = cf*wsave(2*i+2)
        p(i+1)=cf2*(wsave(2*i+2)**2+wsave(2*i+1)**2)
  104 continue
      return
c  105 iw1 = n+n+1
c      do 106 i=1,n
c         wsave(2*i-1) = r(i)
c         wsave(2*i) = 0.
c  106 continue
c      call cfftf (n,wsave,wsave(iw1))
c      cf = 2./float(n)
c       cf2=cf*cf
cc      azero = .5*cf*wsave(1)
c      nm1s2 = (n-1)/2
c      do 107 i=1,nm1s2
cc         a(i) = cf*wsave(2*i+1)
cc         b(i) = -cf*wsave(2*i+2)
c       p(i)=cf2*(wsave(2*i+2)**2+wsave(2*i+1)**2)
c  107 continue
c      return
      end
c----------------------------------------------------------------------
        subroutine ezsc(n,r,b,a,wsave)
        implicit real*4 (a-h,o-z)
c--     modified from ezfftf to return sine and cosine series for a real
c-      variable in a format suitable for ana, this means that first
c       cosine is set to mean (azero here) and first and last sine are 0
c       the sine corresponds to b and the cosine to a
c       we also assume that n is always even (checked by calling routine
c       which is sc in fun2.mar
c-      also note that call to ezffti to set up wsave is made in fun2
      dimension       r(*),wsave(*),b(*),a(*)
c
        if (n.le.2) then
C        type 2,n
2       format(' a sine-cosine fit for',i5,' points is silly !')
        return
        endif
  102 ns2 = n/2
      iw1 = n+1
      do 103 i=1,n
         wsave(i) = r(i)
  103 continue
      call rfftf (n,wsave,wsave(iw1))
      cf = 1./float(n)
      a(1) = .5*cf*wsave(1)     !this was azero
        b(1)=0.0                !first and last sine terms are 0
      a(ns2+1) = .5*cf*wsave(2)
      b(ns2+1) = 0.
      do 104 i=2,ns2
         a(i) = cf*wsave(2*i-1)
         b(i) = cf*wsave(2*i)
  104 continue
      return
      end
c----------------------------------------------------------------------
        subroutine ezdsc(n,r,b,a,wsave)
        implicit real*4 (a-h,o-z)
c--     modified from ezfftf to return sine and cosine series for a real
c-      variable, this is different and more sensible than the sc
c-	routine for multi-dimensional extensions
c-	the lengths are each n/2 and the first value in the sine
c-	array is actually the highest cosine coefficient
c       we also assume that n is always even (checked by calling routine
c       which is dsc in fun2.mar
c-	result is unnormalized to save a little time, divide by
c-	n to convert most to the actual sine/cosine coefficients
c-      also note that call to ezffti to set up wsave is made in fun2
      dimension       r(*),wsave(*),b(*),a(*)
c
        if (n.le.2) then
C        type 2,n
2       format(' a sine-cosine fit for',i5,' points is silly !')
        return
        endif
  102 ns2 = n/2
      iw1 = n+1
      do 103 i=1,n
         wsave(i) = r(i)
  103 continue
      call rfftf (n,wsave,wsave(iw1))
      do 104 i=1,ns2
         a(i) = wsave(2*i-1)
         b(i) = wsave(2*i)
  104 continue
      return
      end
c----------------------------------------------------------------------
      subroutine ezscb (n,r,a,b,wsave)
        implicit real*4 (a-h,o-z)
c--     this is a modified version of ezfftb for ana
c       it is similar but the sine and cosine series are 1 element longer
c       with c(1)=azero and s(1)=0 always assumed
c       also we assume n is even
      dimension       r(n)       ,a(1)       ,b(1)       ,wsave(1)
c
      if (n .gt. 1) go to 101
      r(1) = a(1)
      return
  101 if (n .gt. 2) go to 102
      r(1) = a(1)+a(2)
      r(2) = a(1)-a(2)
      return
  102 ns2 = n/2
      iw1 = n+1
      do 103 i=2,ns2
         r(2*i-1) = .5*a(i)
         r(2*i) = .5*b(i)
  103 continue
      r(1) = a(1)
      r(2) = a(ns2+1)
      call rfftb (n,r,wsave(iw1))
      return
      end
c----------------------------------------------------------------------
      subroutine ezfftb (n,r,azero,a,b,wsave)
c     subroutine ezfftb computes a real perodic sequence from its
c     fourier coefficients (fourier synthesis). the transform is
c     defined below at output parameter r. ezfftb is a simplified
c     version of rfftb. it is not as fast as rfftb since scaling and
c     initialization are computed for each transform. the repeated
c     initialization can be suppressed by removeing the statment
c     ( call ezffti(n,wsave) ) from both ezfftf and ezfftb and inserting
c     it at the appropriate place in your program.
c
c     input parameters
c
c     n       the length of the array r. ezfftb is about twice as fast
c             for even n as it is for odd n. also ezfftb is more
c             efficient when n is a product of small primes.
c
c     azero   the constant fourier coefficient
c
c     a,b     arrays which contain the remaining fourier coefficients
c             these arrays are not destroyed.
c
c             the length of these arrays depends on whether n is even or
c             odd.
c
c             if n is even n/2    locations are required
c             if n is odd (n-1)/2 locations are required
c
c     wsave   a work array whose length depends on whether n is even or
c             odd.
c
c                  if n is even 3.5*n+15 locations are required
c                  if n is odd  6*n+15   locations are required
c
c
c     output parameters
c
c     r       if n is even define kmax=n/2
c             if n is odd  define kmax=(n-1)/2
c
c             then for i=1,...,n
c
c                  r(i)=azero plus the sum from k=1 to k=kmax of
c
c                  a(k)*cos(k*(i-1)*2*pi/n)+b(k)*sin(k*(i-1)*2*pi/n)
c
c     ********************* complex notation **************************
c
c             for j=1,...,n
c
c             r(j) equals the sum from k=-kmax to k=kmax of
c
c                  c(k)*exp(i*k*(j-1)*2*pi/n)
c
c             where
c
c                  c(k) = .5*cmplx(a(k),-b(k))   for k=1,...,kmax
c
c                  c(-k) = conjg(c(k))
c
c                  c(0) = azero
c
c                       and i=sqrt(-1)
c
c     *************** amplitude - phase notation ***********************
c
c             for i=1,...,n
c
c             r(i) equals azero plus the sum from k=1 to k=kmax of
c
c                  alpha(k)*cos(k*(i-1)*2*pi/n+beta(k))
c
c             where
c
c                  alpha(k) = sqrt(a(k)*a(k)+b(k)*b(k))
c
c                  cos(beta(k))=a(k)/alpha(k)
c
c                  sin(beta(k))=-b(k)/alpha(k)
c
      dimension       r(1)       ,a(1)       ,b(1)       ,wsave(1)
c
      if (n .gt. 1) go to 101
      r(1) = azero
      return
  101 if (n .gt. 2) go to 102
      r(1) = azero+a(1)
      r(2) = azero-a(1)
      return
  102 ns2 = n/2
c
c     to supress repeated initialization remove the following statment
c     ( call ezffti(n,wsave)) from both this program and ezfftf and
c     insert it at the appropriate place in your program.
c
      call ezffti (n,wsave)
c
      modn = n-ns2-ns2
      if (modn .ne. 0) go to 104
      iw1 = n+1
      ns2m = ns2-1
      do 103 i=1,ns2m
         r(2*i+1) = .5*a(i)
         r(2*i+2) = .5*b(i)
  103 continue
      r(1) = azero
      r(2) = a(ns2)
      call rfftb (n,r,wsave(iw1))
      return
  104 iw1 = n+n+1
      nm1s2 = (n-1)/2
      do 105 i=1,nm1s2
         ic = n-i
         wsave(2*i+1) = a(i)
         wsave(2*i+2) = -b(i)
         wsave(2*ic+1) = a(i)
         wsave(2*ic+2) = b(i)
  105 continue
      wsave(1) = azero+azero
      wsave(2) = 0.
      call cfftb (n,wsave,wsave(iw1))
      do 106 i=1,n
         r(i) = .5*wsave(2*i-1)
  106 continue
      return
      end
c-----------------------------------------------------------------------------
c-----------------------------------------------------------------------------
c the dp versions follow
      subroutine ezfftid (n,wsave)
        implicit real*8 (a-h,o-z)
c---    8/19/82  this version modified to work only with even n
      dimension       wsave(*)
c
      if (n .le. 2) return
      ns2 = n/2
      iw1 = n+1
      call rfftid (n,wsave(iw1))
      return
      end
      subroutine cfftid (n,wsave)
        implicit real*8 (a-h,o-z)
      dimension       wsave(1)
c
      if (n .eq. 1) return
      iw1 = n+n+1
      iw2 = iw1+n+n
      call cffti1d (n,wsave(iw1),wsave(iw2))
      return
      end
      subroutine cffti1d (n,wa,ifac)
        implicit real*8 (a-h,o-z)
      dimension       wa(1)      ,ifac(1)    ,ntryh(4)
      data ntryh(1),ntryh(2),ntryh(3),ntryh(4)/3,4,2,5/
      nl = n
      nf = 0
      j = 0
  101 j = j+1
      if (j-4) 102,102,103
  102 ntry = ntryh(j)
      go to 104
  103 ntry = ntry+2
  104 nq = nl/ntry
      nr = nl-ntry*nq
      if (nr) 101,105,101
  105 nf = nf+1
      ifac(nf+2) = ntry
      nl = nq
      if (ntry .ne. 2) go to 107
      if (nf .eq. 1) go to 107
      do 106 i=2,nf
         ib = nf-i+2
         ifac(ib+2) = ifac(ib+1)
  106 continue
      ifac(3) = 2
  107 if (nl .ne. 1) go to 104
      ifac(1) = n
      ifac(2) = nf
      tpi = 8.*atan(1.)
      arg1 = tpi/float(n)
      dc = dcos(arg1)
      ds = dsin(arg1)
      wa(1) = dc
      wa(2) = ds
      nt = n+n
      do 108 i=4,nt,2
         wa(i-1) = dc*wa(i-3)-ds*wa(i-2)
         wa(i) = ds*wa(i-3)+dc*wa(i-2)
  108 continue
      return
      end
      subroutine cfftbd (n,c,wsave)
        implicit real*8 (a-h,o-z)
      dimension       c(1)       ,wsave(1)
c
      if (n .eq. 1) return
      iw1 = n+n+1
      iw2 = iw1+n+n
      iw3 = iw2+1
      call cfftb1d (n,c,wsave,wsave(iw1),wsave(iw2))
      return
      end
      subroutine cfftb1d (n,c,ch,wa,ifac)
        implicit real*8 (a-h,o-z)
      dimension       ch(1)      ,c(1)       ,wa(1)      ,ifac(1)
      nf = ifac(2)
      l1 = 1
      do 105 k1=1,nf
         ip = ifac(k1+2)
         l2 = ip*l1
         ido = n/l2
         idot = ido+ido
         idl1 = idot*l1
         if (ip .ne. 4) go to 101
         ix2 = l1+l1
         ix3 = ix2+l1
         call passb4d (idot,l1,idl1,ix2,ix3,c,c,c,ch,ch,wa,wa,wa)
         go to 104
  101    if (ip .ne. 2) go to 102
         call passb2d (idot,l1,idl1,c,c,c,ch,ch,wa)
         go to 104
  102    if (ip .ne. 3) go to 103
         ix2 = l1+l1
         call passb3d (idot,l1,idl1,ix2,c,c,c,ch,ch,wa,wa)
         go to 104
  103    call passbd (idot,ip,l1,idl1,c,c,c,ch,ch,wa)
  104    l1 = l2
  105 continue
      return
      end
      subroutine passb2d (ido,l1,idl1,cc,c1,c2,ch,ch2,wa1)
        implicit real*8 (a-h,o-z)
      dimension       cc(ido,2,l1)           ,c1(ido,l1,2)           ,
     1                c2(idl1,2) ,ch(ido,l1,2)           ,ch2(idl1,2),
     2                wa1(l1,1)
      idot = ido/2
      if (ido .lt. l1) go to 103
      do 102 k=1,l1
         do 101 i=1,ido
            ch(i,k,1) = cc(i,1,k)+cc(i,2,k)
            ch(i,k,2) = cc(i,1,k)-cc(i,2,k)
  101    continue
  102 continue
      go to 106
  103 do 105 i=1,ido
         do 104 k=1,l1
            ch(i,k,1) = cc(i,1,k)+cc(i,2,k)
            ch(i,k,2) = cc(i,1,k)-cc(i,2,k)
  104    continue
  105 continue
  106 do 107 ik=1,idl1
         c2(ik,1) = ch2(ik,1)
  107 continue
      do 108 k=1,l1
         c1(1,k,2) = ch(1,k,2)
         c1(2,k,2) = ch(2,k,2)
  108 continue
      if (ido .eq. 2) return
      if (idot .lt. l1) go to 111
      do 110 k=1,l1
         do 109 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
  109    continue
  110 continue
      return
  111 do 113 i=4,ido,2
         do 112 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
  112    continue
  113 continue
      return
      end
      subroutine passb3d (ido,l1,idl1,ix2,cc,c1,c2,ch,ch2,wa1,wa2)
        implicit real*8 (a-h,o-z)
      dimension       cc(ido,3,l1)           ,c1(ido,l1,3)           ,
     1                c2(idl1,3) ,ch(ido,l1,3)           ,ch2(idl1,3),
     2                wa1(l1,1)  ,wa2(ix2,1)
      data taur,taui /-.5,.866025403784439/
      idot = ido/2
      if (ido .lt. l1) go to 103
      do 102 k=1,l1
         do 101 i=1,ido
            ch(i,k,1) = cc(i,1,k)
            ch(i,k,2) = cc(i,2,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)-cc(i,3,k)
  101    continue
  102 continue
      go to 106
  103 do 105 i=1,ido
         do 104 k=1,l1
            ch(i,k,1) = cc(i,1,k)
            ch(i,k,2) = cc(i,2,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)-cc(i,3,k)
  104    continue
  105 continue
c
  106 do 107 ik=1,idl1
         c2(ik,1) = ch2(ik,1)+ch2(ik,2)
         c2(ik,2) = ch2(ik,1)+taur*ch2(ik,2)
         c2(ik,3) = taui*ch2(ik,3)
  107 continue
      do 108 ik=2,idl1,2
         ch2(ik-1,2) = c2(ik-1,2)-c2(ik,3)
         ch2(ik-1,3) = c2(ik-1,2)+c2(ik,3)
  108 continue
      do 109 ik=2,idl1,2
         ch2(ik,2) = c2(ik,2)+c2(ik-1,3)
         ch2(ik,3) = c2(ik,2)-c2(ik-1,3)
  109 continue
      do 111 j=2,3
         do 110 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  110    continue
  111 continue
      if (ido .eq. 2) return
      if (idot-1 .lt. l1) go to 114
      do 113 k=1,l1
         do 112 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)+
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)-
     1                    wa2(ix2,i-2)*ch(i,k,3)
  112    continue
  113 continue
      return
  114 do 116 i=4,ido,2
         do 115 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)+
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)-
     1                    wa2(ix2,i-2)*ch(i,k,3)
  115    continue
  116 continue
      return
      end
      subroutine passb4d (ido,l1,idl1,ix2,ix3,cc,c1,c2,ch,ch2,wa1,wa2,
     1                   wa3)
        implicit real*8 (a-h,o-z)
      dimension       cc(ido,4,l1)           ,c1(ido,l1,4)           ,
     1                c2(idl1,4) ,ch(ido,l1,4)           ,ch2(idl1,4),
     2                wa1(l1,1)  ,wa2(ix2,1) ,wa3(ix3,1)
      idot = ido/2
c
      if (ido .lt. l1) go to 106
      do 103 k=1,l1
         do 101 i=2,ido,2
            ch(i-1,k,4) = cc(i,4,k)-cc(i,2,k)
  101    continue
         do 102 i=2,ido,2
            ch(i,k,4) = cc(i-1,2,k)-cc(i-1,4,k)
  102    continue
  103 continue
      do 105 k=1,l1
         do 104 i=1,ido
            ch(i,k,2) = cc(i,1,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)+cc(i,4,k)
            ch(i,k,1) = cc(i,1,k)-cc(i,3,k)
  104    continue
  105 continue
      go to 111
  106 do 108 i=2,ido,2
         do 107 k=1,l1
            ch(i-1,k,4) = cc(i,4,k)-cc(i,2,k)
            ch(i,k,4) = cc(i-1,2,k)-cc(i-1,4,k)
  107    continue
  108 continue
      do 110 i=1,ido
         do 109 k=1,l1
            ch(i,k,2) = cc(i,1,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)+cc(i,4,k)
            ch(i,k,1) = cc(i,1,k)-cc(i,3,k)
  109    continue
  110 continue
  111 do 112 ik=1,idl1
         c2(ik,1) = ch2(ik,2)+ch2(ik,3)
  112 continue
      do 113 ik=1,idl1
         ch2(ik,3) = ch2(ik,2)-ch2(ik,3)
  113 continue
      do 114 ik=1,idl1
         ch2(ik,2) = ch2(ik,1)+ch2(ik,4)
  114 continue
      do 115 ik=1,idl1
         ch2(ik,4) = ch2(ik,1)-ch2(ik,4)
  115 continue
      do 117 j=2,4
         do 116 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  116    continue
  117 continue
      if (ido .eq. 2) return
      if (idot .lt. l1) go to 120
      do 119 k=1,l1
         do 118 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)+
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)-
     1                    wa2(ix2,i-2)*ch(i,k,3)
            c1(i,k,4) = wa3(ix3-1,i-2)*ch(i,k,4)+
     1                  wa3(ix3,i-2)*ch(i-1,k,4)
            c1(i-1,k,4) = wa3(ix3-1,i-2)*ch(i-1,k,4)-
     1                    wa3(ix3,i-2)*ch(i,k,4)
  118    continue
  119 continue
      return
  120 do 122 i=4,ido,2
         do 121 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)+wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)-
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)+
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)-
     1                    wa2(ix2,i-2)*ch(i,k,3)
            c1(i,k,4) = wa3(ix3-1,i-2)*ch(i,k,4)+
     1                  wa3(ix3,i-2)*ch(i-1,k,4)
            c1(i-1,k,4) = wa3(ix3-1,i-2)*ch(i-1,k,4)-
     1                    wa3(ix3,i-2)*ch(i,k,4)
  121    continue
  122 continue
      return
      end
      subroutine passbd (ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa)
        implicit real*8 (a-h,o-z)
      dimension       ch(ido,l1,ip)          ,cc(ido,ip,l1)          ,
     1                c1(ido,l1,ip)          ,wa(1)      ,c2(idl1,ip),
     2                ch2(idl1,ip)
      idot = ido/2
      nt = ip*idl1
      ipp2 = ip+2
      ipph = (ip+1)/2
      l1t = l1+l1
c
      if (ido .lt. l1) go to 106
      do 103 j=2,ipph
         jc = ipp2-j
         do 102 k=1,l1
            do 101 i=1,ido
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
  101       continue
  102    continue
  103 continue
      do 105 k=1,l1
         do 104 i=1,ido
            ch(i,k,1) = cc(i,1,k)
  104    continue
  105 continue
      go to 112
  106 do 109 j=2,ipph
         jc = ipp2-j
         do 108 i=1,ido
            do 107 k=1,l1
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
  107       continue
  108    continue
  109 continue
      do 111 i=1,ido
         do 110 k=1,l1
            ch(i,k,1) = cc(i,1,k)
  110    continue
  111 continue
  112 do 113 ik=1,idl1
         c2(ik,1) = ch2(ik,1)
  113 continue
      idj = 0
      do 115 j=2,ipph
         jc = ipp2-j
         idj = idj+idl1
         do 114 ik=1,idl1
            c2(ik,j) = ch2(ik,1)+wa(idj-1)*ch2(ik,2)
            c2(ik,jc) = wa(idj)*ch2(ik,ip)
  114    continue
  115 continue
      do 117 j=2,ipph
         do 116 ik=1,idl1
            c2(ik,1) = c2(ik,1)+ch2(ik,j)
  116    continue
  117 continue
c
      idl = 0
      do 120 l=2,ipph
         lc = ipp2-l
         idl = idl+idl1
         idlj = idl
         do 119 j=3,ipph
            jc = ipp2-j
            idlj = idlj+idl
            if (idlj .gt. nt) idlj = idlj-nt
            war = wa(idlj-1)
            wai = wa(idlj)
            do 118 ik=1,idl1
               c2(ik,l) = c2(ik,l)+war*ch2(ik,j)
               c2(ik,lc) = c2(ik,lc)+wai*ch2(ik,jc)
  118       continue
  119    continue
  120 continue
c
      do 122 j=2,ipph
         jc = ipp2-j
         do 121 ik=2,idl1,2
            ch2(ik-1,j) = c2(ik-1,j)-c2(ik,jc)
            ch2(ik-1,jc) = c2(ik-1,j)+c2(ik,jc)
  121    continue
  122 continue
      do 124 j=2,ipph
         jc = ipp2-j
         do 123 ik=2,idl1,2
            ch2(ik,j) = c2(ik,j)+c2(ik-1,jc)
            ch2(ik,jc) = c2(ik,j)-c2(ik-1,jc)
  123    continue
  124 continue
c
      do 126 j=2,ip
         do 125 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  125    continue
  126 continue
      if (ido .eq. 2) return
      idj = 0
      if (idot .gt. l1) go to 130
      do 129 j=2,ip
         idj = idj+l1t
         idij = 0
         do 128 i=4,ido,2
            idij = idij+idj
            do 127 k=1,l1
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
  127       continue
  128    continue
  129 continue
      return
  130 do 134 j=2,ip
         idj = idj+l1t
         do 133 k=1,l1
            idij = 0
            do 131 i=4,ido,2
               idij = idij+idj
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
  131       continue
            idij = 0
            do 132 i=4,ido,2
               idij = idij+idj
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
  132       continue
  133    continue
  134 continue
      return
      end
      subroutine cfftfd (n,c,wsave)
        implicit real*8 (a-h,o-z)
      dimension       c(1)       ,wsave(1)
c
      if (n .eq. 1) return
      iw1 = n+n+1
      iw2 = iw1+n+n
      call cfftf1d (n,c,wsave,wsave(iw1),wsave(iw2))
      return
      end
      subroutine cfftf1d (n,c,ch,wa,ifac)
        implicit real*8 (a-h,o-z)
      dimension       ch(1)      ,c(1)       ,wa(1)      ,ifac(1)
      nf = ifac(2)
      l1 = 1
      do 105 k1=1,nf
         ip = ifac(k1+2)
         l2 = ip*l1
         ido = n/l2
         idot = ido+ido
         idl1 = idot*l1
         if (ip .ne. 4) go to 101
         ix2 = l1+l1
         ix3 = ix2+l1
         call passf4d (idot,l1,idl1,ix2,ix3,c,c,c,ch,ch,wa,wa,wa)
         go to 104
  101    if (ip .ne. 2) go to 102
         call passf2d (idot,l1,idl1,c,c,c,ch,ch,wa)
         go to 104
  102    if (ip .ne. 3) go to 103
         ix2 = l1+l1
         call passf3d (idot,l1,idl1,ix2,c,c,c,ch,ch,wa,wa)
         go to 104
  103    call passfd (idot,ip,l1,idl1,c,c,c,ch,ch,wa)
  104    l1 = l2
  105 continue
      return
      end
      subroutine passf2d (ido,l1,idl1,cc,c1,c2,ch,ch2,wa1)
        implicit real*8 (a-h,o-z)
      dimension       cc(ido,2,l1)           ,c1(ido,l1,2)           ,
     1                c2(idl1,2) ,ch(ido,l1,2)           ,ch2(idl1,2),
     2                wa1(l1,1)
      idot = ido/2
      if (ido .lt. l1) go to 103
      do 102 k=1,l1
         do 101 i=1,ido
            ch(i,k,1) = cc(i,1,k)+cc(i,2,k)
            ch(i,k,2) = cc(i,1,k)-cc(i,2,k)
  101    continue
  102 continue
      go to 106
  103 do 105 i=1,ido
         do 104 k=1,l1
            ch(i,k,1) = cc(i,1,k)+cc(i,2,k)
            ch(i,k,2) = cc(i,1,k)-cc(i,2,k)
  104    continue
  105 continue
  106 do 107 ik=1,idl1
         c2(ik,1) = ch2(ik,1)
  107 continue
      do 108 k=1,l1
         c1(1,k,2) = ch(1,k,2)
         c1(2,k,2) = ch(2,k,2)
  108 continue
      if (ido .eq. 2) return
      if (idot .lt. l1) go to 111
      do 110 k=1,l1
         do 109 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
  109    continue
  110 continue
      return
  111 do 113 i=4,ido,2
         do 112 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
  112    continue
  113 continue
      return
      end
      subroutine passf3d (ido,l1,idl1,ix2,cc,c1,c2,ch,ch2,wa1,wa2)
        implicit real*8 (a-h,o-z)
      dimension       cc(ido,3,l1)           ,c1(ido,l1,3)           ,
     1                c2(idl1,3) ,ch(ido,l1,3)           ,ch2(idl1,3),
     2                wa1(l1,1)  ,wa2(ix2,1)
      data taur,taui /-.5,-.866025403784439/
      idot = ido/2
      if (ido .lt. l1) go to 103
      do 102 k=1,l1
         do 101 i=1,ido
            ch(i,k,1) = cc(i,1,k)
            ch(i,k,2) = cc(i,2,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)-cc(i,3,k)
  101    continue
  102 continue
      go to 106
  103 do 105 i=1,ido
         do 104 k=1,l1
            ch(i,k,1) = cc(i,1,k)
            ch(i,k,2) = cc(i,2,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)-cc(i,3,k)
  104    continue
  105 continue
c
  106 do 107 ik=1,idl1
         c2(ik,1) = ch2(ik,1)+ch2(ik,2)
         c2(ik,2) = ch2(ik,1)+taur*ch2(ik,2)
         c2(ik,3) = taui*ch2(ik,3)
  107 continue
      do 108 ik=2,idl1,2
         ch2(ik-1,2) = c2(ik-1,2)-c2(ik,3)
         ch2(ik-1,3) = c2(ik-1,2)+c2(ik,3)
  108 continue
      do 109 ik=2,idl1,2
         ch2(ik,2) = c2(ik,2)+c2(ik-1,3)
         ch2(ik,3) = c2(ik,2)-c2(ik-1,3)
  109 continue
      do 111 j=2,3
         do 110 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  110    continue
  111 continue
      if (ido .eq. 2) return
      if (idot-1 .lt. l1) go to 114
      do 113 k=1,l1
         do 112 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)-
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)+
     1                    wa2(ix2,i-2)*ch(i,k,3)
  112    continue
  113 continue
      return
  114 do 116 i=4,ido,2
         do 115 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)-
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)+
     1                    wa2(ix2,i-2)*ch(i,k,3)
  115    continue
  116 continue
      return
      end
      subroutine passf4d (ido,l1,idl1,ix2,ix3,cc,c1,c2,ch,ch2,wa1,wa2,
     1                   wa3)
        implicit real*8 (a-h,o-z)
      dimension       cc(ido,4,l1)           ,c1(ido,l1,4)           ,
     1                c2(idl1,4) ,ch(ido,l1,4)           ,ch2(idl1,4),
     2                wa1(l1,1)  ,wa2(ix2,1) ,wa3(ix3,1)
      idot = ido/2
c
      if (ido .lt. l1) go to 106
      do 103 k=1,l1
         do 101 i=2,ido,2
            ch(i-1,k,4) = cc(i,2,k)-cc(i,4,k)
  101    continue
         do 102 i=2,ido,2
            ch(i,k,4) = cc(i-1,4,k)-cc(i-1,2,k)
  102    continue
  103 continue
      do 105 k=1,l1
         do 104 i=1,ido
            ch(i,k,2) = cc(i,1,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)+cc(i,4,k)
            ch(i,k,1) = cc(i,1,k)-cc(i,3,k)
  104    continue
  105 continue
      go to 111
  106 do 108 i=2,ido,2
         do 107 k=1,l1
            ch(i,k,4) = cc(i-1,4,k)-cc(i-1,2,k)
            ch(i-1,k,4) = cc(i,2,k)-cc(i,4,k)
  107    continue
  108 continue
      do 110 i=1,ido
         do 109 k=1,l1
            ch(i,k,2) = cc(i,1,k)+cc(i,3,k)
            ch(i,k,3) = cc(i,2,k)+cc(i,4,k)
            ch(i,k,1) = cc(i,1,k)-cc(i,3,k)
  109    continue
  110 continue
  111 do 112 ik=1,idl1
         c2(ik,1) = ch2(ik,2)+ch2(ik,3)
  112 continue
      do 113 ik=1,idl1
         ch2(ik,3) = ch2(ik,2)-ch2(ik,3)
  113 continue
      do 114 ik=1,idl1
         ch2(ik,2) = ch2(ik,1)+ch2(ik,4)
  114 continue
      do 115 ik=1,idl1
         ch2(ik,4) = ch2(ik,1)-ch2(ik,4)
  115 continue
      do 117 j=2,4
         do 116 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  116    continue
  117 continue
      if (ido .eq. 2) return
      if (idot .lt. l1) go to 120
      do 119 k=1,l1
         do 118 i=4,ido,2
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)-
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)+
     1                    wa2(ix2,i-2)*ch(i,k,3)
            c1(i,k,4) = wa3(ix3-1,i-2)*ch(i,k,4)-
     1                  wa3(ix3,i-2)*ch(i-1,k,4)
            c1(i-1,k,4) = wa3(ix3-1,i-2)*ch(i-1,k,4)+
     1                    wa3(ix3,i-2)*ch(i,k,4)
  118    continue
  119 continue
      return
  120 do 122 i=4,ido,2
         do 121 k=1,l1
            c1(i,k,2) = wa1(l1-1,i-2)*ch(i,k,2)-wa1(l1,i-2)*ch(i-1,k,2)
            c1(i-1,k,2) = wa1(l1-1,i-2)*ch(i-1,k,2)+
     1                    wa1(l1,i-2)*ch(i,k,2)
            c1(i,k,3) = wa2(ix2-1,i-2)*ch(i,k,3)-
     1                  wa2(ix2,i-2)*ch(i-1,k,3)
            c1(i-1,k,3) = wa2(ix2-1,i-2)*ch(i-1,k,3)+
     1                    wa2(ix2,i-2)*ch(i,k,3)
            c1(i,k,4) = wa3(ix3-1,i-2)*ch(i,k,4)-
     1                  wa3(ix3,i-2)*ch(i-1,k,4)
            c1(i-1,k,4) = wa3(ix3-1,i-2)*ch(i-1,k,4)+
     1                    wa3(ix3,i-2)*ch(i,k,4)
  121    continue
  122 continue
      return
      end
      subroutine passfd (ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa)
        implicit real*8 (a-h,o-z)
      dimension       ch(ido,l1,ip)          ,cc(ido,ip,l1)          ,
     1                c1(ido,l1,ip)          ,wa(1)      ,c2(idl1,ip),
     2                ch2(idl1,ip)
      idot = ido/2
      nt = ip*idl1
      ipp2 = ip+2
      ipph = (ip+1)/2
      l1t = l1+l1
c
      if (ido .lt. l1) go to 106
      do 103 j=2,ipph
         jc = ipp2-j
         do 102 k=1,l1
            do 101 i=1,ido
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
  101       continue
  102    continue
  103 continue
      do 105 k=1,l1
         do 104 i=1,ido
            ch(i,k,1) = cc(i,1,k)
  104    continue
  105 continue
      go to 112
  106 do 109 j=2,ipph
         jc = ipp2-j
         do 108 i=1,ido
            do 107 k=1,l1
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
  107       continue
  108    continue
  109 continue
      do 111 i=1,ido
         do 110 k=1,l1
            ch(i,k,1) = cc(i,1,k)
  110    continue
  111 continue
  112 do 113 ik=1,idl1
         c2(ik,1) = ch2(ik,1)
  113 continue
      idj = 0
      do 115 j=2,ipph
         jc = ipp2-j
         idj = idj+idl1
         do 114 ik=1,idl1
            c2(ik,j) = ch2(ik,1)+wa(idj-1)*ch2(ik,2)
            c2(ik,jc) = -wa(idj)*ch2(ik,ip)
  114    continue
  115 continue
      do 117 j=2,ipph
         do 116 ik=1,idl1
            c2(ik,1) = c2(ik,1)+ch2(ik,j)
  116    continue
  117 continue
c
      idl = 0
      do 120 l=2,ipph
         lc = ipp2-l
         idl = idl+idl1
         idlj = idl
         do 119 j=3,ipph
            jc = ipp2-j
            idlj = idlj+idl
            if (idlj .gt. nt) idlj = idlj-nt
            war = wa(idlj-1)
            wai = wa(idlj)
            do 118 ik=1,idl1
               c2(ik,l) = c2(ik,l)+war*ch2(ik,j)
               c2(ik,lc) = c2(ik,lc)-wai*ch2(ik,jc)
  118       continue
  119    continue
  120 continue
c
      do 122 j=2,ipph
         jc = ipp2-j
         do 121 ik=2,idl1,2
            ch2(ik-1,j) = c2(ik-1,j)-c2(ik,jc)
            ch2(ik-1,jc) = c2(ik-1,j)+c2(ik,jc)
  121    continue
  122 continue
      do 124 j=2,ipph
         jc = ipp2-j
         do 123 ik=2,idl1,2
            ch2(ik,j) = c2(ik,j)+c2(ik-1,jc)
            ch2(ik,jc) = c2(ik,j)-c2(ik-1,jc)
  123    continue
  124 continue
c
      do 126 j=2,ip
         do 125 k=1,l1
            c1(1,k,j) = ch(1,k,j)
            c1(2,k,j) = ch(2,k,j)
  125    continue
  126 continue
      if (ido .eq. 2) return
      idj = 0
      if (idot .gt. l1) go to 130
      do 129 j=2,ip
         idj = idj+l1t
         idij = 0
         do 128 i=4,ido,2
            idij = idij+idj
            do 127 k=1,l1
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)+wa(idij)*ch(i,k,j)
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)-wa(idij)*ch(i-1,k,j)
  127       continue
  128    continue
  129 continue
      return
  130 do 134 j=2,ip
         idj = idj+l1t
         do 133 k=1,l1
            idij = 0
            do 131 i=4,ido,2
               idij = idij+idj
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)+wa(idij)*ch(i,k,j)
  131       continue
            idij = 0
            do 132 i=4,ido,2
               idij = idij+idj
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)-wa(idij)*ch(i-1,k,j)
  132       continue
  133    continue
  134 continue
      return
      end
      subroutine cosqid (n,wsave)
        implicit real*8 (a-h,o-z)
      dimension       wsave(1)
c
      iw1 = n+1
      pih = 2.*atan(1.)
      dt = pih/float(n)
      dc = dcos(dt)
      ds = dsin(dt)
      wsave(1) = dc
      wsk = ds
      do 101 k=2,n
         wsave(k) = dc*wsave(k-1)-ds*wsk
         wsk = ds*wsave(k-1)+dc*wsk
  101 continue
      call rfftid (n,wsave(iw1))
      return
      end
      subroutine cosqfd (n,x,wsave)
        implicit real*8 (a-h,o-z)
      dimension       x(1)       ,wsave(1)
      data tsq2/2.82842712474619/
c
      if (n .gt. 2) go to 101
      tx = x(1)+x(1)
      tsqx = tsq2*x(2)
      x(1) = tx+tsqx
      x(2) = tx-tsqx
      return
  101 iw1 = n+1
      iw2 = n/2+iw1
      call cosqf1d (n,x,wsave,wsave(iw1),wsave(iw2))
      return
      end
      subroutine cosqf1d (n,x,w,w1,xh)
        implicit real*8 (a-h,o-z)
      dimension       x(1)       ,xh(1)      ,w1(1)      ,w(1)
      ns2 = n/2
      nm1 = n-1
      np2 = n+2
      do 101 k=2,ns2
         kc = np2-k
         xh(k) = x(k)+x(kc)
         xh(kc) = x(k)-x(kc)
  101 continue
      xh(ns2+1) = x(ns2+1)+x(ns2+1)
      do 102 k=2,ns2
         kc = np2-k
         x(k) = w(k-1)*xh(kc)+w(kc-1)*xh(k)
         x(kc) = w(k-1)*xh(k)-w(kc-1)*xh(kc)
  102 continue
      x(ns2+1) = w(ns2)*xh(ns2+1)
      call rfftfd (n,x,w1)
      xn = x(2)
      do 103 i=3,nm1,2
         x(i-1) = x(i)+x(i+1)
         x(i) = x(i)-x(i+1)
  103 continue
      x(n) = xn
      return
      end
      subroutine cosqbd (n,x,wsave)
        implicit real*8 (a-h,o-z)
      dimension       x(1)       ,wsave(1)
      data tsq2/2.82842712474619/
c
      if (n .gt. 2) go to 101
      x1 = 4.*(x(1)+x(2))
      x(2) = tsq2*(x(1)-x(2))
      x(1) = x1
      return
  101 iw1 = n+1
      iw2 = n/2+iw1
      call cosqb1d (n,x,wsave,wsave(iw1),wsave(iw2))
      return
      end
      subroutine cosqb1d (n,x,w,w1,xh)
        implicit real*8 (a-h,o-z)
      dimension       x(1)       ,w(1)       ,w1(1)      ,xh(1)
      ns2 = n/2
      nm1 = n-1
      np2 = n+2
      xn = x(n)
      do 101 ii=3,nm1,2
         i = np2-ii
         x(i+1) = x(i-1)-x(i)
         x(i) = x(i)+x(i-1)
  101 continue
      x(1) = x(1)+x(1)
      x(2) = xn+xn
      call rfftbd (n,x,w1)
      do 102 k=2,ns2
         kc = np2-k
         xh(k) = w(k-1)*x(kc)+w(kc-1)*x(k)
         xh(kc) = w(k-1)*x(k)-w(kc-1)*x(kc)
  102 continue
      x(ns2+1) = w(ns2)*(x(ns2+1)+x(ns2+1))
      do 103 k=2,ns2
         kc = np2-k
         x(k) = xh(k)+xh(kc)
         x(kc) = xh(k)-xh(kc)
  103 continue
      x(1) = x(1)+x(1)
      return
      end
      subroutine costid (n,wsave)
        implicit real*8 (a-h,o-z)
      dimension       wsave(1)
c
      if (n .le. 3) return
      nm1 = n-1
      ns2 = nm1/2
      ns2m = ns2-1
      iw1 = ns2+1
      pi = 4.*atan(1.)
      dt = pi/float(nm1)
      dcs = dcos(dt)
      dss = dsin(dt)
      wsave(1) = dss+dss
      wck = dcs+dcs
      if (ns2m .lt. 2) go to 102
      do 101 k=2,ns2m
         wsave(k) = dss*wck+dcs*wsave(k-1)
         wck = dcs*wck-dss*wsave(k-1)
  101 continue
  102 call rfftid (nm1,wsave(iw1))
      return
      end
      subroutine costd (n,x,wsave)
        implicit real*8 (a-h,o-z)
c
      dimension       x(1)       ,wsave(1)
c
      nm1 = n-1
      np1 = n+1
      ns2 = nm1/2
      iw1 = ns2+1
      if (n .gt. 3) go to 101
      x1px3 = x(1)+x(3)
      tx13 = x1px3+x1px3
      fx2 = 4.*x(2)
      x(2) = 2.*(x(1)-x(3))
      x(1) = tx13+fx2
      x(3) = tx13-fx2
      return
  101 c1 = x(1)-x(n)
      do 102 k=2,ns2
         kc = nm1-k
         ks = ns2-k
         c1 = c1+wsave(ks+1)*(x(k)-x(kc+2))
  102 continue
      c1 = c1+c1
      x(1) = x(1)+x(n)
      do 103 k=2,ns2
         kc = np1-k
         t1 = x(k)+x(kc)
         t2 = wsave(k-1)*(x(k)-x(kc))
         x(k) = t1-t2
         x(kc) = t1+t2
  103 continue
      x(ns2+1) = x(ns2+1)+x(ns2+1)
      call rfftfd (nm1,x,wsave(iw1))
      cn = x(2)
      x(2) = c1
      do 104 i=4,nm1,2
         x(i) = x(i)+x(i-2)
  104 continue
      x(n) = cn
      return
      end
      subroutine rfftid (n,wsave)
        implicit real*8 (a-h,o-z)
      dimension       wsave(1)
c
      ns2 = n/2
      nqm = (ns2-1)/2
      tpi = 8.*atan(1.)
      dt = tpi/float(n)
      dc = dcos(dt)
      ds = dsin(dt)
      wsave(1) = dc
      wsave(ns2-1) = ds
      if (nqm .lt. 2) go to 102
      do 101 k=2,nqm
         kc = ns2-k
         wsave(k) = dc*wsave(k-1)-ds*wsave(kc+1)
         wsave(kc) = ds*wsave(k-1)+dc*wsave(kc+1)
  101 continue
  102 iw1 = ns2+1
      call cfftid (ns2,wsave(iw1))
      return
      end
      subroutine rfftfd (n,r,wsave)
        implicit real*8 (a-h,o-z)
      dimension       r(2,1)     ,wsave(1)
c
      if (n .gt. 2) go to 101
      r1 = 2.*(r(1,1)+r(2,1))
      r(2,1) = 2.*(r(1,1)-r(2,1))
      r(1,1) = r1
      return
  101 iw1 = n/2+1
      call rfftf1d (n,r,wsave(iw1),wsave)
      return
      end
      subroutine rfftf1d (n,x,xh,w)
        implicit real*8 (a-h,o-z)
      dimension       x(2,1)     ,xh(2,1)    ,w(1)
      ns2 = n/2
      ns2p2 = ns2+2
      nq = ns2/2
      ipar = ns2-nq-nq
      nqm = nq
      if (ipar .eq. 0) nqm = nqm-1
      nqp = nqm+1
      call cfftfd (ns2,x,xh)
      if (nqp .lt. 2) go to 107
      do 101 k=2,nqp
         kc = ns2p2-k
         xh(1,kc) = x(2,k)+x(2,kc)
         xh(2,kc) = x(1,kc)-x(1,k)
  101 continue
      do 102 k=2,nqp
         kc = ns2p2-k
         xh(1,k) = x(1,k)+x(1,kc)
         xh(2,k) = x(2,k)-x(2,kc)
  102 continue
      do 103 k=2,nqp
         kc = ns2p2-k
         x(1,kc) = w(k-1)*xh(1,kc)+w(kc-1)*xh(2,kc)
         x(2,kc) = w(k-1)*xh(2,kc)-w(kc-1)*xh(1,kc)
  103 continue
      do 104 k=2,nqp
         kc = ns2p2-k
         xh(1,kc) = x(1,kc)
         xh(2,kc) = x(2,kc)
  104 continue
      do 105 k=2,nqp
         kc = ns2p2-k
         x(1,kc) = xh(1,k)-xh(1,kc)
         x(2,kc) = xh(2,k)-xh(2,kc)
  105 continue
      do 106 k=2,nqp
         kc = ns2p2-k
         x(1,k) = xh(1,k)+xh(1,kc)
         x(2,k) = -xh(2,k)-xh(2,kc)
  106 continue
      if (ipar .ne. 0) go to 108
  107 x(1,nqp+1) = x(1,nqp+1)+x(1,nqp+1)
      x(2,nqp+1) = x(2,nqp+1)+x(2,nqp+1)
  108 xhold1 = x(2,1)+x(1,1)
      xhold2 = x(1,1)-x(2,1)
      x(2,1) = xhold2+xhold2
      x(1,1) = xhold1+xhold1
      return
      end
      subroutine rfftbd (n,r,wsave)
        implicit real*8 (a-h,o-z)
      dimension       r(2,1)     ,wsave(1)
c
      if (n .gt. 2) go to 101
      r1 = r(1,1)+r(2,1)
      r(2,1) = r(1,1)-r(2,1)
      r(1,1) = r1
      return
  101 iw1 = n/2+1
      call rfftb1d (n,r,wsave(iw1),wsave)
      return
      end
      subroutine rfftb1d (n,x,xh,w)
        implicit real*8 (a-h,o-z)
      dimension       x(2,1)     ,xh(2,1)    ,w(1)
      ns2 = n/2
      ns2p2 = ns2+2
      nq = ns2/2
      ipar = ns2-nq-nq
      nqm = nq
      if (ipar .eq. 0) nqm = nqm-1
      nqp = nqm+1
      xhold1 = x(1,1)-x(2,1)
      x(1,1) = x(2,1)+x(1,1)
      x(2,1) = xhold1
      if (ipar .ne. 0) go to 101
      x(1,nqp+1) = x(1,nqp+1)+x(1,nqp+1)
      x(2,nqp+1) = x(2,nqp+1)+x(2,nqp+1)
  101 if (nqp .lt. 2) go to 108
      do 102 k=2,nqp
         kc = ns2p2-k
         xh(1,k) = x(1,k)+x(1,kc)
         xh(2,k) = x(2,kc)-x(2,k)
  102 continue
      do 103 k=2,nqp
         kc = ns2p2-k
         xh(1,kc) = x(1,k)-x(1,kc)
         xh(2,kc) = -x(2,k)-x(2,kc)
  103 continue
      do 104 k=2,nqp
         kc = ns2p2-k
         x(1,kc) = xh(1,kc)
         x(2,kc) = xh(2,kc)
  104 continue
      do 105 k=2,nqp
         kc = ns2p2-k
         xh(1,kc) = w(k-1)*x(1,kc)-w(kc-1)*x(2,kc)
         xh(2,kc) = w(k-1)*x(2,kc)+w(kc-1)*x(1,kc)
  105 continue
      do 106 k=2,nqp
         kc = ns2p2-k
         x(1,k) = xh(1,k)-xh(2,kc)
         x(2,k) = xh(2,k)+xh(1,kc)
  106 continue
      do 107 k=2,nqp
         kc = ns2p2-k
         x(1,kc) = xh(1,k)+xh(2,kc)
         x(2,kc) = xh(1,kc)-xh(2,k)
  107 continue
  108 call cfftbd (ns2,x,xh)
      return
      end
      subroutine sinqid (n,wsave)
        implicit real*8 (a-h,o-z)
      dimension       wsave(1)
c
      call cosqid (n,wsave)
      return
      end
      subroutine sinqfd (n,x,wsave)
        implicit real*8 (a-h,o-z)
      dimension       x(1)       ,wsave(1)
c
      ns2 = n/2
      do 101 k=1,ns2
         kc = n-k
         xhold = x(k)
         x(k) = x(kc+1)
         x(kc+1) = xhold
  101 continue
      call cosqfd (n,x,wsave)
      do 102 k=2,n,2
         x(k) = -x(k)
  102 continue
      return
      end
      subroutine sinqbd (n,x,wsave)
        implicit real*8 (a-h,o-z)
      dimension       x(1)       ,wsave(1)
c
      ns2 = n/2
      do 101 k=2,n,2
         x(k) = -x(k)
  101 continue
      call cosqbd (n,x,wsave)
      do 102 k=1,ns2
         kc = n-k
         xhold = x(k)
         x(k) = x(kc+1)
         x(kc+1) = xhold
  102 continue
      return
      end
      subroutine sintid (n,wsave)
        implicit real*8 (a-h,o-z)
      dimension       wsave(1)
c
      if (n .le. 1) return
      np1 = n+1
      ns2 = np1/2
      ns2m = ns2-1
      iw1 = ns2+1
      pi = 4.*atan(1.)
      dt = pi/float(np1)
      dcs = dcos(dt)
      dss = dsin(dt)
      wsave(1) = dss+dss
      wck = dcs+dcs
      if (ns2m .lt. 2) go to 102
      do 101 k=2,ns2m
         wsave(k) = dss*wck+dcs*wsave(k-1)
         wck = dcs*wck-dss*wsave(k-1)
  101 continue
  102 call rfftid (np1,wsave(iw1))
      return
      end
      subroutine sintd (n,x,wsave)
        implicit real*8 (a-h,o-z)
      dimension       x(1)       ,wsave(1)
c
      if (n .gt. 1) go to 101
      x(1) = 4.*x(1)
      return
  101 np1 = n+1
      ns2 = np1/2
      ns2m = ns2-1
      iw1 = ns2+1
      x1 = x(1)
      x(1) = 0.
      do 102 k=1,ns2m
         kc = np1-k
         t1 = x1-x(kc)
         t2 = wsave(k)*(x1+x(kc))
         x1 = x(k+1)
         x(k+1) = t1+t2
         x(kc+1) = t2-t1
  102 continue
      x(ns2+1) = 4.*x1
      call rfftfd (np1,x,wsave(iw1))
      x(1) = .5*x(1)
      do 103 i=3,n,2
         x(i-1) = x(i+1)
         x(i) = x(i)+x(i-2)
  103 continue
      return
      end
        subroutine ezpowerd(n,r,p,wsave)
        implicit real*8 (a-h,o-z)
c     ******************************************************************
c       modified from ezfftf to just return power
      dimension       r(*),wsave(*),p(*)
c
        if (n.le.2) then
C        type 2,n
2       format(' a power spectrum for',i5,' points is silly !')
        return
        endif
  102 ns2 = n/2
c      modn = n-ns2-ns2
c      if (modn .ne. 0) go to 105
      iw1 = n+1
      do 103 i=1,n
         wsave(i) = r(i)
  103 continue
      call rfftfd (n,wsave,wsave(iw1))
      cf = 1./float(n)
        cf2=cf*cf
c      azero = .5*cf*wsave(1)
        p(1)=.25*cf2*wsave(1)*wsave(1)
c      a(ns2) = .5*cf*wsave(2)
c      b(ns2) = 0.
        p(ns2+1)=cf2*.25*wsave(2)*wsave(2)
      ns2m = ns2-1
      do 104 i=1,ns2m
c         a(i) = cf*wsave(2*i+1)
c         b(i) = cf*wsave(2*i+2)
        p(i+1)=cf2*(wsave(2*i+2)**2+wsave(2*i+1)**2)
  104 continue
      return
c  105 iw1 = n+n+1
c      do 106 i=1,n
c         wsave(2*i-1) = r(i)
c         wsave(2*i) = 0.
c  106 continue
c      call cfftfd (n,wsave,wsave(iw1))
c      cf = 2./float(n)
c       cf2=cf*cf
cc      azero = .5*cf*wsave(1)
c      nm1s2 = (n-1)/2
c      do 107 i=1,nm1s2
cc         a(i) = cf*wsave(2*i+1)
cc         b(i) = -cf*wsave(2*i+2)
c       p(i)=cf2*(wsave(2*i+2)**2+wsave(2*i+1)**2)
c  107 continue
c      return
      end
c----------------------------------------------------------------------
        subroutine ezscd(n,r,b,a,wsave)
        implicit real*8 (a-h,o-z)
c     ******************************************************************
c--     modified from ezfftf to return sine and cosine series for a real
c-      variable in a format suitable for ana, this means that first
c       cosine is set to mean (azero here) and first and last sine are 0
c       the sine corresponds to b and the cosine to a
c       we also assume that n is always even (checked by calling routine
c       which is sc in fun2.mar
c-      also note that call to ezffti to set up wsave is made in fun2
      dimension       r(*),wsave(*),b(*),a(*)
c
        if (n.le.2) then
C        type 2,n
2       format(' a sine-cosine fit for',i5,' points is silly !')
        return
        endif
  102 ns2 = n/2
      iw1 = n+1
      do 103 i=1,n
         wsave(i) = r(i)
  103 continue
      call rfftfd (n,wsave,wsave(iw1))
      cf = 1./float(n)
        cf2=cf*cf
      a(1) = .5*cf*wsave(1)     !this was azero
        b(1)=0.0                !first and last sine terms are 0
      a(ns2+1) = .5*cf*wsave(2)
      b(ns2+1) = 0.
      do 104 i=2,ns2
         a(i) = cf*wsave(2*i-1)
         b(i) = cf*wsave(2*i)
  104 continue
      return
      end
c----------------------------------------------------------------------
        subroutine ezdscd(n,r,b,a,wsave)
        implicit real*8 (a-h,o-z)
      dimension       r(*),wsave(*),b(*),a(*)
c
        if (n.le.2) then
C        type 2,n
2       format(' a sine-cosine fit for',i5,' points is silly !')
        return
        endif
  102 ns2 = n/2
      iw1 = n+1
      do 103 i=1,n
         wsave(i) = r(i)
  103 continue
      call rfftfd (n,wsave,wsave(iw1))
      do 104 i=1,ns2
         a(i) = wsave(2*i-1)
         b(i) = wsave(2*i)
  104 continue
      return
      end
c----------------------------------------------------------------------
      subroutine ezscbd(n,r,a,b,wsave)
        implicit real*8 (a-h,o-z)
c--     this is a modified version of ezfftb for ana
c       it is similar but the sine and cosine series are 1 element longer
c       with c(1)=azero and s(1)=0 always assumed
c       also we assume n is even
      dimension       r(n)       ,a(1)       ,b(1)       ,wsave(1)
c
      if (n .gt. 1) go to 101
      r(1) = a(1)
      return
  101 if (n .gt. 2) go to 102
      r(1) = a(1)+a(2)
      r(2) = a(1)-a(2)
      return
  102 ns2 = n/2
      iw1 = n+1
      do 103 i=2,ns2
         r(2*i-1) = .5*a(i)
         r(2*i) = .5*b(i)
  103 continue
      r(1) = a(1)
      r(2) = a(ns2+1)
      call rfftbd (n,r,wsave(iw1))
      return
      end
      subroutine ezfftbd (n,r,azero,a,b,wsave)
        implicit real*8 (a-h,o-z)
      dimension       r(1)       ,a(1)       ,b(1)       ,wsave(1)
c
      if (n .gt. 1) go to 101
      r(1) = azero
      return
  101 if (n .gt. 2) go to 102
      r(1) = azero+a(1)
      r(2) = azero-a(1)
      return
  102 ns2 = n/2
c
c     to supress repeated initialization remove the following statment
c     ( call ezffti(n,wsave)) from both this program and ezfftf and
c     insert it at the appropriate place in your program.
c
      call ezfftid (n,wsave)
c
      modn = n-ns2-ns2
      if (modn .ne. 0) go to 104
      iw1 = n+1
      ns2m = ns2-1
      do 103 i=1,ns2m
         r(2*i+1) = .5*a(i)
         r(2*i+2) = .5*b(i)
  103 continue
      r(1) = azero
      r(2) = a(ns2)
      call rfftbd (n,r,wsave(iw1))
      return
  104 iw1 = n+n+1
      nm1s2 = (n-1)/2
      do 105 i=1,nm1s2
         ic = n-i
         wsave(2*i+1) = a(i)
         wsave(2*i+2) = -b(i)
         wsave(2*ic+1) = a(i)
         wsave(2*ic+2) = b(i)
  105 continue
      wsave(1) = azero+azero
      wsave(2) = 0.
      call cfftbd (n,wsave,wsave(iw1))
      do 106 i=1,n
         r(i) = .5*wsave(2*i-1)
  106 continue
      return
      end
     
