c calculating three dimensional sum for regularized sunset diagram c uses interpolation to do the innermost sum. 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 'headerint' 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 tableloc is the position in the interpolation table to center c the "fit" on, and it uses 2*PTS points either side of central value. c A high order fit (eg 6th or 8th) seems to improve accuracy, so c we use that. integer tableloc,PTS c interpolation values in table real*8 table(0:INTERPVALS),TABLEINT,innersum c dvals real*8 dvals(0:INTERPVALS),error parameter(TABLEINT = 1d0/qreal(interpvals),PTS=4) c 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 c d and ddash are the scaled and unscaled variables used the c interpolation (d actually runs between 0 and infinity, but c we scale it fit within [0,1] real*8 result,factor,mom,d,ddash c defines functions called by program real*8 scaledvar, unscaledvar integer maxloc,InnerCompQ 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=200000000000q0 factor = h**3 c initialize p's and c's (assuming Lambda-> infty for now) do loops1 = 1,PROPS do loops2 = MAXNEG,MAXPOS c(loops1,loops2) = exp(qreal(loops2)*h) c uncomment this piece to include lambda terms. c c c(loops1,loops2)) =c(loops1,loops2) + c mass*mass/(lambda*lambda) if(c(loops1,loops2) .ge.0d0) then p(loops1,loops2) = exp(qreal(loops2)*h - . exp(qreal(loops2)*h))/ . (c(loops1,loops2)*c(loops1,loops2)) else p(loops1,loops2) = 1d0 endif enddo enddo do momvals = 1,40 mom=qreal(momvals)/10. 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 c interpolation table set up here. do loops1 = 0,INTERPVALS-1 loops2 = 0 inner = 1 direct(PROPS)=1 ddash = qreal(loops1)/qreal(INTERPVALS) dvals(loops1) = ddash d = unscaledvar(ddash) table(loops1) = innersum(ddash,mom) enddo c This last case corresponds to have d->infinity in the denominator, c the inner sum is identically zero. table(INTERPVALS)=0d0 dvals(INTERPVALS)=1d0 c now we do the remaining two outer sums finished= 0 clicker =1 result=TOL c this is the main loop do while(finished.eq.0) c We are doing PROPS-1 loops, so deptj is set differently. depth = PROPS-1 running(depth) = 0d0 inner = 1 do while(inner.eq.1) c looks up value from interpolation table, after setting d. d = 1d0/c(1,k(1)) + 1d0/c(2,k(2)) ddash = scaledvar(d) tableloc = int(dreal(INTERPVALS)*ddash) c checks we are not at the end of the range. if((tableloc).le.INTERPVALS) then if ((tableloc).gt. INTERPVALS-PTS ) . tableloc =INTERPVALS-PTS if (tableloc .lt. PTS) tableloc =PTS c does the interpolation call polint(dvals(tableloc-PTS+1), . table(tableloc-PTS+1),PTS*2, . ddash,partial(depth,k(depth)),error) partial(depth,k(depth)) = partial(depth,k(depth))* . p(1,k(1))*p(2,k(2)) else partial(depth,k(depth)) = 0d0 endif c adjusts running total. 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(running(depth)).eq.0) .or. & (abs(k(depth)-startvals(depth)).lt.2)).and. & (ddash.lt.1d0 -dreal(PTS)*TABLEINT ))then c add one to inner most loop k(depth) = k(depth) + direct(depth) c if k is increasing, turn around and work backwards. k(depth) = k(depth) + direct(depth) 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) = 0d0 result = result+running(depth) c if terms no longer contribute to sum, either turn round or finish c write(6,*) k(1),k(2),k(3) c write(6,*) partial(depth,k(depth)),result c write(6,*) partial(depth,k(depth)-direct(depth)) if((partial(depth,k(depth))/result.gt.TOL).or. & (partial(depth,k(depth)).gt. & partial(depth,k(depth)-direct(depth))).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 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 real*8 result include 'headerint' InnerCompQ = partial(depth,k(depth))/result.lt.(TOL).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? c include 'headerint' 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 assumes does not make a trip through the while loop if c test is false to begin with. 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 real*8 function innerterm(term,d,mom) include 'headerint' real*8 expfac,prefac,d real*8 mom integer term prefac = p(3,term)/((d + 1./c(3,term) )**2) expfac = mom*mom/(mass*mass*(d + 1./c(3,term))) if(expfac.gt.1q-6) then innerterm = prefac*(exp(-expfac) - 1q0 + expfac) else innerterm = prefac*expfac*expfac*(.5q0 -expfac/6q0 . + expfac*expfac/24q0 ) endif end c this does the innersum (called when setting up the interpolation c table) real*8 function innersum(ddash,mom) include 'headerint' integer counter,direct,inner,heading real*8 ddash,d,lastterm,innerterm,mom,unscaledvar inner = 1 heading=1 d = unscaledvar(ddash) innersum = 0d0 counter=0 do while(inner.eq.1) lastterm = innerterm(counter,d,mom) innersum = innersum+lastterm if((lastterm/innersum).gt.TOL) then counter = counter+heading else if(heading.eq.1) then counter = -1 heading=-1 else inner = 0 endif endif enddo return end c maps [0,infty) into [0,1] for interpolation. real*8 function scaledvar(uns) real*8 uns uns = uns**(.25) scaledvar = uns/(8d0 + uns) return end c inverse of scaledvar. real*8 function unscaledvar(scl) real*8 scl unscaledvar = (8d0 * scl/(1d0 - scl))**4 return end c does polynomial interpolation. subroutine polint(xvals,yvals,n,x,y,dy) integer n,NMAX real*8 dy,x,y,yvals(n),xvals(n) PARAMETER (NMAX = 10) integer i,m,ns real*8 den,dift,ho,hp,w,c(NMAX),d(NMAX) ns = 1 dif = abs(x-xvals(1)) do i =1,n dift = abs(x-xvals(i)) if (dift.lt.dif) then ns=i dif=dift endif c(i) = yvals(i) d(i) = yvals(i) enddo y= yvals(ns) ns = ns-1 do m=1,n-1 do i=1,n-m ho = xvals(i) -x hp = xvals(i+m) -x w = c(i+1) -d(i) den = ho-hp den = w/den d(i) = hp*den c(i) = ho*den enddo if(2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns) ns=ns-1 endif y=y+dy enddo return end