program ppropag
      implicit none
      integer*4 nat,ia,ix,nxmin,nxmax,nymin,nymax,nzmin,nzmax,
     1iy,iz,natb
      real*8 xx,xy,yy,xz,yz,zz,c,rotx,roty,rotz,pi,shx,shy,shz,
     1tx,ty,tz,rx,ry,rz,cm(3),p(3,3),x,y,z,xn,yn,zn,s,a(3,3),
     1o(3,3),on(3,3),xp,yp,zp,bohr
      real*8,allocatable::r(:,:)
      integer*4,allocatable::q(:)
      character*80 s80
      write(6,600)
 600  format(' Propagate polarizabilities for birefringence models',/,/,
     1       '              Input:  ONE.POL',/,
     1       '                   :  ONE.X',/,
     1       '                   :  PPROPAG.OPT',/,
     1       '             Output:  POL.LST',/,
     1       '                      MORE.X')
      pi=4.0d0*datan(1.0d0)
      bohr=0.529177d0
      open(9,file='ONE.POL')
      read(9,*)xx,xy,yy,xz,yz,zz,xp,yp,zp
      close(9)
      p(1,1)=xx
      p(1,2)=xy
      p(2,1)=xy
      p(2,2)=yy
      p(2,3)=yz
      p(1,3)=xz
      p(3,1)=xz
      p(3,2)=yz
      p(3,3)=zz
      open(9,file='ONE.X')
      read(9,*)
      read(9,*)nat
      allocate(q(nat),r(nat,3))
      do 1 ia=1,nat
  1   read(9,*)q(ia),(r(ia,ix),ix=1,3)
      close(9)
      nxmin=-0
      nxmax= 0
      nymin= 0
      nymax= 0
      nzmin=-5
      nzmax= 5
      rotx=0.0d0
      roty=0.0d0
      rotz=45.0d0
      shx=4.0d0
      shy=4.0d0
      shz=4.0d0
      open(9,file='PPROPAG.OPT')
  100 read(9,90,end=99,err=99)s80
  90  format(a80)
      if(s80(1:5).eq.'NXMIN')read(9,*)nxmin
      if(s80(1:5).eq.'NXMAX')read(9,*)nxmax
      if(s80(1:5).eq.'NYMIN')read(9,*)nymin
      if(s80(1:5).eq.'NYMAX')read(9,*)nymax
      if(s80(1:5).eq.'NZMIN')read(9,*)nzmin
      if(s80(1:5).eq.'NZMAX')read(9,*)nzmax
      if(s80(1:4).eq.'ROTX')read(9,*)rotx
      if(s80(1:4).eq.'ROTY')read(9,*)roty
      if(s80(1:4).eq.'ROTZ')read(9,*)rotz
      if(s80(1:6).eq.'SHIFTX')read(9,*)shx
      if(s80(1:6).eq.'SHIFTY')read(9,*)shy
      if(s80(1:6).eq.'SHIFTZ')read(9,*)shz
      goto 100
  99  close(9)
      rotx=rotx*pi/180.0d0
      roty=roty*pi/180.0d0
      rotz=rotz*pi/180.0d0

      cm=0.0d0
      do 3 ix=1,3
      do 31 ia=1,nat
   31 cm(ix)=cm(ix)+r(ia,ix)
    3 cm(ix)=cm(ix)/dble(nat)

      open(9,file='MORE.X')
      open(8,file='POL.LST')
      natb=nat*(nxmax-nxmin+1)*(nymax-nymin+1)*(nzmax-nzmin+1)
      write(9,900)natb
  900 format('Mropagated molecules',/,i6)
      do 2 ix=nxmin,nxmax
      rx=rotx*dble(ix)
      tx=shx*dble(ix)
      do 2 iy=nymin,nymax
      ry=roty*dble(iy)
      ty=shy*dble(iy)
      do 2 iz=nzmin,nzmax
      rz=rotz*dble(iz)
      tz=shz*dble(iz)
      do 20 ia=1,nat
      x=r(ia,1)-cm(1)
      y=r(ia,2)-cm(2)
      z=r(ia,3)-cm(3)
      s=sin(rx)
      c=cos(rx)
      a=0.0d0
      a(1,1)= 1.0d0
      a(2,2)= c
      a(3,3)= c
      a(2,3)= s
      a(3,2)=-s
      if(ia.eq.1)call rt(p,a,on)
      call rv(x,y,z,a,xn,yn,zn)

      s=sin(ry)
      c=cos(ry)
      a(2,2)= 1.0d0
      a(1,1)= c
      a(3,3)= c
      a(1,3)= s
      a(3,1)=-s
      if(ia.eq.1)call rt(on,a,o)
      call rv(xn,yn,zn,a,x,y,z)

      s=sin(rz)
      c=cos(rz)
      a(3,3)= 1.0d0
      a(1,1)= c
      a(2,2)= c
      a(1,2)= s
      a(2,1)=-s
      if(ia.eq.1)call rt(o,a,on)
      call rv(x,y,z,a,xn,yn,zn)

  20  write(9,908)q(ia),cm(1)+tx+xn,cm(2)+ty+yn,cm(3)+tz+zn
 908  format(i6,3f12.6,' 0 0 0 0 0 0 0 0.0')

      xx=on(1,1)
      xy=on(1,2)
      yy=on(2,2)
      xz=on(1,3)
      yz=on(2,3)
      zz=on(3,3)
 2    write(8,601)xx,xy,yy,xz,yz,zz,xp+tx/bohr,yp+ty/bohr,
     1zp+tz/bohr,nat
601   format(6f10.3,1x,3f10.3,i4)
      close(9)
      close(8)

      end

      subroutine rt(p,a,o)
      implicit none
      real*8 p(3,3),o(3,3),a(3,3),t(3,3)
      integer*4 i,j
      do 201 i=1,3
      do 201 j=1,3
 201  t(i,j)=a(i,1)*p(1,j)+a(i,2)*p(2,j)+a(i,3)*p(3,j)
      do 202 i=1,3
      do 202 j=1,3
 202  o(i,j)=a(j,1)*t(i,1)+a(j,2)*t(i,2)+a(j,3)*t(i,3)
      return
      end

      subroutine rv(x,y,z,a,xn,yn,zn)
      implicit none
      real*8 x,y,z,xn,yn,zn,a(3,3)
      xn=a(1,1)*x+a(1,2)*y+a(1,3)*z
      yn=a(2,1)*x+a(2,2)*y+a(2,3)*z
      zn=a(3,1)*x+a(3,2)*y+a(3,3)*z
      return
      end