program roav
c     Euler rotational averaging
      implicit none
      integer*4 i,j,jj,nat,nt,nf,nc,it,ifi,ic,ig,is,ni,na
      real*8 df,dt,dc,t,f,c,pi,cc(3),zz,dv,y,a(3,3),sw
      character*50 fn,s
      character*10 s10
      real*8,allocatable::r(:,:),rt(:,:),r0(:,:)
      integer*4,allocatable::z(:)

      if(iargc().ne.6)then
       write(6,*)'usage:'
       write(6,*)'roav <file.x> Ntheta Nfi Nchi natmin natmax'
       stop
      endif

      call getarg(1,fn)
      call getarg(2,s)
      read(s,*)nt
      call getarg(3,s)
      read(s,*)nf
      call getarg(4,s)
      read(s,*)nc
      call getarg(5,s)
      read(s,*)ni
      call getarg(6,s)
      read(s,*)na

      open(10,file=fn,status='old')
      read(10,*)
      read(10,*)nat
      allocate (r(3,nat),z(nat),r0(3,nat),rt(3,nat))
      do 1 i=1,nat
1     read(10,*)z(i),(r(j,i),j=1,3)
      close(10)
      write(6,*)nat,' atoms'

c     charge center
      do 3 j=1,3
      cc(j)=0.0d0
      zz=0.0d0
      do 31 i=1,nat
      if(i.ge.ni.and.i.le.na)then
       zz=zz+dble(z(i))
       cc(j)=cc(j)+z(i)*r(j,i)
      endif
31    continue
3     cc(j)=cc(j)/zz

c     charge center coordinates
      do 4 i=1,nat
      if(i.ge.ni.and.i.le.na)then
       do 41 j=1,3
41     r0(j,i)=r(j,i)-cc(j)
      endif
4     continue

      ig=0
      sw=0.0d0
      open(12,file='weights')
      open(13,file='rot.x')
      pi=4.0d0*datan(1.0d0)
      dt=pi/dble(nt)
      df=2.0d0*pi/nf
      dc=2.0d0*pi/nc
      y=180.0d0/pi
      t=0.0d0
      do 2 it=1,nt
      t=t+dt
      f=0.0d0
      do 2 ifi=1,nf
      f=f+df
      c=0.0d0
      do 2 ic=1,nc
      c=c+dc
      dv=sin(t)*dt*df*dc/8.0d0/pi**2
      sw=sw+dv
c     Euler angle transformation matrix:
      a(1,1)= cos(f)*cos(t)*cos(c)-sin(f)*sin(c)
      a(1,2)= sin(f)*cos(t)*cos(c)+cos(f)*sin(c)
      a(2,1)=-cos(f)*cos(t)*sin(c)-sin(f)*cos(c)
      a(2,2)=-sin(f)*cos(t)*sin(c)+cos(f)*cos(c)
      a(1,3)=-sin(t)*cos(c)
      a(3,1)= cos(f)*sin(t)
      a(2,3)= sin(t)*sin(c)
      a(3,2)= sin(f)*sin(t)
      a(3,3)= cos(t)
      ig=ig+1
      do 5 i=1,nat
      if(i.ge.ni.and.i.le.na)then
       do 51 j=1,3
       rt(j,i)=cc(j)
       do 51 jj=1,3
51     rt(j,i)=rt(j,i)+a(j,jj)*r0(jj,i)
      else
       rt(1,i)=r(1,i)
       rt(2,i)=r(2,i)
       rt(3,i)=r(3,i)
      endif
5     continue

      write(s10,100)ig
100   format(i10)
      do 6 is=1,len(s10)
6     if(s10(is:is).ne.' ')goto 7
7     write(13,1111)ig,it,ifi,ic,t*y,f*y,c*y,dv
c     open(11,file=s10(is:len(s10))//'.x')
c     write(11,1111)ig,it,ifi,ic,t*y,f*y,c*y,dv
      write(6,1111)ig,it,ifi,ic,t*y,f*y,c*y,dv
      write(12,1111)ig,it,ifi,ic,t*y,f*y,c*y,dv
1111  format(i10,3i4,3f10.2,1x,g15.6)
c     write(11,*)nat
      write(13,*)nat
      do 8 i=1,nat
c     write(11,1112)z(i),(rt(j,i),j=1,3)
8     write(13,1112)z(i),(rt(j,i),j=1,3)
1112  format(i3,3f12.6,' 0 0 0 0 0 0 0 0.0')
c     close(11)
2     continue
      close(12)
      close(13)
      write(6,*)ig,' geometries in rot.x'
      write(6,*)' weights in weights'
      write(6,60)sw
60    format(' weight sum:',f12.6)
      end