! ----------------------------------------------------------------------
! **  Compose sine functions (inverse sine transform)                 **
! ----------------------------------------------------------------------

   program sine_sqw
!
!  ****  initialize variables  ****
!
      implicit none

      real( 8 ), allocatable:: a( : ), f( : ), g( : )
	  
      real( 8 ):: Tperi
      real( 8 ):: t, dt
	  
      integer::   nmax, mcyc, msam, mtime
      integer::   i, j, n, is
!
!  **** set pi  **** 
!	  
      real( 8 ), parameter:: pi = 4.0d0 * atan(1.0d0)
!
!  **** order n  ****
!
      write ( 6, * ) 'order of Fourier transform (default 1)'
      read  ( 5, * ) nmax
      if ( nmax == 0 ) then
        write ( 6, * ) 'we set the order nmax = 1'
        nmax = 1
      end if
!
!   ****  cycles  ****
!
      write ( 6, * ) '# of cycles for plot'
      read ( 5, * ) mcyc
      if ( mcyc == 0 ) then 
        write ( 6, * ) 'We set mcyc = 10' 
        mcyc = 10
      end if
!
!   ****  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
      if ( nmax > msam / 2 ) then 
        write (6,*) 'We set nmax = msam/2 ' 
        nmax = msam / 2
      end if
!
!    ****  # of time steps  ****
!
      mtime = mcyc * msam
!
!  ****  allocate dimension variables  ****
!
      allocate( a( 0:nmax ), f( 0:mtime ), g( 0:mtime ) )
!
!  ****  set period, dt, coefficeint a_n  ****
!
      Tperi = 2.0d0 * pi
      dt = Tperi / dfloat( msam )

      open (  9, file = 'coef1.dat' )

      do n = 1, nmax
        if ( mod( n, 2 ) == 1 ) then
          a( n ) = 4.0d0 /  ( pi *  dfloat( n ) )
         else
          a( n ) = 0.0d0
        end if
        write ( 9, * ) n, a( n )
      end do
!
!  ****   open file  ****
!
      open ( 10, file = 'sinesw.dat' )
      open ( 11, file = 'sw.dat' )
!
!  **** compose sine functions
!
      do j = 0, mtime
        t = dt * dfloat( j )

        do n = 1, nmax
          f( j ) = f( j ) + a( n ) * sin( 2.0d0 * pi * dfloat( n ) * t / Tperi )
        end do

        write ( 10, * ) t, f( j )
      end do


      j = 0
      t = dt * dfloat( j )
      g( j ) = 0.0d0
      write ( 11 , * ) t, g( j )

      do is = 1, mcyc
	  
        do i = 1, msam / 2 - 1
          j = ( is - 1 ) * msam + i
          t = dt * dfloat( j )
          g( j ) = 1.0d0
          write ( 11, * ) t, g( j )
        end do

        j = ( is - 1 ) * msam + msam / 2
        t = dt * dfloat( j )
        g( j ) = 0.0d0
        write ( 11 , * ) t, g( j )

        do i= msam / 2 + 1, msam - 1
          j = ( is - 1 ) * msam + i
          t = dt * dfloat( j )
          g( j ) = -1.0
          write ( 11, * ) t, g( j )
        end do

        j = is * msam
        t = dt * dfloat( j )
        g( j ) = 0.0d0
        write ( 11 , * ) t, g( j )

      end do

      close ( 10 )
      close ( 11 )

      stop
	  
   end program sine_sqw