program ftsmooth
      implicit none
      integer*4 i,n,np1,iw,ni,na,iargc
      real*8 xmin,xmax,dx,pi,fik,L,dw
      real*8,allocatable::sx(:),sy(:),sw(:),sxw(:)
      character*80 fn,fnn,s
      pi=4.0d0*atan(1.0d0)
      if(iargc().ne.4)then
       write(6,600)
600    format('Usage: ftsmooth <file.prn> Nmin Nmax <newfile.prn>')
       stop
      endif
      call getarg(1,fn)
      call getarg(2,s)
      read(s,*)ni
      call getarg(3,s)
      read(s,*)na
      call getarg(4,fnn)
c     number of points:
      n=np1(fn,-1.d9,1.d9)      
      write(6,*)n,' points'
      allocate(sx(n),sy(n),sw(2*n),sxw(2*n))
c     read  spectrum:
      call ls(sx,sy,fn,-1.d9,1.d9)      
      xmin=sx(1)
      xmax=sx(1)
      do 1 i=1,n
      if(sx(i).lt.xmin)xmin=sx(i)
1     if(sx(i).gt.xmax)xmax=sx(i)
      L=(xmax-xmin)
c     slow FT:
      sw=0.0d0
      do 2 i=1,n
      if(i.eq.1.or.i.eq.n)then
       if(i.eq.1)dx=dabs(sx(i+1)-sx(1))/2.0d0
       if(i.eq.n)dx=dabs(sx(n)-sx(n-1))/2.0d0
      else
       dx=dabs(sx(i+1)-sx(i-1))/2.0d0
      endif
      do 2 iw=1,2*n
      sxw(iw)=dble(iw)/2.0d0
2     sw(iw)=sw(iw)+sy(i)*dx*fik(iw,sx(i),L)
      call wrs('tem-w.prn',2*n,sxw,sw)
c     slow back FT:
      sy=0.0d0
      do 3 iw=1,2*n
      if(iw.ge.ni.and.iw.le.na)then
       dw=1.0d0
       do 21 i=1,n
21     sy(i)=sy(i)+sw(iw)*dw*fik(iw,sx(i),L)
      endif
3     continue
      call wrs(fnn,n,sx,sy)
      end

      subroutine wrs(s,np,x,y)
      implicit none
      character*(*) s
      integer*4 i,np
      real*8 x(*),y(*)
      open(67,file=s)
      do 1 i=1,np
1     write(67,900)x(i),y(i)
900   format(f9.2,g13.5)
      close(67)
      return
      end

      function fik(k,x,L)
      implicit none
      integer*4 k
      real*8 fik,x,L,pi,sf
      pi=4.0d0*atan(1.0d0)
      if(k.eq.1)then
       sf=1.0d0/dsqrt(L)
      else
       sf=dsqrt(2.0d0/L)
      endif
      if(mod(k,2).ne.0)then
       fik=sf*cos(dble(2*((k-1)/2))*pi*x/L)
      else
       fik=sf*sin(dble(2*(k/2))*pi*x/L)
      endif
      return
      end


      function np1(s,wmin,wmax)
      implicit none
      integer*4 np1,n
      real*8 wmin,wmax,x
      character*(*) s
      open(9,file=s)
      n=0
101   read(9,*,end=99,err=99)x
      if(x.ge.wmin.and.x.le.wmax)n=n+1
      goto 101
99    close(9)
      np1=n
      return
      end

      subroutine ls(xs,ys,name,wmin,wmax)
c     loads in the spectrum  within wmin wmax    
      implicit none
      integer*4 i
      real*8 xs(*),ys(*),wmin,wmax,x,y
      character*(*) name
      i=0
      open(9,file=name)
 11   read(9,*,end=99,err=99)x,y
      if(wmin.le.x.and.x.le.wmax)then
       i=i+1       
       xs(i)=x
       ys(i)=y
      endif
      goto 11
  99  close(9)
      return
      end