program sft
      parameter (np0=2000)
      character*80 s80
      dimension xx(np0),yy(np0),si(np0),sm(np0),st(np0),sm2(np0),
     1sm3(np0),start(np0)

      ia=iargc()
      if(ia.lt.2)write(6,600)
600   format('Usage: sft limit nstart',/,
     1       '       usually limit = 2...10',/,
     2       '       nstart = 1')
      call getarg(1,s80)
      read(s80,*)alim
      call getarg(2,s80)
      read(s80,*)nstart
      ns=0
      open(9,file='spec.lst')
1     read(9,900,end=99,err=99)s80
900   format(a80)
      ns=ns+1
      goto 1
99    close(9)
      nst=ns
      write(6,*)ns,' spectra'

      open(9,file='spec.lst')
      do 10 ns=1,nst
      read(9,900)s80
      open(91,file=s80)
      np=0
2     read(91,*,end=88,err=88)x,y
      np=np+1
      if(np.gt.np0)call report('np>np0')
      xx(np)=x
      yy(np)=y
      goto 2
88    close(91)
      if(ns.eq.nstart)then
       do 11 i=1,np
11     start(i)=yy(i)
       write(6,*)ns,' remembered'
      endif
      if(ns.eq.1)then
       write(6,*)np,' points'
       np1=np
       do 3 i=1,np
       si(i)=yy(i)
       sm(i)=yy(i)
       sm2(i)=yy(i)
       sm3(i)=yy(i)
3      st(i)=0.0
      else
       if(np.ne.np1)then
        write(6,900)s80
        call report('different number of points')
       endif
       do 4 i=1,np
       der=(si(i)+sm(i)-sm2(i)-sm3(i))/4.0
       ader=abs(yy(i)-si(i))
       if(der.ne.0.)ader=ader/abs(der)
       if(ader.gt.alim)then
        yave=(si(i)+sm(i))/2.0+der*2.0
        st(i)=st(i)+(yy(i)-yave)*real(ns)/real(nst)
       endif
4      continue
       do 6 i=1,np
       sm3(i)=sm2(i)
       sm2(i)=sm(i)
       sm(i)=si(i)
6      si(i)=yy(i)
      endif
10    continue
      close(9)

      open(40,file='total.prn')
      do 7 i=1,np
7     write(40,400)xx(i),yy(i)-st(i)-start(i)*real(nstart)/real(nst)
400   format(2f18.4)
      close(40)
      stop
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end