mdiffn.f

Go to the documentation of this file.

      subroutine mdiffn(ib,ie,il,dlyr)
c
c     compute the integral: int(pidw)
c
c*************************************************************************
c.....includes the common blocks
      implicit real*8 (a-h,o-z)
      include 'common_all.cmn'
c
      real*8 fiit(4,50,46)
c**************************************************************************
c
      do it=ib,ie
      do ip=1,jpart
      if(ip.eq.1 .or. ip.eq.jpart)then
      do is=1,4
         sumta=0.0d0
         do kk=1,nmum1
            sumtb=0.0d0
            do ll=1,jpart
               if(ip.eq.1)then
                 llp=ll
                 ipc=ip+ll-1
               else
                 llp=jpart-ll+1
                 ipc=ip+ll-1
               endif
               do ic=1,4
                  fiit(ic,kk,llp)=fiic(ic,kk,llp)
               enddo
               if(ipc.gt. jpart)then
                 fiit(3,kk,llp)=-fiic(3,kk,llp)
                 fiit(4,kk,llp)=-fiic(4,kk,llp)
               endif
c
               prod1=0.0d0
               do j=1,4
                  ij=(is-1)*4+j
                  prod1=prod1+(c(kk)*tmsl*ppin(ij,ll,it,kk)+
     1                trsl*ppin(ij+16,ll,it,kk))*fiit(j,kk,llp)/dlyr
               enddo
               if(ll.eq.1 .or. ll.eq.jpart)then
                 sumtb=sumtb+prod1
               else if(is.le.2)then
                 sumtb=sumtb+2.0d0*prod1
               endif
            enddo
            sumta=sumta+sumtb*dmu(kk)
         enddo
         ftmp(is,it,ip)=sumta*ddphi
      enddo
      else
      do is=1,4
         sumta=0.0d0
         do kk=1,nmum1
            sumtb=0.0d0
            do ll=1,nophi
               mnz=ll+ip-2
               mmp=mnz+1-(nophi*(mnz/nophi))
c
               if(ll.le.jpart .and. mmp.le.jpart)then
                 iflg1=0
                 iflg2=0
                 llp=ll
               else if(ll.le.jpart .and. mmp.gt.jpart)then
                 iflg1=1
                 iflg2=0
                 llp=ll
                 mmp=nophi-mmp+2
               else if(ll.gt.jpart .and. mmp.le.jpart)then
                 iflg1=1
                 iflg2=1
                 llp=nophi-ll+2
               else if(ll.gt.jpart .and. mmp.gt.jpart)then
                 iflg1=0
                 iflg2=1
                 llp=nophi-ll+2
                 mmp=nophi-mmp+2
               endif
               do ic=1,4
                  fiit(ic,kk,llp)=fiic(ic,kk,llp)
               enddo
c
               if(iflg1.eq.1)then
                 fiit(3,kk,mmp)=-fiic(3,kk,mmp)
                 fiit(4,kk,mmp)=-fiic(4,kk,mmp)
               endif
               prod1=0.0d0
               do j=1,4
                  ij=(is-1)*4+j
                 prod1=prod1+(c(kk)*tmsl*ppin(ij,llp,it,kk)+
     1             trsl*ppin(ij+16,llp,it,kk))*fiit(j,kk,mmp)/dlyr
               enddo
               if(is.ge.3 .and. iflg2.eq.1)then
                 prod1=-prod1
               endif
               sumtb=sumtb+prod1
            enddo
            sumta=sumta+sumtb*dmu(kk)
         enddo
         ftmp(is,it,ip)=sumta*ddphi
      enddo
      endif
      enddo
      enddo
      return
      end
c**********************************************************************