subroutine denmat(rho)
c     prepare density matrix describing the polarized beam, 
c     following G.G.Ohlsen, Rep.Prog.Phys.35('72)717, taking phi from y-axis
      implicit real*8(a-h,o-z)
      integer*4  s,rank,q
      parameter (s=1,rank=s*2)
      real*8     t(0:rank)
      complex*16 rho(-s:s,-s:s),Y(0:(rank+1)*(rank+2)),fac
      read(5,*) beta,phi,pZ,pZZ
      write(6,'(/A/)') ' beam polarization (following G.G.Ohlsen).'
      write(6,'(2(A,F6.1),2(A,F6.3))')
     &'  beta=',beta,'  phi=',phi,'  pZ=',pZ,'  pZZ=',pZZ
      t(0)=1.D0
      t(1)=pZ *sqrt(3.D0/2.D0)
      t(2)=pZZ*sqrt(1.D0/2.D0)
c
c            1                          k
c     rho = ---- sum   t   [sum [tau   D   ]] , where t    = t delta
c           2s+1    kq' kq'    q    kq  qq'            kq'    k     q'0
c
c            1                     (4*pi)  q*             pi
c         = ---- sum   t  tau  sqrt(----) Y  (-beta,-(phi+--))
c           2s+1    kq  k    kq    (2k+1)  k              2
c                                                  c.f. Ohlsen, p.745, Fig.8
c           (s)                     s-m1
c     <m1|tau  |m2> = sqrt(2s+1) (-)    < s m2 s -m1 | k q >
c            kq
c
      do m1=-s,s
         do m2=-s,s
            rho(m1,m2)=0.D0  ! initialize density matrix
         enddo
      enddo
      pi=atan(1.D0)*4.D0
      call spherical(-beta/180.D0*pi,(phi+90.D0)/180.D0*pi,rank,Y)
      do k=0,rank
         fac0=t(k)*sqrt(4.D0*pi/(2*k+1.D0))/sqrt(2*s+1.D0)
         do q=-k,k
            fac=fac0*Y(k*(k+1)+q)
            do m1=-s,s
               do m2=-s,s
                  rho(m1,m2)=rho(m1,m2)+
     &                 fac*(-1)**(s-m1)*cleb(2*s,2*m2,2*s,-2*m1,2*k,2*q)
               enddo
            enddo
         enddo
      enddo
      end

c     implicit real*8(a-h,o-z)
c     complex*16 rho(-1:1,-1:1)
c     call inirac       ! initialize faclog for cleb, rac, u9
c     call denmat(rho)
c     do i=1,-1,-1
c        write(6,'(3(2H (,F6.3,1H,F6.3,1H)))')(rho(i,j),j=1,-1,-1)
c     enddo
c     end
c     include 'racah.f'
c     include 'sphere.f'