c     calculating three dimensional sum for regularized sunset diagram

c     copyright Richard Easther, Gerry Guralnik and Stephen Hahn (1999)
c
c     You are free to use, copy and modify this code.
c     However, please ask us before distributing a modified version.



c     the "header" file loads the definitions for the common blocks.

      include 'header'

c     finished,inner and clicker are flags used to determine when
c     loops have finished.

      integer finished,inner,clicker

c     Since the sums run from -\infty to +\infty, the actual starting
c     point is arbitrary. We choose it to be the value of index for which
c     the general term was a maximum on the last pass through the same loop.
c     We store the last starting value in startvals.

      integer startvals(MAXPROPS)

c     these variables are used by various loops in the code.

      integer loops1,loops2,momvals
       
c     result is the answer, and factor is the multiplicative prefactor. 
c     In general the sums contain a term like m^2/(4 Pi)^n WHICH WE HAVE
c     OMITTED IN ALL THESE CALCULATIONS!

c     mom is the the value of the momentum for which the diagram is 
c     calculated.

      real*8 result,factor,mom

c     these are functions called by the the summation code.

      integer maxloc,InnerCompQ,OuterCompQ

c     diag(mom) gives returns the general term as a function of the k_i
c     indices.

      real*8 diag

      
c     the value of h is defined here.

      h =  .9q0
     
     
c     mass and cutoff set. Note that this code never actually *uses* the
c     cutoff since we are computing the renormalized contribution.

      mass=1q0
      lambda=4q0
      factor =  h**PROPS
     

c     initialize p's and c's (assuming Lambda-> infty)

      do loops1 = 1,PROPS

        
         do loops2 = MAXNEG,MAXPOS

            c(loops1,loops2) = exp(qreal(loops2)*h) 
     .           + mass*mass/(lambda*lambda)


c     This is an ugly hack to prevent a divide by zero error, since the c()
c     can get very small for extreme values of k. For reasonable values
c     h the code never sees these values, but we should do some proper
c     range checking on the k_i!

            if(c(loops1,loops2) .ge.0q0) then

               p(loops1,loops2) = exp(qreal(loops2)*h - 
     .              exp(qreal(loops2)*h))/
     .              (c(loops1,loops2)*c(loops1,loops2))

            else
               p(loops1,loops2) = 1q0
            endif

         enddo
      enddo


      do momvals = 1,40

      mom=qreal(momvals)/10q0
     

c     initialize variables before using them.
c     running is the running total for the i-th nested sum
c     k(i) are the indices of the sums.
c     direct(i) determines whether the summation index is increasing or
c     decreasing. startvals is where we began the sum from.

      do loops1 = 1,MAXPROPS
         running(loops1)= 0q0
         k(loops1)= 0
         direct(loops1)=  1
         startvals(loops1)= 0
      enddo


      finished= 0
      clicker =1

c     we set this to avoid a divide by zero error on the first run through.
    
      result=diag(mom)
      

c     this is the main loop

      do while(finished.eq.0)

c     the counter being iterated is the innermost loop, so set depth 
c     to the number of  props.  Then begin the loop for the innermost sum.

         depth = PROPS
         running(depth) = 0q0
         inner = 1

         do while(inner.eq.1)

            partial(depth,k(depth)) = diag(mom)
            running(depth) =running(depth)+partial(depth,k(depth))
          

c     decide whether to turn round, or to finish. The test is 
c     handled by a call to InnerCompQ which returns a "yes" if the
c     terms are no longer significant.


            if((InnerCompQ(result).eq.0) .or.
     &           (abs(k(depth)-startvals(depth)).lt.2))then 

c     add one to inner most loop
               k(depth) = k(depth) + direct(depth)

c     if k is increasing, turn around and work backwards.

            else              
               if(direct(depth).eq.1) then
                  direct(depth)=  -1
                  endpos(depth) = k(depth)
                  k(depth) = startvals(depth)-1

c     or finish this loop.
               else
                  direct(depth)=1
                  endneg(depth)=k(depth)

c     locate next starting point
                  startvals(depth)= maxloc(startvals(depth))
                  k(depth) = startvals(depth)
                  inner =0
               endif
            endif
         enddo

c     end of inner loop. 



c     this handles the "outer" loops.

         clicker =1 
         do while((clicker.eq.1).and.(depth.gt.1))

c     depth says which of the sums we are looking (depth=1 is the outermost)

            depth=depth-1
            partial(depth,k(depth)) = running(depth+1)
            running(depth) = running(depth) + running(depth+1)
            running(depth+1) = 0q0
            
            if(depth.eq.props-1)  result = result+running(depth)


c     if terms no longer contribute to sum, either turn round or finish
c     analagous to procedure for innermost sum.
            
            if((OuterCompQ(result).eq.0).or.
     &           (abs(k(depth)-startvals(depth)).lt.2)) then
               
               k(depth) = k(depth) + direct(depth)
               clicker =0
            else              
               if(direct(depth).eq.1) then
                  direct(depth)=  -1
                  endpos(depth) = k(depth)
                  k(depth) = startvals(depth)-1
                  clicker=0
               else
                  direct(depth)=1
                  endneg(depth)=k(depth)
                  
                  if(depth.eq.1) then
                     finished=1
                  else

                     startvals(depth)= maxloc(startvals(depth))
                     k(depth)=startvals(depth)

                  endif
               endif
            endif
         enddo



      enddo

c     When we reach here the sum is finished. Print out the results.

      write(6,*) mom
      write(6,*) running(1) * factor

      enddo
      end


c     this does the aitken extrapolation (as described in the paper)
c     it is called by the innermost sum.

      real*8 function aitken(deep)

      integer deep
      include 'header'

      aitken = running(deep)+partial(deep,k(deep))**2 / 
     &     ( partial(deep,k(deep)-direct(deep))
     &     - partial(deep,k(deep)))


      return
      end

      integer function InnerCompQ(result)

c     when InnerCompQ = 0, successive terms in the inner loop
c     are still making a signficant contribution to the result
c     
c     this loop also uses the aitken delta^2 extrapolation and
c     if this is "within tolerance" it replaces running with the
c     improved result, and terminates the sum.

      real*8 result,aitken

      include 'header'

c     This determines whether the sum is "completed". TOL is a global
c     parameter, but since the round-off error is additive, we demand
c     more sensitivity from the inner sums that we do from the outer 
c     sums. Hence the TOL/(5**PROPS). 
c
c     result is the OVERALL running total.
c     running(depth) is the running total for the innermost sum.
c    
c     We also check that the general term is decreasing, obviously
c
c     This routine is very heuristic, and I strongly suspect it could be
c     improved upon.
c      
      InnerCompQ = partial(depth,k(depth))/result.lt.(TOL/10d0).and. 
     &           (partial(depth,k(depth)).lt.
     &           partial(depth,k(depth)-direct(depth)))

      if((direct(PROPS) .eq.-1 .and. InnerCompQ.eq.0)) then
         if(abs(aitken(props)- 
     &        lastaitken(props))/running(props).lt.TOL) then
            running(PROPS) = aitken(PROPS)
            InnerCompQ = 1
            lastaitken(PROPS) = 0q0
         else
            lastaitken(props) = aitken(props)
         endif         
      endif


      return
      
      end

      integer function OuterCompQ(result)

c     when OuterCompQ = 0, successive terms in the k^th loop
c     are still making a signficant contribution to the result
c
c     This function is analagous to InnerCompQ

      real*8 result

      include 'header'
      



      OuterCompQ = ((partial(depth,k(depth))/result.lt.
     &           TOL*4**(PROPS-depth)).and. 
     &           (partial(depth,k(depth)).lt.
     &           partial(depth,k(depth)-direct(depth))))

  
      return
      
      end


      integer function maxloc(midpt)

c
c     This function locates the maximal value of partial(depth,k_n)
c     in terms of k. It provides the starting value of k for the
c     next trip through the loop, thereby optimizing performance
c
c     This is very crude, but it works.
c     Replace with parabolic version later, for better convergence?


      include 'header'
 
      integer lowend,midpt,hiend


      lowend=endneg(depth)
      hiend=endpos(depth)


c     if endpoints are close, just split the difference and return.

      if((hiend-lowend).lt.4) then
         maxloc= (hiend+lowend)/2
         return
      endif

c     check to make sure lowend,midpt and hiend really do bracket a maximum
c     This assumes that the code does not make a trip through the while loop 
c     if the test is false to begin with.  This may be compiler dependent.


      do while((partial(depth,lowend).gt.partial(depth,midpt)).or.
     &     (partial(depth,hiend).gt.partial(depth,midpt)))

         if(partial(depth,lowend).gt.partial(depth,hiend)) then
            hiend = midpt
            midpt = (lowend+hiend)/2
         else
            lowend = midpt
            midpt = (lowend+hiend)/2
         endif

         if((hiend-lowend).lt.4) then
            maxloc=(hiend+lowend)/2
            return
         endif
      enddo

c     now we know that partial(depth,midpt) is bigger than 
c     partial(depth,lowend) and partial(depth,hiend)

      
      do while((hiend-lowend).gt.8)
      
         if(partial(depth,(midpt+lowend)/2).gt.
     &        partial(depth,midpt)) then
            hiend=midpt
            midpt = (lowend+hiend)/2
         else
            if(partial(depth,(midpt+hiend)/2).gt.
     &           partial(depth,midpt)) then
               
               lowend=midpt
               midpt= (lowend+hiend)/2
            else
               lowend= (lowend+midpt)/2
               hiend= (midpt+hiend)/2
            endif
         endif
      enddo

      maxloc = midpt
      return
      end

c     compute the general term.


      real*8 function diag(mom)


      include 'header'
      real*8 expfac,prefac
      real*8 mom



      prefac =(c(1,k(1))**2*c(2,k(2))**2*c(3,k(3))**2*
     -    c(4,k(4))**2*c(5,k(5))**2*c(6,k(6))**2*
     -    c(7,k(7))**2*p(1,k(1))*p(2,k(2))*p(3,k(3))*
     -    p(4,k(4))*p(5,k(5))*p(6,k(6))*p(7,k(7)))/
     -  (c(4,k(4))*(c(3,k(3))*c(5,k(5))*c(6,k(6)) + 
     -        c(3,k(3))*c(5,k(5))*c(7,k(7)) + 
     -        c(3,k(3))*c(6,k(6))*c(7,k(7)) + 
     -        c(5,k(5))*c(6,k(6))*c(7,k(7))) + 
     -     c(2,k(2))*(c(3,k(3))*c(4,k(4))*c(6,k(6)) + 
     -        c(3,k(3))*c(5,k(5))*c(6,k(6)) + 
     -        c(4,k(4))*c(5,k(5))*c(6,k(6)) + 
     -        c(3,k(3))*c(4,k(4))*c(7,k(7)) + 
     -        c(3,k(3))*c(5,k(5))*c(7,k(7)) + 
     -        c(4,k(4))*c(5,k(5))*c(7,k(7)) + 
     -        c(3,k(3))*c(6,k(6))*c(7,k(7)) + 
     -        c(5,k(5))*c(6,k(6))*c(7,k(7))) + 
     -     c(1,k(1))*(c(3,k(3))*c(4,k(4))*c(6,k(6)) + 
     -        c(3,k(3))*c(5,k(5))*c(6,k(6)) + 
     -        c(3,k(3))*c(4,k(4))*c(7,k(7)) + 
     -        c(3,k(3))*c(5,k(5))*c(7,k(7)) + 
     -        c(3,k(3))*c(6,k(6))*c(7,k(7)) + 
     -        c(4,k(4))*c(6,k(6))*c(7,k(7)) + 
     -        c(5,k(5))*c(6,k(6))*c(7,k(7)) + 
     -        c(2,k(2))*
     -         (c(4,k(4))*c(6,k(6)) + 
     -           c(5,k(5))*c(6,k(6)) + 
     -           c(4,k(4))*c(7,k(7)) + 
     -           c(5,k(5))*c(7,k(7)) + 
     -           c(6,k(6))*c(7,k(7)))))**2

      expfac =(mom**2*c(1,k(1))*
     -    (c(4,k(4))*(c(5,k(5))*c(6,k(6))*c(7,k(7)) + 
     -         c(3,k(3))*
     -          (c(5,k(5))*c(6,k(6)) + 
     -            c(5,k(5))*c(7,k(7)) + 
     -            c(6,k(6))*c(7,k(7)))) + 
     -      c(2,k(2))*
     -       (c(5,k(5))*
     -          (c(4,k(4))*c(6,k(6)) + 
     -            c(4,k(4))*c(7,k(7)) + 
     -            c(6,k(6))*c(7,k(7))) + 
     -         c(3,k(3))*
     -          (c(4,k(4))*c(6,k(6)) + 
     -            c(5,k(5))*c(6,k(6)) + 
     -            c(4,k(4))*c(7,k(7)) + 
     -            c(5,k(5))*c(7,k(7)) + 
     -            c(6,k(6))*c(7,k(7))))))/
     -  (mass**2*(c(4,k(4))*
     -       (c(5,k(5))*c(6,k(6))*c(7,k(7)) + 
     -         c(3,k(3))*
     -          (c(5,k(5))*c(6,k(6)) + 
     -            c(5,k(5))*c(7,k(7)) + 
     -            c(6,k(6))*c(7,k(7)))) + 
     -      c(1,k(1))*
     -       (c(3,k(3))*c(4,k(4))*c(6,k(6)) + 
     -         c(3,k(3))*c(5,k(5))*c(6,k(6)) + 
     -         c(3,k(3))*c(4,k(4))*c(7,k(7)) + 
     -         c(3,k(3))*c(5,k(5))*c(7,k(7)) + 
     -         c(3,k(3))*c(6,k(6))*c(7,k(7)) + 
     -         c(4,k(4))*c(6,k(6))*c(7,k(7)) + 
     -         c(5,k(5))*c(6,k(6))*c(7,k(7)) + 
     -         c(2,k(2))*
     -          (c(4,k(4))*c(6,k(6)) + 
     -            c(5,k(5))*c(6,k(6)) + 
     -            c(4,k(4))*c(7,k(7)) + 
     -            c(5,k(5))*c(7,k(7)) + 
     -            c(6,k(6))*c(7,k(7)))) + 
     -      c(2,k(2))*
     -       (c(5,k(5))*
     -          (c(4,k(4))*c(6,k(6)) + 
     -            c(4,k(4))*c(7,k(7)) + 
     -            c(6,k(6))*c(7,k(7))) + 
     -         c(3,k(3))*
     -          (c(4,k(4))*c(6,k(6)) + 
     -            c(5,k(5))*c(6,k(6)) + 
     -            c(4,k(4))*c(7,k(7)) + 
     -            c(5,k(5))*c(7,k(7)) + 
     -            c(6,k(6))*c(7,k(7))))))



      diag = prefac*exp(-expfac)

      end