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 'afrt_rt2.cmn'
c
      real*8 fiit(4,2*nsz,nph)
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
                           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*dcmu(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))
                           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
                           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*dcmu(kk)
                     enddo
                     ftmp(is,it,ip)=sumta*ddphi
                  enddo
              endif
          enddo
      enddo
      return
      end
c**********************************************************************