      subroutine decay
c     particle-decay calculation [c.f. G.Finkel et al., PRC 19('79)1782.]
c     with a polarized beam of spin-1 particle
      implicit none
      integer*4  i,ithc,nthc,ndmax,ldmax,nthcmx,ndecay
      integer*4  m,m0,ijtr,ijtr0,njtr,j,m1,m2,md,ldx
      real*8     phic,dthc,pi,thc,fac,cleb,w,thcmin,thcmax
c     implicit real*8(a-h,o-z)
      include 'parameter.f'
      integer*4  sa
      parameter  (sa=1)
      complex*16 S0(-sa:sa),S(-sa:sa,-jtrmax:jtrmax,njtrmx)
      complex*16 rho(-sa:sa,-sa:sa)
      parameter  (ndmax=5,ldmax=jtrmax+1,nthcmx=100)
      real*8     IA,IB,IC,sc,jd(ndmax),Ad(ndmax),mud,MC,msc,MA,MB
      real*8     angcor(0:nthcmx)
      integer*4  ld(ndmax),jtr(*),err
      complex*16 Y(0:(ldmax+1)*(ldmax+2)),T1q,T2,T3m1,T3m2,T3q
      save       IC,sc,ndecay,ld,jd,Ad,thcmin,dthc,nthc,pi,rho,S,angcor
      real*8     thr,d2
      save       thr
      complex*16 Y0(0:(ldmax+1)*(ldmax+2))
c
c     readin parameters
      entry decay0(IB,*,*)
      err=0
      read(5,*,end=90000) IC,sc
      read(5,*) ndecay,(ld(i),jd(i),Ad(i),i=1,ndecay)
      if (ndecay.gt.ndmax) call decay9(err,'too large ndecay')
      if (ndecay.eq.1) Ad(1)=1.D0
      ldx=0
      do i=1,ndecay
         if (jd(i).lt.abs(ld(i)-sc).or.jd(i).gt.ld(i)+sc)
     &        call decay9(err,'inconsistent [ld,sc]jd')
         if (IB.lt.abs(IC-jd(i)).or.IB.gt.IC+jd(i))
     &        call decay9(err,'inconsistent [IC,jd]IB')
         if (ld(i).gt.ldx) ldx=ld(i)
      enddo
      read(5,*) phic,thcmin,thcmax,dthc,thr
      nthc=nint((thcmax-thcmin)/dthc)
      if (nthc.gt.nthcmx) call decay9(err,'too many angles for thc')
      write(6,'(/A/)') ' particle decay parameters.'
      write(6,'(2(A,F4.1))')
     &'  IC=',IC,'  sc=',sc
      write(6,'(A,10(A,I2,A,F4.1,A))')
     &'  (ld,jd)=',(' (',ld(i),',',jd(i),')',i=1,ndecay)
      write(6,'(5(2X,A,F5.1),A)') 'phic=',phic,
     &'thcmin=',thcmin,'thcmax=',thcmax,'dthc=',dthc,'thr=',thr,' [deg]'
      pi=atan(1.D0)*4.D0
      call denmat(rho)
      if (err.ne.0) return 1
      return
c
c     clear storage
c
      entry decay1(njtr,jtr)
      do ithc=0,nthc
         angcor(ithc)=0.D0
      enddo
c     clear Samp.
      do ijtr=1,njtr
         j=jtr(ijtr)
         do m=-j,j
            do m1=-sa,sa
               S(m1,m,ijtr)=(0.D0,0.D0)
            enddo
         enddo
      enddo
      return
c
c     store Samp.(ma,m,j)
c
      entry decay2(S0,m0,ijtr0)
      do m1=-sa,sa
         S(m1,m0,ijtr0)=S0(m1)
      enddo
      return
c
c     calculate angular correlation by T3q/T1q, where 
c     T3q = sum <MC,mc|T3|MA,ma>rho(ma,ma')<MC,mc|T3|MA,ma'>*
c           over MC,mc,MA,ma,ma' incoherently, 
c     <MC,mc|T3|MA,ma> = sum <MC,mc|T2|MB><MB|T1|MA,ma>
c                        over MB coherently, 
c     <MB|T1|MA,ma>    = sum sqrt(2j+1)CG(MA,m|MB)S(m,ma,j)
c                        over j coherently (m=MB-MA, thus no summation)
c     <MC,mc|T2|MB>    = sum Ad*CG(md,mc|mud)CG(mud,MC|MB)Y(ld,md)
c                        over d (mud=MB-MC, md=mud-mc, no summation)
c     T1q = sum <MB|T1|MA,ma>rho(ma,ma')<MB|T1|MA,ma'>*
c           over MA,MB,ma,ma'
c         = (2IB+1) sum S(m,ma,j)rho(ma,ma')S(m,ma',j)*
c                   over j,m,ma,ma'
c
      entry decay3(njtr,jtr,IA,IB,w)
      T1q=(0.D0,0.D0)
      fac=2*IB+1.D0
      do ijtr=1,njtr
         j=jtr(ijtr)
         do m=-j,j
            do m1=-sa,sa
               do m2=-sa,sa
                  if (rho(m1,m2).ne.0.D0) ! for saving time...
     &            T1q=T1q+fac*S(m1,m,ijtr)*rho(m1,m2)*
     &                  conjg(S(m2,m,ijtr))
               enddo ! ma'
            enddo ! ma
         enddo ! m
      enddo ! j
      do ithc=0,nthc
         T3q=(0.D0,0.D0)
         thc=(ithc*dthc+thcmin)/180.D0*pi
         call spherical(thc,phic/180.D0*pi,ldx,Y0)
         do j=0,ldx
            do m1=-j,j
               Y(j*(j+1)+m1)=0.D0
               do m2=-j,j
                  Y(j*(j+1)+m1)=Y(j*(j+1)+m1)
     &                         +d2(j,m2,m1,thr/180.D0*pi)*Y0(j*(j+1)+m2)
               enddo
            enddo
         enddo
         do MC=-IC,IC,1.D0
            do msc=-sc,sc,1.D0
               do MA=-IA,IA,1.D0
                  do m1=-sa,sa ! ma
                     do m2=-sa,sa ! ma'
                        if (rho(m1,m2).ne.0.D0) then ! for saving time...
                        T3m1=(0.D0,0.D0)
                        T3m2=(0.D0,0.D0)
                        do MB=-IB,IB,1.D0
                           T2  =(0.D0,0.D0)
                           mud=MB-MC
                           md =nint(mud-msc)
                           do i=1,ndecay
                              if (ld(i).ge.abs(md ).and.
     &                            jd(i).ge.abs(mud)-1.D-3)
     &                            T2=T2+cleb(nint(2*jd(i)),nint(2*mud),
     &                                       nint(2*IC   ),nint(2*MC ),
     &                                       nint(2*IB   ),nint(2*MB ))*
     &                                  cleb(     2*ld(i) ,     2*md  ,
     &                                       nint(2*sc)   ,nint(2*msc),
     &                                       nint(2*jd(i)),nint(2*mud))*
     &                                  Y(ld(i)*(ld(i)+1)+md)*Ad(i)
                           enddo ! ndecay
                           m=nint(MB-MA)
                           do ijtr=1,njtr
                              j=jtr(ijtr)
                              if (j.ge.abs(m)) then
                                 fac=sqrt(2*j+1.D0)*
     &                                cleb(nint(2*IA),nint(2*MA),
     &                                          2*j  ,     2*m,
     &                                     nint(2*IB),nint(2*MB))
                                 T3m1=T3m1+fac*S(m1,m,ijtr)*T2
                                 T3m2=T3m2+fac*S(m2,m,ijtr)*T2
                              endif
                           enddo ! j
                        enddo ! MB
                        T3q=T3q+T3m1*rho(m1,m2)*conjg(T3m2)
                        endif ! if (rho(m1,m2).ne.0.D0)
                     enddo ! ma'
                  enddo ! ma
               enddo ! MA
            enddo ! mc
         enddo ! MC
         angcor(ithc)=angcor(ithc)+abs(T3q/T1q)*w
      enddo ! ithc
      return
c
c     print angular correlation
      entry decay4(w)
      do ithc=0,nthc
         if (w.ne.0.D0) angcor(ithc)=angcor(ithc)/w !non-pntprj, 1/(E12max-min)
         write(6,'(F10.1,1pE12.3)') ithc*dthc,angcor(ithc)
      enddo
      return
90000 return 2
      end

      subroutine decay9(i,s)
      character s*(*)
      i=i+1
      write(6,'(/2A/)') ' !!! ',s
      return
      end

      real*8 function d2(l,m10,m20,th)
      implicit real*8(a-h,o-z)
      real*8    faclog(0:499)
      common/clebma/faclog
      m1=min(m10,m20)
      m2=max(m10,m20)
      d2=0.D0
      do i=0,min(l+m1,l-m2)
         d2=d2
     &     +(-1)**i*cos(th/2.)**(2*l+m1-m2-2*i)*sin(-th/2.)**(2*i-m1+m2)
     &     /exp(faclog(l+m1-i)+faclog(l-m2-i)+faclog(i)+faclog(i-m1+m2))
      enddo
      d2=d2
     &  *exp((faclog(l+m1)+faclog(l-m1)+faclog(l+m2)+faclog(l-m2))/2.)
      if (m20.gt.m10) d2=d2*(-1)**(m20-m10)
      return
      end
