double precision function df(s,ma,mb) * c derivative of the real part of the function f(s,ma,mb) * implicit real*8(a-z) parameter(eps=1.d-6) if(s.lt.(ma-mb)**2) then rplus = dsqrt((ma+mb)**2-s) rmin = dsqrt((ma-mb)**2-s) df = ((mb**2-ma**2)*dlog(mb**2/ma**2)/2d0/s- . ((rmin**2+2.d0*rmin**2*rplus**2/s+rplus**2)/(2.d0*rmin*rplus))* . dlog((rmin+rplus)**2/(4.d0*ma*mb))-1.d0 )/s elseif(s.lt.(ma+mb)**2) then rplus = dsqrt((ma+mb)**2/s-1.d0) rmin = dsqrt(1.d0-(ma-mb)**2/s) df = ((mb**2-ma**2)*dlog(mb**2/ma**2)/2d0/s+ . ((rmin**2+2.d0*rmin**2*rplus**2-rplus**2)/(rmin*rplus))* . datan(rmin/rplus)-1.d0 )/s else rplus = dsqrt(s-(ma+mb)**2) rmin = dsqrt(s-(ma-mb)**2) df = ((mb**2-ma**2)*dlog(mb**2/ma**2)/2d0/s- . ((rmin**2-2.d0*rmin**2*rplus**2/s+rplus**2)/(2.d0*rmin*rplus))* . dlog((rmin+rplus)**2/(4.d0*ma*mb))-1.d0 )/s * y = dsqrt(s-4.d0*ma**2) * z = dsqrt(s) * df = ((y**2-s)*dlog((z+y)/(z-y))/(2.d0*z*y)-1.d0)/s endif end * double precision function g(s,ma,mb) * c imaginary part of the function f(s,ma,mb) * implicit double precision (a-z) pi = 4d0*datan(1d0) g = 0.d0 pm = (dabs(ma)+dabs(mb))**2 mm = (dabs(ma)-dabs(mb))**2 if(s.gt.pm) g = pi*dsqrt((s-pm)*(s-mm))/s end * double precision function p(s,m) *************************************************************** * real part of the 1-loop qed vacuumpolarisation contribution * * from a fermion with mass m * * * * relation with the function f: * * p(s,m) = 1/3 - (1 + 2m**2/s)*f(s,m,m) * *************************************************************** implicit double precision (a-z) if(s.eq.0) then p = 0.d0 else if(s.lt.0) then x = dsqrt(1.d0-4.d0*m**2/s) p = -8.d0/3.d0+x**2- . x*(3.d0-x**2)*dlog(-4.d0*m**2/(s*(1.d0+x)**2))/2.d0 else if(s.lt.4.d0*m**2) then x = dsqrt(4.d0*m**2/s-1.d0) p = -8.d0/3.d0-x**2+x*(3.d0+x**2)*datan(1.d0/x) else x = dsqrt(1.d0-4.d0*m**2/s) p = -8.d0/3.d0+x**2- . x*(3.d0-x**2)*dlog(4.d0*m**2/(s*(1.d0+x)**2))/2.d0 endif end * complex*16 function l2(x) implicit double precision (a-z) parameter(z2=1.64493406684823d0) * statement function (sp(u)=li2(x), u=-log(1-x) ): sp(u) = ((((((-1.d0/10886400.d0)*u**2+1.d0/211680.d0)*u**2- . 1.d0/3600.d0)*u**2+1.d0/36.d0)*u-1.d0/4.d0)*u+1.d0)*u * pi = 4d0*datan(1d0) rel2 = -7.d0/2.d0-2.d0*x-(x+3.d0/2.d0)*dlog(x**2) if(x.gt.1) then u = -dlog(1.d0+1.d0/x) iml2 = -pi*(2.d0*(1.d0+x)**2*u+3.d0+2.d0*x) rel2 = rel2+2.d0*(1.d0+x)**2*(-sp(u)-u*dlog(x)) elseif(x.gt.0.d0) then u = -dlog(1.d0+x) u1 = dlog(x) iml2 = -pi*(2.d0*(1.d0+x)**2*(u+log(x))+3.d0+2.d0*x) rel2 = rel2+2.d0*(1.d0+x)**2*(sp(u)+z2-u1**2/2.d0-u*u1) elseif(x.gt.-1.d0/2.d0) then u = -log(1.d0+x) u1 = log(-x) iml2 = 0.d0 rel2 = rel2+2.d0*(1.d0+x)**2*(sp(u)-2.d0*z2-u1**2/2.d0-u*u1) elseif(x.gt.-2.d0) then u = -log(-x) iml2 = 0.d0 rel2 = rel2+2.d0*(1.d0+x)**2*(-sp(u)-z2-u**2/2.d0) else u = -log(1.d0+1.d0/x) iml2 = 0.d0 rel2 = rel2+2.d0*(1.d0+x)**2*(-sp(u)-u*log(-x)) endif l2=dcmplx(rel2,iml2) end * complex*16 function l3(x) implicit double precision (a-z) pi = 4d0*datan(1d0) if(x.gt.1.d0/4.d0) then sq =sqrt(4.d0*x-1.d0) f =atan(1.d0/sq) iml3=0.d0 rel3=(2.5d0-2.d0*x+(4.d0*x+2.d0)*sq*f- $ 8.d0*x*(x+2.d0)*f**2)/3.d0 elseif(x.gt.0.d0) then sq =sqrt(1.d0-4.d0*x) f =log((1.d0+sq)/(1.d0-sq)) iml3=-pi*((2.d0*x+1.d0)*sq+2.d0*x*(x+2.d0)*f)/3.d0 rel3=(2.5d0-2.d0*x+(2.d0*x+1.d0)*sq*f+ . 2.d0*x*(x+2.d0)*(f**2-pi**2))/3.d0 else sq =sqrt(1.d0-4.d0*x) f =log((sq+1.d0)/(sq-1.d0)) iml3=0.d0 rel3=(2.5d0-2.d0*x+(2.d0*x+1.d0)*sq*f+ $ 2.d0*x*(x+2.d0)*f**2)/3.d0 endif l3=dcmplx(rel3,iml3) end * double precision function f(s,ma,mb) * c real part of the function f(s,ma,mb) * implicit real*8(a-z) parameter(eps=1.d-6) ma2=ma**2 mb2=mb**2 if(abs(s).lt.eps) then f=0.0d0 elseif(dabs(ma).lt.eps) then if(s.gt.mb2+eps) then f=1.d0+(1.d0-mb2/s)*dlog(1./(s/mb2-1.d0)) elseif(s.lt.mb2-eps) then f=1.d0+(1.d0-mb2/s)*dlog(1./(1.d0-s/mb2)) else f=1.d0 endif elseif(dabs(mb).lt.eps) then if(s.gt.ma2+eps) then f=1.d0+(1.d0-ma2/s)*dlog(1./(s/ma2-1.d0)) elseif(s.lt.ma2-eps) then f=1.d0+(1.d0-ma2/s)*dlog(1./(1.d0-s/ma2)) else f=1.d0 endif else if(dabs(dabs(mb)-dabs(ma)).lt.eps) then f=2.d0 else f=1.d0+((ma2-mb2)/s-(ma2+mb2)/(ma2-mb2))*dlog(mb2/ma2)/2d0 endif if(s.ge.(dabs(ma)+dabs(mb))**2) then rplus=dsqrt(s-(dabs(ma)+dabs(mb))**2) rmin =dsqrt(s-(dabs(ma)-dabs(mb))**2) f=f-rplus*rmin*dlog((rplus+rmin)**4/(16.d0*ma2*mb2))/2d0/s elseif(s.lt.(dabs(ma)-dabs(mb))**2) then rplus=dsqrt((dabs(ma)+dabs(mb))**2-s) rmin =dsqrt((dabs(ma)-dabs(mb))**2-s) f=f+rplus*rmin*dlog((rplus+rmin)**4/(16.d0*ma2*mb2))/2d0/s else rplus=dsqrt((dabs(ma)+dabs(mb))**2-s) rmin =dsqrt(s-(dabs(ma)-dabs(mb))**2) f=f-2.d0*rplus*rmin*datan(rmin/rplus)/s endif endif end c********************************************************* * subroutine cfunc(s,mf,m1,m2,m3,c0,c1p,c1m,c20,c2p,c2m) * c definition of the invariant functions in the 3-point integrals c with equal external masses mf. c s = momentum transfer; m1,m2,m3 are the internal masses c c p2 ( = pf) c m2 . c c s . . m3 c c m1 . c p1 ( = -pf) c c c c0 = scalar integral, c1p/m = c1+/-, c2p/m= c2+/- c implicit real*8(a-z) complex*16 c0, c1p,c1m,c20,c2p,c2m,cscal,b0s12,b031,b032, & b132,b131,b1s12,c1,c2 xmf=mf*mf c0=cscal(s,mf,m1,m2,m3) call bquer(s,m1,m2,b0s12,b1s12) call bquer(xmf,m3,m1,b031,b131) call bquer(xmf,m3,m2,b032,b132) c c b0jk := b0(xmf,mj,mk), b1jk:= b1(xmf,mf,mk) c b0s12 := b0(s,m1,m2) c c the c1 functions: c ms=4.d0*xmf-s xm1=m1**2 xm2=m2**2 xm3=m3**2 c1=.25d0*dlog(xm3**2/xm1/xm2)+b0s12-.5d0*(b031+b032) & +(xmf+xm3-xm1/2.d0-xm2/2.d0)*c0 c1p=c1/ms c1m=(dlog(xm2/xm1)/2d0+b031-b032+(xm2-xm1)*c0)/2.d0/s c c the c2 functions: c c2=1.d0+b0s12+(xm1+xm2-2.d0*xm3-2.d0*xmf)*c1p+(xm1-xm2)*c1m & +2.d0*xm3*c0 c20=c2/4.d0 c2= (b131+b132+2.d0*b0s12-.5d0)/4.d0 & +(2.d0*xm3-xm1-xm2+2.d0*xmf)/2.d0*c1p-c20 c2p=c2/ms c2= -(b131+b132-.5d0)/4.d0 & -(xm1-xm2)/2.d0*c1m-c20 c2m=c2/s c return end c c******************************************************** c subroutine bquer(s,m1,m2,b0,b1) c c b0 and b1 are the (finite) invariant functions in the c 2-point integrals. c s = q**2; m1,m2 are the internal masses c implicit real*8(a-z) external f,g complex*16 b0,b1,cf xm1=m1**2 xm2=m2**2 lm=dlog(xm2/xm1)/2d0 if ((s.eq.0d0).and.(dabs(dabs(m1)-dabs(m2)).gt.1d-8) $ .and.(dabs(m1*m2).gt.0d0)) then b0 = dcmplx(1d0+(xm1+xm2)/(xm1-xm2)*lm,0d0) b1 = dcmplx(-1d0/4d0-(xm1+xm2)/(xm1-xm2)/4d0- $ xm1/(xm1-xm2)*lm-xm1*xm2/(xm1-xm2)**2*lm,0d0) goto 20 endif cf=dcmplx(f(s,m1,m2),g(s,m1,m2)) if (dabs(m1).eq.dabs(m2)) goto 10 b0=1.d0-(xm2+xm1)/(xm2-xm1)*lm+cf b1=-.25d0+xm1/(xm2-xm1)*lm+(xm2-xm1-s)/2.d0/s*cf goto 20 10 continue b0=cf c achtung: Aenderung am 14.8.94: factor 1/4 dazugezaehlt! b1=-.5d0*cf+1d0/4d0 20 return end c c******************************************************** c subroutine bquer1(x,m1,m2,b0,b1,p0,p1) c c b0 and b1 are the (finite) inariant functions in the c 2-point integrals, p0 and p1 their derivatives. c real parts only, needed for fermion renormalization. c x = q**2; m1,m2 are the internal masses c implicit real*8(a-z) external f cf=f(x,m1,m2) xm1=m1**2 xm2=m2**2 lm=dlog(xm2/xm1)/2d0 if (m1.eq.m2) goto 10 b0=1.d0-(xm2+xm1)/(xm2-xm1)*lm+cf b1=-.25d0+xm1/(xm2-xm1)*lm+(xm2-xm1-x)/2.d0/x*cf goto 20 10 b0=cf b1=-b0/2.d0+1d0/4d0 20 continue c c calculation of the derivatives: c sm=xm1+xm2 dm=xm2-xm1 sm12=(m1+m2)**2 dm12=(m1-m2)**2 s=dsqrt(dabs(sm12-x)) d=dsqrt(dabs(dm12-x)) klam=(dm*dm/(x*x)-sm/x)/s/d anf=-1.d0/x+dm/(x*x)*lm if (x.lt.dm12) goto 30 if (x.gt.sm12) goto 40 fact=2.d0*datan(d/s) goto 41 30 eps=1.d0 fact=dlog(dabs((s+d)/(s-d))) goto 41 40 eps=-1.d0 fact=-dlog(dabs((s+d)/(s-d))) 41 continue deriv=anf-klam*fact p0=deriv b1p=.5d0-lm-2.d0*b1-b0+(xm2-xm1-x)*deriv p1=b1p/2.d0/x return end c c************************************************************** c * c the scalar vertex integral with equal external masses mf * c * c************************************************************** c complex*16 function cscal(s,mf,m1,m2,m3) 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 xmf=mf*mf if (mf.lt.1d-1) then mfstrich = 1d-1 xmf=mfstrich**2 end if c.........xmf etc. are fermion and boson masses squared xm1=m1*m1 xm2=m2*m2 xm3=m3*m3 c c..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 cscal=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******************************************************************* c complex*16 function cspen(x,sig) implicit real*8(a-y) complex*16 z,cpi,spence,zx pi=4d0*datan(1d0) pi6=pi*pi/6.d0 cpi=dcmplx(0.d0,pi) if (x.lt.1.d0) goto 10 if(x.eq.1.d0) goto 11 lx=dlog(x) x1=1.d0-x lx1=dlog(-x1) z=dcmplx(x1,0.d0) if (sig.gt.0.d0) goto 5 cspen=-spence(z)+pi6-lx*(lx1+cpi) goto 20 5 cspen=-spence(z)+pi6-lx*(lx1-cpi) goto 20 10 zx=dcmplx(x,0.d0) cspen=spence(zx) goto 20 11 cspen=dcmplx(pi6,0.d0) 20 return end c c************************************************************** c complex*16 function cln(x,sig) implicit real*8(a-z) complex*16 cpi pi=4d0*datan(1d0) cpi=dcmplx(0.d0,pi) if (x.gt.0.d0) goto 10 x1=-x if (sig.gt.0.d0) goto 5 cln=dlog(x1)-cpi goto 20 5 cln=dlog(x1)+cpi goto 20 10 y=dlog(x) cln=dcmplx(y,0.d0) 20 return end complex*16 function spence(xx) c the dilogarithm for general complex argument. c not allowed: real(xx) gt 1 with imag(xx)=0. implicit real*8(a-z) integer n complex*16 xx,x,z,d,p,cdlog dimension a(19) pi=4d0*datan(1d0) x=xx xr=dreal(x) xi=dimag(x) if(xr.ne.1.) goto 111 if(xi.eq.0.) goto 20 111 continue c projection into the convergence radius vor=1.d0 p=dcmplx(0.d0,0.d0) r=dreal(x) if (r .le. 0.5d0) goto 1 p=pi*pi/6.d0- cdlog(x)*cdlog(1.d0-x) vor=-1.d0 x=1.d0-x 1 continue b=cdabs(x) if (b .lt. 1.d0) goto 2 p=p - (pi*pi/6.d0+ cdlog(-x)*cdlog(-x)/2.d0)*vor vor=vor*(-1.d0) x=1.d0/x 2 continue c calculation of the spence function a(1)=1.d0 a(2)=-0.5d0 a(3)=1.d0/6.d0 a(5)=-1.d0/30.d0 a(7)=1.d0/42.d0 a(9)=-1.d0/30.d0 a(11)=5.d0/66.d0 a(13)=-691.d0/2730.d0 a(15)=7.d0/6.d0 a(17)=-3617.d0/510.d0 a(19)=43867.d0/798.d0 do 5 n=2,9,1 a(2*n)=0.d0 5 continue z=(-1.d0)*cdlog(1.d0-x) d=dcmplx(a(19),0.d0) do 10 n=1,18,1 d=d*z/(20.d0-n) + a(19-n) 10 continue d=d*z spence=d*vor + p goto 30 20 continue spence=pi*pi/6.d0 30 continue return end