c--------------------------------------------------------------------- c checked version of "integrallib.f" c--------------------------------------------------------------------- c==================================================================== c checked with hollik's routines ; status 16.2.93 ; c==================================================================== c-------------------------------------------------------------------- subroutine fint(s,m1,m2,f) c-------------------------------------------------------------------- c the f-function as specified in : boehm et.al., c fortschr.phys.34 ; and in : hollik, fortschr.phys. 38 c-------------------------------------------------------------------- implicit real*8(a-z) parameter (eps=1d-6) complex*16 f pi = 4d0*datan(1d0) x1 = m1**2 x2 = m2**2 p = s-(dabs(m1)+dabs(m2))**2 m = s-(dabs(m1)-dabs(m2))**2 imagf = 0d0 rest = 0d0 if (abs(s).lt.eps) then realf = 0d0 goto 999 endif c if (x1.lt.eps) then if (x2.lt.eps) then realf = 0d0 else if (s.gt.x2+eps) then realf = 1d0+(1d0-x2/s)*dlog(1d0/(s/x2-1d0)) imagf = pi*(1d0-x2/s) else if (s.lt.x2-eps) then realf = 1d0+(1d0-x2/s)*dlog(1d0/(1d0-s/x2)) else realf = 1d0 endif endif goto 999 else if (x2.lt.eps) then if (x1.lt.eps) then realf = 0d0 else if (s.gt.x1+eps) then realf = 1d0+(1d0-x1/s)*dlog(1d0/(s/x1-1d0)) imagf = pi*(1d0-x1/s) else if (s.lt.x1-eps) then realf = 1d0+(1d0-x1/s)*dlog(1d0/(1d0-s/x1)) else realf = 1d0 endif endif goto 999 endif c if ((abs(s).lt.x1/1d3).or.(abs(s).lt.x2/1d3)) then if (dabs(dabs(m1)-dabs(m2)).lt.eps) then realf = s/6d0/x1 else realf = s/(x1-x2)**2*((x1+x2)/2d0 $ -x1*x2/(x1-x2)*dlog(x1/x2)) endif goto 999 endif c if (dabs(dabs(m1)-dabs(m2)).lt.eps) then rest = 2d0 else rest = 1d0-((x1-x2)/s-(x1+x2)/(x1-x2))*dlog(x1/x2)/2d0 endif c if (m.lt.0d0) then realf = dsqrt(p*m)*dlog((dsqrt(-p)+dsqrt(-m))**2 $ /(4d0*dsqrt(x1*x2)))/s else if (p.lt.0d0) then realf = -2d0*dsqrt(-p*m)*datan(dsqrt(-m/p))/s else realf = -dsqrt(p*m)*dlog((dsqrt(p)+dsqrt(m))**2 $ /(4d0*dsqrt(x1*x2)))/s imagf = dsqrt(p*m)/s*pi endif 999 continue f = dcmplx(rest+realf,imagf) end c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine afunc(m,a0) implicit real*8(a-z) parameter (eps=1d-4) common / dimreg / delta,muedr if (m**2.lt.eps) then a0 = 0d0 else a0 = m**2*(1d0-dlog(m**2/muedr**2)) endif end c--------------------------------------------------------------------- subroutine bfunc(s,m1,m2,b0,b1,b20) c--------------------------------------------------------------------- c b0 and b1 are the b0"quer" and b1"quer"-quantities defined c in : hollik, fortschr.phys. 38 (but with log-terms incl.) c b20 is the g_{mu}.{nu} coefficient of the c tensor-2-point-integral (log-terms incl, too) c--------------------------------------------------------------------- implicit real*8(a-z) parameter (eps=1d-4) complex*16 b0,b1,b20,ff common / dimreg / delta,muedr pi = 4d0*datan(1d0) x1 = m1**2 x2 = m2**2 if (dabs(s).lt.eps) then call bnull(m1,m2,realb0) call beins(m1,m2,realb1) call bzwnull(m1,m2,realb20) b0 = dcmplx(realb0,0d0) b1 = dcmplx(realb1,0d0) b20 = dcmplx(realb20,0d0) goto 999 endif if ((x1.lt.eps).and.(x2.lt.eps)) goto 47 call fint(s,m1,m2,ff) 47 continue if (x1.lt.eps) then if (x2.lt.eps) then b0 = 2d0-dlog(s/muedr**2)+dcmplx(0d0,1d0)*pi b1 = -b0/2d0 else b0 = 1d0+ff-dlog(x2/muedr**2) b1 = -(1d0+(1d0-x2/s)*ff-dlog(x2/muedr**2))/2d0 endif goto 666 endif if (x2.lt.eps) then if (x1.lt.eps) then b0 = 2d0-dlog(s/muedr**2)+dcmplx(0d0,1d0)*pi b1 = -b0/2d0 else b0 = 1d0+ff-dlog(x1/muedr**2) b1 = -(1d0+(1d0+x1/s)*ff-dlog(x1/muedr**2))/2d0 endif goto 666 endif if (dabs(x1-x2).lt.1d3*eps) then b0 = ff b1 = -ff/2d0 else b0 = 1d0-(x1+x2)/(x1-x2)*dlog(x1/x2)/2d0+ff b1 = -1d0/2d0+x1/(x1-x2)*dlog(x1/x2)/2d0-(s+x1-x2)/2d0/s*ff endif b0 = b0 - dlog(x1*x2/muedr**4)/2d0 b1 = b1 + dlog(x2/muedr**2)/2d0 c 666 continue call afunc(m2,a0) b20 = a0/6d0+x1/3d0*b0 $ +(s+x1-x2)/6d0*b1+(x1+x2-s/3d0)/6d0 999 continue end c--------------------------------------------------------------------- subroutine bderiv(x,m1,m2,db0,db1,db20) c--------------------------------------------------------------------- c real parts of d(b0(s,m1,m2))/ds , d(b1(s,m1,m2))/ds c and d(b20(s,m1,m2))/ds c (s -> x in the code) c--------------------------------------------------------------------- implicit real*8(a-z) complex*16 cf,b0,b1,b20 parameter (eps=1d-1) * pi = 4d0*datan(1d0) xm1 = m1**2 xm2 = m2**2 * if (x.lt.eps) then if (dabs(xm1-xm2).lt.eps) then db0 = 1d0/6d0/xm2 db1 = -1d0/12d0/xm2 else if (xm1.lt.eps) then db0 = 1d0/2d0/xm2 db1 = -1d0/6d0/xm2 else if (xm2.lt.eps) then db0 = 1d0/2d0/xm1 db1 = -1d0/6d0/xm1 else db0 = (xm1**2-xm2**2-2d0*xm1*xm2*dlog(xm1/xm2)) $ /2d0/(xm1-xm2)**3 fss0 = (xm1**3-xm2**3+9d0*xm1*xm2*(xm1-xm2) $ -6d0*xm1*xm2*dlog(xm1/xm2)*(xm1+xm2)) $ /3d0/(xm1-xm2)**5 db1 = -db0/2d0-(xm1-xm2)/4d0*fss0 endif goto 11 endif 13 continue if ((x.lt.xm1/1d3).or.(x.lt.xm2/1d3)) then if (dabs(xm1-xm2).lt.eps) then db0 = 1d0/6d0/xm1 else db0 = ((xm1+xm2)/2d0-xm1*xm2/(xm1-xm2)*dlog(xm1/xm2)) $ /(xm1-xm2)**2 endif db1 = -db0/2d0 goto 11 endif c if ((xm1.lt.eps).and.(xm2.lt.eps)) then db0 = -1d0/x db1 = -db0/2d0 goto 11 endif c if (xm2.lt.eps) then if (x.gt.xm1+eps) then deriv = -(1d0+xm1/x*dlog(x/xm1-1d0))/x else if (x.lt.xm1-eps) then deriv = -(1d0+xm1/x*dlog(1d0-x/xm1))/x else deriv = 0d0 endif goto 10 endif c if (xm1.lt.eps) then if (x.gt.xm2+eps) then deriv = -(1d0+xm2/x*dlog(x/xm2-1d0))/x else if (x.lt.xm2-eps) then deriv = -(1d0+xm2/x*dlog(1d0-x/xm2))/x else deriv = 0d0 endif goto 10 endif c sm = xm1+xm2 dm = xm2-xm1 sm12 = (dabs(m1)+dabs(m2))**2 dm12 = (dabs(m1)-dabs(m2))**2 lm = dlog(xm2/xm1)/2d0 * if (x.lt.dm12) then s = dsqrt(sm12-x) d = dsqrt(dm12-x) fact = dlog((d+s)**2/(4d0*dabs(m1)*dabs(m2))) deriv = (dm*lm/x-((d**2+s**2+2d0*d**2*s**2/x) $ /(2d0*s*d))*fact-1d0)/x else if (x.lt.sm12) then if (dabs(xm1-xm2).lt.eps) then if (x.lt.eps) then deriv = 1d0/6d0/xm1 else sx = 4d0*xm1/x deriv = (sx/dsqrt(sx-1d0) $ *datan(1d0/dsqrt(sx-1d0))-1d0)/x endif else * c Achtung: fuer x=sm12 geht deriv -> infty! * s = dsqrt(sm12/x-1d0) d = dsqrt(1d0-dm12/x) fact = datan(d/s) deriv = (dm*lm/x+((d**2-s**2+2d0*d**2*s**2) $ /(s*d))*fact-1d0)/x endif else s = dsqrt(x-sm12) d = dsqrt(x-dm12) fact = dlog((d+s)**2/(4d0*dabs(m1)*dabs(m2))) deriv = (dm*lm/x-((d**2+s**2-2d0*d**2*s**2/x) $ /(2d0*s*d))*fact-1d0)/x endif * 10 continue c-------------------------------------------------------------------- db0 = deriv call fint(x,m1,m2,cf) f = dreal(cf) db1 = -(x+xm1-xm2)/2d0/x*deriv+(xm1-xm2)/2d0/x**2*f c 11 continue c call bfunc(x,m1,m2,b0,b1,b20) db20 = xm1/3d0*db0+dreal(b1)/6d0+(x+xm1-xm2)/6d0*db1 $ -1d0/18d0 c 99 continue end c--------------------------------------------------------------------- subroutine cfuncnew(m,s,m1,m2,m3,c0,c1p,c1m,c20,c2p,c2m) c--------------------------------------------------------------------- c contains the scalar coefficient functions of the c tensor three point integral c--------------------------------------------------------------------- implicit real*8(a-z) complex*16 c0,c1p,c1m,c20,c2p,c2m $ ,b0x32,b0x31,b0s12,b1x32,b1x31,adummy,bdummy x = m**2 x1 = m1**2 x2 = m2**2 x3 = m3**2 c call bfunc(s,m1,m2,b0s12,adummy,bdummy) call bfunc(x,m3,m2,b0x32,b1x32,adummy) call bfunc(x,m3,m1,b0x31,b1x31,adummy) call c0int(m,s,m1,m2,m3,c0) c c1m = ((b0x32+b0x31)/2d0-b0s12 $ -(x+x3-(x1+x2)/2d0)*c0)/(4d0*x-s) c1p = ((b0x32-b0x31)/2d0+(x1-x2)/2d0*c0)/s c20 = b0s12/4d0+x3/2d0*c0-(x1-x2)/4d0*c1p $ -(x1+x2-2d0*(x+x3))/4d0*c1m+1d0/4d0 c2m = (-c20+(b1x32+b1x31)/4d0+b0s12/2d0 $ -(x+x3-(x1+x2)/2d0)*c1m)/(4d0*x-s) c2p = (-c20-(b1x32+b1x31)/4d0+(x1-x2)/2d0*c1p)/s end c--------------------------------------------------------------------- subroutine c0int(mf,s,m1,m2,m3,c0) c************************************************************** c * c the scalar vertex integral with equal external masses mf * c * c------------------------------------------------------------------- c s = momentum transfer; m1,m2,m3 are the internal masses c implicit real*8 (a-y) complex*16 z1,z2,z11,z12,z21,z22,cl1,cl2,cl3,cspen,spence, & int,c0 xmf=mf*mf if (mf.lt.1d-1) then mfstrich = 1d-1 xmf = mfstrich**2 endif c xm's : are fermion and boson masses squared xm1=m1*m1 xm2=m2*m2 xm3=m3*m3 c the t'hooft-veltman parameters a=1.d0 b=xmf/s c=-1.d0 d=xm1-xm2-s e=xm3-xm1-xmf+s f=xm2/s d=d/s e=e/s c discriminante for alpha-equation disc=c*c-4.d0*a*b if (disc .lt. 0.d0) goto 500 al=(-c-dsqrt(disc))/2.d0/b nenner=c+2.d0*al*b c..the first integral............................................. y0=-(d+e*al+2.d0*a+c*al)/nenner y01=y0-1.d0 d1=(c+e)**2-4.d0*b*(a+d+f) x1=-(c+e)/2.d0/b if (d1.gt.0.d0) goto 10 c.......complex zeroes of logarithms sq1=dsqrt(-d1) x2=sq1/2.d0/b z1=dcmplx(x1,x2) z2=dcmplx(x1,-x2) z11=y0/(y0-z1) z12=y01/(y0-z1) z21=y0/(y0-z2) z22=y01/(y0-z2) cl1=spence(z11)-spence(z12)+spence(z21)-spence(z22) goto 15 10 continue c........real zeroes sq1=dsqrt(d1) x2=sq1/2.d0/b y1=x1+x2 y2=x1-x2 sig1= y0/dabs(y0) sig2= y01/dabs(y01) y11=y0/(y0-y1) y12=y01/(y0-y1) y21=y0/(y0-y2) y22=y01/(y0-y2) cl1=cspen(y11,sig1)-cspen(y12,sig2)+cspen(y21,-sig1) & -cspen(y22,-sig2) 15 continue c..the second integral............................................ y0=-(d+e*al)/nenner/(1.d0-al) y01=y0-1.d0 d2=(e+d)**2-4.d0*f*(a+b+c) x1=-(e+d)/2.d0/(a+b+c) if(d2.gt.0.d0) goto 20 c.......complex zeroes of logarithms sq2=dsqrt(-d2) x2=sq2/2.d0/(a+b+c) z1=dcmplx(x1,x2) z2=dcmplx(x1,-x2) z11=y0/(y0-z1) z12=y01/(y0-z1) z21=y0/(y0-z2) z22=y01/(y0-z2) cl2=spence(z11)-spence(z12)+spence(z21)-spence(z22) goto 25 20 continue c........real zeroes x2=dsqrt(d2)/2.d0/(a+b+c) y1=x1+x2 y2=x1-x2 y11=y0/(y0-y1) y12=y01/(y0-y1) y21=y0/(y0-y2) y22=y01/(y0-y2) sig1= y0/dabs(y0) sig2= y01/dabs(y01) cl2=cspen(y11,sig1)-cspen(y12,sig2)+cspen(y21,-sig1) & -cspen(y22,-sig2) 25 continue c..the third integral............................................ y0=(d+e*al)/nenner/al y01=y0-1.d0 d3=d*d-4.d0*a*f x1=-d/2.d0/a if (d3.gt.0.d0) goto 30 c........complex zeroes of logarithms sq3=dsqrt(-d3) x2=sq3/2.d0/a z1=dcmplx(x1,x2) z2=dcmplx(x1,-x2) z11=y0/(y0-z1) z12=y01/(y0-z1) z21=y0/(y0-z2) z22=y01/(y0-z2) cl3=spence(z11)-spence(z12)+spence(z21)-spence(z22) goto 35 30 continue c........real zeroes x2=dsqrt(d3)/2.d0/a y1=x1+x2 y2=x1-x2 31 format(1h ,3e12.4) y11=y0 /(y0-y1) y12=y01/(y0-y1) y21=y0/(y0-y2) y22=y01/(y0-y2) sig1= y0/dabs(y0) sig2= y01/dabs(y01) cl3=cspen(y11,sig1)-cspen(y12,sig2)+cspen(y21,-sig1) & -cspen(y22,-sig2) 35 continue c..summation of the 3 integrals .................................... int = -cl1+cl2-cl3 c0 = int/nenner/s goto 501 500 continue c..error message for complex alpha................................ write(6,21) 21 format(1h ,' i cannot handle a complex alpha') stop 501 return end c-------------------------------------------------------------------- c the x"null" subroutines calculate the b0-, b1-, c0-, c d0- and d20-integrals at zero external momenta c-------------------------------------------------------------------- subroutine bnull(a,b,b0) implicit real*8(a-z) parameter (eps=1d-8) common / dimreg / delta,muedr xa = a**2 xb = b**2 x = xa+xb if (x.lt.eps) then c write(6,*)'all args zero is not allowed for b0-function !' b0 = 0d0 goto 2 endif if (xa*xb.eq.0d0) then zw = 1d0-dlog(x/muedr**2) else zw = -dlog(xa*xb/muedr**4)/2d0 if (dabs(dabs(a)-dabs(b)).gt.eps) then zw = zw+1d0-(xa+xb)/(xa-xb)*dlog(xa/xb)/2d0 endif endif b0 = zw 2 continue end c----------------------------------------------------------------- subroutine beins(a,b,b1) implicit real*8(a-z) parameter (eps=1d-8) common / dimreg / delta,muedr xa = a**2 xb = b**2 x = xa+xb if (x.lt.eps) then c write(6,*)'all args zero is not allowed for b1-function !' b1 = 0d0 goto 2 endif if (xa.eq.0d0) then zw = -(1d0/2d0-dlog(xb/muedr**2))/2d0 else if (xb.eq.0d0) then zw = -(3d0/2d0-dlog(xa/muedr**2))/2d0 else zw = dlog(xb/muedr**2)/2d0 if (dabs(dabs(a)-dabs(b)).gt.eps) then zw = zw -(xa+xb)/(xa-xb)/4d0-1d0/2d0 $ +xa/(xa-xb)*dlog(xa/xb)/2d0 $ +xa*xb/(xa-xb)**2*dlog(xa/xb)/2d0 endif endif b1 = zw 2 continue end c----------------------------------------------------------------- subroutine bzwnull(a,b,b20) implicit real*8(a-z) common / dimreg / delta,muedr xa = a**2 xb = b**2 call bnull(a,b,b0) call beins(a,b,b1) zw = xa*b0/3d0+(xa-xb)*b1/6d0+(xa+xb)/6d0 if (xb.ne.0d0) then zw = zw+xb*(1d0-dlog(xb/muedr**2))/6d0 endif b20 = zw end c----------------------------------------------------------------- subroutine cnull(a,b,c,c0) implicit real*8(a-z) parameter (eps=1d-8) xa = a**2 xb = b**2 xc = c**2 if (dabs(dabs(a)-dabs(b)).lt.eps) then if (dabs(dabs(a)-dabs(c)).lt.eps) then zw = -1d0/2d0/xa else zw = (-1d0+xc/(xa-xc)*dlog(xa/xc))/(xa-xc) endif else call bnull(a,c,bac) call bnull(b,c,bbc) zw = (bac-bbc)/(xa-xb) endif c0 = zw end c----------------------------------------------------------------- subroutine dnull(a,b,c,d,d0) implicit real*8(a-z) integer i,j real*8 m(4) parameter (eps=1d-8) m(1) = a m(2) = b m(3) = c m(4) = d do 30 i = 1,4 do 40 j = i+1,4 if (dabs(dabs(m(i))-dabs(m(j))).lt.eps) then m(i) = m(i)+eps m(j) = m(j)-eps endif 40 continue 30 continue call cnull(a,b,c,c0abc) call cnull(a,b,d,c0abd) zw = (c0abc-c0abd)/(c**2-d**2) d0 = zw end c----------------------------------------------------------------- subroutine dzwnull(a,b,c,d,d20) implicit real*8(a-z) call cnull(b,c,d,c0bcd) call dnull(a,b,c,d,d0abcd) zw = c0bcd+a**2*d0abcd d20 = zw end c-----------------------------------------------------------------