program psmooth
      implicit none
      integer*4 i,n,np1,iargc,it,maxit
      real*8 p,smin,smax,a,one,two,ym,y0,yp
      real*8,allocatable::sx(:),sy(:),sp(:),sn(:)
      character*80 fn,fnn,s
      one=1.0d0
      two=2.0d0
      if(iargc().ne.4)then
       write(6,600)
600    format('Usage: psmooth <file.prn> maxit p <newfile.prn>')
       stop
      endif
      call getarg(1,fn)
      call getarg(2,s)
      read(s,*)maxit
      call getarg(3,s)
      read(s,*)p
      call getarg(4,fnn)
c     number of points:
      n=np1(fn,-1.d9,1.d9)      
      write(6,*)n,' points'
      allocate(sx(n),sy(n),sp(n),sn(n))
c     read  spectrum:
      call ls(sx,sy,fn,-1.d9,1.d9)      
      smin=sy(1)
      smax=sy(1)
      do 1 i=1,n
      if(sy(i).lt.smin)smin=sy(i)
1     if(sy(i).gt.smax)smax=sy(i)
      a=p*(smax-smin)
      sp=sy
      it=0

101   it=it+1
      if(it.gt.maxit)goto 999
c     i=1
      sn(1)=min(0.9d0*sp(1)+0.1d0*sp(2),sy(1))
c     i=n
      sn(n)=min(0.9d0*sp(n)+0.1d0*sp(n-1),sy(n))
      do 2 i=2,n-1
c     sn(i)=min((sp(i)+(sp(i-1)+sp(i+1))/two + a)/two,sy(i))
2     sn(i)=min((sp(i-1)+sp(i+1))/two + a,sy(i))
      sp=sn
      goto 101

999   call wrs('tem-b.prn',n,sx,sp)

      do 3 i=1,n
3     sp(i)=sy(i)-sp(i)
      call wrs(fnn,n,sx,sp)
      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 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