program dft
c      Parameter estimation for periodic signal
c      ******************************************************************
c      mflg=1 Periodic case, mflg2=Non-periodic spatial data
c      N: Total number of data points should be even and equal to $2*N+2$.
c      t1: Start time of the data (usually 0)
c      t2: End time of the data
c      f(i): i=1,2,...,2*N+2. Amplitudes of the signal (values at 
c      equally-spaced times with $f(1) the value at $t=t1$ and $f(2*N+2)$
c      the value at $t=t2$)
c      For a periodic function, program finds amplitudes of cosine and sine 
c      components and the frequencies to machine precision
c      In this implementation, $f(i)$ is being generated thorugh an 
c      example function. Replace this line by your own function or
c      Comment this line out and supply $f(i)$ through dft.dat
c      Redirect output to file: Example  dft.e >> dft.out
c      ******************************************************************
       real*8, allocatable :: A(:),K(:,:),udelt(:), f(:),
     & fdelt(:)
       real*8 t1,t2,Ts,sum1,sum2,pi,c0,c1,t2orig,fdnorm,fden
       integer i,j,k0,m,N,mflg,N0
       open(unit=5,file='dft.dat',status='old')
       read(5,*) mflg,N,t1,t2
       t2orig=t2
       pi=4.0d0*datan(1.0d0)
       N0=N
       Ts=(t2-t1)/(2*N+1)
       allocate(f(2*N+2))
c       read(5,*) (f(i),i=1,2*N+2)
        print*,'        Time                     Values of f'
       do 15 i=1,2*N+2
        f(i)=2.0d0+1.0d0*dcos(2.0d0*Pi*(t1+(i-1)*Ts))+2.0d0*
     & dsin(2.4d0*Pi*(t1+(i-1)*Ts))
c        f(i)=0.1d0+0.3*dcos(50.0d0*pi*(t1+(i-1)*Ts)+7.0d0*pi/18.0d0)+
c     & 0.2d0*dcos(70.7d0*pi*(t1+(i-1)*Ts)+pi/3.0d0)+
c     & 0.7*dcos(71.7d0*pi*(t1+(i-1)*Ts)+4.0*pi/9.0d0)+
c     & dcos(100.0d0*pi*(t1+(i-1)*Ts))+
c     & 0.5*dcos(173.2d0*pi*(t1+(i-1)*Ts)+pi/2.0d0)+
c     & 0.4*dcos(300.0d0*pi*(t1+(i-1)*Ts)+2.0d0*pi/9.0d0)
        print*,i,t1+(i-1)*Ts,f(i)
15     continue
       fden=dsqrt(dot_product(f(1:2*N+2),f(1:2*N+2)))
       if(mflg.eq.1)then
       do 200 N=N0,12,-1
       allocate(K(2*N+2,2*N+2),udelt(2*N+2),A(2*N+2),fdelt(2*N+2))
       t2=t1+(2*N+1)*Ts
       sum1=0.0d0
       do 20 m=2,2*N+1
       sum1=sum1+f(m)
20     continue
       A(2*N+1)=(sum1+(f(1)+f(2*N+2))/2.0d0)/(2*N+1)
       A(2*N+2)=t2-t1
       print*,'T, N ',A(2*N+2),N
       do 40 j=1,N
        sum1=0.0d0
        sum2=0.0d0
         do 30 m=2,2*N+1
         sum1=sum1+f(m)*dcos(2*pi*j*(t1+(m-1)*Ts)/A(2*N+2))
         sum2=sum2+f(m)*dsin(2*pi*j*(t1+(m-1)*Ts)/A(2*N+2))
30       continue
        A(2*j-1)=2.0d0*(sum1+(f(1)*cos(2*pi*j*t1/A(2*N+2))+
     &  f(2*N+2)*cos(2*pi*j*t2/A(2*N+2)))/2.0d0)/(2*N+1)
        A(2*j)=2.0d0*(sum2+(f(1)*sin(2*pi*j*t1/A(2*N+2))+
     &  f(2*N+2)*sin(2*pi*j*t2/A(2*N+2)))/2.0d0)/(2*N+1)
40     continue
       fdnorm=0.0d0
       do 1050 i=1,2*N0+2
       sum1=0.0d0
       do 1040 j=1,N
       sum1=sum1+A(2*j-1)*dcos(2*pi*j/A(2*N+2)*(t1+(i-1)*Ts))+
     & A(2*j)*dsin(2*pi*j/A(2*N+2)*(t1+(i-1)*Ts))
1040   continue
       sum1=sum1+A(2*N+1)
       fdnorm=fdnorm+(sum1-f(i))**2
1050   continue
       fdnorm=dsqrt(fdnorm)/fden
       print*,'fdnorm ',fdnorm
       if(fdnorm.lt.0.4d0)then
       do 80 k0=1,8
        do 60 i=1,2*N+2
         do 50 j=1,N
         K(i,2*j-1)=dcos(2*pi*j*(t1+(i-1)*Ts)/A(2*N+2))
         K(i,2*j)=dsin(2*pi*j*(t1+(i-1)*Ts)/A(2*N+2))
50       continue
         sum1=0.0d0
         sum2=0.0d0
         do 55 m=1,N
         sum1=sum1+2.0d0*pi*m*(t1+(i-1)*Ts)/A(2*N+2)**2*(A(2*m-1)*
     & dsin(2*Pi*m*(t1+(i-1)*Ts)/A(2*N+2))-A(2*m)*
     & dcos(2*Pi*m*(t1+(i-1)*Ts)/A(2*N+2)))
         sum2=sum2+A(2*m-1)*dcos(2*Pi*m*(t1+(i-1)*Ts)/A(2*N+2))+A(2*m)*
     & dsin(2*Pi*m*(t1+(i-1)*Ts)/A(2*N+2))
55       continue
         K(i,2*N+1)=1.0d0
         K(i,2*N+2)=sum1
         fdelt(i)=f(i)-A(2*N+1)-sum2
60      continue
       fdnorm=0.0d0
       do 105 i=1,2*N0+2
       sum1=0.0d0
       do 104 j=1,N
       sum1=sum1+A(2*j-1)*dcos(2*pi*j/A(2*N+2)*(t1+(i-1)*Ts))+
     & A(2*j)*dsin(2*pi*j/A(2*N+2)*(t1+(i-1)*Ts))
104    continue
       sum1=sum1+A(2*N+1)
       fdnorm=fdnorm+(sum1-f(i))**2
105    continue
       fdnorm=dsqrt(fdnorm)/fden
       print*,'fdnorm ',fdnorm
       if(fdnorm.gt.2.0d0)go to 180
       if(fdnorm.le.1.0d-10.and.A(2*N+2).gt.0.0d0.and.
     & A(2*N+2).le.t2orig-t1+1.0d-7)then
       print*,''
       print*,'         c_0                       T '
       print*, A(2*N+1),A(2*N+2)
       print*,''
       print*, '          i    Amplitude (Cos)            Amplitude (Sin
     &)            Frequency (Hz)'
       do 100 i=1,N
       print*,i,(A(2*i-2+j),j=1,2),i/A(2*N+2)
100    continue
       stop
        endif
       call r8mat_fs(2*N+2,K,fdelt,udelt)
        A=A+udelt
80     continue
       endif
180    deallocate(K,udelt,A,fdelt)
200    continue     
c        print*,'Strategy failed to converge'
c        stop
       else
       allocate(K(2*N,2*N),A(2*N),fdelt(2*N))
       c0=f(1)
       c1=f(2*N+2)
       do 400 i=1,2*N
         fdelt(i)=f(i+1)-c0-(c1-c0)*i/(2*N+1)
        do 300 j=1,2*N
         K(i,j)=dsin(pi*i*j/(2*N+1)) 
300     continue
400    continue
       call r8mat_fs(2*N,K,fdelt,A)
        print*,''
        print*,'         c0                        c1 '
        print*,c0,c1
       print*,''
       print*, '          i    Amplitude (Sin)            Frequency (Hz)
     &                Phase'
       do 500 i=1,2*N
       print*,i,A(i),i/(2.0d0*(t2-t1)),-i*pi*t1/(t2-t1)
500    continue
       endif
      stop
      end

      subroutine r8mat_fs ( n, a, b, x )

c*********************************************************************72
c
cc R8MAT_FS factors and solves a system with one right hand side.
c
c  Discussion:
c
c    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
c
c    This routine differs from R8MAT_FSS in two ways:
c    * only one right hand side is allowed;
c    * the input matrix A is not modified.
c
c    This routine uses partial pivoting, but no pivot vector is required.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Modified:
c
c    20 January 2013
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer N, the order of the matrix.
c    N must be positive.
c
c    Input, double precision A(N,N), the coefficient matrix.
c
c    Input, double precision B(N), the right hand side of the linear system.
c
c    Output, double precision X(N), the solution of the linear system.
c
      implicit none

      integer n

      double precision a(n,n)
      double precision a2(n,n)
      double precision b(n)
      integer i
      integer info
      integer ipiv
      integer j
      integer jcol
      integer k
      double precision piv
      double precision temp
      double precision x(n)
c
c  Copy the matrix.
c
      do j = 1, n
        do i = 1, n
          a2(i,j) = a(i,j)
        end do
      end do

      do i = 1, n
        x(i) = b(i)
      end do

      info = 0

      do jcol = 1, n
c
c  Find the maximum element in column I.
c
        piv = abs ( a2(jcol,jcol) )
        ipiv = jcol
        do i = jcol + 1, n
          if ( piv .lt. abs ( a2(i,jcol) ) ) then
            piv = abs ( a2(i,jcol) )
            ipiv = i
          end if
        end do

        if ( piv .eq. 0.0D+00 ) then
          info = jcol
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'R8MAT_FS - Fatal error!'
          write ( *, '(a,i8)' ) '  Zero pivot on step ', info
          stop 1
        end if
c
c  Switch rows JCOL and IPIV, and B.
c
        if ( jcol .ne. ipiv ) then

          do j = 1, n
            temp = a2(jcol,j)
            a2(jcol,j) = a2(ipiv,j)
            a2(ipiv,j) = temp
          end do

          temp = x(jcol)
          x(jcol) = x(ipiv)
          x(ipiv) = temp

        end if
c
c  Scale the pivot row.
c
        do j = jcol + 1, n
          a2(jcol,j) = a2(jcol,j) / a2(jcol,jcol)
        end do
        x(jcol) = x(jcol) / a2(jcol,jcol)
        a2(jcol,jcol) = 1.0D+00
c
c  Use the pivot row to eliminate lower entries in that column.
c
        do i = jcol + 1, n
          if ( a2(i,jcol) .ne. 0.0D+00 ) then
            temp = - a2(i,jcol)
            a2(i,jcol) = 0.0D+00
            do j = jcol + 1, n
              a2(i,j) = a2(i,j) + temp * a2(jcol,j)
            end do
            x(i) = x(i) + temp * x(jcol)
          end if
        end do

      end do
c
c  Back solve.
c
      do k = n, 2, -1
        do i = 1, k - 1
          x(i) = x(i) - a2(i,k) * x(k)
        end do
      end do

      return
      end