! ================================================================================ !
! Programa em FORTRAN para calcular os Harmônicos de Fourier em uma série temporal !
! Enilson Palmeira Cavalanti - Em agosto de 2021.                                  !
! ================================================================================ !
      integer t, n
      real media, sd2m
      dimension x(1000), w(500), fi(500), alfa(500), beta(500)
      dimension amp(500), sd2(500), pvar(500), h(1000,500)  
! Lendo a série do arquivo "serie.dat" (informar o número de elementos da série "T")
      open(2, file='fourier.dat', status='old')
      read(2,*) t
      n=t/2
      write(*,10) t
10    format('Série',5x,'T=',i5)
      do i=1,t 
         read(2,*) x(i)
         write(*,20) i, x(i)
20       format(5x,i5,5x,f8.3)
      end do
      close(2) 
      write(*,*)
      media=0.0
      do i=1,t
         media=media+x(i)
      end do
      media=media/t
      write(*,25) media
25    format(8x,'Média= ',f8.3)
      do j=1,n
         w(j)=2*3.1416*j/t
         alfa(j)=0
         beta(j)=0
      end do
      write(*,30) n
30    format('Harmonicos',5x,'N=',i5,/,6x,'j',8x,'w(j)',4x,'alfa(j)',3x,
     *'beta(j)',2x,'aplit(j)',3x,'fase(j)',4x,'var(j)',3x,'%var(j)')
      sd2m=0.0
      do j=1,n
         do i=1,t
            alfa(j)=alfa(j)+x(i)*cos(w(j)*i)
            beta(j)=beta(j)+x(i)*sin(w(j)*i)
         end do
         if(j .lt. n)then
            alfa(j)=alfa(j)*2/t
            beta(j)=beta(j)*2/t 
         else
            alfa(j)=alfa(j)/t
            beta(j)=0.0
         end if
         amp(j)=sqrt(alfa(j)**2+beta(j)**2)
         fi(j)=atan(beta(j)/alfa(j))
         if(j .lt. n)then
            sd2(j)=(amp(j)**2)/2
         else
            sd2(j)=amp(j)**2
         end if
         sd2m=sd2m+sd2(j)
      end do
      do j=1,n
         pvar(j)=100*sd2(j)/sd2m
         write(*,40) j, w(j), alfa(j), beta(j), amp(j), fi(j),
     *   sd2(j), pvar(j)
40       format(2x,i5,3x,7(1x,f9.4))
      end do
! ********************
      do j=1,n
         do i=1,t
            h(i,j)=alfa(j)*cos(w(j)*i)+beta(j)*sin(w(j)*i)    
         end do
      end do
      write(*,50)
50    format(3x,"N",3x,"Obs.",8x,"Est.")
      do i=1,t
        soma1=0
        soma2=0
        do j=1,n
           soma1=soma1+h(i,j)
        end do
         write(*,60) i, x(i), media+soma1
60       format(2x,i3, f8.3,2x,f8.3)
      end do   
      end

 
