matrx.f

Go to the documentation of this file.

      subroutine matrx
c                                                                       
      implicit real*8 (a-h,o-z)
      include 'afrt_rt2.cmn'
c***********************************************************************
      amusq=cosmu(ii)**2
      amumu=cosmu(ii)*cosmu(kk)
      amups=cosmu(kk)**2
      do l=1,jpart
          copsi=sinmu(ii)*sinmu(kk)+cosmu(ii)*cosmu(kk)*costh(l)
          copsq=copsi**2
          cfisq=costh(l)**2
          copcs=copsi*costh(l)
          sfisq=sinth(l)**2
          if((cosmu(kk)-cosmu(ii)).eq.0) then
              if((l-1).eq.0) then
                  if((ii-nx).lt.0) then
                      iq=4*(ii-1)
                  else
                     iq=(nmum1-ii)*4
                  endif
                  mn=0
                  do ik=1,16
                      jp=((ik-1)/4)*4
                      ip=iq+jp/4+1
                      ir=ik-jp
                      p(ik,l)=pp(ip,ir)
                  enddo
                  call depol(l)
                  if((ii-nx).lt.0) then
                      if((l-1).eq.0) then
                          call mats(jpart+1)
                      endif
                  endif
              else
                 call mats(l)
              endif
          else if((cosmu(ii)+cosmu(kk)).eq.0) then
              if((l-jpart).eq.0) then
                  if((ii-nx).lt.0) then
                      iq=4*(ii-1)
                  else
                      iq=(nmum1-ii)*4
                  endif
                  mn=jpart*32-31
                  do ik=mn,mn+15
                      jp=((ik-mn)/4)*4
                      ip=iq+jp/4+1
                      ir=ik-jp-mn+1
                      p(ik-mn+1,l)=qq(ip,ir)
                  enddo
                  mn=mn-1
                  call depol(l)
                  if((ii-nx).lt.0) then
                      if((l-1).eq.0) then
                          call mats(jpart+1)
                      endif
                  endif
              else
                  call mats(l)
              endif
          else
              call mats(l)
          endif
      enddo
      return                                                            
      end                                                               
c***********************************************************************
c***********************************************************************
 
      subroutine mats(l)                                                
c***********************************************************************
      implicit real*8 (a-h,o-z)
      include 'afrt_rt2.cmn'
c***********************************************************************
      x = cosmu(ii)*cosmu(kk)+sinmu(ii)*sinmu(kk)*costh(l)
      if(x.gt.1.0d0)x=1.0d0
      if(x.lt.-1.0d0)x=-1.0d0
      if(dabs(x) .lt. 1.0e-6)x=0.0
      tf=10.0d0*(dacos(x)/conv)+1.0
      tfm1=idint(tf+0.01d0)
      tfp1=tfm1+1.0d0
      if((l-jpart)>0) then
          mt=1
          mn=jpart*32
      else
          mt=tfm1
          mn=(l-1)*32
      endif   
      call xntpln(tf,tfm1,tfp1,t(mt,1),t(mt+1,1),tmt1)
      call xntpln(tf,tfm1,tfp1,t(mt,2),t(mt+1,2),tmt2)
      call xntpln(tf,tfm1,tfp1,t(mt,3),t(mt+1,3),tmt3)
      call xntpln(tf,tfm1,tfp1,t(mt,4),t(mt+1,4),tmt4)
      p(1,l)=cfisq*tmt1+copsq*tmt2+2.0d0*copcs*tmt3
      p(2,l)=sfisq*(amups*tmt1+amusq*tmt2+2.0d0*amumu*tmt3)
      p(3,l)=-sinth(l)*(cosmu(kk)*costh(l)*tmt1+cosmu(ii)*copsi*tmt2+
     1       (cosmu(kk)*copsi+cosmu(ii)*costh(l))*tmt3)
      p(4,l)=sinth(l)*(cosmu(kk)*copsi-cosmu(ii)*costh(l))*tmt4
      p(5,l)=sfisq*(amups*tmt2+amusq*tmt1+2.0d0*amumu*tmt3)
      p(6,l)=copsq*tmt1+cfisq*tmt2+2.0d0*copcs*tmt3
      pmats=const*0.5d0*(p(1,l)+p(2,l)+p(5,l)+p(6,l))*4.0d0*pi
      p(7,l)=sinth(l)*(cosmu(kk)*costh(l)*tmt2+cosmu(ii)*copsi*tmt1+
     1         (cosmu(kk)*copsi+cosmu(ii)*costh(l))*tmt3)
      p(8,l)=-p(4,l)
      p(9,l)=2.0d0*sinth(l)*(cosmu(ii)*costh(l)*tmt1+cosmu(kk)*copsi*tmt2+
     1        (cosmu(kk)*costh(l)+cosmu(ii)*copsi)*tmt3)
      p(10,l)=-2.0d0*sinth(l)*(cosmu(ii)*costh(l)*tmt2+cosmu(kk)*copsi*tmt1+
     1        (cosmu(kk)*costh(l)+cosmu(ii)*copsi)*tmt3)
      p(11,l)=(copcs-amumu*sfisq)*(tmt1+tmt2)+(cfisq+copsq-sfisq
     1        *(amusq+amups))*tmt3
      p(12,l)=(cfisq-copsq-sfisq*(amusq-amups))*tmt4
      p(13,l)=2.0d0*sinth(l)*(cosmu(ii)*copsi-cosmu(kk)*costh(l))*tmt4
      p(14,l)=-p(13,l)
      p(15,l)=(copsq-cfisq-sfisq*(amusq-amups))*tmt4
      p(16,l)=(copcs+amumu*sfisq)*(tmt1+tmt2)+(cfisq+copsq+sfisq
     1        *(amusq+amups))*tmt3     
      call depol(l)
      return                                                            
      end                                                               
c
c************************************************************************
c***********************************************************************
 
      subroutine depol(l)                                                
c***********************************************************************
      implicit real*8 (a-h,o-z)
      include 'afrt_rt2.cmn'
c***********************************************************************
      p(17,l)=copsq
      p(18,l)=amusq*sfisq
      p(19,l)=-cosmu(ii)*copsi*sinth(l)
      p(20,l)=0.0d0
      p(21,l)=amups*sfisq
      p(22,l)=cfisq
      p(23,l)=cosmu(kk)*costh(l)*sinth(l)
      p(24,l)=0.0d0
      p(25,l)=2.0d0*cosmu(kk)*copsi*sinth(l)
      p(26,l)=-2.0d0*cosmu(ii)*costh(l)*sinth(l)
      p(27,l)=-amumu*sfisq+copcs
      p(28,l)=0.0d0
      p(29,l)=0.0d0
      p(30,l)=0.0d0
      p(31,l)=0.0d0
      p(32,l)=amumu*sfisq+copcs
c     apply molecular depolarization correction
      if(ipol.eq.1)then
          do i=1,3
              do j=1,3
                 k=16+(i-1)*4+j
                 dgm=bgm
                 if(i.eq.3 .or. j.eq.3) dgm=0.0d0
                 p(k,l)=agm*p(k,l)+dgm
              enddo
          enddo
          p(32,l)=p(32,l)+cgm
      endif
      return                                                            
      end                                                               
c
c************************************************************************