program speave
      implicit none
      character*80 fl,fn,fw
      integer*4 nl,i,nf,np,ie,iargc,j
      real*8 yy
      real*8,allocatable::x(:),ave(:),dev(:),ave2(:),adev(:),w(:)

      write(6,6000)
6000  format(/,' Average of spectra listed in a file',/)
      if(iargc().ne.1.and.iargc().ne.2)then
       write(6,*)'Usage: speave <file.lst> <<weigthsfile>>'
       stop
      endif 

      call getarg(1,fl)
c     how many lines in fl:
      nf=nl(fl)
      write(6,*)nf,' files'
      allocate(w(nf))
      if(iargc().eq.2)then
       call getarg(2,fw)
       open(9,file=fw,status='old')
       do 9 i=1,nf
9      read(9,*)w(i)
       close(9)
       write(6,*)'weights read from '//fw
      else
       do 10 i=1,nf
10     w(i)=1.0d0/dble(nf)
       write(6,*)'equal weights'
      endif

c     read first line and allocate:
      open(9,file=fl)
      read(9,80)fn
80    format(a80)
      close(9)
      np=nl(fn)
      write(6,*)np,' points'
      allocate(x(np),ave(np),dev(np),ave2(np),adev(np))
      call vz(x,np)
c     average <S>:
      call vz(ave,np)
c     average <S^2>:
      call vz(ave2,np)
c     quadratic deviation 0.5 x sqrt(<S^2>-<S>^2):
      call vz(dev,np)
c     absolute deviation 0.5 x <|S-<S>|>:
      call vz(adev,np)

      open(9,file=fl)
      do 1 i=1,nf
      read(9,80)fn
      open(10,file=fn)
      do 2 j=1,np
      read(10,*)x(j),yy
      ave( j)=ave( j)+yy*w(i)
2     ave2(j)=ave2(j)+yy**2*w(i)
1     close(10)
      close(9)

      do 3 j=1,np
3     dev( j)=dsqrt(ave2(j)-ave(j)**2)/2.0d0

c     second round for absolute deviation:
      open(9,file=fl)
      do 6 i=1,nf
      read(9,80)fn
      open(10,file=fn)
      do 7 j=1,np
      read(10,*)yy,yy
7     adev(j)=adev(j)+dabs(yy-ave(j))*w(i)
6     close(10)
      close(9)

      do 8 j=1,np
8     adev(j)=adev(j)/2.0d0

      do 4 i=1,len(fl)
      if(fl(i:i).eq.'.')then
       ie=i-1
       goto 5
      endif
4     continue
      ie=len(fl)
5     continue

      call ws(fl(1:ie)//'.ave.prn',np,x,ave)
      call ws(fl(1:ie)//'.dev.prn',np,x,dev)
      call ws(fl(1:ie)//'.adev.prn',np,x,adev)
      end

      subroutine ws(f,n,x,s)
      implicit none
      integer*4 n,i
      real*8 x(*),s(*)
      character*(*) f
      open(9,file=f)
      do 1 i=1,n
 1    write(9,900)x(i),s(i)
 900  format(f12.2,e14.4)
      close(9)
      write(6,*)f//' written'
      return
      end

      function nl(f)
      implicit none
      character*(*) f
      integer*4 nl,n
      open(9,file=f)
      n=0
1     read(9,*,end=99,err=99)
      n=n+1
      goto 1
99    close(9)
      nl=n
      return
      end

      subroutine vz(x,n)
      implicit none
      real*8 x(*)
      integer*4 n,i
      do 1 i=1,n
  1   x(i)=0.0d0
      return
      end