      implicit real*8(a-h,o-z)
      integer ndim, n
      parameter (ndim=10000)
      parameter (maxsqrt=64,max=8192)
      real*8 a(4,0:ndim), time(0:ndim)
      integer ip(0:maxsqrt+1)
      real*8 w(0:max*5/4-1), af(0:max-1), freqr(0:max-1),
     $     freqi(0:max-1)
c
c     read dimension of matrix to be solved
c
 10   continue
      read(5,'(i5)') n
      if( n .lt. 0 ) then
         write(6,'(a)')
     $        'n is not negative. Please retry to enter n.'
         go to 10
      end if
c
      call signal(n, a, time)
c
c     call FFT
c
      call fft(n, a, w, af, freqr, freqi, ip)
      call trans(n, a, time)
c
      end
c======================================================================
      real*8 function wave1(v,t)
      implicit real*8(a-h,o-z)
      real*8 c, pi, v
      parameter(c=2.9979246d10)
c
      pi = atan(1.0d0)*4.0d0
      wave1 = cos( 2.0d0*pi*c*v*t )
      return
      end
c======================================================================
      subroutine signal(n,a,time)
      implicit real*8(a-h,o-z)
      real*8 a(4,*), time(*)
      parameter(dtime=1.0d-15)
c
      do i = 0 , n-1
         time(i) = dble(i)*dtime
         a(1,i) = wave1(1000.0d0,time(i))
         a(2,i) = wave1(1100.0d0,time(i))*2.0d0
         a(3,i) = wave1(1150.0d0,time(i))*0.5d0
         a(4,i) = a(1,i) + a(2,i) + a(3,i)
c
c$$$         write(6,'(5e18.10)') time(i)/dtime, (a(j,i),j=1,4) 
c
      end do
c
      return
      end
c======================================================================
      subroutine fft(n, a, w, af, freqr, freqi, ip)
      implicit real*8(a-h,o-z)
      parameter (zero=0.0d0,one=1.0d0)
      integer ip(*)
      real*8 w(*), af(*), freqr(*), freqi(*)
      real*8 a(4,*)
c
c     check the number of data
c
      k = 1
      do j = 1 , 13
         k = k*2
         if( k .ge. n ) go to 10
      end do
 10   continue
c      write(6,*) 'number of data for FFT = ',k
c
      do inum = 1 , 4
c     **********
c     cosine FFT
c     **********
         do i = 0 , k-1
            if( i .gt. n-1 ) then
               af(i) = zero
            else
               af(i) = a(inum,i)
            end if
         end do
c
         call ddct(k, 1, af, ip, w)
c
         do i = 0 , n-1
            freqr(i) = af(i)
         end do
c     *********
c     sine FFT
c     *********
         do i = 0 , k-1
            if( i .gt. n-1 ) then
               af(i) = zero
            else
               af(i) = a(inum,i)
            end if
         end do
c
         call ddst(n, 1, af, ip, w)
c
         do i = 0 , n-1
            freqi(i) = af(i)
         end do
c
         do i = 0 , n-1
            a(inum,i) = sqrt( ( freqr(i)**2 + freqi(i)**2 ) / n )
         end do
      end do
c
      return
      end
c======================================================================
      subroutine trans(n,a,time)
      implicit real*8(a-h,o-z)
      real*8 a(4,*)
      real*8 time(*), c
      parameter(c=2.9979246d10)
      parameter(dtime=1.0d-15)
      parameter(zero=0.0d0,one=1.0d0,two=2.0d0)
c
      dt = dtime * c
      totalt = dt * (n-1)
      df = one/two/totalt
c***
c      write(6,*) 'dt = ',dt, ' totalt = ',totalt, ' df = ',df
c***
c
      do i = 0 , n-1
         time(i) = df*i
c
         write(6,'(5e18.10)') time(i), (a(j,i),j=1,4) 
c
      end do
c
      return
      end
