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 = .4q0
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=2000000q0
factor = h**3
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)
c If you are using a finite lambda (cut-off) add this term to the c(,)
c above: 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 proceedure 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/(5**PROPS)).and.
& (partial(depth,k(depth)).lt.
& partial(depth,k(depth)-direct(depth))).and.
& (partial(depth,k(depth))/running(depth).lt.TOL*10)
c Do the aitken extrap if k is decreasing.
if((direct(PROPS) .eq.-1 .and. InnerCompQ.eq.0)) then
if(abs(aitken(props)-
& lastaitken(props))/running(props).lt.TOL/10q0) 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/(5**depth)).and.
& (partial(depth,k(depth)).lt.
& partial(depth,k(depth)-direct(depth))).and.
& (partial(depth,k(depth))/running(depth).lt.TOL*10)
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 = (p(1,k(1)+k(3))*p(2,k(2)+k(3))*p(3,k(3)))/
. ((1q0/c(1,k(1)+k(3)) + 1q0/c(2,k(2)+k(3)) +
. 1q0/c(3,k(3)))**2)
expfac = mom*mom/(mass**2*(1q0/c(1,k(1)+k(3)) +
. 1q0/c(2,k(2)+k(3)) + 1q0/c(3,k(3)) ))
c if expfac is very small, use the Taylor series to avoid subtracting
c numbers with a very similar size....
if(expfac.gt.1q-5) then
diag = prefac*(exp(-expfac) - 1q0 + expfac)
else
diag = prefac*expfac*expfac*(.5q0 - expfac/6q0
. + expfac*expfac/24q0 -
. expfac*expfac*expfac/120q0)
endif
end