1      subroutine ana_rdf_init()
2c
3c $Id$
4c
5      implicit none
6c
7#include "ana_common.fh"
8#include "global.fh"
9#include "mafdecls.fh"
10#include "msgids.fh"
11#include "rtdb.fh"
12c
13      if(.not.ma_push_get(mt_int,nrdf*nsel*mwa,'irdf',l_rdf,i_rdf))
14     + call md_abort('Could not allocate irdf',0)
15      print*,'rdf allocated in rdfhdr ',nrdf*nsel*mwa
16c
17      return
18      end
19      subroutine ana_rdfhdr(irdf)
20c
21      implicit none
22c
23#include "ana_common.fh"
24#include "global.fh"
25#include "mafdecls.fh"
26#include "msgids.fh"
27#include "rtdb.fh"
28c
29      integer irdf(nsel,mwa,nrdf)
30c
31      character*255 fname
32      integer i,j,k,lq,l,m
33c
34      fname=filtrj
35      lq=index(filtrj,'?')
36      if(lq.gt.0) then
37      fname=filtrj(1:lq-1)//cnum//filtrj(lq+1:index(filtrj,' ')-1)
38      endif
39      lq=index(fname,'.trj')
40      fname(lq:lq+3)='.rdf'
41c
42      open(unit=lfnrdf,file=fname(1:index(fname,' ')-1),
43     + status='unknown')
44c
45      write(*,3333) fname(1:index(fname,' ')-1)
46 3333 format(/,' Opening rdf file ',a)
47      rewind(lfnrdf)
48c
49      numrdf=0
50c
51      do 1 i=1,nrdf
52      do 2 j=1,mwa
53      do 3 k=1,nsel
54      irdf(k,j,i)=0
55    3 continue
56    2 continue
57    1 continue
58c
59      print*,'rdf set to zero in rdfhdr'
60c
61      return
62      end
63      subroutine ana_rdf(isel,xs,xw,irdf)
64c
65      implicit none
66c
67#include "ana_common.fh"
68#include "global.fh"
69#include "mafdecls.fh"
70#include "msgids.fh"
71#include "rtdb.fh"
72c
73      integer isel(msa),irdf(nsel,mwa,nrdf)
74      real*8 xs(msa,3),xw(mwm,mwa,3)
75      integer i,j,k,l,m
76      real*8 d
77c
78      i=0
79      do 1 j=1,nsa
80      if(isel(j).gt.0) then
81      i=i+1
82      do 2 k=1,nwm
83      do 3 l=1,nwa
84      d=sqrt((xs(j,1)-xw(k,l,1))**2+(xs(j,2)-xw(k,l,2))**2+
85     + (xs(j,3)-xw(k,l,3))**2)
86      m=int(dble(nrdf*rrdf)/d)
87      if(m.le.nrdf) irdf(i,l,m)=irdf(i,l,m)+1
88    3 continue
89    2 continue
90      endif
91    1 continue
92c
93      numrdf=numrdf+1
94c
95      return
96      end
97      subroutine ana_rdfwrt(isel,irdf)
98c
99      implicit none
100c
101#include "ana_common.fh"
102#include "global.fh"
103#include "mafdecls.fh"
104#include "msgids.fh"
105#include "rtdb.fh"
106c
107      integer isel(msa)
108      integer irdf(nsel,mwa,nrdf)
109c
110      integer i,j,k,l,m
111      real*8 d
112c
113      write(*,'(a)') 'ANA_RDFWRT'
114c
115      do 4 i=1,nsel
116      do 5 m=1,nrdf
117      write(*,'(2i5,10i10)') i,m,(irdf(i,l,m),l=1,nwa)
118    5 continue
119    4 continue
120c
121c      i=0
122c      do 1 j=1,nsa
123c      if(isel(j).gt.0) then
124c      i=i+1
125c      do 2 k=1,nrdf
126c      d=dble(k)*rrdf/dble(nrdf)
127c      write(*,1000) d,(irdf(i,l,k),l=1,nwa)
128c 1000 format(f12.6,10i5)
129c    2 continue
130c      endif
131c    1 continue
132c
133      close(unit=lfnrdf,status='keep')
134      write(*,'(a)') ' Closing rdf file '
135c
136      return
137      end
138