!
!  FOURIER.f
!
!  Fixed-Format Fortran Source File
!  Generated by PGI Visual Fortran(R)
!  9/27/2013 9:05:33 AM
!
      program FOURIER

      implicit none
	  integer i,k,Nint,Nmax,Nyquist,Ndati,iscratch,Nup,Ndown
      integer itipo
      parameter(Nmax=10001)

	  real Ak(Nmax/2),Bk(Nmax/2),Amp(Nmax/2),Fase(Nmax/2)
      real t(Nmax),eta(Nmax),four(Nmax),scratch,alt(Nmax),scov
      real period,A0,omega,deltat,pigr,periodo_fisico

      write(*,*)'LEGGI DATI DA FILE (0) OPPURE GENERA DATI (1)?'
      read(*,*)itipo

      if(itipo.eq.0) then

! LEGGO IL FILE DATI
!	    open(10,file='Newyearwave2.txt',status='old',access='sequential')
!	    open(10,file='gradino.dat',status='old',access='sequential')
	    open(10,file='serie_generata.dat',
     >       status='unknown',access='sequential')

	  i=1
99    read(10,*,end=66)iscratch
	  i=i+1
	  go to 99
  66  continue
      Ndati=i-1
	  write(*,*)'Ndati=',Ndati



      Nint=Ndati-1
	  if(Nint+1.gt.Nmax) then
        write(*,*)'Ocio che il dimensionamento e sbagliato'
        stop
      end if
      rewind(10)

      do i=1,Nint+1
!        read(10,*)iscratch,scratch,eta(i),t(i)
        read(10,*)t(i),eta(i)
      end do
	  close(10)

	  elseif(itipo.eq.1) then
		omega=2.0
        pigr=acos(-1.)
        Periodo_fisico=2*pigr/omega

	    period=Periodo_fisico*3
        deltat=0.01
        Nint=int(period/deltat)
        Ndati=Nint+1
        write(*,*)Ndati

	    open(25,file='serie_generata.dat',
     >       status='unknown',access='sequential')

      do i=1,Ndati
        t(i)=deltat*(i-1)
        eta(i)=5.0*cos(omega*t(i))+3.0*sin(omega*t(i))
        write(25,*)t(i),eta(i)
      end do
	  close(25)


	  end if

! FINE LETTURA DATI

! CALCOLO VARIABILI DI LAVORO
	  period=t(Nint+1)-t(1)
      omega=2*acos(-1.)/period
	  write(*,*)'Periodo della REGISTRAZIONE= ',period

! CALCOLO A0
      A0=0.
      do i=1,Nint
        A0=A0+(eta(i+1)+eta(i))
      end do
      A0=A0/Nint

!CALCOLO Ak
      Nyquist=Nint/2
      do k=1,Nyquist
        Ak(k)=0.e0
        Bk(k)=0.e0
        do i=1,Nint
          Ak(k)=Ak(k)+(eta(i+1)*cos(k*omega*t(i+1))+
     >                 eta(i)*cos(k*omega*t(i)))
          Bk(k)=Bk(k)+(eta(i+1)*sin(k*omega*t(i+1))+
     >                 eta(i)*sin(k*omega*t(i)))
        end do
        Ak(k)=Ak(k)/Nint
        Bk(k)=Bk(k)/Nint
		Amp(k)=sqrt(Ak(k)**2+Bk(k)**2)
		Fase(k)=atan2(Bk(k),Ak(k))
      end do

! stampa su file
	  open(111,file='Fourier.out',status='unknown',access='sequential')
      write(111,'I4,1x,5(F15.5,1x)')0,omega*0,A0/2.,0.,A0/2.,0.
	  do k=1,Nyquist
        write(111,'I4,1x,5(F15.5,1x)')k,omega*k,Ak(k),Bk(k),
     >                                Amp(k),Fase(k)
	  end do
	  close(111)

	  go to 100

  55  write(*,*)'Dammi Nup'
	  read(*,*) Nup
	  write(*,*)'Dammi Ndown'
	  read(*,*) Ndown

	  open(112,file='Serie_Fourier.out',
     >     status='unknown',access='sequential')
	  do i=1,Nint+1
	    four(i)=A0/2.
        do k=Ndown, Nup   !Nyquist  !/10
!          write(*,*)k
	      four(i)=four(i)+Ak(k)*cos(k*omega*t(i))+
     >                    Bk(k)*sin(k*omega*t(i))
        end do
        write(112,'3(F15.5,1x)')t(i),four(i),eta(i)
	  end do
	  close(112)

	  go to 55

 100  continue

	  call finestra_mobile(Ndati,t,eta,period)

      stop
      end