program specave
c     average of N spectra from SP.LST
c     can also calculate integral of average spectrum in a given interval:
c     ok: 1    int  s dx
c         2    int |s| dx
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,n,i,ii,is,ie,imin,iargc
      real*8,allocatable::sx(:),sa(:),sav(:),p(:),e(:),s2(:)
      real*8 sy,x,y,ai,a,b,dx,aiv,er,aia,kT,emin,q,rms,eoa
      logical boltzmann
      character*80 s80,fn
      character*1 ok
      write(6,601)
601   format(' usage: specave N c a b',/,
     1' N: number of spectra, if missed acerage all in SP.LST',/,
     1'    if N<0 use Boltzmann factors for energies in E.LST',/,
     1' if c<>n: integrate from a to b, if c=1 integrate S',/,
     1'                                 is c=2 integrate |S|',/)

      call readopt(boltzmann,kT,fn)
      if(iargc().gt.0)then
       call getarg(1,s80)
       read(s80,*)n
      else
       n=0
       open(9,file=fn)
2002   read(9,*,end=2003,err=2003)
       n=n+1
       goto 2002
2003   close(9)
       write(6,*)n,' spectra'
      endif



      boltzmann=n.lt.0
      n=abs(n)
      ok='n'
      if(iargc().gt.1)then
       call getarg(2,ok)
       call getarg(3,s80)
       read(s80,*)a
       call getarg(4,s80)
       read(s80,*)b
      endif

      allocate(p(n),e(n))
      if(boltzmann)then
       open(9,file='E.LST',status='old')
       do 9 i=1,n
  9    read(9,*)e(i)
       close(9)
       emin=e(1)
       imin=1
       do 10 i=1,n
       if(e(i).lt.emin)then
        imin=i
        emin=e(i)
       endif
10     continue
       write(6,600)n,imin,emin
600    format(' Boltzmann weighting requested',/,
     1 i6,' energies from E.LST',/,
     2 ' minimum #',i6,':',e12.4,' kcal/mol')
       q=0.0d0
       do 11 i=1,n
11     q=q+exp(-(e(i)-emin)/kT)
       do 12 i=1,n
12     p(i)=exp(-(e(i)-emin)/kT)/q
       do 13 i=1,n
       write(6,602)i,p(i)
13     if(mod(i,5).eq.0)write(6,*)
602    format(i6,f7.4,$)
       write(6,*)
      else
       p=1.0d0
      endif
      
      ns=0
      ii=0
      open(9,file=fn)
992   read(9,900,end=918,err=918)s80
900   format(a80)
      ns=ns+1
      if(ns.eq.1)then
       open(10,file=s80)
       np=0
993    read(10,*,end=919,err=919)x,y
       np=np+1
       goto 993
919    close(10)
       allocate(sx(np),sa(np),sav(np),s2(np))
       sa=0.0d0
       sav=0.0d0
       s2=0.0d0
      endif

      open(10,file=s80)
      do 1 i=1,np
      read(10,*)sx(i),sy
      if(ns.le.n)then
       sa(i)=sa(i)+sy*p(ns)
       s2(i)=s2(i)+sy**2*p(ns)
      endif
1     sav(i)=sav(i)+sy*p(ns)
      close(10)
      if(ns.le.n)ii=ii+1

      goto 992
918   close(9)

      if(.not.boltzmann)then
       do 2 i=1,np
       sa(i)=sa(i)/dble(ii)
       s2(i)=s2(i)/dble(ii)
2      sav(i)=sav(i)/dble(ns)
      endif

      write(s80,'(i80)')ii
      do 5 is=1,len(s80)
5     if(s80(is:is).ne.' ')goto 6
6     do 7 ie=len(s80),1,-1
7     if(s80(ie:ie).ne.' ')goto 8
8     open(10,file='AVE'//s80(is:ie)//'.PRN')
      if(ii.eq.n)then
c      to show rms deviation:
       open(11,file='rms1.prn')
       open(12,file='rms2.prn')
c      to show error of the average:
       open(13,file='ae1.prn')
       open(14,file='ae2.prn')
      endif
      ai =0.0d0
      aia=0.0d0
      aiv=0.0d0
      do 3 i=1,np
      if(ok.ne.'n')then
       if(sx(i).ge.a.and.sx(i).le.b)then
        if(i.gt.1)then
         dx=dabs(sx(i)-sx(i-1))
        else
         dx=dabs(sx(i+1)-sx(i))
        endif
        aia=aia+dabs(sa (i)-sav(i))*dx
        aiv=aiv+dabs(sav(i)       )*dx
        if(ok.eq.'1')then
         ai=ai+sa(i)*dx
        else
         if(ok.eq.'2')then
          ai=ai+dabs(sa(i))*dx
         else
          write(6,*)'unknown option',ok
          stop
         endif
        endif
       endif
      endif
      if(ii.eq.n)then
c      RMS deviation:
       rms=dsqrt(s2(i)-sa(i)**2)
c      error of average:
       eoa=rms/dsqrt(dble(n))
       write(11,100)sx(i),sa(i)-rms
       write(12,100)sx(i),sa(i)+rms
       write(13,100)sx(i),sa(i)-eoa
       write(14,100)sx(i),sa(i)+eoa
      endif
3     write(10,100)sx(i),sa(i)
100   format(f10.4,' ',g12.6)
      if(ii.eq.n)then
       close(11)
       close(12)
       close(13)
       close(14)
      endif
      er=0.0d0
      if(aiv.ne.0.0d0)er=aia/aiv*100.0d0
      close(10)

      write(6,6000)ii,ai,er,ns,np
6000  format(i5,' spectra averaged in AVExxx.PRN ',g12.5,f6.1,2i6)

      end

      subroutine readopt(boltzmann,kT,fn)
      implicit none
      logical boltzmann,lex
      real*8 kelvin,kT
      character*80 fn
      character*4 key
      kelvin=300.0d0
      fn='SP.LST'
      inquire(file='SPECAVE.OPT',exist=lex)
      if(lex)then
       open(9,file='SPECAVE.OPT')
 1     read(9,90,end=99,err=99)key
 90    format(a4)
       if(key.eq.'BOLT')read(9,*)boltzmann
       if(key.eq.'KELV')read(9,*)kelvin
       if(key.eq.'FILE')read(9,*)fn
       goto 1
 99    close(9)
      endif
c     kt in kcal/mol:
      write(6,600)kelvin
600   format(' T = ',f10.1,' K')
      kT=kelvin/503.228d0
      return
      end