***************************************************************         
*      Programme de calcul de la transformée de Fourier       *         
*      d'un profil de raie d'absorption afin de déterminer    *         
*       la vitesse de rotation projetée v*sin(i)              *  
*      Programme écrit par Y. Nazé & G. Rauw                  *
*                                                             *         
***************************************************************         

*     1st version: April 2007

*     Déclaration des variables 
*     ========================= 
      implicit none
      real*8 lam(10000),flux(10000),dv,fl(32768),flu(65536)
      integer i,j,k,l,n,np,nzero,index,kk
      character*20 name      

*     Ouverture des fichiers données
*     ==============================
      write(6,*)
     %'Enter the name of the file with the normalised line profile:'
      read(5,'(a)')name
      open(unit=8,file=name)
      write(6,*)
     %'Name of the file with the results of the Fourier method?'
      read(5,'(a)')name
      open(unit=9,file=name)

      do 1 i=1,10000
        read(8,*,end=20)kk,lam(i),flux(i)
        flux(i)=1.0-flux(i)
C     Profil normalisé ramené à 0.0
C
        n=i
1     continue
20    write(6,*)'Number of data points read in:',n

*     Préparation des données pour la FFT: nombre de points (np) = puissance 
*     de 2
*     ======================================================================

      np=INT(DLOG(n*1.00D+00)/DLOG(2.0D+00))+6
      np=2**np
      if(np.le.32768) then 
        write(6,*)'The number of data points in the FFT is ',np
      else
        write(6,*)'The number of data points in the FFT exceeds 32768.'
        write(6,*)'Please reduce the number of points in the original
     %  data file or adapt the FORTRAN code. The programme
     %  will stop here'
        goto 1002
      endif
C     nzero = nombre de zéros à ajouter pour avoir un nombre de points 
C     égal à une puissance de 2; dv = pas en longueur d'onde

      nzero=int((np-n)/2)
      dv=(lam(n)-lam(1))/(n-1)

C     Ajout des zéros.

       do 21 k=1,nzero
        fl(k)=0.0D+00
21     continue
       do 22 k=nzero+n+1,np
        fl(k)=0.0D+00
22     continue
       do 23 k=nzero+1,nzero+n
        fl(k)=flux(k-nzero)
23     continue

C      Partie réelle et imaginaire (nulle) de la série de données
       do 28 k=1,np
         flu(2*k-1)=flux(k)
         flu(2*k)=0.0D+00
28     continue

      
      write(*,*)dv

*     FFT avec lemme de Danielson-Lanczos (Numerical Recipes)
*     =======================================================
      write(*,*)'Start of FFT routine'
      call four1(flu,np,1)
      write(*,*)'FFT done'
      do 12 i=2,np/2
      write(9,*)DLOG10(1.D+00*(i-1)/(dv*np)),0.5*DLOG10(flu(2*i)*
     %flu(2*i)+flu(2*i-1)*flu(2*i-1)),
     %(flu(2*i)*flu(2*i)+flu(2*i-1)*flu(2*i-1))**0.5
12    continue
1002  STOP                                                              
      END  
                        
c     SOUS-ROUTINE: fourier
c     ---------------------

      subroutine four1(data,nn,isign)

      integer nn,j,i,isign,istep,m,mmax,n
      real*8 wr,wi,wpr,wpi,wtemp,theta
      real*8 tempi,tempr
      real*8 data(2*nn)

      n = 2*nn
      j = 1
    
      do i=1,n,2
         if (j .gt. i) then
            tempr = data(j)
            tempi = data(j+1)
            data(j) = data(i)
            data(j+1) = data(i+1)
            data(i) = tempr
            data(i+1) = tempi
         end if
         m = nn
30       if ((m .ge. 2) .and. (j .gt. m)) then
            j = j - m
            m = m/2
            go to 30
         end if
         j = j + m
      end do

      mmax = 2

31    if (n .gt. mmax) then
         istep = 2*mmax
         theta = 6.28318530717959D0/(isign*mmax)
         wpr = -2.D0*dsin(0.5D0*theta)**2
         wpi = dsin(theta)
         wr = 1.D0
         wi =  0.D0
         do m=1,mmax,2
            do i=m,n,istep
               j = i + mmax
               tempr = wr*data(j)-wi*data(j+1)
               tempi = wr*data(j+1)+wi*data(j)
               data(j) = data(i) - tempr
               data(j+1) = data(i+1) - tempi
               data(i) = data(i) + tempr
               data(i+1) = data(i+1) + tempi
            end do
            wtemp = wr
            wr = wr * wpr - wi * wpi + wr
            wi = wi * wpr + wtemp * wpi +wi
         end do
         mmax = istep
         go to 31
      end if   

      return
      end   
                             
                                
                                                             

