      program fft

      implicit none

      integer nmax
      parameter (nmax = 65536)

      integer choice,i1,i,j,n
      double precision pi,hbc,w1,f1,w2,f2,dw,wi,f(1:nmax+1)
      double precision dummy,dt,t1,t2,ti
      double precision l,wmax,sumrul(0:1)
      character*65 infile,outfile
      character*4 eunit

      data pi /3.141592653589793d0/

C nuclear physics version
      data hbc /197.32696d0/
      data eunit /' MeV'/

C cold atoms version
C      data hbc /1.d0/
C      data eunit /'    '/

      write (*,*) 'select transformation:'
      write (*,*) '1 -> from F(w) to f(t) (sine transform)'
      write (*,*) '2 -> from f(t) to F(w) (sine transform)'
      write (*,*) '3 -> from f(t) to F(w) (cosine transform)'
      choice = 2
C begin special  version
      read (*,*) choice
      write (*,*) choice
C end special version

      if (choice.eq.1) then

        write (*,*) 'file containing spectrum F(w): '
        read (*,*) infile
        open (15,file=infile,status='old')

        read (15,*,end=100) w1,f1
        read (15,*,end=100) w2,f2
        dw = w2-w1
        i1 = 1+nint(w1/dw)
        if (abs((i1-1)*dw-w1).gt.1.d-4*dw) then
          write (*,*) 'w = ',w1,' is not a multiple of dw = ',dw
          stop
        end if
        if (i1.ge.nmax)
     1    stop 'lowest frequency (w1) too high.'
        do j = 1,i1-1
          f(j) = 0.d0
        end do
        f(i1) = f1
        f(i1+1) = f2
        do i = i1+2,nmax
          read (15,*,end=100) wi,f(i)
          if (abs((i-1)*dw-wi).gt.1.d-4*dw) then
            write (*,*) 'w = ',wi,' is not a multiple of dw = ',dw
            stop
          end if
        end do
        read (15,*,end=100) wi,dummy
        stop 'input file too long!'
 100    continue
        close (15)
        i = i-1
        if (i.lt.2) stop 'input file too short!'
        write (*,*) i-i1+1,' lines read.'
        n = 2**int(log(dble(i))/log(2.d0))
        if (n.lt.i) n = 2*n
        write (*,*) 'padding with zeros to n = ',n
        do j = i+1,n
          f(j) = 0.d0
        end do

        write (*,*) 'select desired resolution:'
        do i = 0,nint(log(dble(nmax/n))/log(2.d0))
          write (*,*) i,' -> dt = ',pi/(n*2**i*dw)
        end do
        read (*,*) i
        do j = n+1,n*2**i
          f(j) = 0.d0
        end do
        n = n*2**i

        if (f(1).ne.0.d0) then
          f(1) = 0.d0
          write (*,*) 'warning: need to put F(0) to 0'
        end if

        call sinft(f,n)

        write (*,*) 'file for output data f(t): '
        read (*,*) outfile
        open (15,file=outfile,status='new')

        dt = pi/(n*dw)
        do i = 0,n-1
          write (15,'(2(e14.7,:,1x))') i*dt,-2/pi*dw*f(i+1)
        end do
        close (15)

      else if (choice.eq.2) then

        write (*,*) 'file containing the data f(t): '
        infile = 'eld'
C begin special version
        read (*,*) infile
        write (*,*) infile
C end special version
        open (15,file=infile,status='old')
        read (15,*,end=200) t1,f1
        read (15,*,end=200) t2,f2
        dt = t2-t1
        i1 = 1+nint(t1/dt)
        if (abs((i1-1)*dt-t1).gt.1.d-4*dt) then
          write (*,*) 't = ',t1,' is not a multiple of dt = ',dt
          stop
        end if
        if (i1.ge.nmax)
     1    stop 'first time (t1) too large.'
        do j = 1,i1-1
          f(j) = 0.d0
        end do
        f(i1) = f1
        f(i1+1) = f2
        do i = i1+2,nmax
          read (15,*,end=200) ti,f(i)
          if (abs((i-1)*dt-ti).gt.1.d-4*dt) then
            write (*,*) 't = ',ti,' is not a multiple of dt = ',dt
            stop
          end if
        end do
        read (15,*,end=200) ti,dummy
        stop 'input file too long!'
 200    continue
        close (15)
        i = i-1
        if (i.lt.2) stop 'input file too short!'
        write (*,*) i-i1+1,' lines read.'
        n = 2**int(log(dble(i))/log(2.d0))
        if (n.lt.i) n = 2*n
        write (*,*) 'padding with zeros to n = ',n
        do j = i+1,n
          f(j) = 0.d0
        end do

C begin special version
        write (*,*) 'width of the Lorentzian in ',eunit
        read (*,*) l
        l = l/(2*hbc/(dt*i))
        write (*,*) 'width of the Lorentzian in units of ',
     1     real(2*hbc/(dt*i)),eunit
        write (*,*) '(none = 0, recommended = 2) ?'
        l = 2.d0
        read (*,*) l
        write (*,*) l
C end special version
        do j = 1,i
          f(j) = f(j)*exp(-l*j/i)
        end do

        write (*,*) 'select desired resolution:'
        do i = 0,nint(log(dble(nmax/n))/log(2.d0))
          write (*,*) i,' -> dE = ',real(hbc*pi/(n*2**i*dt)),eunit
        end do
        i = i-1
C begin special version
        read (*,*) i
        write (*,*) i
C end special version
        do j = n+1,n*2**i
          f(j) = 0.d0
        end do
        n = n*2**i
        dw = pi/(n*dt)

        if (f(1).ne.0.d0) then
          f(1) = 0.d0
          write (*,*) 'warning: need to put f(0) to 0'
        end if

        call sinft(f,n)

        if (eunit.eq.'    ') then
          write (*,*) 'max. frequency to be written (max. ',
     1      real(hbc*(n-1)*dw),'):'
        else
          write (*,*) 'max. frequency in',eunit,
     1      ' to be written (max. ',real(hbc*(n-1)*dw),'):'
        end if
        wmax = 40.d0
C begin special version
        read (*,*) wmax
        write (*,*) wmax
C end special version 
        wmax = wmax/hbc

        write (*,*) 'output file'
        outfile = 'eld.fft'
C begin special version
        read (*,*) outfile
        write (*,*) outfile
C end special version
        open (15,file=outfile)

        do j = 0,1
          sumrul(j) = 0.d0
        end do
        do i = 0,n-1
          do j = 0,1
            sumrul(j) = 
     1        sumrul(j)+hbc*dw*(i*hbc*dw)**j*(-dt*f(i+1))
          end do
          if (i*dw.le.wmax)
     1      write (15,'(4(e14.7,:,1x))') i*dw*hbc,-dt*f(i+1),
     2        sumrul
        end do
        close (15)
        write (*,*) 'sumrules:'
        write (*,*) 'S_0 =  ',sumrul(0)
        write (*,*) 'S_1 =  ',sumrul(1)

      else if (choice.eq.3) then

        write (*,*) 'file containing the data f(t): '
        read (*,*) infile
        open (15,file=infile,status='old')

        read (15,*,end=400) t1,f1
        read (15,*,end=400) t2,f2
        dt = t2-t1
        i1 = 1+nint(t1/dt)
        if (abs((i1-1)*dt-t1).gt.1.d-4*dt) then
          write (*,*) 't = ',t1,' is not a multiple of dt = ',dt
          stop
        end if
        if (i1.ge.nmax)
     1    stop 'first time (t1) too large.'
        do j = 1,i1-1
          f(j) = 0.d0
        end do
        f(i1) = f1
        f(i1+1) = f2
        do i = i1+2,nmax
          read (15,*,end=400) ti,f(i)
          if (abs((i-1)*dt-ti).gt.1.d-4*dt) then
            write (*,*) 't = ',ti,' is not a multiple of dt = ',dt
            stop
          end if
        end do
        read (15,*,end=400) ti,dummy
        stop 'input file too long!'
 400    continue
        close (15)
        i = i-1
        if (i.lt.2) stop 'input file too short!'
        write (*,*) i-i1+1,' lines read.'
        n = 2**int(log(dble(i))/log(2.d0))
        if (n.lt.i) n = 2*n
        write (*,*) 'padding with zeros to n = ',n
        do j = i+1,n
          f(j) = 0.d0
        end do

        write (*,*) 'relative width of the Lorentzian'
        write (*,*) '(none = 0, recommended = 2) ?'
        read (*,*) l
        do j = 1,i
          f(j) = f(j)*exp(-l*j/i)
        end do

        write (*,*) 'select desired resolution:'
        do i = 0,nint(log(dble(nmax/n))/log(2.d0))
          write (*,*) i,' -> dw = ',pi/(n*2**i*dt)
        end do
        read (*,*) i
        do j = n+1,n*2**i
          f(j) = 0.d0
        end do
        n = n*2**i

        f(n+1) = 0.d0

        call cosft1(f,n)

        dw = pi/(n*dt)

        write (*,*) 'max. frequency (energy) to be written (max. ',
     1    real(n*dw),'):'
        read (*,*) wmax

        write (*,*) 'file for output spectrum F(w): '
        read (*,*) outfile
        open (15,file=outfile,status='new')

        do i = 0,n
          write (15,'(2(e14.7,:,1x))') i*dw,dt*f(i+1)
          if (i*dw.gt.wmax) goto 500
        end do
 500    continue
        close (15)

        sumrul = 0.d0
        do i = 0,n
          sumrul = sumrul+dw*(i*dw)*dt*f(i+1)
        end do
        write (*,*) 'Energy weighted sumrule = ',sumrul 

      end if
      end

      SUBROUTINE sinft(y,n)
      INTEGER n
      DOUBLE PRECISION y(n)
CU    USES realft
      INTEGER j
      DOUBLE PRECISION sum,y1,y2
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      theta=3.141592653589793d0/dble(n)
      wr=1.0d0
      wi=0.0d0
      wpr=-2.0d0*sin(0.5d0*theta)**2
      wpi=sin(theta)
      y(1)=0.0d0
      do 11 j=1,n/2
        wtemp=wr
        wr=wr*wpr-wi*wpi+wr
        wi=wi*wpr+wtemp*wpi+wi
        y1=wi*(y(j+1)+y(n-j+1))
        y2=0.5d0*(y(j+1)-y(n-j+1))
        y(j+1)=y1+y2
        y(n-j+1)=y1-y2
11    continue
      call realft(y,n,+1)
      sum=0.0d0
      y(1)=0.5d0*y(1)
      y(2)=0.0
      do 12 j=1,n-1,2
        sum=sum+y(j)
        y(j)=y(j+1)
        y(j+1)=sum
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software .

      SUBROUTINE realft(data,n,isign)
      INTEGER isign,n
      DOUBLE PRECISION data(n)
CU    USES four1
      INTEGER i,i1,i2,i3,i4,n2p3
      DOUBLE PRECISION c1,c2,h1i,h1r,h2i,h2r,wis,wrs
      DOUBLE PRECISION 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*sin(0.5d0*theta)**2
      wpi=sin(theta)
      wr=1.0d0+wpr
      wi=wpi
      n2p3=n+3
      do 11 i=2,n/4
        i1=2*i-1
        i2=i1+1
        i3=n2p3-i2
        i4=i3+1
        wrs=sngl(wr)
        wis=sngl(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
11    continue
      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
C  (C) Copr. 1986-92 Numerical Recipes Software .

      SUBROUTINE four1(data,nn,isign)
      INTEGER isign,nn
      DOUBLE PRECISION data(2*nn)
      INTEGER i,istep,j,m,mmax,n
      DOUBLE PRECISION tempi,tempr
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      n=2*nn
      j=1
      do 11 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=n/2
1       if ((m.ge.2).and.(j.gt.m)) then
          j=j-m
          m=m/2
        goto 1
        endif
        j=j+m
11    continue
      mmax=2
2     if (n.gt.mmax) then
        istep=2*mmax
        theta=6.28318530717959d0/(isign*mmax)
        wpr=-2.d0*sin(0.5d0*theta)**2
        wpi=sin(theta)
        wr=1.d0
        wi=0.d0
        do 13 m=1,mmax,2
          do 12 i=m,n,istep
            j=i+mmax
            tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
            tempi=sngl(wr)*data(j+1)+sngl(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
12        continue
          wtemp=wr
          wr=wr*wpr-wi*wpi+wr
          wi=wi*wpr+wtemp*wpi+wi
13      continue
        mmax=istep
      goto 2
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software .

      SUBROUTINE cosft1(y,n)
      INTEGER n
      DOUBLE PRECISION y(n+1)
CU    USES realft
      INTEGER j
      DOUBLE PRECISION sum,y1,y2
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      theta=3.141592653589793d0/n
      wr=1.0d0
      wi=0.0d0
      wpr=-2.0d0*sin(0.5d0*theta)**2
      wpi=sin(theta)
      sum=0.5d0*(y(1)-y(n+1))
      y(1)=0.5d0*(y(1)+y(n+1))
      do 11 j=1,n/2-1
        wtemp=wr
        wr=wr*wpr-wi*wpi+wr
        wi=wi*wpr+wtemp*wpi+wi
        y1=0.5d0*(y(j+1)+y(n-j+1))
        y2=(y(j+1)-y(n-j+1))
        y(j+1)=y1-wi*y2
        y(n-j+1)=y1+wi*y2
        sum=sum+wr*y2
11    continue
      call realft(y,n,+1)
      y(n+1)=y(2)
      y(2)=sum
      do 12 j=4,n,2
        sum=sum+y(j)
        y(j)=sum
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software .

