* Dget.F * the four-point tensor coefficients * this file is part of LoopTools * last modified 18 Oct 01 th #include "defs.h" integer function Dget(p1, p2, p3, p4, p1p2, p2p3, + m1, m2, m3, m4) implicit none double precision p1, p2, p3, p4, p1p2, p2p3 double precision m1, m2, m3, m4 integer cachelookup external cachelookup, Dijkl double precision para(10) double precision Dcache(2) common /Dbase/ Dcache * The following data statement initializes two *pointers* to NULL * (see cache.c). This is sneaky but ok because IEEE says 0D0 is * represented as 8 0-bytes. data Dcache /0, 0/ para(1) = p1 para(2) = p2 para(3) = p3 para(4) = p4 para(5) = p1p2 para(6) = p2p3 para(7) = m1 para(8) = m2 para(9) = m3 para(10) = m4 Dget = cachelookup(para, Dcache, Dijkl, 10, 46) end ************************************************************************ subroutine Dijkl(P, D) implicit none double precision P(10) double complex D(46),CP(13) integer ier double precision mw,gw,mw2 double complex mwc, mwcbis common /resonant/ mwc integer Cget double complex D0, prova, prova2, prova3 external Cget, D0 double complex Ccache(1) common /Cbase/ Ccache XREAL f1, f2, f3, det3 XREAL M11, M12, M13, M22, M23, M33 double complex s1, s2, s3 double complex c1, c2, c3, c4, c5 integer C234, C134, C124, C123 double complex D0resonant,D0_ir external D0resonant, D0_ir #ifndef LANGUAGE_FORTRAN_90 double complex lc1, lc2, lc3 lc1() = dcmplx( + M11*dble(s1) + M12*dble(s2) + M13*dble(s3), + M11*dimag(s1) + M12*dimag(s2) + M13*dimag(s3)) lc2() = dcmplx( + M12*dble(s1) + M22*dble(s2) + M23*dble(s3), + M12*dimag(s1) + M22*dimag(s2) + M23*dimag(s3)) lc3() = dcmplx( + M13*dble(s1) + M23*dble(s2) + M33*dble(s3), + M13*dimag(s1) + M23*dimag(s2) + M33*dimag(s3)) #endif mw = dreal(mwc) gw = -dimag(mwc)*2d0 mw2=mw*mw c mw = 80.35d0 c gw = 2.08699d0 f1 = .5D0*(QEXT(P(1)) + QEXT(P(5)) - P(2)) f2 = .5D0*(QEXT(P(1)) + QEXT(P(4)) - P(6)) f3 = .5D0*(QEXT(P(5)) + QEXT(P(4)) - P(3)) M22 = QEXT(P(1))*P(4) - f2*f2 M23 = f1*f2 - P(1)*f3 det3 = 2*(P(5)*M22 - P(4)*f1*f1 + f3*(f1*f2 + M23)) M11 = (QEXT(P(4))*P(5) - f3*f3)/det3 M12 = (f2*f3 - P(4)*f1)/det3 M13 = (f1*f3 - P(5)*f2)/det3 M22 = M22/det3 M23 = M23/det3 M33 = (QEXT(P(1))*P(5) - f1*f1)/det3 f1 = QEXT(P(1)) + QEXT(P(7)) - P(8) f2 = QEXT(P(5)) + QEXT(P(7)) - P(9) f3 = QEXT(P(4)) + QEXT(P(7)) - P(10) C234 = Cget(P(2), P(3), P(6), P(8), P(9), P(10)) C134 = Cget(P(5), P(3), P(4), P(7), P(9), P(10)) C124 = Cget(P(1), P(6), P(4), P(7), P(8), P(10)) C123 = Cget(P(1), P(2), P(5), P(7), P(8), P(9)) c1 = Ccache(C234 + cc0) + Ccache(C234 + cc1) + + Ccache(C234 + cc2) c2 = Ccache(C234 + cc12) + Ccache(C234 + cc22) + + Ccache(C234 + cc2) c3 = 2*(Ccache(C234 + cc1) + c2) + Ccache(C234 + cc11) - + Ccache(C234 + cc22) + Ccache(C234 + cc0) c4 = Ccache(C234 + cc1) + Ccache(C234 + cc11) + + Ccache(C234 + cc12) c5 = Ccache(C234 + cc111) + Ccache(C234 + cc112) + + Ccache(C234 + cc11) if (P(1).eq.P(7).and.P(6).eq.P(10).and.P(8).eq.0d0.and. + P(9).eq.mw2) then CP(1) = dcmplx(P(7),0d0) CP(2) = dcmplx(P(8),0d0) CP(3) = dcmplx(P(9),gw ) CP(4) = dcmplx(P(10),0d0) CP(5) = dcmplx(P(1),0d0) CP(6) = dcmplx(P(2),0d0) CP(7) = dcmplx(P(3),0d0) CP(8) = dcmplx(P(4),0d0) CP(9) = dcmplx(P(5),0d0) CP(10)= dcmplx(P(6),0d0) CP(11)= dcmplx(0d0 ,0d0) CP(12)= dcmplx(0d0 ,0d0) CP(13)= dcmplx(0d0 ,0d0) ier = 0 c call ffcd0(prova,CP,ier) c mwcbis = dcmplx(mw,-0d0) mwcbis = mwc ! 100 GeV^2 is for sure above any threshold! if (P(2).lt.100d0) mwcbis=dcmplx(dreal(mwc),0d0) prova = D0_ir(CP(5),CP(7),CP(9),CP(10),CP(6),CP(8),mwcbis) mwcbis = mwc c prova2 = D0resonant(P(2),P(4),P(1),P(6),mw,gw) c prova3 = D0(P(1), P(2), P(3), P(4), P(5), P(6), c + P(7), P(8), P(9), P(10)) D(dd0) = prova c print*,'gw=',mwcbis c print*,'D0 compl=',prova c print*,'D0 reson=',prova2 c print*,'D0 reale=',prova3 else if (P(5).eq.P(7).and.P(3).eq.P(10).and.P(9).eq.0d0.and. + P(8).eq.mw2) then c print*,'infattisi 2' CP(1) = dcmplx(P(7),0d0) CP(2) = dcmplx(P(8),gw) CP(3) = dcmplx(P(9),0d0 ) CP(4) = dcmplx(P(10),0d0) CP(5) = dcmplx(P(1),0d0) CP(6) = dcmplx(P(2),0d0) CP(7) = dcmplx(P(3),0d0) CP(8) = dcmplx(P(4),0d0) CP(9) = dcmplx(P(5),0d0) CP(10)= dcmplx(P(6),0d0) CP(11)= dcmplx(0d0 ,0d0) CP(12)= dcmplx(0d0 ,0d0) CP(13)= dcmplx(0d0 ,0d0) ier = 0 c mwcbis = dcmplx(mw,-0d0) mwcbis = mwc if (P(2).lt.100d0) mwcbis=dcmplx(dreal(mwc),0d0) prova = D0_ir(CP(5),CP(10),CP(9),CP(7),CP(6),CP(8),mwcbis ) mwcbis = mwc c prova2 = D0resonant(P(2),P(4),P(3),P(5),mw,gw) c prova3 = D0(P(1), P(2), P(3), P(4), P(5), P(6), c + P(7), P(8), P(9), P(10)) c print*,'gw=',gw c print*,'D0 compl=',prova c print*,'D0 reson=',prova2 c print*,'D0 reale=',prova3 D(dd0) = prova else D(dd0) = D0(P(1), P(2), P(3), P(4), P(5), P(6), + P(7), P(8), P(9), P(10)) * print*,'bos vero =',D(dd0) end if s1 = Ccache(C134 + cc0) - Ccache(C234 + cc0) - f1*D(dd0) s2 = Ccache(C124 + cc0) - Ccache(C234 + cc0) - f2*D(dd0) s3 = Ccache(C123 + cc0) - Ccache(C234 + cc0) - f3*D(dd0) c print*,'inizio Ccache' c print*,Ccache(C134+cc0) c print*,Ccache(C124+cc0) c print*,Ccache(C123+cc0) c print*,Ccache(C234+cc0) c print*,'fine Ccache' c print*,'s1=',s1 c print*,'s2=',s2 c print*,'s3=',s3 D(dd1) = lc1() D(dd2) = lc2() D(dd3) = lc3() D(dd00) = P(7)*D(dd0) + .5D0*(Ccache(C234 + cc0) + + D(dd1)*f1 + D(dd2)*f2 + D(dd3)*f3) s1 = c1 - f1*D(dd1) - 2*D(dd00) s2 = Ccache(C124 + cc1) + c1 - f2*D(dd1) s3 = Ccache(C123 + cc1) + c1 - f3*D(dd1) D(dd11) = lc1() s1 = Ccache(C134 + cc1) - Ccache(C234 + cc1) - f1*D(dd2) s2 = -(Ccache(C234 + cc1) + f2*D(dd2)) - 2*D(dd00) s3 = Ccache(C123 + cc2) - Ccache(C234 + cc1) - f3*D(dd2) D(dd12) = lc1() D(dd22) = lc2() s1 = Ccache(C134 + cc2) - Ccache(C234 + cc2) - f1*D(dd3) s2 = Ccache(C124 + cc2) - Ccache(C234 + cc2) - f2*D(dd3) s3 = -(Ccache(C234 + cc2) + f3*D(dd3)) - 2*D(dd00) D(dd13) = lc1() D(dd23) = lc2() D(dd33) = lc3() s1 = Ccache(C134 + cc00) - Ccache(C234 + cc00) - f1*D(dd00) s2 = Ccache(C124 + cc00) - Ccache(C234 + cc00) - f2*D(dd00) s3 = Ccache(C123 + cc00) - Ccache(C234 + cc00) - f3*D(dd00) D(dd001) = lc1() D(dd002) = lc2() D(dd003) = lc3() s1 = -f1*D(dd11) - c3 - 4*D(dd001) s2 = Ccache(C124 + cc11) - f2*D(dd11) - c3 s3 = Ccache(C123 + cc11) - f3*D(dd11) - c3 D(dd111) = lc1() D(dd112) = lc2() D(dd113) = lc3() s1 = Ccache(C134 + cc11) - Ccache(C234 + cc11) - f1*D(dd22) s2 = -(Ccache(C234 + cc11) + f2*D(dd22)) - 4*D(dd002) s3 = Ccache(C123 + cc22) - Ccache(C234 + cc11) - f3*D(dd22) D(dd122) = lc1() D(dd222) = lc2() D(dd223) = lc3() s1 = Ccache(C134 + cc22) - Ccache(C234 + cc22) - f1*D(dd33) s2 = Ccache(C124 + cc22) - Ccache(C234 + cc22) - f2*D(dd33) s3 = -(Ccache(C234 + cc22) + f3*D(dd33)) - 4*D(dd003) D(dd133) = lc1() D(dd233) = lc2() D(dd333) = lc3() s1 = c2 - f1*D(dd13) - 2*D(dd003) s2 = Ccache(C124 + cc12) + c2 - f2*D(dd13) s3 = c2 - f3*D(dd13) - 2*D(dd001) D(dd123) = lc2() D(dd0000) = (P(7)*D(dd00) + .5D0*(Ccache(C234 + cc00) + + f1*D(dd001) + f2*D(dd002) + f3*D(dd003)))/3D0 + 1/36D0 D(dd0011) = (P(7)*D(dd11) + .5D0*(c3 + f1*D(dd111) + + f2*D(dd112) + f3*D(dd113)))/3D0 D(dd0012) = (P(7)*D(dd12) - .5D0*(c4 - f1*D(dd112) - + f2*D(dd122) - f3*D(dd123)))/3D0 D(dd0013) = (P(7)*D(dd13) - .5D0*(c2 - f1*D(dd113) - + f2*D(dd123) - f3*D(dd133)))/3D0 D(dd0022) = (P(7)*D(dd22) + .5D0*(Ccache(C234 + cc11)+ + f1*D(dd122) + f2*D(dd222) + f3*D(dd223)))/3D0 D(dd0023) = (P(7)*D(dd23) + .5D0*(Ccache(C234 + cc12)+ + f1*D(dd123) + f2*D(dd223) + f3*D(dd233)))/3D0 D(dd0033) = (P(7)*D(dd33) + .5D0*(Ccache(C234 + cc22)+ + f1*D(dd133) + f2*D(dd233) + f3*D(dd333)))/3D0 c3 = c2 + Ccache(C234 + cc112) + Ccache(C234 + cc122) c2 = c2 + c3 - Ccache(C234 + cc2) + Ccache(C234 + cc122) + + Ccache(C234 + cc222) c4 = c4 + c1 + c2 + c3 + c4 + c5 + Ccache(C234 + cc12) s1 = c4 - f1*D(dd111) - 6*D(dd0011) s2 = Ccache(C124 + cc111) + c4 - f2*D(dd111) s3 = Ccache(C123 + cc111) + c4 - f3*D(dd111) D(dd1111) = lc1() D(dd1112) = lc2() D(dd1113) = lc3() s1 = -(c2 + f1*D(dd113)) - 4*D(dd0013) s2 = Ccache(C124 + cc112) - c2 - f2*D(dd113) s3 = -(c2 + f3*D(dd113)) - 2*D(dd0011) D(dd1123) = lc2() D(dd1133) = lc3() s1 = c5 - f1*D(dd122) - 2*D(dd0022) s2 = c5 - f2*D(dd122) - 4*D(dd0012) s3 = Ccache(C123 + cc122) + c5 - f3*D(dd122) D(dd1122) = lc1() D(dd1223) = lc3() s1 = Ccache(C134 + cc111) - Ccache(C234 + cc111) - f1*D(dd222) s2 = -(Ccache(C234 + cc111) + f2*D(dd222)) - 6*D(dd0022) s3 = Ccache(C123 + cc222) - Ccache(C234 + cc111) - f3*D(dd222) D(dd1222) = lc1() D(dd2222) = lc2() D(dd2223) = lc3() s1 = Ccache(C134 + cc122) - Ccache(C234 + cc122) - f1*D(dd233) s2 = -(Ccache(C234 + cc122) + f2*D(dd233)) - 2*D(dd0033) s3 = -(Ccache(C234 + cc122) + f3*D(dd233)) - 4*D(dd0023) D(dd1233) = lc1() D(dd2233) = lc2() s1 = Ccache(C134 + cc222) - Ccache(C234 + cc222) - f1*D(dd333) s2 = Ccache(C124 + cc222) - Ccache(C234 + cc222) - f2*D(dd333) s3 = -(Ccache(C234 + cc222) + f3*D(dd333)) - 6*D(dd0033) D(dd1333) = lc1() D(dd2333) = lc2() D(dd3333) = lc3() #ifdef LANGUAGE_FORTRAN_90 contains double complex function lc1() lc1 = dcmplx( + M11*dble(s1) + M12*dble(s2) + M13*dble(s3), + M11*dimag(s1) + M12*dimag(s2) + M13*dimag(s3)) end function lc1 double complex function lc2() lc2 = dcmplx( + M12*dble(s1) + M22*dble(s2) + M23*dble(s3), + M12*dimag(s1) + M22*dimag(s2) + M23*dimag(s3)) end function lc2 double complex function lc3() lc3 = dcmplx( + M13*dble(s1) + M23*dble(s2) + M33*dble(s3), + M13*dimag(s1) + M23*dimag(s2) + M33*dimag(s3)) end function lc3 #endif end double complex function D0resonant(s,r,m1,m2,mw,gw) implicit none double precision s,r,m1,m2 double precision MW,GW,lambda2,pi2,getlambda double complex cli2,dummy integer ier,ipi12 parameter (pi2=3.1415926535897932384626433832795029d0**2) call ffzxdl(cli2,ipi12,dummy,(mw**2+r)/r,1,ier) lambda2 = getlambda() D0resonant = 1d0/(r*(mw*mw-dcmplx(0d0,0d0*gw*mw)-s))* - (log(mw/sqrt(m1))**2+log(mw/sqrt(m2))**2+ - 2*log(-sqrt(m1)*sqrt(m2)/r)* - log((mw*mw-dcmplx(0d0,gw*mw)-s)/ - mw/sqrt(lambda2))+ - pi2/3d0+dble(cli2)+ipi12*pi2/12d0) end function D0resonant !*********************************************************************** function D0_ir(cp1, cp2, cp3, cp4, cp1p2, cp2p3, m3c) implicit none double complex D0_ir double complex cp1, cp2, cp3, cp4, cp1p2, cp2p3,gzmz double precision p1, p2, p3, p4, p1p2, p2p3, m3 double precision d, delta double complex m3c, cm3, m1_, m4_, ieps double complex xs, x2, x3, y, c, f double complex logxs, logx2, logx3, log1x2, log1x3, logy double precision mgamma2, pi, getlambda double complex ln, spence, bdK, addeps, cln external ln, spence, bdK, addeps, cln character*1 boson common/vectorboson/boson integer ifirst,k common/d0irifirst/ifirst data ifirst /1/ if (ifirst.eq.1) then do k = 1,50 c print*,'ATTENTION IN THE DO_ir !!' enddo ifirst = 0 endif pi = 3.1415926535897932384626433832795029d0 ieps = dcmplx(0d0,1d-20) mgamma2 = getlambda() delta = dble(mgamma2) m3 = dble(m3c) cm3 = m3**2 gzmz = dcmplx(0.d0,-2.d0*m3*imag(m3c)) c print*,-2.d0*imag(m3c),gzmz c print*,m3c,m3,cm3 p1 = dble(cp1) p2 = dble(cp2) p3 = dble(cp3) p4 = dble(cp4) p1p2 = dble(cp1p2) p2p3 = dble(cp2p3) m1_ = sqrt(p1) m4_ = sqrt(p4) c print*,'m1m4=',p1,p2,p3,p4,p1p2,p2p3 d = p2p3 - (m1_ - m4_)**2 f = .5D0/m1_/m4_/(p1p2 - cm3) cccATTENZIONE!! if (boson.eq.'Z') f = .5D0/m1_/m4_/(p1p2 - cm3+gzmz) if(abs(d).gt.0d0) then xs = bdK(p2p3, m1_, m4_) logxs = log(xs) f = f*2D0*xs/(1D0 - xs**2) endif ! massless case if (m3.eq.(0.D0,0.D0)) then !write(*,*) ' massless' if(p1.eq.p2 .and. p3.eq.p4) then D0_ir = 2D0*f*ln(-delta/p1p2, 1D0) if(abs(d).gt.0d0) D0_ir = -logxs*D0_ir return endif if(p1.eq.p2) then !write(*,*) 'p1=p2',1D0/(p2p3-p4)/p1p2 !write(*,*) 'p1=p2',-m1_*m4_/(p2p3-p4) !write(*,*) 'p1=p2',cln(-m1_*m4_/(p2p3-p4),1D0) !write(*,*) 'p1=p2',2D0*ln((p3-p4)/p1p2,1D0) D0_ir = 1D0/(p2p3-p4)/p1p2* + (cln(-m1_*m4_/(p2p3-p4),1D0)* + (cln(-m1_*m4_/(p2p3-p4),1D0)+log(delta/p4)+ + 2D0*ln((p3-p4)/p1p2,1D0))+ + pi**2/6D0) print*,'p1eqp2' return endif if(p3.eq.p4) then print*, 'p3=p4' D0_ir = 1D0/(p2p3-p1)/p1p2* + (cln(-m4_*m1_/(p2p3-p1),1D0)* + (cln(-m4_*m1_/(p2p3-p1),1D0)+log(delta/p1)+ + 2D0*ln((p2-p1)/p1p2,1D0))+ + pi**2/6D0) return endif y = m1_/m4_*(p3 - p4 + IEPS)/ (p2 - p1 + IEPS) logy = log(y) c = cln(delta/m1_/m4_, 0D0) + + ln((p2 - p1)/p1p2, p1 - p2) + ln((p3 - p4)/p1p2, p4 - p3) if(abs(d).gt.0d0) then D0_ir = f*(pi**2/6D0 +logxs*(-.5D0*logxs+2D0*log(1D0-xs**2)-c)+ + spence(xs**2, 0D0) + .5D0*logy**2 -spence(xs/y, 0D0) - + (logxs + log(1D0/y))*log(1D0 - xs/y) -spence(xs*y, 0D0) - + (logxs + logy)*log(1D0 - xs*y)) print*, ' massless',D0_ir return endif D0_ir = f*(c - 2D0 - (1D0 + y)/(1D0 - y)*logy) return endif ! massive case c print*, ' massive' x2 = bdK(p2, m1_, dcmplx(m3)) x3 = bdK(p3, m4_, dcmplx(m3)) logx2 = log(x2) logx3 = log(x3) log1x3 = log(1D0/x3) c c = cln(m3*sqrt(delta)/(m3c**2 - p1p2), 1D0) c = log(m3*sqrt(delta)/(m3c**2 - p1p2)) if(abs(d).gt.0d0) then log1x2 = log(1D0/x2) D0_ir = f*(.5D0*pi**2 +2D0*log(xs)*(log(1D0 - xs**2) - c) + + spence(xs**2, 0D0)+logx2**2 + logx3**2 -spence(xs/x2/x3, 0D0) - + (logxs+log1x2+log1x3)*log(1D0-xs/x2/x3) -spence(xs*x2/x3, 0D0) - + (logxs + logx2 +log1x3)*log(1D0-xs*x2/x3)-spence(xs/x2*x3, 0D0) - + (logxs+log1x2+logx3)*log(1D0 - xs/x2*x3) -spence(xs*x2*x3, 0D0) - + (logxs + logx2 + logx3)*log(1D0 - xs*x2*x3)) return endif D0_ir = f*(2D0*c - (1D0 + x2/x3)/(1D0 - x2/x3)*(logx2 + log1x3) - + (1D0 + x2*x3)/(1D0 - x2*x3)*(logx2 + logx3) - 2D0) end function D0_ir !********************************************************************** function addeps(k) implicit none double complex addeps double complex k addeps = k*dcmplx(1D0, -sign(1d-8, dreal(k))) end function addeps !********************************************************************* subroutine roots(p, m1, m2, x1, x2, y1, y2, r) implicit none double precision p, m1, m2 double complex x1, x2, y1, y2, r double precision q,eps eps = 1d-20 r = sqrt(dcmplx(p*(p - 2D0*(m1 + m2)) + (m1 - m2)**2)) q = p + m1 - m2 x1 = (q + r)/2D0/p x2 = (q - r)/2D0/p if(abs(x2) > abs(x1)) then x1 = m1/p/x2 else if(abs(x1) > abs(x2)) then x2 = m1/p/x1 endif x1 = x1 + dcmplx(0D0, abs(p*x1)/p*EPS) x2 = x2 + dcmplx(0D0, -abs(p*x2)/p*EPS) q = p - m1 + m2 y2 = (q + r)/2D0/p y1 = (q - r)/2D0/p if(abs(y2) > abs(y1)) then y1 = m2/p/y2 else if(abs(y1) > abs(y2)) then y2 = m2/p/y1 endif y1 = y1 + dcmplx(0D0, -abs(p*y1)/p*EPS) y2 = y2 + dcmplx(0D0, abs(p*y2)/p*EPS) end subroutine roots !******************************************************************* function fpv(n, x, y) implicit none double complex fpv integer n double complex x, y integer m double complex xm double precision calacc calacc = 1d-12 if(abs(x) < 10D0) then if(n.eq.0) then fpv = -log(-y/x) else if(x.eq.dcmplx(0D0,0d0)) then fpv = -1D0/n else fpv = 0D0 xm = 1D0 do m = 0, n - 1 fpv = fpv - xm/(n - m) xm = xm*x enddo fpv = fpv - xm*log(-y/x) endif else fpv = 0D0 xm = 1D0 do m = 1, 30 xm = xm/x fpv = fpv + xm/(m + n) if(abs(xm/fpv) < CALACC**2) return enddo endif end function fpv !******************************************************************** function yfpv(n, x, y) implicit none double complex yfpv integer n double complex x, y double complex fpv external fpv if(abs(y).eq.0D0) then yfpv = 0D0 else yfpv = y*fpv(n, x, y) endif end function yfpv !********************************************************************* function xlogx(x) implicit none double complex xlogx double complex x if(abs(x).eq.0D0) then xlogx = 0D0 else xlogx = x*log(x) endif end function xlogx !*********************************************************************** function ln(x, isig) implicit none double complex ln double precision x, isig, pi pi = 3.1415926535897932384626433832795029d0 if(x > 0D0) then ln = log(x) else ln = log(-x) +dcmplx(0D0, sign(pi, isig)) endif end function ln !********************************************************************* function cln(z, isig) implicit none double complex cln double complex z double precision isig, pi pi = 3.1415926535897932384626433832795029d0 if(dimag(z).eq.0D0 .and. dble(z).lt.0D0) then #ifdef WARNINGS if(isig.eq.0D0) print *, "cln: argument on cut" #endif cln = log(-z) +dcmplx(0D0, sign(pi, isig)) else cln = log(z) endif end function cln !********************************************************************* function spence(z, isig) implicit none double complex spence double complex z double precision isig,pi double complex z1, ieps double precision az1 double complex li2series, cln external li2series, cln pi = 3.1415926535897932384626433832795029d0 ieps = dcmplx(0d0,1d-20) z1 = 1D0 - z az1 = abs(z1) #ifdef WARNINGS if(isig.eq.0D0 .and.dimag(z).eq.0D0 .and. abs(dble(z1)) < + CALACC) print *, "spence: argument on cut" #endif if(az1 < 1D-15) then spence = pi**2/6D0 else if(dble(z) < .5D0) then if(abs(z) < 1D0) then spence = li2series(z, isig) else spence = -pi**2/6D0 -.5D0*cln(-z, -isig)**2-li2series(1D0/z,-isig) endif else if(az1 < 1D0) then spence = pi**2/6D0 -cln(z, isig)*cln(z1, -isig)-li2series(z1,-isig) else spence = pi**2/3D0 + + .5D0*cln(-z1, -isig)**2 - cln(z, isig)*cln(z1, -isig) + + li2series(1D0/z1, isig) endif endif end function spence !********************************************************************* function li2series(z, isig) implicit none double complex li2series double complex z double precision isig double complex xm, x2, new integer j double complex cln external cln ! these are the even-n Bernoulli numbers, already divided by (n + 1)! ! as in Table[BernoulliB[n]/(n + 1)!, {n, 2, 50, 2}] double precision b(25) data b / 0.02777777777777777777777777777777777777777778774D0, & -0.000277777777777777777777777777777777777777777778D0, & 4.72411186696900982615268329554043839758125472D-6, & -9.18577307466196355085243974132863021751910641D-8, & 1.89788699889709990720091730192740293750394761D-9, & -4.06476164514422552680590938629196667454705711D-11, & 8.92169102045645255521798731675274885151428361D-13, & -1.993929586072107568723644347793789705630694749D-14, & 4.51898002961991819165047655285559322839681901D-16, & -1.035651761218124701448341154221865666596091238D-17, & 2.39521862102618674574028374300098038167894899D-19, & -5.58178587432500933628307450562541990556705462D-21, & 1.309150755418321285812307399186592301749849833D-22, & -3.087419802426740293242279764866462431595565203D-24, & 7.31597565270220342035790560925214859103339899D-26, & -1.740845657234000740989055147759702545340841422D-27, & 4.15763564461389971961789962077522667348825413D-29, & -9.96214848828462210319400670245583884985485196D-31, & 2.394034424896165300521167987893749562934279156D-32, & -5.76834735536739008429179316187765424407233225D-34, & 1.393179479647007977827886603911548331732410612D-35, & -3.372121965485089470468473635254930958979742891D-37, & 8.17820877756210262176477721487283426787618937D-39, & -1.987010831152385925564820669234786567541858996D-40, & 4.83577851804055089628705937311537820769430091D-42 / xm = -cln(1D0 - z, -isig) x2 = xm**2 li2series = xm - x2/4D0 do j = 1, 25 xm = xm*x2 new = li2series + xm*b(j) if(new == li2series) return li2series = new enddo #ifdef WARNINGS print *, "li2series: bad convergence" #endif end function li2series !********************************************************************* function contspence(z1, z2) implicit none double complex contspence, z1, z2, z12, ris double complex spence external spence z12 = z1*z2 !ris = spence(1D0-z12,0D0)+log(1D0-z12)*(log(z12)-Log(z1)-log(z2)) !write(*,*) 'contspence=',ris contspence=spence(1D0-z12,0D0)+log(1D0-z12)*(log(z12)-Log(z1)-log(z2)) end function contspence !********************************************************************** function cspence(z1, z2, im1, im2) implicit none double complex cspence double complex z1, z2 double precision im1, im2, pi double complex cln, spence, ieps integer eta external cln, spence, eta double complex z12 double precision im12 integer etas pi = 3.1415926535897932384626433832795029d0 ieps = dcmplx(0d0,1d-20) if (cdabs(z1)<1D-15.or.cdabs(z2)<1D-15) then cspence = pi**2/6D0 return end if z12 = z1*z2 im12 = im2*sign(1D0, dble(z1)) if(dble(z12) > .5D0) then cspence = spence(1D0 - z12, 0D0) etas = eta(z1, z2, im1, im2, im12) if(abs(etas).gt.0) cspence = cspence +etas*dcmplx(0D0, 2D0*pi)* + cln(1D0 - z12, -im12) else if(abs(z12) < 1D-4) then cspence = pi**2/6D0 -spence(z12, 0D0) + (cln(z1, im1) + + cln(z2, im2))*z12* + (1D0 + z12*(.5D0 + z12*(1D0/3D0 + z12/4D0))) else cspence = pi**2/6D0 -spence(z12, 0D0) - + (cln(z1, im1) + cln(z2, im2))*cln(1D0 - z12, 0D0) endif end function cspence !********************************************************************* function eta(c1, c2, im1x, im2x, im12x) implicit none integer eta double complex c1, c2 double precision im1x, im2x, im12x double precision im1, im2, im12 im1 = dimag(c1) if(im1 .eq. 0D0) im1 = im1x im2 = dimag(c2) if(im2 .eq. 0D0) im2 = im2x im12 = dimag(c1*c2) if(im12 .eq. 0D0) im12 = im12x if(im1 < 0D0 .and. im2 < 0D0 .and. im12 > 0D0) then eta = 1 else if(im1 > 0D0 .and. im2 > 0D0 .and. im12 < 0D0) then eta = -1 else eta = 0 #ifdef WARNINGS if(.not. (im2 .eq. 0D0 .and. dble(c2) > 0D0 .or. + im1 .eq. 0D0 .and. dble(c1) > 0D0) .and. + (im1 .eq. 0D0 .and. dble(c1) < 0D0 .or. + im2 .eq. 0D0 .and. dble(c2) < 0D0 .or. + im12 .eq. 0D0 .and. dble(c1*c2) < 0D0) + ) print *, 'eta not defined' #endif endif end function eta !******************************************************************** function eta_n(c1, c2) implicit none integer eta_n double complex c1, c2 integer eta external eta eta_n = eta(c1, c2, 0D0, 0D0, 0D0) end function eta_n !********************************************************************* function eta_tilde(c1, c2, im1x, im2x) implicit none integer eta_tilde double complex c1, c2 double precision im1x, im2x double precision im1, im2 integer eta external eta im1 = dimag(c1) if(im1 .eq. 0D0) im1 = im1x im2 = dimag(c2) if(abs(im2).gt.0D0) then eta_tilde = eta(c1, c2, im1x, 0D0, 0D0) else if(dble(c2) > 0D0) then eta_tilde = 0 else if(im1 > 0D0 .and. im2x > 0D0) then eta_tilde = -1 else if(im1 < 0D0 .and. im2x < 0D0) then eta_tilde = 1 else eta_tilde = 0 #ifdef WARNINGS if(im1 .eq. 0D0 .and. dble(c1) < 0D0 .or. + im2x .eq. 0D0 .and. dble(c1*c2) < 0D0 + ) print *, 'eta_tilde not defined' #endif endif end function eta_tilde !********************************************************************** function bdK(x, m1, m2) ! this is actually -K from the Beenakker/Denner paper for D0_ir implicit none double complex bdK double complex m1, m2, d double precision x double complex zz,ieps ieps = dcmplx(0d0,1d-20) d = x - (m1 - m2)**2 !write(*,*) 'bdk=',x,m1,m2,d if(d .eq. dcmplx(0D0,0D0)) then bdK = dcmplx(1D0,0D0) !write(*,*) 'bdK1=',bdK else zz = sqrt(1D0 - 4D0*m1*m2/(d + IEPS)) bdK = (zz - 1D0)/(zz + 1D0) c print*, 'bdK2=',zz,sqrt(1D0 - 4D0*m1*m2/(d - IEPS)) endif end function bdK