program intrdf
      parameter (nat0=100,np0=200)
      real noise,noiset,ni
      dimension ai(nat0),yi(4,nat0,np0),x(np0),noise(nat0,np0),
     1signal(nat0,np0),noiset(np0),signalt(np0)
      character*80 s80,fn
      character*6 s6
      logical lex
c     number of points in RDF distributions:
      np=100
c     number of sets to average:
      ns=4
c     number of atoms in each set
      nat=19
c     integration limits (A):
      d1=3.0
      d2=20.0
      inquire(file='INTRDF.OPT',exist=lex)
      if(lex)then
       open(9,file='INTRDF.OPT')
1001   read(9,80,end=99,err=99)s6
       if(s6(1:2).eq.'NP')read(9,*)np
       if(s6(1:2).eq.'NS')read(9,*)ns
       if(s6(1:3).eq.'NAT')read(9,*)nat
       if(s6(1:2).eq.'d1')read(9,*)d1
       if(s6(1:2).eq.'d2')read(9,*)d2
       goto 1001
99     close(9)
      endif
      write(6,*)np,ns,nat,d1,d2
      if(np.gt.np0)then
       write(6,*)'too many points'
       stop
      endif
      if(nat.gt.nat0)then
       write(6,*)'too many atoms'
       stop
      endif
      do 3 ia=1,nat
3     ai(ia)=0.

      open(20,file='NAMES.LST')
      do 1 is=1,4
c     is=1:"RR:"
c        2 "RS:"
c        3 "SS:"
c        4 "SR:"
      do 111 ia=1,nat
      read(20,200)fn
      write(6,200)fn
200   format(a80)
      open(21,file=fn)
      do 1111 ip=1,np
1111  read(21,*)x(ip),yi(is,ia,ip)
111   close(21)
1     continue
      close(20)


      write(6,*)'averaging Fs and aromate'
      do 20 is=1,4
      do 20 ip=1,np
      t=(yi(is,9,ip)+yi(is,10,ip)+yi(is,11,ip))/3.
      yi(is,9,ip)=t
      yi(is,10,ip)=t
      yi(is,11,ip)=t
      t=(yi(is,3,ip)+yi(is,14,ip))/2.
      yi(is,3,ip)=t
      yi(is,14,ip)=t
      t=(yi(is,4,ip)+yi(is,15,ip))/2.
      yi(is,4,ip)=t
      yi(is,15,ip)=t
      t=(yi(is,1,ip)+yi(is,16,ip))/2.
      yi(is,1,ip)=t
      yi(is,16,ip)=t
      t=(yi(is,2,ip)+yi(is,17,ip))/2.
      yi(is,2,ip)=t
20    yi(is,17,ip)=t

      do 23 ia=1,nat
      write(s6,80)ia
      do 24 istart=1,len(s6)
24    if(s6(istart:istart).ne.' ')goto 25
25    open(8,file=s6(istart:6)//'rr.prn')
      do 21 ip=1,np
21    write(8,81)x(ip),(yi(1,ia,ip)+yi(3,ia,ip))/2.
      close(8)
      open(8,file=s6(istart:6)//'rs.prn')
      do 22 ip=1,np
22    write(8,81)x(ip),(yi(2,ia,ip)+yi(4,ia,ip))/2.
23    close(8)

      do 12 ip=1,np
      signalt(ip)=0.0
12    noiset(ip)=0.0
      do 18 ia=1,nat
      do 18 ip=1,np
c     SS+RR-RS-SR:
      signal(ia,ip)=yi(1,ia,ip)+yi(3,ia,ip)-yi(2,ia,ip)-yi(4,ia,ip)
c     (|SS-RR|+|RS-SR|)/2:
      noise(ia,ip)=(abs(yi(1,ia,ip)-yi(3,ia,ip))+
     1             abs(yi(2,ia,ip)-yi(4,ia,ip)))/2.0
      noiset(ip)=noiset(ip)+noise(ia,ip)
      signalt(ip)=signalt(ip)+abs(signal(ia,ip))
18    if(x(ip).ge.d1.and.x(ip).le.d2)ai(ia)=ai(ia)+signal(ia,ip)**2

      si=0.0
      ni=0.0
      do 13 ip=1,np
      if(x(ip).ge.d1.and.x(ip).le.d2)si=si+signalt(ip)
13    if(x(ip).ge.d1.and.x(ip).le.d2)ni=ni+noiset(ip)
      si=si*abs(x(1)-x(np))/real(np)
      ni=ni*abs(x(1)-x(np))/real(np)
      write(6,609)si,ni
609   format('Signal total integral: ',f12.6,/,
     1       'Noise total integral : ',f12.6)

      open(7,file='signalt.prn')
      open(8,file='noiset.prn')
      do 14 ip=1,np
      write(8,81)x(ip),noiset(ip)
14    write(7,81)x(ip),signalt(ip)
      close(7)
      close(8)

      do 19 ia=1,nat
      write(s6,80)ia
 80   format(i6)
      do 15 istart=1,len(s6)
15    if(s6(istart:istart).ne.' ')goto 16
16    open(8,file=s6(istart:6)//'s.prn')
      open(9,file=s6(istart:6)//'n.prn')
      do 19 ip=1,np
      write(8,81)x(ip),signal(ia,ip)
19    write(9,81)x(ip),noise(ia,ip)
81    format(f12.4,f12.6)
      close(8)
      close(9)

      am=ai(1)
      do 4 ia=1,nat
4     if(ai(ia).gt.am)am=ai(ia)
      write(6,*)'Relative Integrals (%):'
      do 2 ia=1,nat
2     write(6,600)ia,ai(ia)/am*100.
600   format(i4,f10.2)
      stop
      end