mdiffs.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), fiiz(4,2*nsz,jpart,nph)
c**************************************************************************
c
      do kk=1,nmum1
          do ll=1,jpart
              do is=1,4
                  fiit(is,kk,ll)=fiic(is,kk,ll)
              enddo
          enddo
          do ll=jpart+1,nophi+1
              fiit(1,kk,ll)=fiic(1,kk,nophi-ll+2)
              fiit(2,kk,ll)=fiic(2,kk,nophi-ll+2)
              fiit(3,kk,ll)=-fiic(3,kk,nophi-ll+2)
              fiit(4,kk,ll)=-fiic(4,kk,nophi-ll+2)
          enddo
      enddo
      do kk=1,nmum1
          do ip=1,jpart
              do ll=1,nophi+1
                  mmp=ll+ip-1
                  if(ll.gt.jpart) then
                      mmp=nophi*(1+(mmp-1)/nophi)-mmp+2
                  endif
                  do is=1,4
                      fiiz(is,kk,ip,ll)=fiit(is,kk,mmp)
                  enddo
              enddo
          enddo
      enddo
      do ii=1,nmum1
          do kk=1,nmum1
              do ll=jpart+1,nophi+1
                 do is=1,4
                     do j=1,4
                        ij=(is-1)*4+j
                        if (is.lt.3) then
                            ppin(ij,ll,ii,kk)=ppin(ij,nophi-ll+2,ii,kk)
                            ppin(ij+16,ll,ii,kk)=ppin(ij+16,nophi-ll+2,ii,kk)
                        else
                            ppin(ij,ll,ii,kk)=-ppin(ij,nophi-ll+2,ii,kk)
                            ppin(ij+16,ll,ii,kk)=-ppin(ij+16,nophi-ll+2,ii,kk)
                        endif
                     enddo
                 enddo
              enddo
          enddo
      enddo
 
      do it=ib,ie
          do ip=1,jpart
              do is=1,4
                 sumta=0.0d0
                 do kk=1,nmum1
                    sumtb=0.0d0
                    do ll=1,nophi
                       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))*fiiz(j,kk,ip,ll)/dlyr
                       enddo
                       sumtb=sumtb+prod1
                    enddo
                    sumta=sumta+sumtb*dcmu(kk)
                 enddo
                 ftmp(is,it,ip)=sumta*ddphi
              enddo
          enddo
      enddo
      return
      end
c**********************************************************************