***************************************************************         
*      Programme de calcul de la transformée de Fourier       *         
*    de la fonction d'élargissement rotationnel d'une raie    *         
*      Programme écrit par Y. Nazé & G. Rauw                  *
*                                                             *         
***************************************************************         

*     1st version: April 2007

*     Déclaration des variables et lecture des données
*     ================================================
      implicit none
      character*40 name
      real*8 eps,lambda0,vsini,lambda,gg(131072),clum
      real*8 vitesse,parenthese,pi,machin,machin2
      real*8 sigma(65537),amplitude(65537),paslambda
      integer ii,nfin,newnbre,nreal,nimagin,nbre,machin3,jj

      clum=299792.458d0
      pi=6.28318530717959d0/2.d0

c     rest wavelength ?
      write(*,*) 'longueur d onde centrale ?'
      read(*,*) lambda0

c assombrissement centre-bord pour les etoiles massives (B-V=-0.3)
      if (lambda0.le.4500) then
          eps=0.32d0
      elseif (lambda0.le.5500) then
          eps=0.28d0
      elseif (lambda0.le.6500) then
          eps=0.24d0
      else
          eps=0.16d0
      endif

c     vsini ?
      write(*,*) 'v sin(i) ?'
      read(*,*) vsini

*     Calcul de la fonction d'élargissement rotationnel
*     =================================================
c     resultat ?
      write(*,*) 'nom du fichier de sortie ?'
      read(*,*) name
      open(unit=11,file=name)
      write(11,*) 'profil de rotation theorique pour vsini=',vsini
     &            ,' et raie a ',lambda0,eps

      nfin=4096
      nbre=512
      paslambda = 2.d0*vsini/(nbre-1)
      machin3 = int((nfin-nbre)/2)
      do ii=1,machin3
         jj=ii-machin3
         vitesse=-vsini+(jj-1)*paslambda
         lambda=lambda0*(1.d0+vitesse/clum)
         gg(ii)=0.d0
         write(11,*) lambda,vitesse,gg(ii)
      enddo
      do ii=machin3+1,machin3+nbre
         jj=ii-machin3
         vitesse=-vsini+(jj-1)*paslambda
         lambda=lambda0*(1.d0+vitesse/clum)
         vitesse=vitesse/vsini
         parenthese=1.d0-(vitesse*vitesse)
         gg(ii)=2.d0*(1.d0-eps)*dsqrt(parenthese)+pi*eps*parenthese/2.d0
         gg(ii)=gg(ii)/(pi*(1.d0-eps/3.d0)*vsini)
         write(11,*) lambda,vitesse,gg(ii)
      enddo
      do ii=machin3+nbre+1,nfin
         jj=ii-machin3
         vitesse=-vsini+(jj-1)*paslambda
         lambda=lambda0*(1.d0+vitesse/clum)
         gg(ii)=0.d0
         write(11,*) lambda,vitesse,gg(ii)
      enddo

      paslambda=paslambda*lambda0/clum
      write(*,*) nfin,' pas de ',paslambda,'angstroms'

      write(11,*) 'amplitude de la FFT du profil de rotation theorique'

*     Calcul de la transformée de Fourier
*     ==================================
c calcul de la fourier du profil
      call realft(gg,nfin,+1)
      newnbre = nfin/2+1
      amplitude(1)=dabs(gg(1))
      sigma(1)=0.d0
      amplitude(newnbre)=dabs(gg(2))
      sigma(newnbre)=1.d0/(2.d0*paslambda)
      newnbre = nfin/2-1
      do ii=1,newnbre
         nreal=2*ii+1
         nimagin=2*ii+2
         amplitude(ii+1)=gg(nreal)*gg(nreal)+
     &                  gg(nimagin)*gg(nimagin)
         amplitude(ii+1)=dsqrt(amplitude(ii+1))
         sigma(ii+1)=(1.d0*ii)/(nfin*paslambda)
      enddo
      do ii=1,newnbre+2
         write(11,*) sigma(ii),amplitude(ii)
      enddo
      end


      SUBROUTINE realft(data,n,isign)
c Calculates the Fourier transform of a set of n real-valued data points. 
c Replaces this data (which is stored in array data(1:n)) by the positive 
c frequency half of its complex Fourier transform. The real-valued first 
c and last components of the complex transform are returned as elements 
c data(1) and data(2), respectively. n must be a power of 2. This routine
c also calculates the inverse transform of a complex data array if it is 
c the transform of real data. (Result in this case must be multiplied by 
c 2/n, and isign must be -1)

      INTEGER isign,n
      real*8 data(n)
      INTEGER i,i1,i2,i3,i4,n2p3
      real*8 c1,c2,h1i,h1r,h2i,h2r,wis,wrs
      real*8 theta,wi,wpi,wpr,wr,wtemp 
      theta=3.141592653589793d0/dble(n/2) 
      c1=0.5d0
      if (isign.eq.1) then
        c2=-0.5d0
        call four1(data,n/2,+1)
      else
        c2=0.5d0
        theta=-theta
      endif
      wpr=-2.0d0*dsin(0.5d0*theta)**2
      wpi=dsin(theta)
      wr=1.0d0+wpr
      wi=wpi
      n2p3=n+3
      do i=2,n/4 
        i1=2*i-1
        i2=i1+1
        i3=n2p3-i2
        i4=i3+1
        wrs=wr
        wis=wi
        h1r=c1*(data(i1)+data(i3)) 
        h1i=c1*(data(i2)-data(i4))
        h2r=-c2*(data(i2)+data(i4))
        h2i=c2*(data(i1)-data(i3))
        data(i1)=h1r+wrs*h2r-wis*h2i 
        data(i2)=h1i+wrs*h2i+wis*h2r
        data(i3)=h1r-wrs*h2r+wis*h2i
        data(i4)=-h1i+wrs*h2i+wis*h2r
        wtemp=wr
        wr=wr*wpr-wi*wpi+wr
        wi=wi*wpr+wtemp*wpi+wi
      enddo
      if (isign.eq.1) then
        h1r=data(1)
        data(1)=h1r+data(2)
        data(2)=h1r-data(2) 
      else
        h1r=data(1)
        data(1)=c1*(h1r+data(2))
        data(2)=c1*(h1r-data(2))
        call four1(data,n/2,-1) 
      endif
      return
      END

      SUBROUTINE four1(data,nn,isign)
c Replaces data(1:2*nn) by its discrete Fourier transform, if isign is 
c input as 1; or replaces data(1:2*nn) by nn times its inverse discrete 
c Fourier transform, if isign is input as -1.
c data is a complex array of length nn or, equivalently, a real array of 
c length 2*nn. nn MUST be an integer power of 2 (this is not checked for!).

      INTEGER isign,nn
      real*8 data(2*nn)
      INTEGER i,istep,j,m,mmax,n
      real*8 tempi,tempr
      real*8 theta,wi,wpi,wpr,wr,wtemp 
      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
        endif
        m=nn
1       if ((m.ge.2).and.(j.gt.m)) then
           j=j-m
           m=m/2
           goto 1
        endif
        j=j+m
      enddo
      mmax=2 
2     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
          enddo 
          wtemp=wr
          wr=wr*wpr-wi*wpi+wr
          wi=wi*wpr+wtemp*wpi+wi
        enddo
        mmax=istep
        goto 2 
      endif 
      return
      END
