program specave2
c     average of N spectra from SP.LST
c     can also calculate integral of average spectrum in a given interval:
c     also normalized spectral error er = int |s-save| / int save|
c     s ... n point average
c     save ... all point average
      implicit none
      integer ns,np,i,j,iargc
      real*8,allocatable::x(:),sav(:),sr(:),sc(:),dd(:),ee(:)
      real*8 d,pi,spi,p,ef,dx
      character*80 s80
      character*8  s8 
      pi=4.0d0*datan(1.0d0)
      spi=dsqrt(pi)

      write(6,*)' Spectral average with convolution'
      if(iargc().lt.1)then
       write(6,*)' usage: specave2 d'
       stop
      endif

      call getarg(1,s80)
      read(s80,*)d

c     number of points from the first spectrum
      call snp(np)

      allocate(x(np),sav(np),sr(np),sc(np),ee(np+1),dd(np))
      sav=0.0d0

c     read first spectrum's scale
      call xnp(np,x)
      dx=dabs(x(1)-x(2))

c     precalculate exp
      ee=0.0d0
      do 22 j=1,np
      p=((x(1)-x(j))/d)**2
      ef=0.0d0
      if(p.lt.20.0d0)ef=exp(-p)
22    ee(j)=ef/spi/d*dx

      ns=0
      open(9,file='SP.LST')
992   read(9,900,end=918,err=918)s80
900   format(a80)
      ns=ns+1
c     read raw spectrum:
      open(10,file=s80)
      do 1 i=1,np
1     read(10,*)x(i),sr(i)
      close(10)
c     make convoluted spectrum:
c     sc(i) = int sr(j) exp -(xi-xj)^2/d^2
      sc=0.0d0
      do 11 i=1,np
      do 11 j=1,np
11    sc(i)=sc(i)+sr(j)*ee(abs(i-j)+1)
c     difference of raw and convoluted
      do 13 i=1,np
13    dd(i)=sr(i)-sc(i)
c     record difference of raw and convoluted
      write(s8,800)ns
800   format(i8)
      do 5 i=1,len(s8)
5     if(s8(i:i).ne.' ')goto 6
6     call ws(np,x,dd,'_'//s8(i:len(s8))//'.prn')
c     add difference of raw and convoluted
      do 12 i=1,np
12    sav(i)=sav(i)+dd(i)

      goto 992
918   close(9)

      do 2 i=1,np
2     sav(i)=sav(i)/dble(ns)

      call ws(np,x,sav,'AVE.PRN')

      write(6,*)ns,' spectra, AVE.PRN'

      end

      subroutine xnp(np,x)
      real*8 x(*)
      integer*4 np,i
      character*80 s80
      open(9,file='SP.LST')
      read(9,900)s80
900   format(a80)
      close(9)
      open(10,file=s80)
      do 1 i=1,np
1     read(10,*)x(i)
      close(10)
      return
      end

      subroutine snp(np)
      real*8 x,y
      integer*4 np
      character*80 s80
      open(9,file='SP.LST')
      read(9,900)s80
900   format(a80)
      close(9)
      open(10,file=s80)
      np=0
993   read(10,*,end=919,err=919)x,y
      np=np+1
      goto 993
919   close(10)
      return
      end

      subroutine ws(np,x,y,s)
      integer*4 np,i
      real*8 x(*),y(*)
      character*(*) s
      open(13,file=s)
      do 1 i=1,np
 1    write(13,1313)x(i),y(i)
 1313 format(f12.2,e12.4)
      close(13)
      return
      end