! ----------------------------------------------------------------------
! *  Discrete Fourier series (sine transform)                          *
! ----------------------------------------------------------------------
   program dft_w
!
!  ****  initialize variables  ****
!
      implicit none

      real( 8 ), allocatable:: a( : ), f( : ), t( : )
	  
      real( 8 ):: Tperi
      real( 8 ):: dt
	  
      integer::   nmax, msam
      integer::   j, n
!
!  **** set pi  **** 
!
      real( 8 ), parameter:: pi = 4.0d0 * atan( 1.0d0 )
!
!   ****  sampling freqency for one period  ****
!
      write ( 6, * ) '# of sampling for one period'
      read  ( 5, * ) msam
      if ( msam == 0 ) then 
        write ( 6, * ) 'We set msam = 100' 
        msam = 100
      end if
!
!  ****  set max order of Fourier coefficeint  ****
!
      nmax = msam / 2
      write ( 6, * ) 'nmax is ', nmax
!
!  ****  allocate dimension variables  ****
!
      allocate( a( 0:nmax ), f( 0:msam ), t( 0:msam ) )
!
!  ****   read sine.dat  ****
!
      open ( 10, file = 'sw.dat' )
      open ( 11, file = 'four.dat' )

      Tperi = 2.0d0 * pi
      dt   = Tperi / dfloat( msam )
      mtime = msam

      do j = 0, mtime
        read  ( 10, * ) t( j ), f( j )
      end do 
!
!  ****  set a_n = 0 for summation  ****
!
      do n = 1,nmax
        a( n )  = 0.0
      end do
!
!  ****  Fourier transform  ****
!	  
      do n = 1, nmax
        do j = 0, msam - 1
          a( n ) = a( n ) + f( j ) * sin( 2.0d0 * pi * dfloat( n ) * t( j ) / Tperi ) * dt
        end do
        a( n ) = a( n ) * 2.0d0 / Tperi
        write ( 11, * ) n, a( n )
      end do

      close ( 10 )
      close ( 11 )

      stop
	  
   end program dft_w