* block data init implicit real*8 (a-h,o-z) * common/param/eps,pi,pis,qdelta common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 common/fmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy common/fmass2/em2,rmm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * *-----Official set of data input * data eps/1.d-37/,qdelta/0.d0/ data gf/8.2476172696d-6/,zm/91.1867d0/, # em/0.51099906d-3/,pi/3.141592653589793238462643d0/, # zg/2.4948d0/,tm/1.7771d0/, # bqm/4.7d0/,cqm/1.55d0/,alphai/137.0359895d0/ data rmm/0.10565839d0/,uqm/0.041d0/,dqm/0.041d0/,sqm/0.15d0/ data rnm/1.d-10/,dmy/0.d0/ * end * *---------------------------------------------------------------------- * implicit real*8 (a-h,o-z) character*1 op,opas * parameter(n=2,nt=1,lwa=(n*(3*n+13))/2,lwat=(nt*(3*nt+13))/2) parameter(nout=99) * common/opt/op common/g2/g2z common/ohm/ohm common/mix/alst,alstz common/param/eps,pi,pis,qdelta common/iaz/alzr,alzi,alzir,alzii common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 common/fmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy common/fmass2/em2,rmm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * dimension zs(2),zl(2) dimension p3qz(2),fzz(2) dimension pggfz(2),pgglqz(2),pggbz(2) dimension st(nt),fvect(nt),wat(lwat) dimension fvecw(n),waw(lwa),waz(lwa),sw(n),fvecz(n),sz(n) * external c05nbf,c02ajf,x02ajf,fsw,fsz,fsh,f06ejf * read'(a)',opas if(opas.eq.'y') then read*,wm,tqm,als else if(opas.eq.'n') then read*,wm,tqm als= 1.d-6 endif * op= 'p' pis= pi*pi em2= em*em rmm2= rmm*rmm tm2= tm*tm rnm2= rnm*rnm uqm2= uqm*uqm dqm2= dqm*dqm cqm2= cqm*cqm sqm2= sqm*sqm bqm2= bqm*bqm dmy2= dmy*dmy * zm2= zm*zm zm3= zm2*zm wm2= wm*wm * wmc= wm2*wm wg= 1.5d0*gf*wmc/pi fwr= wm2 fwi= -wm*wg * tqm2= tqm*tqm * cw2= wm2/zm2 sw2= 1.d0-cw2 * *-----QCD is initialized * if(opas.eq.'y') then nf= 5 if(tqm.gt.zm) then rsqcd= tqm else rsqcd= zm endif alst= tralphas(zm,rsqcd,als,nf) else alst= als endif alstz= als * *-----m_H is fixed * gmh= 100.d0 st(1)= gmh*gmh tol= sqrt(x02ajf()) ifail= 1 call c05nbf(fsh,nt,st,fvect,tol,wat,lwat,ifail) print 1,ifail bhm2= st(1) bhm= sqrt(bhm2) ohm= bhm * print*,'-------------------------------------------------------' print 20,wm,bhm,tqm,als 20 format(' M_W (GeV) = ',e20.5,/ # ' m_H (GeV) = ',e20.5,/ # ' m_t (GeV) = ',e20.5,/ # ' alpha_s = ',e20.5,//) print*,'-------------------------------------------------------' * *-----alpha(zm) * pzr= -zm2 pzi= 0.d0 call pself(pzr,pzi,pggfz,pggbz,pgglqz,pggnpz) call vbself0(pggf0,pggb0,fw0) x= alphai-0.25d0*(pggfz(1)+pggbz(1)-pggf0-pggb0+pggnpz)/pi y= -0.25d0*(pggfz(2)+pggbz(2)+pgglqz(2))/pi alzir= x alzii= y alzr= x/(x*x+y*y) alzi= -y/(x*x+y*y) print 3,x,y 3 format(/' 1/alpha(zm) ',2e20.7,/) * *-----Re(1/g^2(zm)) * apis= 16.d0*pis call vbself0(pggf0,pggb0,fw0) p2zr= -zm2 p2zi= 0.d0 call vbself(p2zr,p2zi,p3qz,fzz) bx= 1.d0/8.d0/gf/zm2+(fw0-fzz(1))/zm2/apis xih= -p3qz(2)/apis ac= 4.d0*pi*alzr bc= -1.d0-8.d0*pi*xih*alzi cc= -4.d0*pi*xih*xih*alzr+bx ifail= 0 call c02ajf(ac,bc,cc,zs,zl,ifail) g2z= zs(1) * *-----Iteration for sw starts * print*,' Iteration for sw starts ' print*,' ' dswi= 1.d-1 twm= wm+0.01d0 twm2= twm*twm swi0= twm*wg*(1.d0+2.d0/3.d0*als/pi) it= 1 11 print 6,it sw(1)= twm2 sw(2)= -swi0*(1.d0-(it-1)*dswi) tol= sqrt(x02ajf()) ifail= 1 call c05nbf(fsw,n,sw,fvecw,tol,waw,lwa,ifail) print 1,ifail if(ifail.eq.0) then wmu= sqrt(sw(1)) wlg= -sw(2)/wmu print 2,wm,wg print 7,wmu,wlg else if(it.lt.20) then it= it+1 go to 11 else print*,' Failure' endif endif print*,' ' print*,'---------------------------------------------------------' print*,' ' * *-----Iteration for sz starts * print*,' Iteration for sz starts ' print*,' ' dszi= 1.d-1 szi0= zm*zg it= 1 10 print 6,it sz(1)= zm2 sz(2)= -szi0*(1.d0-(it-1)*dszi) tol= sqrt(x02ajf()) ifail= 1 call c05nbf(fsz,n,sz,fvecz,tol,waz,lwa,ifail) print 1,ifail if(ifail.eq.0) then zmu= sqrt(sz(1)) zlg= -sz(2)/zmu print 2,zm,zg print 7,zmu,zlg fnorm= f06ejf(n,fvecz,1) else if(it.lt.20) then it= it+1 go to 10 else print*,' Failure' endif endif * 1 format(/' On exit IFAIL = ',i2,/) 2 format(/' Input ',2e20.5,//) 7 format(/' Output ',2e20.5,//) 6 format(' Iteration = ',i3,/) * stop end * *-------------------------------------------------------------------------- * subroutine fsh(nt,st,fvect,iflag) implicit real*8(a-h,o-z) * common/param/eps,pi,pis,qdelta common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 common/fmass2/em2,rmm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * dimension zs(2),zl(2) dimension st(nt),fvect(nt) dimension pggfz(2),pgglqz(2),p3qz(2),fzz(2),p3qw(2),fw(2), # fzw(2),pggbz(2) * bhm2= st(1) bhm= sqrt(bhm2) * *-----alpha(zm) * pzr= -zm2 pzi= 0.d0 call pself(pzr,pzi,pggfz,pggbz,pgglqz,pggnpz) call vbself0(pggf0,pggb0,fw0) x= alphai-0.25d0*(pggfz(1)+pggbz(1)-pggf0-pggb0+pggnpz)/pi y= -0.25d0*(pggfz(2)+pggbz(2)+pgglqz(2))/pi alzr= x/(x*x+y*y) alzi= -y/(x*x+y*y) * *-----Re(1/g^2(zm)) * apis= 16.d0*pis p2zr= -zm2 p2zi= 0.d0 call vbself(p2zr,p2zi,p3qz,fzz) bx= 1.d0/8.d0/gf/zm2+(fw0-fzz(1))/zm2/apis xih= -p3qz(2)/apis ac= 4.d0*pi*alzr bc= -1.d0-8.d0*pi*xih*alzi cc= -4.d0*pi*xih*xih*alzr+bx ifail= 0 call c02ajf(ac,bc,cc,zs,zl,ifail) g2z= zs(1) * p2r= -wm2 p2i= 0.d0 call wself(p2r,p2i,fw) call vbself(p2r,p2i,p3qw,fzw) * fvect(1)= 1.d0+0.5d0*gf/pis*((fw0-fw(1))+p3qw(1)*wm2)- # 8.d0*gf*wm2*(g2z+p3qz(1)/apis) * return end * *-------------------------------------------------------------------------- * subroutine fsw(n,sw,fvecw,iflag) implicit real*8(a-h,o-z) * common/param/eps,pi,pis,qdelta common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 * dimension sw(n),fvecw(n) dimension p3qw(2),fzw(2),fw(2),p3q(2),fz(2),fcw(2) * swr= sw(1) swi= sw(2) wmu= sqrt(swr) wlg= -swi/wmu print 2,wmu,wlg 2 format(1x,2e20.7,//) * call vbself0(pggf0,pggb0,fw0) p2r= -wm2 p2i= 0.d0 call vbself(p2r,p2i,p3qw,fzw) call wself(p2r,p2i,fw) p2r= -swr p2i= -swi call vbself(p2r,p2i,p3q,fz) call wself(p2r,p2i,fcw) * x0= 1.d0/8.d0/gf/wm2 xr= 1.d0+0.5d0*gf*wm2/pis*((fw0-fw(1))/wm2+p3qw(1)-p3q(1)) xi= -0.5d0*gf*wm2/pis*p3q(2) xr= x0*xr xi= x0*xi * fvecw(1)= gf*(swr*xr-swi*xi+(fcw(1)-fw0)/16.d0/pis) # -1.d0/8.d0 fvecw(2)= swr*xi+swi*xr+fcw(2)/16.d0/pis * return end * *--------------------------------------------------------------------------- * subroutine fsz(n,sz,fvecz,iflag) implicit real*8(a-h,o-z) * common/g2/g2z common/param/eps,pi,pis,qdelta common/iaz/alzr,alzi,alzir,alzii common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 * dimension sz(n),fvecz(n) dimension p3qz(2),fzz(2),p3q(2),fz(2) dimension pggfz(2),pgglqz(2),pggf(2),pggb(2),pggbz(2),pgglq(2) dimension evl(2) * szr= sz(1) szi= sz(2) if(szr.ge.0.d0) then zmu= sqrt(szr) zlg= -szi/zmu print 2,zmu,zlg else print 3,szr iflag= -1 endif 2 format(1x,2e20.7,//) 3 format(1x,'*',1e20.7,//) * apis= 16.d0*pis call vbself0(pggf0,pggb0,fw0) pzr= -zm2 pzi= 0.d0 call pself(pzr,pzi,pggfz,pggbz,pgglqz,pggnpz) call vbself(pzr,pzi,p3qz,fzz) p2r= -szr p2i= -szi call pself(p2r,p2i,pggf,pggb,pgglq,pggnp) call vbself(p2r,p2i,p3q,fz) do i=1,2 evl(i)= pggfz(i)+pggbz(i)-pggf(i)-pggb(i)+pgglqz(i) # -pgglq(i) enddo axr= alzir+0.25d0*evl(1)/pi axi= alzii+0.25d0*evl(2)/pi alr= axr/(axr*axr+axi*axi) ali= -axi/(axr*axr+axi*axi) * xr= g2z+(p3qz(1)-p3q(1))/apis xi= -p3q(2)/apis * yr= xr-4.d0*pi*(alr*(xr*xr-xi*xi)-2.d0*ali*xr*xi) yi= xi-4.d0*pi*(2.d0*alr*xr*xi+ali*(xr*xr-xi*xi)) * fvecz(1)= 8.d0*gf*(szr*yr-szi*yi+(fz(1)-fw0)/apis)-1.d0 fvecz(2)= szr*yi+szi*yr+fz(2)/apis * return end * *-----pself-------------------------------------------------------- * subroutine pself(p2r,p2i,pggf,pggb,pgglq,pggnp) implicit real*8 (a-h,o-z) * common/mix/alst,alstz common/param/eps,pi,pis,qdelta common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 common/fmass2/em2,rmm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * dimension bfl(2),bfh(2) dimension b0l(2),b1l(2),b21l(2) dimension b0h(2),b1h(2),b21h(2) dimension b0w(2),b1w(2),b21w(2) dimension pggf(2),pggb(2),pgglq(2) dimension cmpt(2),cmpz(2),clnt(2),clnz(2) dimension v1(2),a1(2),df1(2),v1p(2),a1p(2) * data z3/1.20205690315959428540d0/ * call sbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l) call sbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h) call sbff(p2r,p2i,wm2,wm2,b0w,b1w,b21w) * do i=1,2 bfl(i)= 2.d0*b21l(i)-b0l(i) bfh(i)= 2.d0*b21h(i)-b0h(i) enddo * sllc= 1.d0+3.d0/4.d0/alphai/pi * p2m2= p2r*p2r+p2i*p2i sc2= tqm2 aexpt= alst/pi aexpz= alstz/pi rm= tqm2/sc2 rl= log(rm) bx= -rl-4.d0*z3+55.d0/12.d0 rvr= -0.25d0*p2r/tqm2 rvi= -0.25d0*(p2i-eps)/tqm2 cmpt(1)= p2r/sc2 cmpt(2)= (p2i-eps)/sc2 call cqlnx(cmpt,clnt) cmpz(1)= p2r/zm2 cmpz(2)= (p2i-eps)/zm2 call cqlnx(cmpz,clnz) * if(p2i.eq.0.d0) then add= 0.d0 else add= -2.d0*pi endif fbr= -z3+55.d0/48.d0-0.25d0*clnt(1) fbi= -0.25d0*(clnt(2)+add) flqr= -z3+55.d0/48.d0-0.25d0*clnz(1) flqi= -0.25d0*(clnz(2)+add) call alals(p2r,ap2i,v1,a1,df1,v1p,a1p) * do i=1,2 pggf(i)= 12.d0*bfl(i)*sllc+16.d0/3.d0*bfh(i) pgglq(i)= 44.d0/3.d0*bfl(i) enddo * cmbr= rvr*bx+v1(1) cmbi= rvi*bx+v1(2) pggf(1)= pggf(1)+64.d0/9.d0*aexpt*tqm2/p2m2*(p2r*cmbr+p2i*cmbi) pggf(2)= pggf(2)+64.d0/9.d0*aexpt*tqm2/p2m2*(-p2i*cmbr+p2r*cmbi) pgglq(1)= pgglq(1)+16.d0/9.d0*aexpt*fbr-160.d0/9.d0*aexpz*flqr pgglq(2)= pgglq(2)+16.d0/9.d0*aexpt*fbi-160.d0/9.d0*aexpz*flqi * pggb(1)= 2.d0/3.d0-12.d0*b21w(1)+7.d0*b0w(1) pggb(2)= -12.d0*b21w(2)+7.d0*b0w(2) * p2z= -zm2 dmp= 0.d0 call rbff(p2z,dmp,em2,em2,b0eez,b1eez,b21eez) call rbff(p2z,dmp,rmm2,rmm2,b0mmz,b1mmz,b21mmz) call rbff(p2z,dmp,tm2,tm2,b0tauz,b1tauz,b21tauz) call rbff0(em2,em2,b0ee0,b1ee0,b21ee0) call rbff0(rmm2,rmm2,b0mm0,b1mm0,b21mm0) call rbff0(tm2,tm2,b0tau0,b1tau0,b21tau0) bfez= 2.d0*b21eez-b0eez bfmz= 2.d0*b21mmz-b0mmz bftauz= 2.d0*b21tauz-b0tauz bfe0= 2.d0*b21ee0-b0ee0 bfm0= 2.d0*b21mm0-b0mm0 bftau0= 2.d0*b21tau0-b0tau0 pggl= 4.d0*(bfez+bfmz+bftauz)*sllc pggl0= 4.d0*(bfe0+bfm0+bftau0)*sllc b= 0.00299d0 c= 1.d0 an= -b*log(1.d0+c*zm2)-1.d0/alphai/4.d0/pi*(pggl-pggl0)+ # 1.d0-128.896d0/alphai ap2x= abs(p2r) px= sqrt(ap2x) if(px.lt.0.3d0) then a= 0.d0 b= 0.00835d0 c= 1.d0 else if(px.lt.3.d0) then a= 0.d0 b= 0.00238d0 c= 3.927d0 else if(px.lt.100.d0) then a= 0.00165d0 b= 0.00299d0 c= 1.d0 else if(px.gt.100.d0) then a= 0.00221d0 b= 0.00293d0 c= 1.d0 endif * da= an-0.00165d0 a= a+da pggnp= 4.d0*pi*alphai*(a+b*log(1.d0+c*ap2x)) * return end * *-----vbself0----------------------------------------------------------- * calcola le funzioni pggf pggb fw ad impulso scambiato NULLO * (pggf0, pggb0, fw0) * subroutine vbself0(pggf0,pggb0,fw0) implicit real*8(a-h,o-z) * common/mix/alst,alstz common/param/eps,pi,pis,qdelta common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 common/fmass2/em2,rmm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * data z3/1.20205690315959428540d0/ * call rbff0(em2,em2,b0ee0,b1ee0,b21ee0) call rbff0(rmm2,rmm2,b0mm0,b1mm0,b21mm0) call rbff0(tm2,tm2,b0tau0,b1tau0,b21tau0) call rbff0(tqm2,tqm2,b0tt0,b1tt0,b21tt0) * bfe= 2.d0*b21ee0-b0ee0 bfm= 2.d0*b21mm0-b0mm0 bftau= 2.d0*b21tau0-b0tau0 bft= 2.d0*b21tt0-b0tt0 * sllc= 1.d0+3.d0/4.d0/alphai/pi * sc2= tqm2 aexpt= alst/pi rm= tqm2/sc2 rl= log(rm) bx= -rl-4.d0*z3+55.d0/12.d0 by= 3.d0*rl*rl-11.d0/2.d0*rl+6.d0*z3+pis/2.d0- # 11.d0/8.d0 f10= -3.d0/2.d0*z3-pis/12.d0+23.d0/16.d0 * pggf0= 4.d0*(bfe+bfm+bftau)*sllc+16.d0/3.d0*bft- # 16.d0/9.d0*aexpt*(bx+4.d0*z3-5.d0/6.d0) * call rbff0(wm2,wm2,b0ww0,b1ww0,b21ww0) pggb0= 2.d0/3.d0-12.d0*b21ww0+7.d0*b0ww0 * call rbff0(tqm2,dmy2,b0tb0,b1tb0,b21tb0) fw0f= -3.d0*tqm2*(b0tb0+b1tb0) * call rbff0(wm2,dmy2,b0wp0,b1wp0,b21wp0) call rbff0(zm2,wm2,b0zw0,b1zw0,b21zw0) call rbff0(wm2,bhm2,b0wh0,b1wh0,b21wh0) * s0ww0b= 9.d0/2.d0*(zm2-wm2)*b1zw0+0.25d0*(13.d0*zm2- # 21.0*wm2)*b0zw0+4.d0*wm2*b0ww0 s1ww0b= 2.d0*(wm2-zm2)*(2.d0*b1zw0+b0zw0)+2.d0*wm2* # (2.d0*b1wp0+b0wp0) osh= 0.5d0*(wm2-bhm2)*b1wh0+0.25d0*(5.d0*wm2-bhm2)*b0wh0 fw0b= s0ww0b+sw2*s1ww0b+osh * dg= 2.d0*(3.d0/sw2+(7.d0-4.d0*sw2)/4.d0/sw2**2*log(cw2)) dsww= sw2*dg fw0= fw0b+dsww+fw0f+4.d0*aexpt*tqm2*(0.25d0*by+f10) * return end * *-----vbself------------------------------------------------------------ * calcola le funzioni p3q e fz a impulso arbitrario * (con contributi fermionici e bosonici) * subroutine vbself(p2r,p2i,p3q,fz) implicit real*8(a-h,o-z) * common/mix/alst,alstz common/param/eps,pi,pis,qdelta common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 common/fmass2/em2,rmm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * dimension p3q(2),fz(2) dimension bfl(2),bfh(2) dimension b0l(2),b1l(2),b21l(2) dimension b0h(2),b1h(2),b21h(2) dimension b0w(2),b1w(2),b21w(2) dimension b0zh(2),b1zh(2),b21zh(2) dimension p3qf(2),p3qb(2),fzb(2),fzf(2) dimension cmpt(2),cmpz(2),clnt(2),clnz(2) dimension v1(2),a1(2),f1(2),v1p(2),a1p(2) * data z3/1.20205690315959428540d0/ * p2m2= p2r*p2r+p2i*p2i * call sbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l) call sbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h) * do i=1,2 bfl(i)= 2.d0*b21l(i)-b0l(i) bfh(i)= 2.d0*b21h(i)-b0h(i) enddo * do i=1,2 p3qf(i)= 10.d0*bfl(i)+2.d0*bfh(i) enddo * if(p2i.eq.0.d0) then add= 0.d0 else add= -2.d0*pi endif sc2= tqm2 aexpt= alst/pi aexpz= alstz/pi rm= tqm2/sc2 rl= log(rm) bx= -rl-4.d0*z3+55.d0/12.d0 by= 3.d0*rl*rl-11.d0/2.d0*rl+6.d0*z3+pis/2.d0-11.d0/8.d0 rvr= -0.25d0*p2r/tqm2 rvi= -0.25d0*(p2i-eps)/tqm2 xvr= -p2r/tqm2 xvi= -(p2i-eps)/tqm2 * cmpt(1)= p2r/sc2 cmpt(2)= (p2i-eps)/sc2 call cqlnx(cmpt,clnt) cmpz(1)= p2r/zm2 cmpz(2)= (p2i-eps)/zm2 call cqlnx(cmpz,clnz) * fbr= -z3+55.d0/48.d0-0.25d0*clnt(1) fbpr= -z3+43.d0/48.d0-0.25d0*clnt(1) flqr= -z3+55.d0/48.d0-0.25d0*clnz(1) flqpr= -z3+43.d0/48.d0-0.25d0*clnz(1) fbi= -0.25d0*(clnt(2)+add) fbpi= -0.25d0*(clnt(2)+add) flqi= -0.25d0*(clnz(2)+add) flqpi= -0.25d0*(clnz(2)+add) call alals(p2r,p2i,v1,a1,f1,v1p,a1p) * cmbr= 8.d0/3.d0*aexpt*tqm2*(rvr*bx+v1(1)) cmbi= 8.d0/3.d0*aexpt*tqm2*(rvi*bx+v1(2)) p3qf(1)= p3qf(1)+1.d0/p2m2*(p2r*cmbr+p2i*cmbi)-4.d0/3.d0*aexpt* # fbr-8.d0*aexpz*flqr p3qf(2)= p3qf(2)+1.d0/p2m2*(p2r*cmbi-p2i*cmbr)-4.d0/3.d0*aexpt* # fbi-8.d0*aexpz*flqi * bfsr= bfh(1)-bfl(1) bfsi= bfh(2)-bfl(2) fzf(1)= -0.5d0*(p2r*bfsr-p2i*bfsi)-3.d0/2.d0*tqm2*b0h(1) fzf(2)= -0.5d0*(p2r*bfsi+p2i*bfsr)-3.d0/2.d0*tqm2*b0h(2) * fzf(1)= fzf(1)+aexpt*tqm2*(-2.d0/3.d0*rvr*bx+by-5.d0/3.d0*v1(1)+ # a1(1))-2.d0/3.d0*aexpt*(p2r*fbr-p2i*fbi) fzf(2)= fzf(2)+aexpt*tqm2*(-2.d0/3.d0*rvi*bx-5.d0/3.d0*v1(2)+ # a1(2))-2.d0/3.d0*aexpt*(p2r*fbi+p2i*fbr) * call rbff0(wm2,wm2,b0ww0,b1ww0,b21ww0) call sbff(p2r,p2i,wm2,wm2,b0w,b1w,b21w) call sbff(p2r,p2i,zm2,bhm2,b0zh,b1zh,b21zh) * pzzbr= -9.d0*b21w(1)+25.d0/4.d0*b0w(1)+2.d0/3.d0 pzzhr= -b21zh(1)-b1zh(1)-0.25d0*b0zh(1) pzzbr= pzzbr+pzzhr pzzbi= -9.d0*b21w(2)+25.d0/4.d0*b0w(2) pzzhi= -b21zh(2)-b1zh(2)-0.25d0*b0zh(2) pzzbi= pzzbi+pzzhi sigzzbr= -2.d0*wm2*b0w(1)+4.d0*wm2*b0ww0 sigzzhr= zm2*(0.5d0*b1zh(1)+5.d0/4.d0*b0zh(1))- # 0.5d0*bhm2*(b1zh(1)+0.5d0*b0zh(1)) sigzzbr= sigzzbr+sigzzhr sigzzbi= -2.d0*wm2*b0w(2) sigzzhi= zm2*(0.5d0*b1zh(2)+5.d0/4.d0*b0zh(2))- # 0.5d0*bhm2*(b1zh(2)+0.5d0*b0zh(2)) sigzzbi= sigzzbi+sigzzhi * s33br= p2r*pzzbr-p2i*pzzbi+sigzzbr s33bi= p2r*pzzbi+p2i*pzzbr+sigzzbi * s3gbr= p2r*(-10.d0*b21w(1)+13.d0/2.d0*b0w(1)+2.d0/3.d0)- # p2i*(-10.d0*b21w(2)+13.d0/2.d0*b0w(2))-2.d0*wm2* # (b0w(1)-b0ww0) s3gbi= p2i*(-10.d0*b21w(1)+13.d0/2.d0*b0w(1)+2.d0/3.d0)+ # p2r*(-10.d0*b21w(2)+13.d0/2.d0*b0w(2))-2.d0*wm2*b0w(2) p3qb(1)= 1.d0/p2m2*(p2r*s3gbr+p2i*s3gbi) p3qb(2)= 1.d0/p2m2*(p2r*s3gbi-p2i*s3gbr) fzb(1)= s33br-s3gbr fzb(2)= s33bi-s3gbi * do i=1,2 p3q(i)= p3qf(i)+p3qb(i) fz(i)= fzf(i)+fzb(i) enddo * return end * *-----wself------------------------------------------------------------ * calcola fw (a p**2 arbitrario) * subroutine wself(p2r,p2i,fw) implicit real*8(a-h,o-z) * common/mix/alst,alstz common/param/eps,pi,pis,qdelta common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 common/fmass2/em2,rmm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * dimension fw(2),fwb(2),fwf(2) dimension b0c(2),b1c(2),b21c(2) dimension b0l(2),b1l(2),b21l(2) dimension b0n(2),b1n(2),b21n(2) dimension b0w(2),b1w(2),b21w(2) dimension b0wp(2),b1wp(2),b21wp(2) dimension b0zw(2),b1zw(2),b21zw(2) dimension b0wh(2),b1wh(2),b21wh(2) dimension cmpt(2),cmpz(2),clnt(2),clnz(2) dimension v1(2),a1(2),f1(2),v1p(2),a1p(2) * data z3/1.20205690315959428540d0/ * call sbff(p2r,p2i,tqm2,dmy2,b0c,b1c,b21c) call sbff(p2r,p2i,tqm2,tqm2,b0n,b1n,b21n) call sbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l) * auxr= 3.d0*(b21c(1)+b1c(1))-2.d0*b21n(1)+b0n(1)-b21l(1)- # 3.d0*b1l(1)-b0l(1) auxi= 3.d0*(b21c(2)+b1c(2))-2.d0*b21n(2)+b0n(2)-b21l(2)- # 3.d0*b1l(2)-b0l(2) * fwf(1)= 2.d0*(p2r*auxr-p2i*auxi)-3.d0*tqm2*(b1c(1)+b0c(1)) fwf(2)= 2.d0*(p2r*auxi+p2i*auxr)-3.d0*tqm2*(b1c(2)+b0c(2)) * sc2= tqm2 aexpt= alst/pi aexpz= alstz/pi rm= tqm2/sc2 rl= log(rm) bx= -rl-4.d0*z3+55.d0/12.d0 by= 3.d0*rl*rl-11.d0/2.d0*rl+6.d0*z3+pis/2.d0-11.d0/8.d0 rvr= -0.25d0*p2r/tqm2 rvi= -0.25d0*(p2i-eps)/tqm2 xvr= -p2r/tqm2 xvi= -(p2i-eps)/tqm2 * if(p2i.eq.0.d0) then add= 0.d0 else add= -2.d0*pi endif cmpt(1)= p2r/sc2 cmpt(2)= (p2i-eps)/sc2 call cqlnx(cmpt,clnt) cmpz(1)= p2r/zm2 cmpz(2)= (p2i-eps)/zm2 call cqlnx(cmpz,clnz) * fbr= -z3+55.d0/48.d0-0.25d0*clnt(1) flqr= -z3+55.d0/48.d0-0.25d0*clnz(1) fbi= -0.25d0*(clnt(2)+add) flqi= -0.25d0*(clnz(2)+add) call alals(p2r,p2i,v1,a1,f1,v1p,a1p) * fwf(1)= fwf(1)+aexpt*tqm2*(xvr*bx+by+4.d0*f1(1)-8.d0/3.d0* # (rvr*bx+v1(1)))+4.d0/3.d0*aexpt*(p2r*fbr-p2i*fbi) fwf(2)= fwf(2)+aexpt*tqm2*(xvi*bx+4.d0*f1(2)-8.d0/3.d0* # (rvi*bx+v1(2)))+4.d0/3.d0*aexpt*(p2r*fbi+p2i*fbr) * call rbff0(wm2,wm2,b0ww0,b1ww0,b21ww0) call sbff(p2r,p2i,wm2,dmy2,b0wp,b1wp,b21wp) call sbff(p2r,p2i,wm2,wm2,b0w,b1w,b21w) call sbff(p2r,p2i,zm2,wm2,b0zw,b1zw,b21zw) call sbff(p2r,p2i,wm2,bhm2,b0wh,b1wh,b21wh) * s0wwr= 9.d0/2.d0*(zm2-wm2)*b1zw(1)+(13.d0*zm2-21.d0*wm2)/4.d0* # b0zw(1)+0.5d0*(wm2-bhm2)*b1wh(1)+0.25d0* # (5.d0*wm2-bhm2)*b0wh(1)+4.d0*wm2*b0ww0 s0wwi= 9.d0/2.d0*(zm2-wm2)*b1zw(2)+(13.d0*zm2-21.d0*wm2)/4.d0* # b0zw(2)+0.5d0*(wm2-bhm2)*b1wh(2)+0.25d0* # (5.d0*wm2-bhm2)*b0wh(2) * s1wwr= 2.d0*(wm2-zm2)*(2.d0*b1zw(1)+b0zw(1))+2.d0*wm2*(2.d0* # b1wp(1)+b0wp(1)) s1wwi= 2.d0*(wm2-zm2)*(2.d0*b1zw(2)+b0zw(2))+2.d0*wm2*(2.d0* # b1wp(2)+b0wp(2)) * p0wwr= 2.d0/3.d0-9.d0*b21zw(1)-9.d0*b1zw(1)+7.d0/4.d0*b0zw(1)- # b21wh(1)-b1wh(1)-0.25d0*b0wh(1) p0wwi= -9.d0*b21zw(2)-9.d0*b1zw(2)+7.d0/4.d0*b0zw(2)- # b21wh(2)-b1wh(2)-0.25d0*b0wh(2) * p1wwr= 8.d0*b21zw(1)-2.d0*b0zw(1)+8.d0*b1zw(1)-8.d0*b21wp(1)- # 8.d0*b1wp(1)+2.d0*b0wp(1) p1wwi= 8.d0*b21zw(2)-2.d0*b0zw(2)+8.d0*b1zw(2)-8.d0*b21wp(2)- # 8.d0*b1wp(2)+2.d0*b0wp(2) * pwwr= p0wwr+sw2*p1wwr swwr= s0wwr+sw2*s1wwr pwwi= p0wwi+sw2*p1wwi swwi= s0wwi+sw2*s1wwi swwbr= p2r*pwwr-p2i*pwwi+swwr swwbi= p2r*pwwi+p2i*pwwr+swwi * s3gbr= p2r*(-10.d0*b21w(1)+13.d0/2.d0*b0w(1)+2.d0/3.d0)- # p2i*(-10.d0*b21w(2)+13.d0/2.d0*b0w(2))-2.d0*wm2* # (b0w(1)-b0ww0) s3gbi= p2i*(-10.d0*b21w(1)+13.d0/2.d0*b0w(1)+2.d0/3.d0)+ # p2r*(-10.d0*b21w(2)+13.d0/2.d0*b0w(2))-2.d0*wm2*b0w(2) * fwb(1)= swwbr-s3gbr fwb(2)= swwbi-s3gbi * do i=1,2 fw(i)= fwb(i)+fwf(i) enddo * return end * *-----rbff0----------------------------------------------------------------- * the scalar form factors b0,b1,b21 at zero momentum * subroutine rbff0(rm12,rm22,b0,b1,b21) implicit real*8(a-h,o-z) * common/param/eps,pi,pis,qdelta common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 * dimension arg(2),cln(2),fr(3) * if(rm12.eq.0.d0) then arm12= 1.d-20 else arm12= rm12 endif if(rm22.eq.0.d0) then arm22= 1.d-20 else arm22= rm22 endif if (arm12.eq.arm22) then fact= qdelta-log(arm12) b0= fact b1= -0.5d0*fact b21= fact/3.d0 return else n= 3 yr= arm12/(arm12-arm22) omyr= -arm22/(arm12-arm22) yi= -eps/(arm12-arm22) call rcg(n,yr,yi,omyr,fr) arg(1)= omyr arg(2)= -yi call cqlnx(arg,cln) f1r= fr(1)-cln(1) f2r= 2.d0*fr(2)-cln(1) f3r= 3.d0*fr(3)-cln(1) fact= qdelta-log(arm22) b0= fact-f1r b1= 0.5d0*(-fact+f2r) b21= 1.d0/3.d0*(fact-f3r) return endif end * * *-------------------------sbff------------------------------------* * computes the one-loop two-point form factors /i*pi^2 * * b0,b1,b21 * * !!! b22 must be computed separately !!! * * input parameters are p^2,m1^2,m2^2 * * p^2 can be the complex pole * * m1=m2=0 * * m2=0 * * m1=m2 * *-----------------------------------------------------------------* * subroutine sbff(p2r,p2i,rm12,rm22,b0,b1,b21) implicit real*8(a-h,o-z) character*1 op * common/opt/op common/param/eps,pi,pis,qdelta common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 * dimension cmp2(2) dimension clnmp2(2) dimension b0(2),b1(2),b21(2) * p2m= p2r*p2r+p2i*p2i * if(rm12.eq.0.d0.and.rm22.eq.0.d0) then if(p2i.eq.0.d0) then cmp2(1)= p2r cmp2(2)= -eps call cqlnx(cmp2,clnmp2) b0(1)= qdelta+2.d0-clnmp2(1) b0(2)= -clnmp2(2) else cmp2(1)= p2r cmp2(2)= p2i call cqlnx(cmp2,clnmp2) b0(1)= qdelta+2.d0-clnmp2(1) b0(2)= -clnmp2(2)+2.d0*pi endif b1(1)= -0.5d0*b0(1) b1(2)= -0.5d0*b0(2) b21(1)= 1.d0/18.d0+1.d0/3.d0*b0(1) b21(2)= 1.d0/3.d0*b0(2) else if(rm22.eq.0.d0.and.rm12.ne.0.d0) then a0= rm12*(-qdelta-1.d0+log(rm12)) if(p2i.eq.0.d0) then cmp2(1)= 1.d0+p2r/rm12 cmp2(2)= -eps/rm12 call cqlnx(cmp2,clnmp2) b0(1)= qdelta+2.d0-log(rm12)-(1.d0+rm12/p2r)*clnmp2(1) b0(2)= -(1.d0+rm12/p2r)*clnmp2(2) b1(1)= 0.5d0*(a0+(rm12-p2r)*b0(1))/p2r b1(2)= 0.5d0*(rm12-p2r)*b0(2)/p2r b21(1)= (3.d0*rm12+p2r)/18.d0/p2r+(rm12-p2r)*a0/3.d0/ # (p2r*p2r)+(1.d0-rm12/p2r*(1.d0-rm12/p2r))* # b0(1)/3.d0 b21(2)= (1.d0-rm12/p2r*(1.d0-rm12/p2r))*b0(2)/3.d0 else cmp2(1)= 1.d0+p2r/rm12 cmp2(2)= p2i/rm12 call cqlnx(cmp2,clnmp2) xr= p2r/p2m xi= -p2i/p2m fctr= 1.d0+rm12*xr fcti= rm12*xi b0(1)= qdelta+2.d0-log rg(2).eq.0.d0) then ' if (arg(1).gt.0.d0) then  teta= 0.d0  else- teta= pia endif res(2)= teta else! tnteta= aarg(2)/aarg(1)i teta= atan(tnteta) sr= arg(1)/aarg(1) si= arg(2)/aarg(2) if (sr.gt.0.d0) then res(2)= si*tetaf else" res(2)= si*(pi-teta) endifb endif * return end2*2(*--------------------------------------*(* computes ln(1-x) *(* usually |x| << 1 *(*--------------------------------------**)' subroutine cqlnomx(arg,omarg,res), implicit real*8(a-h,o-z)*2< dimension arg(2),omarg(2),res(2),ares(2),ct(10),sn(10)*0 zr= arg(1) zi= arg(2) zm2= zr*zr+zi*zi zm= sqrt(zm2), if (zm.lt.1.d-7) then  ct(1)= zr/zm sn(1)= zi/zm do n=2,1020 ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1)0 sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1) enddo1 ares(1)= ct(10)/10.d0  ares(2)= sn(10)/10.d0( do k=9,1,-1a) ares(1)= ares(1)*zm+ct(k)/k)) ares(2)= ares(2)*zm+sn(k)/k) enddo  ares(1)= -ares(1)*zm ares(2)= -ares(2)*zm else call cqlnx(omarg,ares) endif0 do i= 1,2  res(i)= ares(i)6 enddo1 return*v endp*t7*-----------------------------------------------------*/7* computes the function g_n(x) *07*-----------------------------------------------------*d*i) subroutine cg(n,zr,zi,omzr,gfr,gfi)  implicit real*8(a-h,o-z)*l$ common/param/eps,pi,pis,qdelta= common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2r*e dimension gfr(n),gfi(n)2B dimension a(2,4),ra(16),z(2),omz(2),oz(2),clnomz(2),clnoz(2)3 dimension ca(2,10),aux(2),zp(2),ct(16),sn(16)a* + data ra/-0.2q+0,0.666666666666667d-1,> # -0.952380952380952d-2,-0.396825396825397d-3,> # 0.317460317460317d-3,-0.132275132275132d-4,= # -0.962000962000962d-5,0.105218855218855d-5,/> # 0.266488361726450d-6,-0.488745528428000d-7,= # -0.675397500794000d-8,0.190720263471000d-8, ? # 0.153663007690000d-9,-0.679697905790000d-10,2? # -0.293683556000000d-11,0.228836696000000d-11/w*z z(1)= zr z(2)= zi omz(1)= omzr omz(2)= -zi ) if (zr.eq.0.d0.and.zi.eq.0.d0) thenm do k=1,n gfr(k)= -1.d0/k/k gfi(k)= 0.d0  enddo2 return. else if (zr.eq.1.d0.and.zi.eq.0.d0) then a(1,1)= -1.d0 a(2,1)= pi0 do j=2,4.7 a(1,j)= ((j-1.d0)*a(1,j-1)-1.d0/j)/j1. a(2,j)= (j-1.d0)/j*a(2,j-1) enddo  do k=1,n gfr(k)= a(1,k) gfi(k)= a(2,k) enddow return else 2 zmod2= zr*zr+zi*zi zmod= sqrt(zmod2)1 call cqlnx(omz,clnomz)* * |z| < 4 * . if (zmod.lt.4.d0) then2 oz(1)= -z(1)) oz(2)= -z(2) *! call cqlnx(oz,clnoz)mF ca(1,1)= omz(1)*clnomz(1)-omz(2)*clnomz(2)+z(1)*clnoz(1)-( # z(2)*clnoz(2)-1.d0F ca(2,1)= omz(1)*clnomz(2)+omz(2)*clnomz(1)+z(1)*clnoz(2)+# # z(2)*clnoz(1)w if (n.eq.1) then gfr(1)= ca(1,1)= gfi(1)= ca(2,1)w return else do j= 2,n. jm1= j-1E ca(1,j)= ((j-1.d0)*(z(1)*ca(1,jm1)-z(2)*ca(2,jm1))+ F # omz(1)*clnomz(1)-omz(2)*clnomz(2)-1.d0/j)/jE ca(2,j)= ((j-1.d0)*(z(1)*ca(2,jm1)+z(2)*ca(1,jm1))+w? # omz(1)*clnomz(2)+omz(2)*clnomz(1))/j3 enddo2 do k=1,n# gfr(k)= ca(1,k)e# gfi(k)= ca(2,k)- enddo- return endif * * |z| > 4*r else' aux(1)= (-zr*omzr+zi**2)/zmod2f aux(2)= zi/zmod2  call cqlnx(aux,zp)& zpm2= zp(1)*zp(1)+zp(2)*zp(2) zpm= sqrt(zpm2) ct(1)= zp(1)/zpmm sn(1)= zp(2)/zpms do k=2,16,, ct(k)= ct(1)*ct(k-1)-sn(1)*sn(k-1), sn(k)= sn(1)*ct(k-1)+ct(1)*sn(k-1) enddo . ca(1,4)= ra(16)*ct(16)*zpm+ra(15)*ct(15). ca(2,4)= ra(16)*sn(16)*zpm+ra(15)*sn(15) do j=14,1,-1* ca(1,4)= ca(1,4)*zpm+ra(j)*ct(j)* ca(2,4)= ca(2,4)*zpm+ra(j)*sn(j) enddo0 ca(1,4)= ca(1,4)*zpm ca(2,4)= ca(2,4)*zpm do j= 3,1,-1 jp1= j+1C ca(1,j)= ((ca(1,jp1)+1.d0/jp1)*z(1)+ca(2,jp1)*z(2))/zmod2 C ca(2,j)= (ca(2,jp1)*z(1)-(ca(1,jp1)+1.d0/jp1)*z(2))/zmod2 enddom do k=1,n ' gfr(k)= (ca(1,k)+clnomz(1))/k ' gfi(k)= (ca(2,k)+clnomz(2))/k2 enddo( return endif3 endif end *l7*-----------------------------------------------------*d7* computes the function re g_n(x) * 7*-----------------------------------------------------*-*-& subroutine rcg(n,zr,zi,omzr,gfr) implicit real*8(a-h,o-z)*o$ common/param/eps,pi,pis,qdelta= common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 *  dimension gfr(n)@ dimension a(4),ra(16),z(2),omz(2),oz(2),clnomz(2),clnoz(2)3 dimension ca(2,10),aux(2),zp(2),ct(16),sn(16) * + data ra/-0.2q+0,0.666666666666667d-1, > # -0.952380952380952d-2,-0.396825396825397d-3,> # 0.317460317460317d-3,-0.132275132275132d-4,= # -0.962000962000962d-5,0.105218855218855d-5, > # 0.266488361726450d-6,-0.488745528428000d-7,= # -0.675397500794000d-8,0.190720263471000d-8, ? # 0.153663007690000d-9,-0.679697905790000d-10,*? # -0.293683556000000d-11,0.228836696000000d-11/ *m z(1)= zr z(2)= zi omz(1)= omzr omz(2)= -zii) if (zr.eq.0.d0.and.zi.eq.0.d0) then2 do k=1,n gfr(k)= -1.d0/k/k enddo return. else if (zr.eq.1.d0.and.zi.eq.0.d0) then a(1)= -1.d0 do j=2,4 3 a(j)= ((j-1.d0)*a(j-1)-1.d0/j)/j  enddo, do k=1,n gfr(k)= a(k)  enddo2 return else = zmod2= zr*zr+zi*zi zmod= sqrt(zmod2), call cqlnx(omz,clnomz)* * |z| < 4 * l if (zmod.lt.4.d0) then  oz(1)= -z(1)b oz(2)= -z(2) )! call cqlnx(oz,clnoz)3F ca(1,1)= omz(1)*clnomz(1)-omz(2)*clnomz(2)+z(1)*clnoz(1)-( # z(2)*clnoz(2)-1.d0F ca(2,1)= omz(1)*clnomz(2)+omz(2)*clnomz(1)+z(1)*clnoz(2)+# # z(2)*clnoz(1)m if (n.eq.1) then gfr(1)= ca(1,1)  return else do j= 2,n  jm1= j-1E ca(1,j)= ((j-1.d0)*(z(1)*ca(1,jm1)-z(2)*ca(2,jm1))+ F # omz(1)*clnomz(1)-omz(2)*clnomz(2)-1.d0/j)/jE ca(2,j)= ((j-1.d0)*(z(1)*ca(2,jm1)+z(2)*ca(1,jm1))+m? # omz(1)*clnomz(2)+omz(2)*clnomz(1))/j  enddop do k=1,n# gfr(k)= ca(1,k)2 enddo  return endifl*l* |z| > 4* else' aux(1)= (-zr*omzr+zi**2)/zmod2  aux(2)= zi/zmod2  call cqlnx(aux,zp) & zpm2= zp(1)*zp(1)+zp(2)*zp(2) zpm= sqrt(zpm2) ct(1)= zp(1)/zpmc sn(1)= zp(2)/zpm  do k=2,16/ ct(k)= ct(1)*ct(k-1)-sn(1)*sn(k-1)p/ sn(k)= sn(1)*ct(k-1)+ct(1)*sn(k-1)p enddo1 ca(1,4)= ra(16)*ct(16)*zpm+ra(15)*ct(15)x1 ca(2,4)= ra(16)*sn(16)*zpm+ra(15)*sn(15)1 do j=14,1,-1(- ca(1,4)= ca(1,4)*zpm+ra(j)*ct(j)*- ca(2,4)= ca(2,4)*zpm+ra(j)*sn(j)( enddo ca(1,4)= ca(1,4)*zpm* ca(2,4)= ca(2,4)*zpm0 do j= 3,1,-1d jp1= j+1 F ca(1,j)= ((ca(1,jp1)+1.d0/jp1)*z(1)+ca(2,jp1)*z(2))/zmod2F ca(2,j)= (ca(2,jp1)*z(1)-(ca(1,jp1)+1.d0/jp1)*z(2))/zmod2 enddo do k=1,n * gfr(k)= (ca(1,k)+clnomz(1))/k enddo returnr endif) endif1 end)*2*------------------------------------------------*2* computes e(i)= xmi(i,j)*d(j) *2* i,j=1,2 *2*------------------------------------------------c*2" subroutine multi2(ea,xmi,da) implicit real*8(a-h,o-z)*,( dimension ea(2,2),da(2,2),xmi(2,2)*- do j= 1,2( do i= 1,27 ea(j,i)= xmi(i,1)*da(j,1)+xmi(i,2)*da(j,2) enddo enddol return end *+2*------------------------------------------------*2* computes e(i)= xmi(i,j)*d(j) *2* i,j=1,3 *2*------------------------------------------------**m" subroutine multi3(ea,xmi,da) implicit real*8(a-h,o-z)*2( dimension ea(2,3),da(2,3),xmi(3,3)*1 do j= 1,22 do i= 1,38 ea(j,i)= xmi(i,1)*da(j,1)+xmi(i,2)*da(j,2)+& # xmi(i,3)*da(j,3) enddo enddo return end * 3*-------------------------------------------------*13* computes the roots of the quadratic form *)3* -p2*x^2+(p2+m2^2-m1^2)*x+m1^2 = 0 *23*-------------------------------------------------*p*)9 subroutine roots(p2r,p2i,rm12,rm22,rpr,rpi,rmr,rmi,b# # omrpr,omrmr)  implicit real*8(a-h,o-z)* $ common/param/eps,pi,pis,qdelta= common/bpar/wm,zm,wm2,zm2,zg,gf,alphai,bhm,bhm2,cw2,sw2 *  dimension zsm(2),zlg(2) *= external c02ahf * - ar= -p2r ai= -p2i br= p2r+rm22-rm12 bi= p2ix cr= rm12 ci= -eps*  ifail= 02 call c02ahf(ar,ai,br,bi,cr,ci,zsm,zlg,ifail)*  rmr= zsm(1)  omrmr= 1.d0-rmr. rmi= zsm(2)3 rpr= zlg(1)# omrpr= 1.d0-rprb rpi= zlg(2)1* return end *aJ*-------------------------------------------------------------------------*- real*8 function tralphas(rs0,rs,als,nf)2 implicit real*8(a-h,o-z)*2$ common/param/eps,pi,pis,qdelta< common/fmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy*l! dimension b0(5),b1(5),b2(5) *27*-----limits for lambda_5 are 1 mev < lambda_5 < 10 gev *m if(als.eq.0.d0) then tralphas= 0.d0a else x1= 0.001d0 x2= 10.0d0- xacc= 1.d-12 - qcdl= tqcdlam(nf,als,rs0,x1,x2,xacc) pqcdl5= qcdlt do i=1,52, b0(i)= (11.d0-2.d0/3.d0*i)/4.d0/ b1(i)= (102.d0-38.d0/3.d0*i)/16.d0 3 b2(i)= 0.5d0*(2857.d0-i*(5033.d0/9.d0-r* # 325.d0/27.d0*i))/64.d0 enddo*  if(rs.lt.bqm) then+ rat= bqm/qcdl  rl= 2.d0*log(rat)* rll= log(rl) rb= log(b0(5)/b0(4))( fac= b1(5)/b0(5)-b1(4)/b0(4)) facp= b2(5)/b0(5)-b2(4)/b0(4))A fac2= b1(5)*b1(5)/b0(5)/b0(5)-b1(4)*b1(4)/b0(4)/b0(4)d9 rhs= (b0(5)-b0(4))*rl+fac*rll-b1(4)/b0(4)*rb+ = # b1(5)/b0(5)/b0(5)*fac*rll/rl+1.d0/b0(5)/rl*(& # fac2-facp-7.d0/72.d0) rhs= rhs/b0(4) ratl2= exp(rhs) ( qcdl= qcdl/sqrt(ratl2)  pqcdl4= qcdl nfe= nf-1  if(rs.lt.cqm) then rat= cqm/qcdl rl= 2.d0*log(rat) rll= log(rl)r# rb= log(b0(4)/b0(3))0+ fac= b1(4)/b0(4)-b1(3)/b0(3)c, facp= b2(4)/b0(4)-b2(3)/b0(3)D fac2= b1(4)*b1(4)/b0(4)/b0(4)-b1(3)*b1(3)/b0(3)/b0(3)< rhs= (b0(4)-b0(3))*rl+fac*rll-b1(3)/b0(3)*rb+@ # b1(4)/b0(4)/b0(4)*fac*rll/rl+1.d0/b0(4)/rl*() # fac2-facp-7.d0/72.d0)) rhs= rhs/b0(3)m ratl2= exp(rhs)+ qcdl= qcdl/sqrt(ratl2)  pqcdl3= qcdls nfe= nf-2 endif elsea nfe= nf endif*)# qcdb0= 11.d0-2.d0/3.d0*nfe % qcdb1= 102.d0-38.d0/3.d0*nfe*0 qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nfe+% # 325.d0/27.d0*nfe*nfe)# qcda= 2.d0*log(rs/qcdl)*2@ tralphas= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*A # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)-< # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0))* endif- return endd*rM*-----QCDLAM----------------------------------------------------------------- /*-----COMPUTES LAMBA^NF_MSBAR FROM ALPHA_S(RS0))*3 real*8 function tqcdlam(nf,als,rs,x1,x2,xacc)r implicit real*8(a-h,o-z)*( parameter (jmax=50) *## fmid= tqcdscale(nf,als,rs,x2)1 f= tqcdscale(nf,als,rs,x1) if (f*fmid.ge.0.d0) then7 print*,'root must be bracketed for bisection'm print 1,als 1 1 format(/' error detected by tqcdlam ',/r8 # ' current value of alpha_s = ',e20.5) stop endif) if (f.lt.0.d0) then- tqcdlam= x1- dx= x2-x1 else tqcdlam= x2n dx= x1-x2 endif  do j=1,jmaxi dx= dx*0.5d0 xmid= tqcdlam+dx) fmid= tqcdscale(nf,als,rs,xmid) ) if (fmid.le.0.d0) tqcdlam= xmido4 if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return enddos! pause 'too many bisections'  end *-K*-----QCDSCALE-------------------------------------------------------------i/*-----COMPUTES LAMBA^NF_MSBAR FROM ALPHA_S(RS0)m*c, real*8 function tqcdscale(nf,als,rs,x) implicit real*8(a-h,o-z)*m$ common/param/eps,pi,pis,qdelta*f qcdb0= 11.d0-2.d0/3.d0*nf ! qcdb1= 102.d0-38.d0/3.d0*nfe? qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nf+325.d0/27.d0*nf*nf)  qcda= 2.d0*log(rs/x)C tqcdscale= als-(4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*dA # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)-2= # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0)))3** return end *pC*-----ALALS--------------------------------------------------------mD* THE O(ALPHA*ALPHA_S) CORRECTIONS TO SELF-ENERGIES ARE COMPUTED* 7 subroutine alals(op2r,op2i,ov1,oa1,of1,ov1p,oa1p)  implicit real*16(a-h,p-z)f implicit real*8(o)*- common/ohm/ohm) common/tqparam/qpi,qpis,qeps,qdeltaaC common/tqmasses2/qem2,qmm2,qtm2,qnm2,quqm2,qdqm2,qcqm2,qsqm2,)0 # qbqm2,qtqm2,qzm2,qwm2,qhm2*= dimension r(2),x(2),xo(2),umx(2),a(2),b(2),b2(2),ab(2),pA # rmr(2),rumr(2),rp(2),rm(2),f(2),h(2),g(2),rm2(2),-@ # rm4(2),ri(2),r2i(2),xi(2),x2i(2),rar(2),rrar(2),4 # xp(2),umr(2),trmr(2),trumr(2),rpp(2)4 dimension ov1(2),oa1(2),of1(2),ov1p(2),oa1p(2)*6' data z3/1.20205690315959428540q0/i*aA*-----this subroutine uses the wrong metric, thus the sign of p^2 * must be changede* (' qpi= 3.141592653589793238462643q0  qpis= qpi*qpi. qdelta= 0.q0 qeps= 1.q-90*  qem= 0.51099906q-3 qmm= 0.10565839q0  qnm= 1.q-10e qtm= 1.7841q0e quqm= 0.041q0t qdqm= 0.041q0i qcqm= 1.5q0t qsqm= 0.15q0 qbqm= 4.7q0  qtqm= 175.6q0  qzm= 91.1867q0 qwm= 80.48q0 qhm= ohm*1.d15*1.q-15a*  qem2= qem*qemn qmm2= qmm*qmmi qtm2= qtm*qtmn qnm2= qnm*qnmr quqm2= quqm*quqm qdqm2= qdqm*qdqm qcqm2= qcqm*qcqm qsqm2= qsqm*qsqm qbqm2= qbqm*qbqm qtqm2= qtqm*qtqm qzm2= qzm*qzm  qwm2= qwm*qwm  qhm2= qhm*qhmu* p2r= op2r*1.d15*1.q-15 p2i= op2i*1.d15*1.q-15 z2= qpis/6.q0  r(1)= -0.25q0*p2r/qtqm2 $ r(2)= -0.25q0*(p2i-qeps)/qtqm2 x(1)= -p2r/qtqm2 x(2)= -(p2i-qeps)/qtqm2l do i=1,2 xo(i)= -x(i) enddo umx(1)= 1.q0-x(1)) umx(2)= -x(2)( umr(1)= 1.q0-r(1)  umr(2)= -r(2)i*r call tcqlnx(xo,a)* call tcqlnomx(x,umx,b) b2(1)= b(1)*b(1)-b(2)*b(2) b2(2)= 2.q0*b(1)*b(2)1 ab(1)= a(1)*b(1)-a(2)*b(2) ab(2)= a(1)*b(2)+a(2)*b(1)*) if(r(1).gt.0.q0) then=& rmr(1)= 0.5q0*r(2)/sqrt(r(1)) rmr(2)= -sqrt(r(1)) else if(r(1).lt.0.q0) then rmr(1)= sqrt(-r(1))( rmr(2)= -0.5q0*r(2)/sqrt(-r(1)) endif  if(umr(1).gt.0.q0) then rumr(1)= sqrt(umr(1))* rumr(2)= -0.5q0*r(2)/sqrt(umr(1))" else if(umr(1).lt.0.q0) then* rumr(1)= 0.5q0*r(2)/sqrt(-umr(1)) rumr(2)= -sqrt(-umr(1)) endif  do i=1,2 rp(i)= rumr(i)+rmr(i)- rm(i)= rumr(i)-rmr(i)- enddo  do i=1,2 trmr(i)= rp(i)-rm(i) trumr(i)= rp(i)+rm(i)- enddo- call tcqlnx(rp,f)  call tcqlnx(trumr,h) call tcqlnx(trmr,g)r**% rm2(1)= rm(1)*rm(1)-rm(2)*rm(2)p rm2(2)= 2.q0*rm(1)*rm(2) arm= rm(1) arm2= arm*arm, arm4= arm2*arm2o zrm= rm(2)/rm(1) zrm2= zrm*zrm6 zrm4= zrm2*zrm2 , brm= rm(2) brm2= brm*brm2 brm4= brm2*brm2 urm= rm(1)/rm(2) urm2= urm*urm1 urm4= urm2*urm2 8 if(abs(zrm).lt.1.q0) then-+ rm4(1)= arm4*(1.q0-6.q0*zrm2+zrm4)0$ else if(abs(urm).lt.1.q0) then+ rm4(1)= brm4*(1.q0-6.q0*urm2+urm4)/ endif % rm4(2)= 4.q0*rm(1)*rm(2)*rm2(1)8 r2= r(1)*r(1)+r(2)*r(2)0 ri(1)= r(1)/r2 ri(2)= -r(2)/r2  rr2= r(1)*r(1)-r(2)*r(2)+ r22= rr2*rr2+4.q0*r(1)*r(1)*r(2)*r(2)  r2i(1)= rr2/r22-! r2i(2)= -2.q0*r(1)*r(2)/r22  x2= x(1)*x(1)+x(2)*x(2)  xi(1)= x(1)/x2 xi(2)= -x(2)/x2( xx2= x(1)*x(1)-x(2)*x(2)+ x22= xx2*xx2+4.q0*x(1)*x(1)*x(2)*x(2). x2i(1)= xx2/x22i! x2i(2)= -2.q0*x(1)*x(2)/x22 *  rar(1)= 1.q0-r(1)/r2 rar(2)= r(2)/r2h if(rar(1).gt.0.q0) then0 rrar(1)= sqrt(rar(1))+ rrar(2)= 0.5q0*rar(2)/sqrt(rar(1))1" else if(rar(1).lt.0.q0) then, rrar(1)= 0.5q0*rar(2)/sqrt(-rar(1)) rrar(2)= sqrt(-rar(1)) endif* , xp2= (1.q0-x(1))*(1.q0-x(1))+x(2)*x(2) xp(1)= (1.q0-x(1))/xp2 xp(2)= x(2)/xp2 - rpp2= (1.q0-r(1))*(1.q0-r(1))+r(2)*r(2)* rpp(1)= (1.q0-r(1))/rpp2 rpp(2)= r(2)/rpp2l*o ar1r= rm2(1) ar1i= rm2(2) umar1r= 1.q0-rm2(1)o ar2r= rm4(1) ar2i= rm4(2) umar2r= 1.q0-rm4(1)z ar3r= xp(1)  ar3i= xp(2)c umar3r= 1.q0-xp(1)2 call tspence(ar1r,ar1i,umar1r,cli2sr,cli2si)2 call tspence(ar2r,ar2i,umar2r,cli2fr,cli2fi)2 call tspence(ar3r,ar3i,umar3r,cli2xr,cli2xi)) call tcli3(ar1r,ar1i,cli3sr,cli3si)z) call tcli3(ar2r,ar2i,cli3fr,cli3fi) ) call tcli3(ar3r,ar3i,cli3xr,cli3xi)f*) comb3r= 2.q0*cli3sr-cli3fr comb3i= 2.q0*cli3si-cli3fi comb2r= cli2sr-cli2fr  comb2i= cli2si-cli2fi  f2r= f(1)*f(1)-f(2)*f(2) f2i= 2.q0*f(1)*f(2)#*  r1p= r(1)+1.5q0*+ aux1r= -f(1)+g(1)/3.q0+2.q0/3.q0*h(1) + aux1i= -f(2)+g(2)/3.q0+2.q0/3.q0*h(2)z+ aux2r= -3.q0*f(1)+2.q0*g(1)+4.q0*h(1)o+ aux2i= -3.q0*f(2)+2.q0*g(2)+4.q0*h(2) ) aux3r= comb2r+f(1)*aux2r-f(2)*aux2i ( aux3pr= -2.q0*(r1p*f(1)-r(2)*f(2))) aux3i= comb2i+f(1)*aux2i+f(2)*aux2r ( aux3pi= -2.q0*(r1p*f(2)+r(2)*f(1))) aux4r= comb2r+f(1)*aux2r-f(2)*aux2i 3 aux4pr= -2.q0*((r(1)-3.q0+0.25q0*ri(1))*f(1)- ' # (r(2)+0.25q0*ri(2))*f(2))n) aux4i= comb2i+f(1)*aux2i+f(2)*aux2r)3 aux4pi= -2.q0*((r(1)-3.q0+0.25q0*ri(1))*f(2)+1' # (r(2)+0.25q0*ri(2))*f(1)) , aux5r= r(1)-1.q0/6.q0-7.q0/48.q0*ri(1)" aux5i= r(2)-7.q0/48.q0*ri(2): aux6r= r(1)-11.q0/12.q0+5.q0/48.q0*ri(1)+1.q0/32.q0* # r2i(1))4 aux6i= r(2)+5.q0/48.q0*ri(2)+1.q0/32.q0*r2i(2)*o% ex0vr= 4.q0*(r(1)-0.25q0*ri(1))z% ex0vi= 4.q0*(r(2)-0.25q0*ri(2)),*z* ex0ar= 4.q0*(r(1)-1.5q0+0.5q0*ri(1))$ ex0ai= 4.q0*(r(2)+0.5q0*ri(2))*$ ex0xr= x(1)-1.5q0+0.5q0*x2i(1) ex0xi= x(2)+0.5q0*x2i(2)*d7 ex1r= comb3r+8.q0/3.q0*(f(1)*comb2r-f(2)*comb2i)+a& # 4.q0*(f2r*aux1r-f2i*aux1i)7 ex1i= comb3i+8.q0/3.q0*(f(1)*comb2i+f(2)*comb2r)+(& # 4.q0*(f2r*aux1i+f2i*aux1r)*)< ex2r= 8.q0/3.q0*((r(1)+0.5q0)*aux3r-r(2)*aux3i)+aux3pr< ex2i= 8.q0/3.q0*((r(1)+0.5q0)*aux3i+r(2)*aux3r)+aux3pi*-; ex3r= -8.q0*(aux5r*f2r-aux5i*f2i)+13.q0/6.q0+z3*ri(1)n0 ex3i= -8.q0*(aux5r*f2i+aux5i*f2r)+z3*ri(2)*-; ex4r= 8.q0/3.q0*((r(1)-1.q0)*aux4r-r(2)*aux4i)+aux4pro; ex4i= 8.q0/3.q0*((r(1)-1.q0)*aux4i+r(2)*aux4r)+aux4pir*e( ex5r= -8.q0*(aux6r*f2r-aux6i*f2i)+5 # 13.q0/6.q0-3.q0*z2+(0.25q0-2.q0*z3)*ri(1)s9 ex5i= -8.q0*(aux6r*f2i+aux6i*f2r)+(0.25q0-2.q0*z3)*l # ri(2) * A ex6r= cli3xr+2.q0/3.q0*(b(1)*cli2xr-b(2)*cli2xi)-1.q0/6.q0*01 # (b2(1)*(a(1)-b(1))-b2(2)*(a(2)-b(2)))9A ex6i= cli3xi+2.q0/3.q0*(b(1)*cli2xi+b(2)*cli2xr)-1.q0/6.q0*61 # (b2(1)*(a(2)-b(2))+b2(2)*(a(1)-b(1))) * ? ex7r= 1.q0/3.q0*((x(1)+0.5q0-0.5q0*xi(1))*(cli2xr-ab(1))-0@ # (x(2)-0.5q0*xi(2))*(cli2xi-ab(2)))+1.q0/3.q0*(b2(1)*: # (x(1)-1.q0/8.q0-xi(1)+5.q0/8.q0*x2i(1))-b2(2)*8 # (x(2)-xi(2)+5.q0/8.q0*x2i(2)))-0.25q0*(b(1)*> # (x(1)-5.q0/2.q0+2.q0/3.q0*xi(1)+5.q0/6.q0*x2i(1))-9 # b(2)*(x(2)+2.q0/3.q0*xi(2)+5.q0/6.q0*x2i(2))))? ex7i= 1.q0/3.q0*((x(1)+0.5q0-0.5q0*xi(1))*(cli2xi-ab(2))+k@ # (x(2)-0.5q0*xi(2))*(cli2xr-ab(1)))+1.q0/3.q0*(b2(1)*@ # (x(2)-xi(2)+5.q0/8.q0*x2i(2))+b2(2)*(x(1)-1.q0/8.q0-B # xi(1)+5.q0/8.q0*x2i(1)))-0.25q0*(b(1)*(x(2)+2.q0/3.q0*C # xi(2)+5.q0/6.q0*x2i(2))+b(2)*(x(1)-5.q0/2.q0+2.q0/3.q0*k$ # xi(1)+5.q0/6.q0*x2i(1)))* F ex8r= -3.q0/4.q0*z2+13.q0/12.q0-5.q0/24.q0*xi(1)-0.5q0*z3*x2i(1)- ex8i= -5.q0/24.q0*xi(2)-0.5q0*z3*x2i(2)* ? v1r= ex0vr*ex1r-ex0vi*ex1i+rrar(1)*ex2r-rrar(2)*ex2i+ex3r ? v1i= ex0vr*ex1i+ex0vi*ex1r+rrar(1)*ex2i+rrar(2)*ex2r+ex3i? a1r= ex0ar*ex1r-ex0ai*ex1i+rrar(1)*ex4r-rrar(2)*ex4i+ex5rz? a1i= ex0ar*ex1i+ex0ai*ex1r+rrar(1)*ex4i+rrar(2)*ex4r+ex5ic* f1x= ex0xr*ex6r-ex0xi*ex6i+ex7r+ex8r+ f1xi= ex0xr*ex6i+ex0xi*ex6r+ex7i+ex8i)*z' ex0vpr= 4.q0*(1.q0+0.25q0*r2i(1))  ex0vpi= r2i(2)& ex0apr= 4.q0*(1.q0-0.5q0*r2i(1)) ex0api= -2.q0*r2i(2)* * bux1r= -comb2r-f(1)*aux2r+f(2)*aux2i* bux1i= -comb2i-f(1)*aux2i-f(2)*aux2r0 bux2r= 2.q0/3.q0*(ri(1)*bux1r-ri(2)*bux1i)0 bux2i= 2.q0/3.q0*(ri(1)*bux1i+ri(2)*bux1r)/ bux3r= bux2r-(1.q0-5.q0/6.q0*ri(1))*f(1)-! # 5.q0/6.q0*ri(2)*f(2))/ bux3i= bux2i-(1.q0-5.q0/6.q0*ri(1))*f(2)+ ' # 5.q0/6.q0*ri(2)*f(1) m4 bux4r= f2r*(4.q0/3.q0*ri(1)-11.q0/6.q0*r2i(1)-: # 4.q0*rpp(1))-f2i*(4.q0/3.q0*ri(2)-11.q0/6.q0*; # r2i(2)-4.q0*rpp(2))-1.q0-1.5q0*ri(1)-z3*r2i(1) 4 bux4i= f2r*(4.q0/3.q0*ri(2)-11.q0/6.q0*r2i(2)-: # 4.q0*rpp(2))+f2i*(4.q0/3.q0*ri(1)-11.q0/6.q0*6 # r2i(1)-4.q0*rpp(1))-1.5q0*ri(2)-z3*r2i(2)1 cux2r= 4.q0/3.q0*(-ri(1)*bux1r+ri(2)*bux1i)1 cux2i= 4.q0/3.q0*(-ri(1)*bux1i-ri(2)*bux1r) = cux3r= cux2r-(1.q0+13.q0/6.q0*ri(1)-0.5q0*r2i(1))*f(1)+)1 # (13.q0/6.q0*ri(2)-0.5q0*r2i(2))*f(2) = cux3i= cux2i-(1.q0+13.q0/6.q0*ri(1)-0.5q0*r2i(1))*f(2)-) # (13.q0/6.q0*ri(2)-0.5q0*r2i(2))*f(1) cux4r= f2r*ri(1)-f2i*ri(2) cux4i= f2r*ri(2)+f2i*ri(1) cux5r= cux4r*(-20.q0/3.q0+13.q0/6.q0*ri(1)+0.5q0*r2i(1))- # cux4i*(13.q0/6.q0*ri(2)+0.5q0*r2i(2))-1.q0+3.q0*ri(1)+ # (2.q0*z3-0.5q0)*r2i(1) cux5i= cux4r*(13.q0/6.q0*ri(2)+0.5q0*r2i(2))+ # cux4i*(20.q0/3.q0+13.q0/6.q0*ri(1)+0.5q0*r2i(1))+ # 3.q0*ri(2)+(2.q0*z3-0.5q0)*r2i(2) * v1pr= ex0vpr*ex1r-ex0vpi*ex1i+2.q0*(rrar(1)*bux3r- # rrar(2)*bux3i)+bux4r a1pr= ex0apr*ex1r-ex0api*ex1i+2.q0*(rrar(1)*cux3r- # rrar(2)*cux3i)+cux5r v1pi= ex0vpr*ex1i+ex0vpi*ex1r+2.q0*(rrar(1)*bux3i+ # rrar(2)*bux3r)+bux4i a1pi= ex0apr*ex1i+ex0api*ex1r+2.q0*(rrar(1)*cux3i+ # rrar(2)*cux3r)+cux5i * ov1(1)= v1r ov1(2)= v1i oa1(1)= a1r oa1(2)= a1i of1(1)= f1x of1(2)= f1xi ov1p(1)= v1pr ov1p(2)= v1pi oa1p(1)= a1pr oa1p(2)= a1pi * return end * *-----cli3----------------------------------------------------------- * computes li_3(x) for complex x * subroutine tcli3(xr,xi,cli3r,cli3i) implicit real*16(a-h,o-z) * common/tqparam/qpi,qpis,qeps,qdelta * dimension b(0:14),bf(0:14) dimension x(2),y(2),addx(2),ox(2),clnx(2),par(2),ct(15),sn(15), # res(2),u1(2),u2(2),clny(2),omy(2),clnomy(2),addy(2), # par1(2),par2(2),ct1(15),sn1(15),ct2(15),sn2(15), # res1(2),res2(2),t(2),resa(2),resb(2),clnt(2), # res3(2),res4(2),clnomt(2),addt(2),addt2(2),omt(2), # omu1(2),omu2(2) * data b/1.q0,-0.75q0,0.236111111111111111111111111111111q0, # -3.472222222222222222222222222222222q-2, # 6.481481481481481481481481481481482q-4, # 4.861111111111111111111111111111111q-4, # -2.393550012597631645250692869740488q-5, # -1.062925170068027210884353741496599q-5, # 7.794784580498866213151927437641723q-7, # 2.526087595532039976484420928865373q-7, # -2.359163915200471237027273583310139q-8, # -6.168132746415574698402981231264060q-9, # 6.824456748981078267312315451125495q-10, # 1.524285616929084572552216019859487q-10, # -1.916909414174054295837274763110831q-11/ data z3/1.20205690315959428540q0/ * z2= qpis/6.q0 do n=0,14 bf(n)= b(n)/(n+1.q0) enddo * x(1)= xr x(2)= xi * xm2= x(1)*x(1)+x(2)*x(2) xm= sqrt(xm2) * *-----the modulus of x is checked * xtst= xm-1.q0 if(xtst.le.1.q-33) then y(1)= x(1) y(2)= x(2) addx(1)= 0.q0 addx(2)= 0.q0 else if(xm.gt.1.q-33) then y(1)= x(1)/xm2 y(2)= -x(2)/xm2 ox(1)= -x(1) ox(2)= -x(2) call tcqlnx(ox,clnx) rlnx2= clnx(1)*clnx(1) ailnx2= clnx(2)*clnx(2) addx(1)= -clnx(1)*(z2+1.q0/6.q0*(rlnx2-3.q0*ailnx2)) addx(2)= -clnx(2)*(z2+1.q0/6.q0*(3.q0*rlnx2-ailnx2)) endif * *-----once x --> y, |y|<1 the sign of re(y) is checked * if re(y)>0 a transformation is required for re(y)>1/2 * y2r= y(1)*y(1)-y(2)*y(2) if(y(1).ge.0.q0.or.y2r.lt.0.q0) then ytst= y(1)-0.5q0 if(ytst.le.1.q-33) then * *-----li_3(y) is computed * omy(1)= 1.q0-y(1) omy(2)= -y(2) call tcqlnomx(y,omy,par) pr= -par(1) pi= -par(2) p2= pr*pr+pi*pi pm= sqrt(p2) ct(1)= pr/pm sn(1)= pi/pm do n=2,15 ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1) sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1) enddo res(1)= pm*(bf(0)*ct(1)+pm*(bf(1)*ct(2)+pm* # (bf(2)*ct(3)+pm*(bf(3)*ct(4)+pm* # (bf(4)*ct(5)+pm*(bf(5)*ct(6)+pm* # (bf(6)*ct(7)+pm*(bf(7)*ct(8)+pm* # (bf(8)*ct(9)+pm*(bf(9)*ct(10)+pm* # (bf(10)*ct(11)+pm*(bf(11)*ct(12)+pm* # (bf(12)*ct(13)+pm*(bf(13)*ct(14)+pm* # (bf(14)*ct(15)))))))))))))))) res(2)= pm*(bf(0)*sn(1)+pm*(bf(1)*sn(2)+pm* # (bf(2)*sn(3)+pm*(bf(3)*sn(4)+pm* # (bf(4)*sn(5)+pm*(bf(5)*sn(6)+pm* # (bf(6)*sn(7)+pm*(bf(7)*sn(8)+pm* # (bf(8)*sn(9)+pm*(bf(9)*sn(10)+pm* # (bf(10)*sn(11)+pm*(bf(11)*sn(12)+pm* # (bf(12)*sn(13)+pm*(bf(13)*sn(14)+pm* # (bf(14)*sn(15)))))))))))))))) cli3r= res(1)+addx(1) cli3i= res(2)+addx(2) return else if(ytst.gt.1.q-33) then ym2= y(1)*y(1)+y(2)*y(2) u1(1)= 1.q0-y(1)/ym2 u1(2)= y(2)/ym2 u2(1)= 1.q0-y(1) u2(2)= -y(2) call tcqlnx(y,clny) omy(1)= 1.q0-y(1) omy(2)= -y(2) call tcqlnomx(y,omy,clnomy) addy(1)= z3+z2*clny(1)+1.q0/6.d0*clny(1)* # (clny(1)*clny(1)-3.q0*clny(2)*clny(2))- # 0.5q0*clnomy(1)*(clny(1)*clny(1)-clny(2)* # clny(2))+clny(1)*clny(2)*clnomy(2) addy(2)= z2*clny(2)+1.q0/6.q0*clny(2)*(3.q0* # clny(1)*clny(1)-clny(2)*clny(2))-0.5q0* # clnomy(2)*(clny(1)*clny(1)-clny(2)*clny(2))- # clny(1)*clnomy(1)*clny(2) * *-----li_3(1-1/y) is computed * omu1(1)= 1.q0-u1(1) omu1(2)= -u1(2) call tcqlnomx(u1,omu1,par1) pr1= -par1(1) pi1= -par1(2) p12= pr1*pr1+pi1*pi1 pm1= sqrt(p12) ct1(1)= pr1/pm1 sn1(1)= pi1/pm1 do n=2,15 ct1(n)= ct1(1)*ct1(n-1)-sn1(1)*sn1(n-1) sn1(n)= sn1(1)*ct1(n-1)+ct1(1)*sn1(n-1) enddo res1(1)= pm1*(bf(0)*ct1(1)+pm1*(bf(1)*ct1(2)+pm1* # (bf(2)*ct1(3)+pm1*(bf(3)*ct1(4)+pm1* # (bf(4)*ct1(5)+pm1*(bf(5)*ct1(6)+pm1* # (bf(6)*ct1(7)+pm1*(bf(7)*ct1(8)+pm1* # (bf(8)*ct1(9)+pm1*(bf(9)*ct1(10)+pm1* # (bf(10)*ct1(11)+pm1*(bf(11)*ct1(12)+pm1* # (bf(12)*ct1(13)+pm1*(bf(13)*ct1(14)+pm1* # (bf(14)*ct1(15)))))))))))))))) res1(2)= pm1*(bf(0)*sn1(1)+pm1*(bf(1)*sn1(2)+pm1* # (bf(2)*sn1(3)+pm1*(bf(3)*sn1(4)+pm1* # (bf(4)*sn1(5)+pm1*(bf(5)*sn1(6)+pm1* # (bf(6)*sn1(7)+pm1*(bf(7)*sn1(8)+pm1* # (bf(8)*sn1(9)+pm1*(bf(9)*sn1(10)+pm1* # (bf(10)*sn1(11)+pm1*(bf(11)*sn1(12)+pm1* # (bf(12)*sn1(13)+pm1*(bf(13)*sn1(14)+pm1* # (bf(14)*sn1(15)))))))))))))))) * *-----li_3(1-y) is computed * omu2(1)= 1.q0-u2(1) omu2(2)= -u2(2) call tcqlnomx(u2,omu2,par2) pr2= -par2(1) pi2= -par2(2) p22= pr2*pr2+pi2*pi2 pm2= sqrt(p22) ct2(1)= pr2/pm2 sn2(1)= pi2/pm2 do n=2,15 ct2(n)= ct2(1)*ct2(n-1)-sn2(1)*sn2(n-1) sn2(n)= sn2(1)*ct2(n-1)+ct2(1)*sn2(n-1) enddo res2(1)= pm2*(bf(0)*ct2(1)+pm2*(bf(1)*ct2(2)+pm2* # (bf(2)*ct2(3)+pm2*(bf(3)*ct2(4)+pm2* # (bf(4)*ct2(5)+pm2*(bf(5)*ct2(6)+pm2* # (bf(6)*ct2(7)+pm2*(bf(7)*ct2(8)+pm2* # (bf(8)*ct2(9)+pm2*(bf(9)*ct2(10)+pm2* # (bf(10)*ct2(11)+pm2*(bf(11)*ct2(12)+pm2* # (bf(12)*ct2(13)+pm2*(bf(13)*ct2(14)+pm2* # (bf(14)*ct2(15)))))))))))))))) res2(2)= pm2*(bf(0)*sn2(1)+pm2*(bf(1)*sn2(2)+pm2* # (bf(2)*sn2(3)+pm2*(bf(3)*sn2(4)+pm2* # (bf(4)*sn2(5)+pm2*(bf(5)*sn2(6)+pm2* # (bf(6)*sn2(7)+pm2*(bf(7)*sn2(8)+pm2* # (bf(8)*sn2(9)+pm2*(bf(9)*sn2(10)+pm2* # (bf(10)*sn2(11)+pm2*(bf(11)*sn2(12)+pm2* # (bf(12)*sn2(13)+pm2*(bf(13)*sn2(14)+pm2* # (bf(14)*sn2(15)))))))))))))))) cli3r= -res1(1)-res2(1)+addx(1)+addy(1) cli3i= -res1(2)-res2(2)+addx(2)+addy(2) return endif * *-----if re(y)<0 a transformation is required in terms of t = -y * and of t^2 * else if(y(1).lt.0.q0) then * *-----first t * t(1)= -y(1) t(2)= -y(2) if(t(1).le.0.5q0) then * *-----li_3(t) is computed * omt(1)= 1.q0-t(1) omt(2)= -t(2) call tcqlnomx(t,omt,par) pr= -par(1) pi= -par(2) p2= pr*pr+pi*pi pm= sqrt(p2) ct(1)= pr/pm sn(1)= pi/pm do n=2,15 ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1) sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1) enddo resa(1)= pm*(bf(0)*ct(1)+pm*(bf(1)*ct(2)+pm* # (bf(2)*ct(3)+pm*(bf(3)*ct(4)+pm* # (bf(4)*ct(5)+pm*(bf(5)*ct(6)+pm* # (bf(6)*ct(7)+pm*(bf(7)*ct(8)+pm* # (bf(8)*ct(9)+pm*(bf(9)*ct(10)+pm* # (bf(10)*ct(11)+pm*(bf(11)*ct(12)+pm* # (bf(12)*ct(13)+pm*(bf(13)*ct(14)+pm* # (bf(14)*ct(15)))))))))))))))) resa(2)= pm*(bf(0)*sn(1)+pm*(bf(1)*sn(2)+pm* # (bf(2)*sn(3)+pm*(bf(3)*sn(4)+pm* # (bf(4)*sn(5)+pm*(bf(5)*sn(6)+pm* # (bf(6)*sn(7)+pm*(bf(7)*sn(8)+pm* # (bf(8)*sn(9)+pm*(bf(9)*sn(10)+pm* # (bf(10)*sn(11)+pm*(bf(11)*sn(12)+pm* # (bf(12)*sn(13)+pm*(bf(13)*sn(14)+pm* # (bf(14)*sn(15)))))))))))))))) else if(t(1).gt.0.5q0) then tm2= t(1)*t(1)+t(2)*t(2) u1(1)= 1.q0-t(1)/tm2 u1(2)= t(2)/tm2 u2(1)= 1.q0-t(1) u2(2)= -t(2) call tcqlnx(t,clnt) omt(1)= 1.q0-t(1) omt(2)= -t(2) call tcqlnomx(t,omt,clnomt) addt(1)= z3+z2*clnt(1)+1.q0/6.d0*clnt(1)* # (clnt(1)*clnt(1)-3.q0*clnt(2)*clnt(2))- # 0.5q0*clnomt(1)*(clnt(1)*clnt(1)-clnt(2)* # clnt(2))+clnt(1)*clnt(2)*clnomt(2) addt(2)= z2*clnt(2)+1.q0/6.q0*clnt(2)*(3.q0* # clnt(1)*clnt(1)-clnt(2)*clnt(2))-0.5q0* # clnomt(2)*(clnt(1)*clnt(1)-clnt(2)*clnt(2))- # clnt(1)*clnomt(1)*clnt(2) * *-----li3(1-1/t) is computed * omu1(1)= 1.q0-u1(1) omu1(2)= -u1(2) call tcqlnomx(u1,omu1,par1) pr1= -par1(1) pi1= -par1(2) p12= pr1*pr1+pi1*pi1 pm1= sqrt(p12) ct1(1)= pr1/pm1 sn1(1)= pi1/pm1 do n=2,15 ct1(n)= ct1(1)*ct1(n-1)-sn1(1)*sn1(n-1) sn1(n)= sn1(1)*ct1(n-1)+ct1(1)*sn1(n-1) enddo res1(1)= pm1*(bf(0)*ct1(1)+pm1*(bf(1)*ct1(2)+pm1* # (bf(2)*ct1(3)+pm1*(bf(3)*ct1(4)+pm1* # (bf(4)*ct1(5)+pm1*(bf(5)*ct1(6)+pm1* # (bf(6)*ct1(7)+pm1*(bf(7)*ct1(8)+pm1* # (bf(8)*ct1(9)+pm1*(bf(9)*ct1(10)+pm1* # (bf(10)*ct1(11)+pm1*(bf(11)*ct1(12)+pm1* # (bf(12)*ct1(13)+pm1*(bf(13)*ct1(14)+pm1* # (bf(14)*ct1(15)))))))))))))))) res1(2)= pm1*(bf(0)*sn1(1)+pm1*(bf(1)*sn1(2)+pm1* # (bf(2)*sn1(3)+pm1*(bf(3)*sn1(4)+pm1* # (bf(4)*sn1(5)+pm1*(bf(5)*sn1(6)+pm1* # (bf(6)*sn1(7)+pm1*(bf(7)*sn1(8)+pm1* # (bf(8)*sn1(9)+pm1*(bf(9)*sn1(10)+pm1* # (bf(10)*sn1(11)+pm1*(bf(11)*sn1(12)+pm1* # (bf(12)*sn1(13)+pm1*(bf(13)*sn1(14)+pm1* # (bf(14)*sn1(15)))))))))))))))) * *-----li3(1-t) is computed * omu2(1)= 1.q0-u2(1) omu2(2)= -u2(2) call tcqlnomx(u2,omu2,par2) pr2= -par2(1) pi2= -par2(2) p22= pr2*pr2+pi2*pi2 pm2= sqrt(p22) ct2(1)= pr2/pm2 sn2(1)= pi2/pm2 do n=2,15 ct2(n)= ct2(1)*ct2(n-1)-sn2(1)*sn2(n-1) sn2(n)= sn2(1)*ct2(n-1)+ct2(1)*sn2(n-1) enddo res2(1)= pm2*(bf(0)*ct2(1)+pm2*(bf(1)*ct2(2)+pm2* # (bf(2)*ct2(3)+pm2*(bf(3)*ct2(4)+pm2* # (bf(4)*ct2(5)+pm2*(bf(5)*ct2(6)+pm2* # (bf(6)*ct2(7)+pm2*(bf(7)*ct2(8)+pm2* # (bf(8)*ct2(9)+pm2*(bf(9)*ct2(10)+pm2* # (bf(10)*ct2(11)+pm2*(bf(11)*ct2(12)+pm2* # (bf(12)*ct2(13)+pm2*(bf(13)*ct2(14)+pm2* # (bf(14)*ct2(15)))))))))))))))) res2(2)= pm2*(bf(0)*sn2(1)+pm2*(bf(1)*sn2(2)+pm2* # (bf(2)*sn2(3)+pm2*(bf(3)*sn2(4)+pm2* # (bf(4)*sn2(5)+pm2*(bf(5)*sn2(6)+pm2* # (bf(6)*sn2(7)+pm2*(bf(7)*sn2(8)+pm2* # (bf(8)*sn2(9)+pm2*(bf(9)*sn2(10)+pm2* # (bf(10)*sn2(11)+pm2*(bf(11)*sn2(12)+pm2* # (bf(12)*sn2(13)+pm2*(bf(13)*sn2(14)+pm2* # (bf(14)*sn2(15)))))))))))))))) resa(1)= -res1(1)-res2(1)+addt(1) resa(2)= -res1(2)-res2(2)+addt(2) endif * *-----then t^2 * t(1)= y(1)*y(1)-y(2)*y(2) t(2)= 2.q0*y(1)*y(2) if(t(1).le.0.5q0) then * *-----li_3(t^2) is computed * omt(1)= 1.q0-t(1) omt(2)= -t(2) call tcqlnomx(t,omt,par) pr= -par(1) pi= -par(2) p2= pr*pr+pi*pi pm= sqrt(p2) ct(1)= pr/pm sn(1)= pi/pm do n=2,15 ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1) sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1) enddo resb(1)= pm*(bf(0)*ct(1)+pm*(bf(1)*ct(2)+pm* # (bf(2)*ct(3)+pm*(bf(3)*ct(4)+pm* # (bf(4)*ct(5)+pm*(bf(5)*ct(6)+pm* # (bf(6)*ct(7)+pm*(bf(7)*ct(8)+pm* # (bf(8)*ct(9)+pm*(bf(9)*ct(10)+pm* # (bf(10)*ct(11)+pm*(bf(11)*ct(12)+pm* # (bf(12)*ct(13)+pm*(bf(13)*ct(14)+pm* # (bf(14)*ct(15)))))))))))))))) resb(2)= pm*(bf(0)*sn(1)+pm*(bf(1)*sn(2)+pm* # (bf(2)*sn(3)+pm*(bf(3)*sn(4)+pm* # (bf(4)*sn(5)+pm*(bf(5)*sn(6)+pm* # (bf(6)*sn(7)+pm*(bf(7)*sn(8)+pm* # (bf(8)*sn(9)+pm*(bf(9)*sn(10)+pm* # (bf(10)*sn(11)+pm*(bf(11)*sn(12)+pm* # (bf(12)*sn(13)+pm*(bf(13)*sn(14)+pm* # (bf(14)*sn(15)))))))))))))))) else if(t(1).gt.0.5q0) then tm2= t(1)*t(1)+t(2)*t(2) u1(1)= 1.q0-t(1)/tm2 u1(2)= t(2)/tm2 u2(1)= 1.q0-t(1) u2(2)= -t(2) call tcqlnx(t,clnt) omt(1)= 1.q0-t(1) omt(2)= -t(2) call tcqlnomx(t,omt,clnomt) addt2(1)= z3+z2*clnt(1)+1.q0/6.d0*clnt(1)*