* *-------------------------------------------------------------------------- * subroutine zztoregion(n,x,j,a,b) implicit real*8(a-h,o-z) * dimension x(n) * if(j.gt.1) then do i=1,j-1 if(x(i).lt.0.d0.or.x(i).gt.1.d0) then stop endif enddo endif * a= 0.d0 b= 1.d0 * return end * *------------------------------------------------------------------- * subroutine zztocorrqcd(scal,zm,als,bqm2,cqm2,rbm2,rcm2,rsm2) implicit real*8(a-h,o-z) * common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * data rz3/1.20205690315959428540d0/,nf/5/ * scal2= scal*scal bqm= sqrt(bqm2) cqm= sqrt(cqm2) * *-----computes the running charm mass *-----computes the running bottom mass * scal1g= 1.d0 als1g= zztoralphas(zm,scal1g,als,nf)/pi alsc4= zztoralphas(zm,cqm,als,nf)/pi alsc42= alsc4*alsc4 als1g2= als1g*als1g alsb= zztoralphas(zm,bqm,als,nf)/pi alsb2= alsb*alsb alsz= zztoralphas(zm,scal,als,nf)/pi alsz2= alsz*alsz * *-----first the running of the s-quark mass * up to c/b-threshold * fn= 3.d0 b03= (11.d0-2.d0/3.d0*fn)/4.d0 b13= (102.d0-38.d0/3.d0*fn)/16.d0 b23= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0* # fn*fn)/64.d0 g03= 1.d0 g13= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0 g23= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*rz3)*fn- # 140.d0/81.d0*fn*fn)/64.d0 rex3= g03/b03 b03s= b03*b03 b03c= b03s*b03 b13s= b13*b13 sm1g= 0.189d0 alsdf= alsc4-als1g cfm13= g13/b03-b13*g03/b03s cfm23= g23/b03-b13*g13/b03s-b23*g03/b03s+ # b13s*g03/b03c rsmc= sm1g*(alsc4/als1g)**rex3*(1.d0+cfm13* # alsdf+0.5d0*cfm13*cfm13*alsdf*alsdf+0.5d0*cfm23* # (alsc42-als1g2)) * fn= 4.d0 b04= (11.d0-2.d0/3.d0*fn)/4.d0 b14= (102.d0-38.d0/3.d0*fn)/16.d0 b24= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0* # fn*fn)/64.d0 g04= 1.d0 g14= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0 g24= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*rz3)*fn- # 140.d0/81.d0*fn*fn)/64.d0 rex4= g04/b04 b04s= b04*b04 b04c= b04s*b04 b14s= b14*b14 alsdf= alsb-alsc4 cfm14= g14/b04-b14*g04/b04s cfm24= g24/b04-b14*g14/b04s-b24*g04/b04s+ # b14s*g04/b04c rsmb= rsmc*(alsb/alsc4)**rex4*(1.d0+cfm14* # alsdf+0.5d0*cfm14*cfm14*alsdf*alsdf+0.5d0*cfm24* # (alsb2-alsc42)) * *-----c quark mass at c-threshold * zero= 0.d0 x1= 0.5d0 x2= 2.0d0 xacc= 1.d-12 fn= 4.d0 cmm4= zztorrunm(x1,x2,xacc,cqm,alsc4,rsmc,zero,fn) * *-----c quark mass at b-threshold * cmb= cmm4*(alsb/alsc4)**rex4*(1.d0+cfm14* # (alsb-alsc4)+0.5d0*cfm14*cfm14* # (alsb-alsc4)**2+0.5d0*cfm24*(alsb2-alsc42)) * *-----running c mass * fn= 5.d0 b05= (11.d0-2.d0/3.d0*fn)/4.d0 b15= (102.d0-38.d0/3.d0*fn)/16.d0 b25= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0* # fn*fn)/64.d0 g05= 1.d0 g15= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0 g25= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*rz3)*fn- # 140.d0/81.d0*fn*fn)/64.d0 rex5= g05/b05 b05s= b05*b05 b05c= b05s*b05 b15s= b15*b15 cfm15= g15/b05-b15*g05/b05s cfm25= g25/b05-b15*g15/b05s-b25*g05/b05s+b15s*g05/b05c rcm= cmb*(alsz/alsb)**rex5*(1.d0+cfm15* # (alsz-alsb)+0.5d0*cfm15*cfm15* # (alsz-alsb)**2+0.5d0*cfm25*(alsz2-alsb2)) rcm2= rcm*rcm * *-----running s mass * rsm= rsmb*(alsz/alsb)**rex5*(1.d0+cfm15* # (alsz-alsb)+0.5d0*cfm15*cfm15* # (alsz-alsb)**2+0.5d0*cfm25*(alsz2-alsb2)) rsm2= rsm*rsm * *-----b quark mass at b-threshold * x1= 0.5d0 x2= 6.0d0 xacc= 1.d-12 fn= 5.d0 bmm5= zztorrunm(x1,x2,xacc,bqm,alsb,cmb,rsmb,fn) * *-----running b mass * alsdf= alsz-alsb rbm= bmm5*(alsz/alsb)**rex5*(1.d0+cfm15* # alsdf+0.5d0*cfm15*cfm15*alsdf*alsdf+0.5d0*cfm25* # (alsz2-alsb2)) rbm2= rbm*rbm * return end * *------------------------------------------------------------------- * real*8 function zztoralphas(rs0,rs,als,nf) implicit real*8(a-h,o-z) * common/zztmass/em,rmm,tm,uqm,dqm,cqm,sqm,bqm common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * dimension b0(5),b1(5),b2(5) * *-----limits for lambda-5 are 1 mev < lambda_5 < 10 gev * if(als.eq.0.d0) then zztoralphas= 0.d0 else x1= 0.001d0 x2= 10.0d0 xacc= 1.d-12 qcdl= zztoqcdlam(nf,als,rs0,x1,x2,xacc) pqcdl5= qcdl do i=1,5 b0(i)= (11.d0-2.d0/3.d0*i)/4.d0 b1(i)= (102.d0-38.d0/3.d0*i)/16.d0 b2(i)= 0.5d0*(2857.d0-i*(5033.d0/9.d0- # 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) fac2= b1(5)*b1(5)/b0(5)/b0(5)-b1(4)*b1(4)/b0(4)/b0(4) 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) rb= log(b0(4)/b0(3)) fac= b1(4)/b0(4)-b1(3)/b0(3) facp= b2(4)/b0(4)-b2(3)/b0(3) 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) ratl2= exp(rhs) qcdl= qcdl/sqrt(ratl2) pqcdl3= qcdl nfe= nf-2 endif else nfe= nf endif * qcdb0= 11.d0-2.d0/3.d0*nfe qcdb1= 102.d0-38.d0/3.d0*nfe qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nfe+ # 325.d0/27.d0*nfe*nfe) qcda= 2.d0*log(rs/qcdl) * zztoralphas= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda* # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)- # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0)) * endif * return end * *-------------------------------------------------------------------- * real*8 function zztorals(qcdl,rs,nf) implicit real*8(a-h,o-z) * common/zztmass/em,rmm,tm,uqm,dqm,cqm,sqm,bqm common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * dimension b0(5),b1(5),b2(5) * pqcdl5= qcdl do i=1,5 b0(i)= (11.d0-2.d0/3.d0*i)/4.d0 b1(i)= (102.d0-38.d0/3.d0*i)/16.d0 b2(i)= 0.5d0*(2857.d0-i*(5033.d0/9.d0- # 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) fac2= b1(5)*b1(5)/b0(5)/b0(5)-b1(4)*b1(4)/b0(4)/b0(4) 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) rb= log(b0(4)/b0(3)) fac= b1(4)/b0(4)-b1(3)/b0(3) facp= b2(4)/b0(4)-b2(3)/b0(3) 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) ratl2= exp(rhs) qcdl= qcdl/sqrt(ratl2) pqcdl3= qcdl nfe= nf-2 endif else nfe= nf endif * qcdb0= 11.d0-2.d0/3.d0*nfe qcdb1= 102.d0-38.d0/3.d0*nfe qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nfe+ # 325.d0/27.d0*nfe*nfe) qcda= 2.d0*log(rs/qcdl) * zztorals= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda* # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)- # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0)) * return end * *---------------------------------------------------------- * real*8 function zztorrunm(x1,x2,xacc,qm,als,rm1,rm2,fn) implicit real*8(a-h,o-z) * parameter (jmax=50) * fmid= zztoqcdmass(qm,als,rm1,rm2,fn,x2) f= zztoqcdmass(qm,als,rm1,rm2,fn,x1) if (f*fmid.ge.0.d0) then print*,'root must be bracketed for bisection' print 1,qm 1 format(/' error detected by rrunm ',/ # ' current value of quark mass = ',e20.5) stop endif if (f.lt.0.d0) then zztorrunm= x1 dx= x2-x1 else zztorrunm= x2 dx= x1-x2 endif do j=1,jmax dx= dx*0.5d0 xmid= zztorrunm+dx fmid= zztoqcdmass(qm,als,rm1,rm2,fn,xmid) if (fmid.le.0.d0) zztorrunm= xmid if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return enddo pause 'too many bisections' end * *-----zztoqcdmass------------------------------------------------------ * real*8 function zztoqcdmass(qm,als,rm1,rm2,fn,x) implicit real*8(a-h,o-z) * common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * rz2= pis/6.d0 rz3= 1.20205690315959428540d0 ra4= 0.5174790617d0 rz5= 1.03692775514336992633d0 rln= 2.d0*log(qm/x) rlns= rln*rln r1= rm1/x r2= rm2/x delta0= 3.d0/4.d0*rz2-3.d0/8.d0 delta1= r1*(pis/8.d0+r1*(-0.597d0+0.230d0*r1)) delta2= r2*(pis/8.d0+r2*(-0.597d0+0.230d0*r2)) * rhs= 1.d0+als*(4.d0/3.d0+rln+als*(3817.d0/288.d0-8.d0/3.d0+ # 2.d0/3.d0*(2.d0+log(2.d0))*rz2-rz3/6.d0-fn/3.d0*(rz2+ # 71.d0/48.d0)+4.d0/3.d0*(delta0+delta1+delta2)+(173.d0/ # 24.d0-13.d0/36.d0*fn)*rln+(15.d0/8.d0-fn/12.d0)*rlns)) zztoqcdmass= qm-x*rhs * return end * *--------------------------------------------------------------------- * real*8 function zztoqcdlam(nf,als,rs,x1,x2,xacc) implicit real*8(a-h,o-z) * parameter (jmax=50,nout=21) * fmid= zztoqcdscale(nf,als,rs,x2) f= zztoqcdscale(nf,als,rs,x1) if (f*fmid.ge.0.d0) then print*,'root must be bracketed for bisection' print 1,als,x1,x2 1 format(/' error detected by zztoqcdlam ',/ # ' current value of alpha_s = ',e20.5,/ # ' interval ',2e20.5) stop endif if (f.lt.0.d0) then zztoqcdlam= x1 dx= x2-x1 else zztoqcdlam= x2 dx= x1-x2 endif do j=1,jmax dx= dx*0.5d0 xmid= zztoqcdlam+dx fmid= zztoqcdscale(nf,als,rs,xmid) if (fmid.le.0.d0) zztoqcdlam= xmid if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return enddo pause 'too many bisections' end * *------------------------------------------------------------------ * real*8 function zztoqcdscale(nf,als,rs,x) implicit real*8(a-h,o-z) * common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * qcdb0= 11.d0-2.d0/3.d0*nf qcdb1= 102.d0-38.d0/3.d0*nf qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nf+325.d0/27.d0*nf*nf) qcda= 2.d0*log(rs/x) zztoqcdscale= als-(4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda* # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)- # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0))) * return end * *------------------------------------------------------------- * subroutine zztofst(nt,st,fvect,iflag) implicit real*8(a-h,o-z) character*1 opqcd character*2 oval * common/zztomu/sclmu common/zztswq/opqcd common/zztalstz/alstm,alszm common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai common/zztmixz/v1z,a1z,df1z,v1pz,a1pz,v1iz,a1iz,f1iz common/zztmixw/v1w,a1w,df1w,v1pw,a1pw,v1iw,a1iw,f1iw common/zztfmass2/em2,rmm2,tm2,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),fwz(2) * tqm2= st(1) tqm= sqrt(tqm2) if(opqcd.eq.'y') then nf= 5 stqm= sclmu*tqm szm= sclmu*zm alstm= zztoralphas(zm,stqm,alsz,nf) if(sclmu.eq.1.d0) then alszm= alsz else alszm= zztoralphas(zm,szm,alsz,nf) endif pwr= -wm2 pwi= 0.d0 call zztoalals(pwr,pwi,v1w,a1w,df1w,v1pw,a1pw,v1iw,a1iw,f1iw) pzr= -zm2 pzi= 0.d0 call zztoalals(pzr,pzi,v1z,a1z,df1z,v1pz,a1pz,v1iz,a1iz,f1iz) endif * *-----alpha(zm) * pzr= -zm2 pzi= 0.d0 oval= 'rz' call zztopoleg(oval,pzr,pzi,pggfz,pgglqz,pggnpz) call zztopolez0(pggf0,fw0) x= alphai-0.25d0*(pggfz(1)-pggf0+pggnpz)/pi y= -0.25d0*(pggfz(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 oval= 'rz' call zztopole(oval,p2zr,p2zi,p3qz,fzz,fwz) 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 oval= 'rw' call zztopole(oval,p2r,p2i,p3qw,fzw,fw) * fvect(1)= 1.d0+0.5d0*gf*wm2/pis*((fw0-fw(1))/wm2+p3qw(1))- # 8.d0*gf*wm2*(g2z+p3qz(1)/apis) * return end * *-------------------------------------------------------------------- * subroutine zztopoleg(oval,p2r,p2i,pggf,pgglq,pggnp) implicit real*8(a-h,o-z) character*1 opqcd character*2 oval * data z3/1.20205690315959428540d0/ * common/zztomu/sclmu common/zztswq/opqcd common/zzted/eps,delta common/zztalstz/alstm,alszm common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai common/zztmixr/v1r,a1r,df1r,v1pr,a1pr,v1ir,a1ir,f1ir common/zztmixz/v1z,a1z,df1z,v1pz,a1pz,v1iz,a1iz,f1iz common/zztmixw/v1w,a1w,df1w,v1pw,a1pw,v1iw,a1iw,f1iw common/zztmixsz/v1sz,a1sz,df1sz,v1psz,a1psz,v1isz,a1isz,f1isz common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2 * dimension pggf(2),pgglq(2) dimension r(2),fb(2),flq(2) dimension bfl(2),bfh(2),bflq(2) dimension b0e(2),b1e(2),b21e(2) dimension b0m(2),b1m(2),b21m(2) dimension b0t(2),b1t(2),b21t(2) dimension b0h(2),b1h(2),b21h(2) dimension b0lq(2),b1lq(2),b21lq(2) * call zztosbff(p2r,p2i,em2,em2,b0e,b1e,b21e) call zztosbff(p2r,p2i,rmm2,rmm2,b0m,b1m,b21m) call zztosbff(p2r,p2i,tm2,tm2,b0t,b1t,b21t) call zztosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h) call zztosbff(p2r,p2i,uqm2,uqm2,b0lq,b1lq,b21lq) * do i=1,2 bfl(i)= 2.d0*b21e(i)-b0e(i)+ # 2.d0*b21m(i)-b0m(i)+ # 2.d0*b21t(i)-b0t(i) bfh(i)= 2.d0*b21h(i)-b0h(i) bflq(i)= 2.d0*b21lq(i)-b0lq(i) enddo * do i=1,2 pggf(i)= 4.d0*bfl(i)+16.d0/3.d0*bfh(i) pgglq(i)= 44.d0/3.d0*bflq(i) enddo * if(opqcd.eq.'y') then qs2= sclmu*sclmu*tqm2 aexpt= alstm/pi/3.d0 aexp= alszm/pi/3.d0 rm= tqm2/qs2 rl= log(rm) bx= -rl-4.d0*z3+55.d0/12.d0 rvr= -0.25d0*p2r/tqm2 xvr= -p2r/tqm2 if(p2i.eq.0.d0) then rvi= 0.25d0*eps/tqm2 xvi= eps/tqm2 else rvi= -0.25d0*p2i/tqm2 xvi= -p2i/tqm2 endif * r(1)= p2r/(sclmu*sclmu)/tqm2 if(p2i.eq.0.d0) then r(2)= -eps/(sclmu*sclmu)/tqm2 else r(2)= p2i/(sclmu*sclmu)/tqm2 endif call zztocqlnx(r,fb) if(p2i.ne.0.d0) then fb(2)= fb(2)-2.d0*pi endif r(1)= p2r/(sclmu*sclmu)/zm2 if(p2i.eq.0.d0) then r(2)= -eps/(sclmu*sclmu)/zm2 else r(2)= p2i/(sclmu*sclmu)/zm2 endif call zztocqlnx(r,flq) if(p2i.ne.0.d0) then flq(2)= flq(2)-2.d0*pi endif fb(1)= -z3+55.d0/48.d0-0.25d0*fb(1) fb(2)= -0.25d0*fb(2) flq(1)= -z3+55.d0/48.d0-0.25d0*flq(1) flq(2)= -0.25d0*flq(2) p2m2= p2r*p2r+p2i*p2i if(oval.eq.'rw') then pggf(1)= pggf(1)+64.d0/3.d0*aexpt*tqm2/p2m2* # (p2r*(rvr*bx+v1w)+p2i*(rvi*bx+v1iw)) pggf(2)= pggf(2)+64.d0/3.d0*aexpt*tqm2/p2m2* # (p2r*(rvi*bx+v1iw)-p2i*(rvr*bx+v1w)) else if(oval.eq.'rz') then pggf(1)= pggf(1)+64.d0/3.d0*aexpt*tqm2/p2m2* # (p2r*(rvr*bx+v1z)+p2i*(rvi*bx+v1iz)) pggf(2)= pggf(2)+64.d0/3.d0*aexpt*tqm2/p2m2* # (p2r*(rvi*bx+v1iz)-p2i*(rvr*bx+v1z)) else if(oval.eq.'pz') then pggf(1)= pggf(1)+64.d0/3.d0*aexpt*tqm2/p2m2* # (p2r*(rvr*bx+v1sz)+p2i*(rvi*bx+v1isz)) pggf(2)= pggf(2)+64.d0/3.d0*aexpt*tqm2/p2m2* # (p2r*(rvi*bx+v1isz)-p2i*(rvr*bx+v1sz)) else if(oval.eq.'ru') then pggf(1)= pggf(1)+64.d0/3.d0*aexpt*tqm2/p2m2* # (p2r*(rvr*bx+v1r)+p2i*(rvi*bx+v1ir)) pggf(2)= pggf(2)+64.d0/3.d0*aexpt*tqm2/p2m2* # (p2r*(rvi*bx+v1ir)-p2i*(rvr*bx+v1r)) endif do i= 1,2 pgglq(i)= pgglq(i)-16.d0*aexpt*fb(i)-160.d0/3.d0* # aexp*flq(i) enddo endif * 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 * pggnp= 4.d0*pi*alphai*(a+b*log(1.d0+c*ap2x)) * return end * *------------------------------------------------------------------------- * subroutine zztopolez0(pggf0,fw0) implicit real*8(a-h,o-z) character*1 opqcd * data z3/1.20205690315959428540d0/ * common/zztomu/sclmu common/zztswq/opqcd common/zzted/eps,delta common/zztalstz/alstm,alszm common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2 * call zztorbff0(em2,em2,b0ee0,b1ee0,b21ee0) call zztorbff0(rmm2,rmm2,b0mm0,b1mm0,b21mm0) call zztorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0) call zztorbff0(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 * pggf0= 4.d0*(bfe+bfm+bftau)+16.d0/3.d0*bft * call zztorbff0(tqm2,dmy2,b0tb0,b1tb0,b21tb0) * fw0= -3.d0*tqm2*(b0tb0+b1tb0) * if(opqcd.eq.'y') then qs2= sclmu*sclmu*tqm2 aexps= alstm/pi/3.d0 rm= tqm2/qs2 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= pggf0-16.d0/3.d0*aexps*(bx+4.d0*z3-5.d0/6.d0) fw0= fw0+12.d0*aexps*tqm2*(0.25d0*by+f10) endif * return end * *------------------------------------------------------------------- * subroutine zztopole(oval,p2r,p2i,p3q,fz,fw) implicit real*8(a-h,o-z) character*1 opqcd character*2 oval * data z3/1.20205690315959428540d0/ data pi/3.141592653589793238462643d0/ * common/zztomu/sclmu common/zztswq/opqcd common/zzted/eps,delta common/zztalstz/alstm,alszm common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztmixr/v1r,a1r,df1r,v1pr,a1pr,v1ir,a1ir,f1ir common/zztmixz/v1z,a1z,df1z,v1pz,a1pz,v1iz,a1iz,f1iz common/zztmixw/v1w,a1w,df1w,v1pw,a1pw,v1iw,a1iw,f1iw common/zztmixsw/v1sw,a1sw,df1sw,v1psw,a1psw,v1isw,a1isw,f1isw common/zztmixsz/v1sz,a1sz,df1sz,v1psz,a1psz,v1isz,a1isz,f1isz common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2 * dimension bfl(2),bfh(2) dimension r(2),fb(2),flq(2) dimension p3q(2),fz(2),fw(2) dimension b0l(2),b1l(2),b21l(2) dimension b0h(2),b1h(2),b21h(2) dimension b0c(2),b1c(2),b21c(2) * pis= pi*pi * call zztosbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l) call zztosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h) call zztosbff(p2r,p2i,tqm2,dmy2,b0c,b1c,b21c) * do i=1,2 bfl(i)= 2.d0*b21l(i)-b0l(i) bfh(i)= 2.d0*b21h(i)-b0h(i) enddo bfsr= bfh(1)-bfl(1) bfsi= bfh(2)-bfl(2) * do i=1,2 p3q(i)= 10.d0*bfl(i)+2.d0*bfh(i) enddo * fz(1)= -0.5d0*(p2r*bfsr-p2i*bfsi)-3.d0/2.d0*tqm2*b0h(1) fz(2)= -0.5d0*(p2r*bfsi+p2i*bfsr)-3.d0/2.d0*tqm2*b0h(2) * auxr= 3.d0*(b21c(1)+b1c(1))-2.d0*b21h(1)+b0h(1)-b21l(1)- # 3.d0*b1l(1)-b0l(1) auxi= 3.d0*(b21c(2)+b1c(2))-2.d0*b21h(2)+b0h(2)-b21l(2)- # 3.d0*b1l(2)-b0l(2) * fw(1)= 2.d0*(p2r*auxr-p2i*auxi)-3.d0*tqm2*(b1c(1)+b0c(1)) fw(2)= 2.d0*(p2r*auxi+p2i*auxr)-3.d0*tqm2*(b1c(2)+b0c(2)) * if(opqcd.eq.'y') then qs2= sclmu*sclmu*tqm2 aexps= alstm/pi/3.d0 aexpsz= alszm/pi/3.d0 rm= tqm2/qs2 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 xvr= -p2r/tqm2 if(p2i.eq.0.d0) then rvi= 0.25d0*eps/tqm2 xvi= eps/tqm2 else rvi= -0.25d0*p2i/tqm2 xvi= -p2i/tqm2 endif * sclmus= sclmu*sclmu r(1)= p2r/tqm2/sclmus if(p2i.eq.0.d0) then r(2)= -eps/tqm2/sclmus else r(2)= p2i/tqm2/sclmus endif call zztocqlnx(r,fb) if(p2i.ne.0.d0) then fb(2)= fb(2)-2.d0*pi endif if(oval.eq.'rz') then v1rd= v1z v1id= v1iz a1rd= a1z a1id= a1iz f1rd= df1z f1id= f1iz else if(oval.eq.'rw') then v1rd= v1w v1id= v1iw a1rd= a1w a1id= a1iw f1rd= df1w f1id= f1iw else if(oval.eq.'pw') then v1rd= v1sw v1id= v1isw a1rd= a1sw a1id= a1isw f1rd= df1sw f1id= f1isw else if(oval.eq.'pz') then v1rd= v1sz v1id= v1isz a1rd= a1sz a1id= a1isz f1rd= df1sz f1id= f1isz else if(oval.eq.'ru') then v1rd= v1r v1id= v1ir a1rd= a1r a1id= a1ir f1rd= df1r f1id= f1ir endif r(1)= p2r/zm2/sclmus if(p2i.eq.0.d0) then r(2)= -eps/zm2/sclmus else r(2)= p2i/zm2/sclmus endif call zztocqlnx(r,flq) if(p2i.ne.0.d0) then flq(2)= flq(2)-2.d0*pi endif fb(1)= -z3+55.d0/48.d0-0.25d0*fb(1) fb(2)= -0.25d0*fb(2) flq(1)= -z3+55.d0/48.d0-0.25d0*flq(1) flq(2)= -0.25d0*flq(2) * p2m2= p2r*p2r+p2i*p2i p3q(1)= p3q(1)-24.d0*aexpsz*flq(1)-12.d0/3.d0*aexps*fb(1)+ # 24.d0/3.d0*aexps*tqm2/p2m2*(p2r*(rvr*bx+v1rd)+ # p2i*(rvi*bx+v1id)) p3q(2)= p3q(2)-24.d0*aexpsz*flq(2)-12.d0/3.d0*aexps*fb(2)+ # 24.d0/3.d0*aexps*tqm2/p2m2*(p2r*(rvi*bx+v1id)- # p2i*(rvr*bx+v1rd)) fw(1)= fw(1)+12.d0*tqm2*aexps*(0.25d0*(xvr*bx+by)+f1rd)- # 24.d0/3.d0*tqm2*aexps*(rvr*bx+v1rd)+12.d0/3.d0* # aexps*(p2r*fb(1)-p2i*fb(2)) fw(2)= fw(2)+12.d0*tqm2*aexps*(0.25d0*xvi*bx+f1id)- # 24.d0/3.d0*tqm2*aexps*(rvi*bx+v1id)+12.d0/3.d0* # aexps*(p2r*fb(2)+p2i*fb(1)) fz(1)= fz(1)+aexps*tqm2*(-5.d0*(rvr*bx+v1rd)+3.d0* # (rvr*bx+by+a1rd))-2.d0*aexps*(p2r*fb(1)-p2i*fb(2)) fz(2)= fz(2)+aexps*tqm2*(-5.d0*(rvi*bx+v1id)+3.d0* # (rvi*bx+a1id))-2.d0*aexps*(p2r*fb(2)+p2i*fb(1)) endif * return end * *------------------------------------------------------------------ * subroutine zztofsw(n,sw,fvecw,iflag) implicit real*8(a-h,o-z) character*1 opqcd character*2 oval * common/zztomu/sclmu common/zztswq/opqcd common/zztalstz/alstm,alszm common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai common/zztmixw/v1w,a1w,df1w,v1pw,a1pw,v1iw,a1iw,f1iw common/zztmixsw/v1sw,a1sw,df1sw,v1psw,a1psw,v1isw,a1isw,f1isw common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2 * dimension sw(n),fvecw(n) dimension p3qw(2),fzw(2),fw(2),p3q(2),fz(2),fcw(2) * wm2= wm*wm swr= sw(1) swi= sw(2) wmu= sqrt(swr) wlg= -swi/wmu * if(opqcd.eq.'y') then nf= 5 tqm= sqrt(tqm2) stqm= sclmu*tqm szm= sclmu*zm alstm= zztoralphas(zm,stqm,alsz,nf) if(sclmu.eq.1.d0) then alszm= alsz else alszm= zztoralphas(zm,szm,alsz,nf) endif pwr= -wm2 pwi= 0.d0 call zztoalals(pwr,pwi,v1w,a1w,df1w,v1pw,a1pw,v1iw,a1iw,f1iw) pwr= -swr pwi= -swi call zztoalals(pwr,pwi,v1sw,a1sw,df1sw,v1psw,a1psw, # v1isw,a1isw,f1isw) endif * call zztopolez0(pggf0,fw0) p2r= -wm2 p2i= 0.d0 oval= 'rw' call zztopole(oval,p2r,p2i,p3qw,fzw,fw) p2r= -swr p2i= -swi oval= 'pw' call zztopole(oval,p2r,p2i,p3q,fz,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 zztofsz(n,sz,fvecz,iflag) implicit real*8(a-h,o-z) character*1 opqcd character*2 oval * common/zztafsz/g2z common/zztomu/sclmu common/zztswq/opqcd common/zztalstz/alstm,alszm common/zztiaz/alzr,alzi,alzir,alzii common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai common/zztmixz/v1z,a1z,df1z,v1pz,a1pz,v1iz,a1iz,f1iz common/zztmixsz/v1sz,a1sz,df1sz,v1psz,a1psz,v1isz,a1isz,f1isz common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2 * dimension sz(n),fvecz(n) dimension p3qz(2),fzz(2),p3q(2),fz(2),fwz(2),fw(2) dimension pggfz(2),pgglqz(2),pggf(2),pgglq(2),evl(2) * szr= sz(1) szi= sz(2) * if(opqcd.eq.'y') then nf= 5 tqm= sqrt(tqm2) stqm= sclmu*tqm szm= sclmu*zm alstm= zztoralphas(zm,stqm,alsz,nf) if(sclmu.eq.1.d0) then alszm= alsz else alszm= zztoralphas(zm,szm,alsz,nf) endif pzr= -zm2 pzi= 0.d0 call zztoalals(pzr,pzi,v1z,a1z,df1z,v1pz,a1pz,v1iz,a1iz,f1iz) pzr= -szr pzi= -szi call zztoalals(pzr,pzi,v1sz,a1sz,df1sz,v1psz,a1psz, # v1isz,a1isz,f1isz) endif apis= 16.d0*pis call zztopolez0(pggf0,fw0) pzr= -zm2 pzi= 0.d0 oval= 'rz' call zztopoleg(oval,pzr,pzi,pggfz,pgglqz,pggnpz) call zztopole(oval,pzr,pzi,p3qz,fzz,fwz) p2r= -szr p2i= -szi oval= 'pz' call zztopoleg(oval,p2r,p2i,pggf,pgglq,pggnp) call zztopole(oval,p2r,p2i,p3q,fz,fw) do i=1,2 evl(i)= pggfz(i)-pggf(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 * *------------------------------------------------------------------------- * subroutine zztosbff(p2r,p2i,rm12,rm22,b0,b1,b21) implicit real*8(a-h,o-z) * common/zzted/eps,delta common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * dimension cmp2(2) dimension clnmp2(2) dimension b0(2),b1(2),b21(2) * if(rm12.ne.0.d0.and.rm12.eq.rm22.and. # (p2r.gt.0.d0.and.sqrt(p2r).lt.(2.d-5*sqrt(rm12)))) then call zztorbff0(rm12,rm22,ab0,ab1,ab21) b0(1)= ab0-p2r/6.d0/rm12 b1(1)= ab1+p2r/12.d0/rm12 b21(1)= ab21-p2r/20.d0/rm12 b0(2)= 0.d0 b1(2)= 0.d0 b21(2)= 0.d0 return endif * p2m= p2r*p2r+p2i*p2i * if(rm12.eq.0.d0.and.rm22.eq.0.d0) then if(p2i.eq.0) then cmp2(1)= p2r cmp2(2)= -eps call zztoclnx(cmp2,clnmp2) b0(1)= delta+2.d0-clnmp2(1) b0(2)= -clnmp2(2) else cmp2(1)= p2r cmp2(2)= p2i call zztoclnx(cmp2,clnmp2) b0(1)= delta+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*(-delta-1.d0+log(rm12)) if(p2i.eq.0) then cmp2(1)= 1.d0+p2r/rm12 cmp2(2)= -eps/rm12 call zztoclnx(cmp2,clnmp2) b0(1)= delta+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 zztoclnx(cmp2,clnmp2) xr= p2r/p2m xi= -p2i/p2m fctr= 1.d0+rm12*xr fcti= rm12*xi b0(1)= delta+2.d0-log(rm12)-fctr*clnmp2(1)+fcti*clnmp2(2) b0(2)= -fctr*clnmp2(2)-fcti*clnmp2(1) if(p2r.lt.(-rm12)) then b0(1)= b0(1)+2.d0*pi*rm12*p2i/p2m b0(2)= b0(2)+2.d0*pi*(1.d0+rm12*p2r/p2m) endif x2r= xr*xr-xi*xi x2i= 2.d0*xr*xi b1(1)= 0.5d0*a0*xr-0.5d0*b0(1)+0.5d0*rm12*(xr*b0(1)- # xi*b0(2)) b1(2)= 0.5d0*a0*xi-0.5d0*b0(2)+0.5d0*rm12*(xr*b0(2)+ # xi*b0(1)) b21(1)= 3.d0/18.d0*rm12*xr+1.d0/18.d0+a0/3.d0* # (rm12*x2r-xr)+1.d0/3.d0*((1.d0-rm12*xr+ # rm12*rm12*x2r)*b0(1)-rm12*(-xi+rm12*x2i)*b0(2)) b21(2)= 3.d0/18.d0*rm12*xi+a0/3.d0* # (rm12*x2i-xi)+1.d0/3.d0*((1.d0-rm12*xr+ # rm12*rm12*x2r)*b0(2)+rm12*(-xi+rm12*x2i)*b0(1)) endif else if(rm12.ne.0.d0.and.rm12.eq.rm22) then a0= rm12*(-delta-1.d0+log(rm12)) if(p2i.eq.0) then p2d= p2r*p2r+eps*eps b2r= 1.d0+4.d0*rm12*p2r/p2d b2i= 4.d0*rm12*eps/p2d call a02aaf(b2r,b2i,br,bi) cmp2(1)= (br*br+bi*bi-1.d0)/((br-1.d0)*(br-1.d0)+bi*bi) cmp2(2)= -2.d0*bi/((br-1.d0)*(br-1.d0)+bi*bi) call zztoclnx(cmp2,clnmp2) b0(1)= delta+2.d0-log(rm12)-(br*clnmp2(1)- # bi*clnmp2(2)) b0(2)= -(br*clnmp2(2)+bi*clnmp2(1)) b1(1)= -0.5d0*b0(1) b1(2)= -0.5d0*b0(2) b21(1)= (6.d0*rm12+p2r)/18.d0/p2r+a0/3.d0/p2r+ # 1.d0/3.d0*(1.d0+rm12/p2r)*b0(1) b21(2)= 1.d0/3.d0*(1.d0+rm12/p2r)*b0(2) else xr= p2r/p2m xi= -p2i/p2m b2r= 1.d0+4.d0*rm12*xr b2i= 4.d0*rm12*xi call a02aaf(b2r,b2i,br,bi) cmp2(1)= (br*br+bi*bi-1.d0)/((br-1.d0)*(br-1.d0)+bi*bi) cmp2(2)= -2.d0*bi/((br-1.d0)*(br-1.d0)+bi*bi) call zztoclnx(cmp2,clnmp2) b0(1)= delta+2.d0-log(rm12)-(br*clnmp2(1)- # bi*clnmp2(2)) b0(2)= -(br*clnmp2(2)+bi*clnmp2(1)) if(p2r.lt.(-4.d0*rm12)) then b0(1)= b0(1)-2.d0*pi*bi b0(2)= b0(2)+2.d0*pi*br endif x2r= xr*xr-xi*xi x2i= 2.d0*xr*xi b1(1)= -0.5d0*b0(1) b1(2)= -0.5d0*b0(2) b21(1)= rm12*xr/3.d0+1.d0/18.d0+a0/3.d0*xr+ # 1.d0/3.d0*((1.d0+rm12*xr)*b0(1)-rm12*xi* # b0(2)) b21(2)= rm12*xi/3.d0+a0/3.d0*xi+1.d0/3.d0*((1.d0+ # rm12*xr)*b0(2)+rm12*xi*b0(1)) endif else a01= rm12*(-delta-1.d0+log(rm12)) a02= rm22*(-delta-1.d0+log(rm22)) if(p2i.eq.0) then b2r= (p2r-rm22+rm12)*(p2r-rm22+rm12)- # 4.d0*rm12*rm22-eps*eps b2i= -2.d0*(p2r+rm22+rm12)*eps call a02aaf(b2r,b2i,br,bi) cmp2(1)= 0.5d0*(-p2r-rm22-rm12+br)/sqrt(rm22*rm12) cmp2(2)= 0.5d0*(eps+bi)/sqrt(rm22*rm12) call zztoclnx(cmp2,clnmp2) crr= (br*clnmp2(1)-bi*clnmp2(2))/p2r cri= (br*clnmp2(2)+bi*clnmp2(1))/p2r b0(1)= delta-crr-0.5d0*log(rm22*rm12)+0.5d0*(rm22-rm12)/ # p2r*log(rm22/rm12)+2.d0 b0(2)= -cri b1(1)= 0.5d0/p2r*(a01-a02+(rm12-rm22-p2r)*b0(1)) b1(2)= 0.5d0/p2r*(rm12-rm22-p2r)*b0(2) b21(1)= (3.d0*(rm12+rm22)+p2r)/18.d0/p2r+ # (rm12-rm22-p2r)*a01/3.d0/(p2r*p2r)+ # (rm22-rm12+2.d0*p2r)*a02/3.d0/(p2r*p2r)+ # (1.d0+1.d0/p2r*(2.d0*rm22-rm12+(rm22-rm12)**2/ # p2r))*b0(1) b21(2)= (1.d0+1.d0/p2r*(2.d0*rm22-rm12+(rm22-rm12)**2/ # p2r))*b0(2) else xr= p2r/p2m xi= -p2i/p2m x2r= xr*xr-xi*xi x2i= 2.d0*xr*xi b2r= (p2r-rm22+rm12)*(p2r-rm22+rm12)- # 4.d0*rm12*rm22-p2i*p2i b2i= 2.d0*(p2r+rm22+rm12)*p2i call a02aaf(b2r,b2i,br,bi) cmp2(1)= 0.5d0*(-p2r-rm22-rm12+br)/sqrt(rm22*rm12) cmp2(2)= 0.5d0*(-p2i+bi)/sqrt(rm22*rm12) call zztoclnx(cmp2,clnmp2) crr= (br*clnmp2(1)-bi*clnmp2(2))*xr- # (br*clnmp2(2)+bi*clnmp2(1))*xi cri= (br*clnmp2(1)-bi*clnmp2(2))*xi+ # (br*clnmp2(2)+bi*clnmp2(1))*xr b0(1)= delta+crr-0.5d0*log(rm22*rm12)+0.5d0*(rm22-rm12)/ # p2r*log(rm22/rm12)+2.d0 b0(2)= cri tsti= -(sqrt(rm12)+sqrt(rm22))*(sqrt(rm12)+sqrt(rm22)) if(p2r.lt.tsti) then b0(1)= b0(1)-2.d0*pi*(br*xi+bi*xr) b0(2)= b0(2)+2.d0*pi*(br*xr-bi*xi) endif b1(1)= 0.5d0*xr*(a01-a02+(rm12-rm22)*b0(1))- # 0.5d0*xi*(rm12-rm22)*b0(2)-0.5d0*b0(1) b1(2)= 0.5d0*xi*(a01-a02+(rm12-rm22)*b0(1))+ # 0.5d0*xr*(rm12-rm22)*b0(2)-0.5d0*b0(2) b21(1)= (rm12+rm22)/6.d0*xr+1.d0/18.d0+ # (rm12-rm22)/3.d0*x2r*a01-a01/3.d0*xr+ # (rm22-rm12)/3.d0*a02*x2r+2.d0/3.d0*a02*xr+ # b0(1)/3.d0+(2.d0*rm22-rm12)/3.d0*(xr*b0(1)- # xi*b0(2))+(rm22-rm12)**2/3.d0*(x2r*b0(1)-x2i*b0(2)) b21(2)= (rm12+rm22)/6.d0*xi+ # (rm12-rm22)/3.d0*x2i*a01-a01/3.d0*xi+ # (rm22-rm12)/3.d0*a02*x2i+2.d0/3.d0*a02*xi+ # b0(2)/3.d0+(2.d0*rm22-rm12)/3.d0*(xr*b0(2)+ # xi*b0(1))+(rm22-rm12)**2/3.d0*(x2r*b0(2)+x2i*b0(1)) endif endif * return end * *------------------------------------------------------------------- * subroutine zztorbff0(rm12,rm22,b0,b1,b21) implicit real*8(a-h,o-z) * common/zzted/eps,delta common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * 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= delta-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 zztorcg(n,yr,yi,omyr,fr) arg(1)= omyr arg(2)= -yi call zztoclnx(arg,cln) f1r= fr(1)-cln(1) f2r= 2.d0*fr(2)-cln(1) f3r= 3.d0*fr(3)-cln(1) fact= delta-log(arm22) b0= fact-f1r b1= 0.5d0*(-fact+f2r) b21= 1.d0/3.d0*(fact-f3r) return endif end * *---------------------------------------------------------- * subroutine zztoclnx(arg,res) implicit real*8(a-h,o-z) * common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * dimension arg(2),aarg(2),res(2) * do i= 1,2 aarg(i)= abs(arg(i)) enddo azm2= (arg(1))**2+(arg(2))**2 azm= sqrt(azm2) res(1)= log(azm) if (arg(1).eq.0.d0) then if (arg(2).gt.0.d0) then teta= pi/2.d0 else teta= -pi/2.d0 endif res(2)= teta return else if (arg(2).eq.0.d0) then if (arg(1).gt.0.d0) then teta= 0.d0 else teta= pi endif res(2)= teta return else tnteta= aarg(2)/aarg(1) teta= atan(tnteta) sr= arg(1)/aarg(1) si= arg(2)/aarg(2) if (sr.gt.0.d0) then res(2)= si*teta else res(2)= si*(pi-teta) endif return endif end * *--------------------------------------------------------- * subroutine zztoclnomx(arg,omarg,res) implicit real*8(a-h,o-z) * dimension arg(2),omarg(2),res(2),ares(2),ct(10),sn(10) * zr= arg(1) zi= arg(2) azm2= zr*zr+zi*zi azm= sqrt(azm2) if (azm.lt.1.d-7) then ct(1)= zr/azm sn(1)= zi/azm do n=2,10 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 ares(1)= ct(10)/10.d0 ares(2)= sn(10)/10.d0 do k=9,1,-1 ares(1)= ares(1)*azm+ct(k)/k ares(2)= ares(2)*azm+sn(k)/k enddo ares(1)= -ares(1)*azm ares(2)= -ares(2)*azm else call zztoclnx(omarg,ares) endif do i= 1,2 res(i)= ares(i) enddo * return end * *----------------------------------------------------- * subroutine zztorcg(n,zr,zi,omzr,gfr) implicit real*8(a-h,o-z) * common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * dimension gfr(n) dimension a(4),ra(16),z(2),omz(2),oz(2),clnomz(2),clnoz(2) dimension ca(2,10),aux(2),zp(2),ct(16),sn(16) data ra/-0.2d0,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/ * z(1)= zr z(2)= zi omz(1)= omzr omz(2)= -zi if (zr.eq.0.d0.and.zi.eq.0.d0) then 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 a(j)= ((j-1.d0)*a(j-1)-1.d0/j)/j enddo do k=1,n gfr(k)= a(k) enddo return else zmod2= zr*zr+zi*zi zmod= sqrt(zmod2) call zztoclnx(omz,clnomz) * * |z| < 4 * if (zmod.lt.4.d0) then oz(1)= -z(1) oz(2)= -z(2) call zztoclnx(oz,clnoz) ca(1,1)= omz(1)*clnomz(1)-omz(2)*clnomz(2)+z(1)*clnoz(1)- # z(2)*clnoz(2)-1.d0 ca(2,1)= omz(1)*clnomz(2)+omz(2)*clnomz(1)+z(1)*clnoz(2)+ # z(2)*clnoz(1) if (n.eq.1) then gfr(1)= ca(1,1) return else do j= 2,n jm1= j-1 ca(1,j)= ((j-1.d0)*(z(1)*ca(1,jm1)-z(2)*ca(2,jm1))+ # omz(1)*clnomz(1)-omz(2)*clnomz(2)-1.d0/j)/j ca(2,j)= ((j-1.d0)*(z(1)*ca(2,jm1)+z(2)*ca(1,jm1))+ # omz(1)*clnomz(2)+omz(2)*clnomz(1))/j enddo do k=1,n gfr(k)= ca(1,k) enddo return endif * * |z| > 4 * else aux(1)= (-zr*omzr+zi**2)/zmod2 aux(2)= zi/zmod2 call zztoclnx(aux,zp) zpm2= zp(1)*zp(1)+zp(2)*zp(2) zpm= sqrt(zpm2) ct(1)= zp(1)/zpm sn(1)= zp(2)/zpm 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) enddo ca(1,4)= ca(1,4)*zpm ca(2,4)= ca(2,4)*zpm do j= 3,1,-1 jp1= j+1 ca(1,j)= ((ca(1,jp1)+1.d0/jp1)*z(1)+ca(2,jp1)*z(2))/zmod2 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 return endif endif end * *---------------------------------------------------------- * subroutine zztocfl(p2r,p2i,p3qf,s3qb,pgg,s33,sww,r1,r2,r3) implicit real*8(a-h,o-z) * common/zzted/eps,delta common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2 * dimension bfl(2),bfh(2),co(2) dimension b0l(2),b1l(2),b21l(2) dimension b0h(2),b1h(2),b21h(2) dimension b0c(2),b1c(2),b21c(2) dimension b0w(2),b1w(2),b21w(2) dimension b0z(2),b1z(2),b21z(2) dimension b0zh(2),b1zh(2),b21zh(2) dimension b0zw(2),b1zw(2),b21zw(2) dimension b0pw(2),b1pw(2),b21pw(2) dimension b0wh(2),b1wh(2),b21wh(2) dimension pggf(2),pggb(2),p3qf(2),s3qb(2),s3q(2),pgg(2),r1(2), # r2(2),bcl(2),bch(2),s33f(2),s33(2),sr33(2),swwf(2), # p33(2),pwwb(2),srwwb(2),sww(2),r3(2) * co(1)= 1.d0 co(2)= 0.d0 * call zztosbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l) call zztosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h) call zztosbff(p2r,p2i,tqm2,dmy2,b0c,b1c,b21c) * 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 pggf(i)= 80.d0/3.d0*bfl(i)+16.d0/3.d0*bfh(i) enddo * do i=1,2 p3qf(i)= 10.d0*bfl(i)+2.d0*bfh(i) enddo * call zztosbff(p2r,p2i,wm2,wm2,b0w,b1w,b21w) call zztosbff(p2r,p2i,zm2,zm2,b0z,b1z,b21z) call zztosbff(p2r,p2i,zm2,hm2,b0zh,b1zh,b21zh) call zztosbff(p2r,p2i,zm2,wm2,b0zw,b1zw,b21zw) call zztosbff(p2r,p2i,wm2,hm2,b0wh,b1wh,b21wh) call zztosbff(p2r,p2i,dmy2,wm2,b0pw,b1pw,b21pw) * s3qb(1)= p2r*(2.d0/3.d0-10.d0*b21w(1)+13.d0/2.d0*b0w(1))- # p2i*(-10.d0*b21w(2)+13.d0/2.d0*b0w(2))- # 2.d0*wm2*b0w(1)+2.d0*wm2*(delta-log(wm2)) s3qb(2)= p2i*(2.d0/3.d0-10.d0*b21w(1)+13.d0/2.d0*b0w(1))+ # p2r*(-10.d0*b21w(2)+13.d0/2.d0*b0w(2))- # 2.d0*wm2*b0w(2) * pggb(1)= 2.d0/3.d0-12.d0*b21w(1)+7.d0*b0w(1) pggb(2)= -12.d0*b21w(2)+7.d0*b0w(2) * s3q(1)= p2r*p3qf(1)-p2i*p3qf(2)+s3qb(1) s3q(2)= p2r*p3qf(2)+p2i*p3qf(1)+s3qb(2) * do i=1,2 pgg(i)= pggf(i)+pggb(i) enddo * r1(1)= s3q(1)-5.d0/46.d0*(p2r*pgg(1)-p2i*pgg(2)) r1(2)= s3q(2)-5.d0/46.d0*(p2r*pgg(2)+p2i*pgg(1)) * do i=1,2 p33(i)= 2.d0/3.d0*co(i)-9.d0*b21w(i)+25.d0/4.d0*b0w(i)- # b21z(i)-b1z(i)-1.d0/4.d0*b0z(i) sr33(i)= -2.d0*wm2*b0w(i)+1.d0/2.d0*(zm2-hm2)*b1zh(i)+ # 1.d0/4.d0*(5.d0*zm2-hm2)*b0zh(i) enddo * s33f(1)= p2r*p3qf(1)-p2i*p3qf(2)-1/2*(p2r*(bfh(1)-bfl(1))- # p2i*(bfh(2)-bfl(2)))-3.d0/2.d0*tqm2*b0h(1) s33f(2)= p2r*p3qf(2)+p2i*p3qf(1)-1/2*(p2r*(bfh(2)-bfl(2))+ # p2i*(bfh(1)-bfl(1)))-3.d0/2.d0*tqm2*b0h(2) s33(1)= p2r*p33(1)-p2i*p33(2)+sr33(1)+s33f(1)+4.d0*wm2* # (delta-log(wm2)) s33(2)= p2r*p33(2)+p2i*p33(1)+sr33(2)+s33f(2) * r2(1)= s33(1)-5.d0/46.d0*(p2r*pgg(1)-p2i*pgg(2)) r2(2)= s33(2)-5.d0/46.d0*(p2r*pgg(2)+p2i*pgg(1)) * do i=1,2 bcl(i)= b21l(i)+b0l(i) bch(i)= b21c(i)+b0c(i) enddo * swwf(1)= p2r*(12.d0*bfl(1)+6.d0*(bch(1)-bcl(1)))- # p2i*(12.d0*bfl(2)+6.d0*(bch(2)-bcl(2)))- # 3.d0*tqm2*(b1c(1)+b0c(1)) swwf(2)= p2r*(12.d0*bfl(2)+6.d0*(bch(2)-bcl(2)))+ # p2i*(12.d0*bfl(1)+6.d0*(bch(1)-bcl(1)))- # 3.d0*tqm2*(b1c(2)+b0c(2)) * do i=1,2 pwwb(i)= 2.d0/3.d0*co(i)-9.d0*(b21zw(i)+b1zw(i))+7.d0/4.d0* # b0zw(i)-b21wh(i)-b1wh(i)-1.d0/4.d0*b0wh(i)+sth2*( # 8.d0*b21zw(i)-2.d0*b0zw(i)+8.d0*b1zw(i)- # 8.d0*b21pw(i)-8.d0*b1pw(i)+2.d0*b0pw(i)) srwwb(i)= 9.d0/2.d0*(zm2-wm2)*b1zw(i)+1.d0/4.d0*(13.d0*zm2- # 21.d0*wm2)*b0zw(i)+1.d0/2.d0*(wm2-hm2)*b1wh(i)+ # 1.d0/4.d0*(5.d0*wm2-hm2)*b0wh(i)+sth2*( # 2.d0*(wm2-zm2)*(2.d0*b1zw(i)+b0zw(i))-2.d0*wm2*( # 2.d0*b1pw(i)+b0pw(i))) enddo sww(1)= p2r*pwwb(1)-p2i*pwwb(2)+srwwb(1)+swwf(1)+4.d0*wm2* # (delta-log(wm2)) sww(2)= p2r*pwwb(2)+p2i*pwwb(1)+srwwb(2)+swwf(2) * r3(1)= sww(1)-5.d0/46.d0*(p2r*pgg(1)-p2i*pgg(2)) r3(2)= sww(2)-5.d0/46.d0*(p2r*pgg(2)+p2i*pgg(1)) * return end * *---------------------------------------------------------- * subroutine zztohadr5(e,der,eder) ****************************************************************** * * * subroutine for the evaluation of the light hadron * * contributions to delta_r and delta_g * * using fits to the * * qed vacuum polarization from e^+ e^- data * * * * f. jegerlehner, paul scherrer institute, ch-5232 villigen * * * * reference: f. jegerlehner, z. phys. c32 (1986) 195 * * h. burkhardt et al., z. phys. c42 (1989) 497 * * s. eidelman, f. jegerlehner, z. phys. c (1995) * * * ****************************************************************** * version: 24/02/1995 * * notation: e energy ( momentum transfer ): e>0 timelike , e<0 spacelike * st2 is sin^2(theta); st2= 0.2322 is the reference value * the routine returns the hadronic contribution of 5 flavors (u,d,s,c,b) * to der= delta_r with hadronic error errder * and deg= delta_g with hadronic error errdeg * the effective value of the fine structure constant alphaqed at energy * e is alphaqed(e)= alphaqed(0)/(1-delta_r) ,similarly for the su(2) * coupling alphasu2(e)= alphasu2(0)/(1-delta_g), where delta_r(g) is the * sum of leptonic, hadronic contributions (top to be added). * * this program does not yet know how to compute delta r and delta g for * energies in the ranges |e|>1tev and 2m_pi < e < 40(13) gev !!!!!!!!! * implicit real*8(a-h,o-z) parameter(nf=9,ns=4) * dimension res(ns) dimension ae(nf,ns),eu(nf),eo(nf),rm1(nf) dimension c1(nf,ns),c2(nf,ns),c3(nf,ns),c4(nf,ns),rl1(nf,ns) * do i=1,nf eu(i)= 0.d0 eo(i)= 0.d0 rm1(i)= 0.d0 do j=1,ns ae(i,j)= 0.d0 c1(i,j)= 0.d0 c2(i,j)= 0.d0 c3(i,j)= 0.d0 c4(i,j)= 0.d0 enddo enddo * * #1# delta_r * fit parameters spacelike -1000 to -200 gev eu(1) = -1.d3 eo(1) = -2.d2 rm1(1) = -1.d3 c1(1,1)= 4.2069394d-02 c2(1,1)= 2.9253566d-03 c3(1,1)= -6.7782454d-04 c4(1,1)= 9.3214130d-06 * fit parameters spacelike -200 to -20 gev eu(2) = -2.d2 eo(2) = -2.d1 rm1(2) = -1.d2 c1(2,1)= 2.8526291d-02 c2(2,1)= 2.9520725d-03 c3(2,1)= -2.7906310d-03 c4(2,1)= 6.4174528d-05 * fit parameters spacelike -20 to -2 gev eu(3) = -2.d1 eo(3) = -2.d0 rm1(3) = -2.d1 rl1(3,1)= 9.3055d-03 c1(3,1)= 2.8668314d-03 c2(3,1)= 0.3514608d0 c3(3,1)= 0.5496359d0 c4(3,1)= 1.9892334d-04 ae(3,1)= 3.d0 * fit parameters spacelike -2 to 0.25 gev eu(4) = -2.d0 eo(4) = 0.25d0 rm1(4) = -2.d0 rl1(4,1)= 9.3055d-03 c1(4,1)= 2.2694240d-03 c2(4,1)= 8.073429d0 c3(4,1)= 0.1636393d0 c4(4,1)= -3.3545541d-05 ae(4,1)= 2.d0 * fit parameters timelike 0.25 to 2 gev eu(5) = 0.25d0 eo(5) = 2.d0 * fit parameters timelike 2 to 40 gev eu(6) = 2.d0 eo(6) = 4.d1 * fit parameters timelike 40 to 80 gev eu(7) = 4.d1 eo(7) = 8.d1 rm1(7) = 8.d1 c1(7,1)= 2.7266588d-02 c2(7,1)= 2.9285045d-03 c3(7,1)= -4.7720564d-03 c4(7,1)= 7.7295507d-04 * fit parameters timelike 80 to 250 gev eu(8) = 8.d1 eo(8) = 2.5d2 rm1(8) = 91.18880d0 c1(8,1)= 2.8039809d-02 c2(8,1)= 2.9373798d-03 c3(8,1)= -2.8432352d-03 c4(8,1)= -5.2537734d-04 * fit parameters timelike 250 to 1000 gev eu(9) = 2.5d2 eo(9) = 1.d3 rm1(9) = 1.d3 c1(9,1)= 4.2092260d-02 c2(9,1)= 2.9233438d-03 c3(9,1)= -3.2966913d-04 c4(9,1)= 3.4324117d-07 * #2# delta_g * fit parameters spacelike -1000 to -200 gev c1(1,2)= 8.6415343d-02 c2(1,2)= 6.0127582d-03 c3(1,2)= -6.7379221d-04 c4(1,2)= 9.0877611d-06 * fit parameters spacelike -200 to -20 gev c1(2,2)= 5.8580618d-02 c2(2,2)= 6.0678599d-03 c3(2,2)= -2.4153464d-03 c4(2,2)= 6.1934326d-05 * fit parameters spacelike -20 to -2 gev rl1(3,2)= 1.9954d-02 c1(3,2)= 5.7231588d-03 c2(3,2)= 0.3588257d0 c3(3,2)= 0.5532265d0 c4(3,2)= 6.0730567d-04 ae(3,2)= 3.d0 * fit parameters spacelike -2 to 0.25 gev rl1(4,2)= 1.9954d-02 c1(4,2)= 4.8065037d-03 c2(4,2)= 8.255167d0 c3(4,2)= 0.1599882d0 c4(4,2)= -1.8624817d-04 ae(3,2)= 2.d0 * fit parameters timelike 40 to 80 gev c1(7,2)= 5.5985276d-02 c2(7,2)= 6.0203830d-03 c3(7,2)= -5.0066952d-03 c4(7,2)= 7.1363564d-04 * fit parameter timelike 80 to 250 gev c1(8,2)= 5.7575710d-02 c2(8,2)= 6.0372148d-03 c3(8,2)= -3.4556778d-03 c4(8,2)= -4.9574347d-04 * fit parameters timelike 250 to 1000 gev c1(9,2)= 8.6462371d-02 c2(9,2)= 6.0088057d-03 c3(9,2)= -3.3235471d-04 c4(9,2)= 5.9021050d-07 * #3# error delta_r * fit parameters spacelike -1000 to -200 gev c1(1,3)= 6.3289929d-04 c2(1,3)= 3.3592437d-06 c3(1,3)= 0.d0 c4(1,3)= 0.d0 * fit parameters spacelike -200 to -20 gev c1(2,3)= 6.2759849d-04 c2(2,3)= -1.0816625d-06 c3(2,3)= 5.050189d0 c4(2,3)= -9.6505374d-02 ae(2,3)= 1.d0 * fit parameters spacelike -20 to -2 gev rl1(3,3)= 2.0243d-04 c1(3,3)= 1.0147886d-04 c2(3,3)= 1.819327d0 c3(3,3)= -0.1174904d0 c4(3,3)= -1.2404939d-04 ae(3,3)= 3.d0 * fit parameters spacelike -2 to 0.25 gev rl1(4,3)= 2.0243d-04 c1(4,3)= -7.1368617d-05 c2(4,3)= 9.980347d-04 c3(4,3)= 1.669151d0 c4(4,3)= 3.5645600d-05 ae(4,3)= 2.d0 * fit parameters timelike 40 to 80 gev c1(7,3)= 6.4947648d-04 c2(7,3)= 4.9386853d-07 c3(7,3)= -55.22332d0 c4(7,3)= 26.13011d0 * fit parameters timelike 80 to 250 gev c1(8,3)= 6.4265809d-04 c2(8,3)= -2.8453374d-07 c3(8,3)= -23.38172d0 c4(8,3)= -6.251794d0 * fit parameter timelike 250 to 1000 gev c1(9,3)= 6.3369947d-04 c2(9,3)= -2.0898329d-07 c3(9,3)= 0.d0 c4(9,3)= 0.d0 * #4# error delta_g * fit parameters spacelike -1000 to -200 gev c1(1,4)= 1.2999176d-03 c2(1,4)= 7.4505529d-06 c3(1,4)= 0.d0 c4(1,4)= 0.d0 * fit parameters spacelike -200 to -20 gev c1(2,4)= 1.2883141d-03 c2(2,4)= -1.3790827d-06 c3(2,4)= 8.056159d0 c4(2,4)= -0.1536313d0 ae(2,4)= 1.d0 * fit parameters spacelike -20 to -2 gev rl1(3,4)= 4.3408d-04 c1(3,4)= 2.0489733d-04 c2(3,4)= 2.065011d0 c3(3,4)= -0.6172962d0 c4(3,4)= -2.5603661d-04 ae(3,4)= 3.d0 * fit parameters spacelike -2 to 0.25 gev rl1(4,4)= 4.3408d-04 c1(4,4)= -1.5095409d-04 c2(4,4)= 9.9847501d-04 c3(4,4)= 1.636659d0 c4(4,4)= 7.5892596d-05 ae(4,4)= 2.d0 * fit parameters timelike 40 to 80 gev c1(7,4)= 1.3335156d-03 c2(7,4)= 2.2939612d-07 c3(7,4)= -246.4966d0 c4(7,4)= 114.9956d0 * fit parameters timelike 80 to 250 gev c1(8,4)= 1.3196438d-03 c2(8,4)= 2.8937683d-09 c3(8,4)= 5449.778d0 c4(8,4)= 930.3875d0 * fit parameters timelike 250 to 1000 gev c1(9,4)= 1.3016918d-03 c2(9,4)= -3.6027674d-07 c3(9,4)= 0.d0 c4(9,4)= 0.d0 * se= 654.d0/643.d0 ! rescaling error to published version 1995 s= e*e der= 0.d0 deg= 0.d0 eder= 0.d0 edeg= 0.d0 if((e.gt.1.e3).or.(e.lt.-1.e3)) goto 100 if((e.lt.4.e1).and.(e.gt.0.25d0)) goto 100 i= 1 do while (e.ge.eo(i)) i= i+1 enddo if(e.eq.1.e3) i= 9 if(e.eq.0.d0 ) goto 100 s0= sign(1.d0,rm1(i))*rm1(i)**2 s = sign(1.d0,e)*e*e x1= s0/s xi= 1.d0/x1 x2= x1*x1 if(ae(i,1).le.0.d0) then do j= 1,4 xlar= xi+ae(i,j)*exp(-xi) xlog= log(xlar) res(j)= c1(i,j)+c2(i,j)*(xlog+c3(i,j)*(x1-1.d0)+ # c4(i,j)*(x2-1.0)) enddo else if (ae(i,1).eq.2.d0) then hx= xi*xi do j= 1,2 fx= 1.d0-c2(i,j)*s gx= c3(i,j)*s/(c3(i,j)-s) xx= log(abs(fx))+c2(i,j)*gx res(j)= c1(i,j)*xx-rl1(i,j)*gx+c4(i,j)*hx enddo do j= 3,4 u= abs(s) gx= -c3(i,j)*u/(c3(i,j)+u) xx= xi**3/(sqrt(abs(xi))**5+c2(i,j)) res(j)= c1(i,j)*xx-rl1(i,j)*gx+c4(i,j)*hx enddo else if (ae(i,1).eq.3.0) then hx= xi do j= 1,4 fx= 1.d0-c2(i,j)*s gx= c3(i,j)*s/(c3(i,j)-s) xx= log(abs(fx))+c2(i,j)*gx res(j)= c1(i,j)*xx-rl1(i,j)*gx+c4(i,j)*hx enddo endif der= res(1) eder= res(3)*se 100 return end * *----------------------------------------------------- * real*8 function zztoli2(xr,xi) implicit real*8(a,b,d-h,o-z) implicit complex*16 (c) * common/zzted/eps,delta common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai * data b0/1.d0/,b1/-0.25d0/,b2/0.277777777777778d-1/, # b4/-0.277777777777778d-3/,b6/0.472411186696901d-5/ * c0= cmplx(0.d0) pis= pi*pi pis6= pis/6.d0 cx= cmplx(xr,xi) comx= 1.d0-cx if (xr.lt.0.d0) then cy= 1.d0-cx sign1= -1.d0 clnx= log(cx) clnomx= log(comx) cadd1= pis6-clnx*clnomx else cy= cx sign1= 1.d0 cadd1= c0 endif comy= 1.d0-cy ym2= conjg(cy)*cy ym= sqrt(ym2) if (ym.gt.1.d0) then cz= conjg(cy)/ym2 sign2= -1.d0 coy= -cy clnoy= log(coy) cadd2= -pis6-0.5d0*clnoy*clnoy else cz= cy sign2= 1.d0 cadd2= c0 endif comz= 1.d0-cz zr= real(cz) if (zr.gt.0.5d0) then ct= 1.d0-cz comt= 1.d0-ct sign3= -1.d0 clnz= log(cz) clnomz= log(comz) cadd3= pis6-clnz*clnomz else ct= cz comt= 1.d0-ct sign3= 1.d0 cadd3= c0 endif cpar= log(comt) cpar2= cpar*cpar * cres= -cpar*(b0-cpar*(b1-cpar*(b2+cpar2*(b4+b6*cpar2)))) * cli2= sign1*(sign2*(sign3*cres+cadd3)+cadd2)+cadd1 * zztoli2= real(cli2) * return end * subroutine zztocflag(oflag,oval) implicit real*8(a-h,o-z) * character*(*) oflag,oval character*4 ochan character*1 out,ofl,otop,oqed,onu,oqedf,oqcdf,oscheme * common/zztn/onu common/zztout/out common/zztch/ochan common/zztoqed/oqed common/zztqedfsr/oqedf common/zztqedfsr/oqcdf common/zztofl/ofl,otop common/zztosch/oscheme * if(oflag.eq.'onu') then onu= oval else if(oflag.eq.'out') then out= oval else if(oflag.eq.'ochan') then ochan= oval else if(oflag.eq.'oqed') then oqed= oval else if(oflag.eq.'oqedf') then oqedf= oval else if(oflag.eq.'oqcdf') then oqcdf= oval else if(oflag.eq.'ofl') then ofl= oval else if(oflag.eq.'otop') then otop= oval else if(oflag.eq.'oscheme') then oscheme= oval endif * return end * *----------------------------------------------------------------------- * subroutine zztoalals(p2r,p2i,v1,a1,f1,v1p,a1p,v1i,a1i,f1i) implicit real*8(a-h,o-z) * common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2 dimension r(2),x(2),xo(2),umx(2),a(2),b(2),b2(2),ab(2), # 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), # xp(2),umr(2),trmr(2),trumr(2),rpp(2),bet(2),xmo(2), # rr(2),rmo(2),rrmo(2),rhop(2),rhom(2),ga(2),ch(2), # alp(2) * data z3/1.20205690315959428540d0/ data pi/3.141592653589793238462643d0/ * *-----this subroutine uses the wrong metric, thus the sign of p^2 * must be changed * qeps= 1.d-30 pis= pi*pi * z2= pis/6.d0 r(1)= -0.25d0*p2r/tqm2 x(1)= -p2r/tqm2 if(p2i.eq.0.d0) then r(2)= 0.25d0*qeps/tqm2 x(2)= qeps/tqm2 else r(2)= -0.25d0*p2i/tqm2 x(2)= -p2i/tqm2 endif * *-----this subroutine runs at complex momenta on the first Riemann sheet * as long as Re(s_p) < m^2_top * do i=1,2 xo(i)= -x(i) enddo umx(1)= 1.d0-x(1) umx(2)= -x(2) umr(1)= 1.d0-r(1) umr(2)= -r(2) rmo(1)= r(1)-1.d0 rmo(2)= r(2) xmo(1)= x(1)-1.d0 xmo(2)= x(2) * call zztocqlnx(xo,a) call zztocqlnomx(x,umx,b) call zztocqlnx(xmo,bet) call zztocqlnx(x,alp) b2(1)= b(1)*b(1)-b(2)*b(2) b2(2)= 2.d0*b(1)*b(2) ab(1)= a(1)*b(1)-a(2)*b(2) ab(2)= a(1)*b(2)+a(2)*b(1) * if(r(1).gt.1.d0) then if(p2i.ne.0.d0) then rmod= sqrt(r(1)*r(1)+r(2)*r(2)) rr(1)= 0.5d0*(rmod+r(1)) rr(1)= sqrt(rr(1)) rr(2)= 0.5d0*(rmod-r(1)) rr(2)= r(2)/abs(r(2))*sqrt(rr(2)) else rr(1)= sqrt(r(1)) rr(2)= qeps endif if(p2i.ne.0.d0) then rmomod= sqrt(rmo(1)*rmo(1)+rmo(2)*rmo(2)) rrmo(1)= 0.5d0*(rmomod+rmo(1)) rrmo(1)= sqrt(rrmo(1)) rrmo(2)= 0.5d0*(rmomod-rmo(1)) rrmo(2)= rmo(2)/abs(rmo(2))*sqrt(rrmo(2)) else rrmo(1)= sqrt(rmo(1)) rrmo(2)= -qeps endif do i=1,2 rhop(i)= 2.d0*rr(i) rhom(i)= 2.d0*rrmo(i) enddo call zztocqlnx(rhop,ga) call zztocqlnx(rhom,ch) endif * if(p2i.ne.0.d0) then rmod= sqrt(r(1)*r(1)+r(2)*r(2)) rmr(1)= 0.5d0*(rmod-r(1)) rmr(1)= sqrt(rmr(1)) rmr(2)= 0.5d0*(rmod+r(1)) rmr(2)= -r(2)/abs(r(2))*sqrt(rmr(2)) else if(r(1).gt.0.d0) then rmr(1)= 0.5d0*r(2)/sqrt(r(1)) rmr(2)= -sqrt(r(1)) else if(r(1).lt.0.d0) then rmr(1)= sqrt(-r(1)) rmr(2)= -0.5d0*r(2)/sqrt(-r(1)) endif endif * if(p2i.ne.0.d0) then umrmod= sqrt(umr(1)*umr(1)+umr(2)*umr(2)) rumr(1)= 0.5d0*(umrmod+umr(1)) rumr(1)= sqrt(rumr(1)) rumr(2)= 0.5d0*(umrmod-umr(1)) rumr(2)= umr(2)/abs(umr(2))*sqrt(rumr(2)) else if(umr(1).gt.0.d0) then rumr(1)= sqrt(umr(1)) rumr(2)= -0.5d0*r(2)/sqrt(umr(1)) else if(umr(1).lt.0.d0) then rumr(1)= 0.5d0*r(2)/sqrt(-umr(1)) rumr(2)= -sqrt(-umr(1)) endif 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 zztocqlnx(rp,f) call zztocqlnx(trumr,h) call zztocqlnx(trmr,g) * rm2(1)= rm(1)*rm(1)-rm(2)*rm(2) rm2(2)= 2.d0*rm(1)*rm(2) arm= rm(1) arm2= arm*arm arm4= arm2*arm2 zrm= rm(2)/rm(1) zrm2= zrm*zrm zrm4= zrm2*zrm2 brm= rm(2) brm2= brm*brm brm4= brm2*brm2 urm= rm(1)/rm(2) urm2= urm*urm urm4= urm2*urm2 if(p2i.ne.0.d0) then rm4(1)= rm(1)**4+rm(2)**4-6.d0*(rm(1)*rm(2))**2 else if(abs(zrm).lt.1.d0) then rm4(1)= arm4*(1.d0-6.d0*zrm2+zrm4) else if(abs(urm).lt.1.d0) then rm4(1)= brm4*(1.d0-6.d0*urm2+urm4) endif endif rm4(2)= 4.d0*rm(1)*rm(2)*rm2(1) r2= r(1)*r(1)+r(2)*r(2) ri(1)= r(1)/r2 ri(2)= -r(2)/r2 rr2= r(1)*r(1)-r(2)*r(2) r22= rr2*rr2+4.d0*r(1)*r(1)*r(2)*r(2) r2i(1)= rr2/r22 r2i(2)= -2.d0*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.d0*x(1)*x(1)*x(2)*x(2) x2i(1)= xx2/x22 x2i(2)= -2.d0*x(1)*x(2)/x22 * rar(1)= 1.d0-r(1)/r2 rar(2)= r(2)/r2 if(p2i.ne.0.d0) then rarmod= sqrt(rar(1)*rar(1)+rar(2)*rar(2)) rrar(1)= 0.5d0*(rarmod+rar(1)) rrar(1)= sqrt(rrar(1)) rrar(2)= 0.5d0*(rarmod-rar(1)) rrar(2)= rar(2)/abs(rar(2))*sqrt(rrar(2)) else if(rar(1).gt.0.d0) then rrar(1)= sqrt(rar(1)) rrar(2)= 0.5d0*rar(2)/sqrt(rar(1)) else if(rar(1).lt.0.d0) then rrar(1)= 0.5d0*rar(2)/sqrt(-rar(1)) rrar(2)= sqrt(-rar(1)) endif endif * xp2= (1.d0-x(1))*(1.d0-x(1))+x(2)*x(2) xp(1)= (1.d0-x(1))/xp2 xp(2)= x(2)/xp2 rpp2= (1.d0-r(1))*(1.d0-r(1))+r(2)*r(2) rpp(1)= (1.d0-r(1))/rpp2 rpp(2)= r(2)/rpp2 * ar1r= rm2(1) ar1i= rm2(2) umar1r= 1.d0-rm2(1) ar2r= rm4(1) ar2i= rm4(2) umar2r= 1.d0-rm4(1) ar3r= xp(1) ar3i= xp(2) umar3r= 1.d0-xp(1) call zztospence(ar1r,ar1i,umar1r,cli2sr,cli2si) call zztospence(ar2r,ar2i,umar2r,cli2fr,cli2fi) call zztospence(ar3r,ar3i,umar3r,cli2xr,cli2xi) call zztocli3(ar1r,ar1i,cli3sr,cli3si) call zztocli3(ar2r,ar2i,cli3fr,cli3fi) call zztocli3(ar3r,ar3i,cli3xr,cli3xi) * comb3r= 2.d0*cli3sr-cli3fr comb3i= 2.d0*cli3si-cli3fi comb2r= cli2sr-cli2fr comb2i= cli2si-cli2fi f2r= f(1)*f(1)-f(2)*f(2) f2i= 2.d0*f(1)*f(2) * r1p= r(1)+1.5d0 aux1r= -f(1)+g(1)/3.d0+2.d0/3.d0*h(1) aux1i= -f(2)+g(2)/3.d0+2.d0/3.d0*h(2) aux2r= -3.d0*f(1)+2.d0*g(1)+4.d0*h(1) aux2i= -3.d0*f(2)+2.d0*g(2)+4.d0*h(2) aux3r= comb2r+f(1)*aux2r-f(2)*aux2i aux3pr= -2.d0*(r1p*f(1)-r(2)*f(2)) aux3i= comb2i+f(1)*aux2i+f(2)*aux2r aux3pi= -2.d0*(r1p*f(2)+r(2)*f(1)) aux4r= comb2r+f(1)*aux2r-f(2)*aux2i aux4pr= -2.d0*((r(1)-3.d0+0.25d0*ri(1))*f(1)- # (r(2)+0.25d0*ri(2))*f(2)) aux4i= comb2i+f(1)*aux2i+f(2)*aux2r aux4pi= -2.d0*((r(1)-3.d0+0.25d0*ri(1))*f(2)+ # (r(2)+0.25d0*ri(2))*f(1)) aux5r= r(1)-1.d0/6.d0-7.d0/48.d0*ri(1) aux5i= r(2)-7.d0/48.d0*ri(2) aux6r= r(1)-11.d0/12.d0+5.d0/48.d0*ri(1)+1.d0/32.d0* # r2i(1) aux6i= r(2)+5.d0/48.d0*ri(2)+1.d0/32.d0*r2i(2) * ex0vr= 4.d0*(r(1)-0.25d0*ri(1)) ex0vi= 4.d0*(r(2)-0.25d0*ri(2)) * ex0ar= 4.d0*(r(1)-1.5d0+0.5d0*ri(1)) ex0ai= 4.d0*(r(2)+0.5d0*ri(2)) * ex0xr= x(1)-1.5d0+0.5d0*x2i(1) ex0xi= x(2)+0.5d0*x2i(2) * ex1r= comb3r+8.d0/3.d0*(f(1)*comb2r-f(2)*comb2i)+ # 4.d0*(f2r*aux1r-f2i*aux1i) ex1i= comb3i+8.d0/3.d0*(f(1)*comb2i+f(2)*comb2r)+ # 4.d0*(f2r*aux1i+f2i*aux1r) * ex2r= 8.d0/3.d0*((r(1)+0.5d0)*aux3r-r(2)*aux3i)+aux3pr ex2i= 8.d0/3.d0*((r(1)+0.5d0)*aux3i+r(2)*aux3r)+aux3pi * ex3r= -8.d0*(aux5r*f2r-aux5i*f2i)+13.d0/6.d0+z3*ri(1) ex3i= -8.d0*(aux5r*f2i+aux5i*f2r)+z3*ri(2) * ex4r= 8.d0/3.d0*((r(1)-1.d0)*aux4r-r(2)*aux4i)+aux4pr ex4i= 8.d0/3.d0*((r(1)-1.d0)*aux4i+r(2)*aux4r)+aux4pi * ex5r= -8.d0*(aux6r*f2r-aux6i*f2i)+ # 13.d0/6.d0-3.d0*z2+(0.25d0-2.d0*z3)*ri(1) ex5i= -8.d0*(aux6r*f2i+aux6i*f2r)+(0.25d0-2.d0*z3)*ri(2) * ex6r= cli3xr+2.d0/3.d0*(b(1)*cli2xr-b(2)*cli2xi)-1.d0/6.d0* # (b2(1)*(a(1)-b(1))-b2(2)*(a(2)-b(2))) ex6i= cli3xi+2.d0/3.d0*(b(1)*cli2xi+b(2)*cli2xr)-1.d0/6.d0* # (b2(1)*(a(2)-b(2))+b2(2)*(a(1)-b(1))) * ex7r= 1.d0/3.d0*((x(1)+0.5d0-0.5d0*xi(1))*(cli2xr-ab(1))- # (x(2)-0.5d0*xi(2))*(cli2xi-ab(2)))+1.d0/3.d0*(b2(1)* # (x(1)-1.d0/8.d0-xi(1)+5.d0/8.d0*x2i(1))-b2(2)* # (x(2)-xi(2)+5.d0/8.d0*x2i(2)))-0.25d0*(b(1)* # (x(1)-5.d0/2.d0+2.d0/3.d0*xi(1)+5.d0/6.d0*x2i(1))- # b(2)*(x(2)+2.d0/3.d0*xi(2)+5.d0/6.d0*x2i(2))) ex7i= 1.d0/3.d0*((x(1)+0.5d0-0.5d0*xi(1))*(cli2xi-ab(2))+ # (x(2)-0.5d0*xi(2))*(cli2xr-ab(1)))+1.d0/3.d0*(b2(1)* # (x(2)-xi(2)+5.d0/8.d0*x2i(2))+b2(2)*(x(1)-1.d0/8.d0- # xi(1)+5.d0/8.d0*x2i(1)))-0.25d0*(b(1)*(x(2)+2.d0/3.d0* # xi(2)+5.d0/6.d0*x2i(2))+b(2)*(x(1)-5.d0/2.d0+2.d0/3.d0* # xi(1)+5.d0/6.d0*x2i(1))) * ex8r= -3.d0/4.d0*z2+13.d0/12.d0-5.d0/24.d0*xi(1)-0.5d0*z3*x2i(1) ex8i= -5.d0/24.d0*xi(2)-0.5d0*z3*x2i(2) * v1r= ex0vr*ex1r-ex0vi*ex1i+rrar(1)*ex2r-rrar(2)*ex2i+ex3r a1r= ex0ar*ex1r-ex0ai*ex1i+rrar(1)*ex4r-rrar(2)*ex4i+ex5r f1x= ex0xr*ex6r-ex0xi*ex6i+ex7r+ex8r f1xi= ex0xr*ex6i+ex0xi*ex6r+ex7i+ex8i * ex0vpr= 4.d0*(1.d0+0.25d0*r2i(1)) ex0vpi= r2i(2) ex0apr= 4.d0*(1.d0-0.5d0*r2i(1)) ex0api= -2.d0*r2i(2) * bux1r= -comb2r-f(1)*aux2r+f(2)*aux2i bux1i= -comb2i-f(1)*aux2i-f(2)*aux2r bux2r= 2.d0/3.d0*(ri(1)*bux1r-ri(2)*bux1i) bux2i= 2.d0/3.d0*(ri(1)*bux1i+ri(2)*bux1r) bux3r= bux2r-(1.d0-5.d0/6.d0*ri(1))*f(1)- # 5.d0/6.d0*ri(2)*f(2) bux3i= bux2i-(1.d0-5.d0/6.d0*ri(1))*f(2)+ # 5.d0/6.d0*ri(2)*f(1) bux4r= f2r*(4.d0/3.d0*ri(1)-11.d0/6.d0*r2i(1)- # 4.d0*rpp(1))-f2i*(4.d0/3.d0*ri(2)-11.d0/6.d0* # r2i(2)-4.d0*rpp(2))-1.d0-1.5d0*ri(1)-z3*r2i(1) bux4i= f2r*(4.d0/3.d0*ri(2)-11.d0/6.d0*r2i(2)- # 4.d0*rpp(2))+f2i*(4.d0/3.d0*ri(1)-11.d0/6.d0* # r2i(1)-4.d0*rpp(1))-1.5d0*ri(2)-z3*r2i(2) cux2r= 4.d0/3.d0*(-ri(1)*bux1r+ri(2)*bux1i) cux2i= 4.d0/3.d0*(-ri(1)*bux1i-ri(2)*bux1r) cux3r= cux2r-(1.d0+13.d0/6.d0*ri(1)-0.5d0*r2i(1))*f(1)+ # (13.d0/6.d0*ri(2)-0.5d0*r2i(2))*f(2) cux3i= cux2i-(1.d0+13.d0/6.d0*ri(1)-0.5d0*r2i(1))*f(2)- # (13.d0/6.d0*ri(2)-0.5d0*r2i(2))*f(1) cux4r= f2r*ri(1)-f2i*ri(2) cux4i= f2r*ri(2)+f2i*ri(1) cux5r= cux4r*(-20.d0/3.d0+13.d0/6.d0*ri(1)+0.5d0*r2i(1))- # cux4i*(13.d0/6.d0*ri(2)+0.5d0*r2i(2))-1.d0+3.d0*ri(1)+ # (2.d0*z3-0.5d0)*r2i(1) cux5i= cux4r*(13.d0/6.d0*ri(2)+0.5d0*r2i(2))+ # cux4i*(20.d0/3.d0+13.d0/6.d0*ri(1)+0.5d0*r2i(1))+ # 3.d0*ri(2)+(2.d0*z3-0.5d0)*r2i(2) * v1pr= ex0vpr*ex1r-ex0vpi*ex1i+2.d0*(rrar(1)*bux3r- # rrar(2)*bux3i)+bux4r a1pr= ex0apr*ex1r-ex0api*ex1i+2.d0*(rrar(1)*cux3r- # rrar(2)*cux3i)+cux5r * tlog= log(2.d0) if((r(1)+1.d0).lt.(-1.d4)) then v1= -2.d0*r(1)*g(1)-6.d0*g(1)+1.d0/r(1)*(-3.d0*g(1)*g(1)- # 5.d0/4.d0*g(1)+z3+5.d0/144.d0) v1i= 0.d0 a1= -2.d0*r(1)*g(1)+3.d0*(-2.d0*g(1)*g(1)+g(1)-z2)+1.d0/ # r(1)*(3.d0*g(1)*g(1)+11.d0/4.d0*g(1)-2.d0*z3+ # 293.d0/144.d0) a1i= 0.d0 else if(abs(r(1)).lt.1.d-3) then v1= (4.d0*z3-5.d0/6.d0)*r(1)+328.d0/81.d0*r(1)*r(1) v1i= 0.d0 a1= 3.d0*(-2.d0*z3-z2+7.d0/4.d0)+(4.d0*z3-49.d0/18.d0)* # r(1)+689.d0/405.d0*r(1)*r(1) a1i= 0.d0 else if(r(1).lt.1.d0.and.((1.d0-r(1)).lt.1.d-3)) then v1= -12.d0*z2*h(1)-13.d0/2.d0*z3+3.d0*z2*(-2.d0*tlog+ # 11.d0/4.d0)+13.d0/6.d0+8.d0*pi*sqrt(1.d0-r(1))+ # (1.d0-r(1))*(20.d0*z2*h(1)+27.d0/2.d0*z3+z2*(10.d0* # tlog-43.d0/4.d0)-3.d0/2.d0) v1i= 0.d0 a1= -2.d0*z3-3.d0/8.d0*z2+29.d0/12.d0+(1.d0-r(1))*( # 8.d0*z2*h(1)+3.d0*z3+2.d0*z2*(2.d0*tlog-5.d0)- # 3.d0/2.d0) a1i= 0.d0 else if(r(1).gt.1.d0.and.((r(1)-1.d0).lt.1.d-3)) then v1= -12.d0*z2*ch(1)-13.d0/2.d0*z3+3.d0*z2*(-2.d0*tlog+ # 11.d0/4.d0)+13.d0/6.d0+(r(1)-1.d0)*(-20.d0*z2*ch(1)- # 27.d0/2.d0*z3+z2*(-10.d0*tlog+43.d0/4.d0)+ # 3.d0/2.d0) v1i= pi*(6.d0*z2-8.d0*sqrt(r(1)-1.d0)+10.d0*z2*(r(1)-1.d0)) a1= -2.d0*z3-3.d0/8.d0*z2+29.d0/12.d0+(r(1)-1.d0)*(-8.d0* # z2*ch(1)-3.d0*z3+2.d0*z2*(-2.d0*tlog+5.d0)+ # 3.d0/2.d0) a1i= 4.d0*pi*z2*(r(1)-1.d0) else if(r(1).gt.1.d4) then v1= -2.d0*r(1)*ga(1)-6.d0*ga(1)+1.d0/r(1)*(-3.d0*ga(1)* # ga(1)-5.d0/4.d0*ga(1)+z3+9.d0/2.d0*z2+5.d0/144.d0) v1i= pi*(r(1)+3.d0+1.d0/r(1)*(3.d0*ga(1)+5.d0/8.d0)) a1= -2.d0*r(1)*ga(1)+3.d0*(-2.d0*ga(1)*ga(1)+ga(1)+2.d0*z2)+ # 1.d0/r(1)*(3.d0*ga(1)*ga(1)+11.d0/4.d0*ga(1)-2.d0*z3- # 9.d0/2.d0*z2+293.d0/144.d0) a1i= pi*(r(1)+6.d0*ga(1)-3.d0/2.d0+1.d0/r(1)*(-3.d0*ga(1)- # 11.d0/8.d0)) else v1= v1r a1= a1r v1i= ex0vr*ex1i+ex0vi*ex1r+rrar(1)*ex2i+rrar(2)*ex2r+ex3i a1i= ex0ar*ex1i+ex0ai*ex1r+rrar(1)*ex4i+rrar(2)*ex4r+ex5i endif * if((x(1)+1.d0).lt.(-1.d4)) then f1= -x(1)*a(1)/4.d0-3.d0/4.d0*(a(1)*a(1)/2.d0+a(1)/2.d0+ # z2)+1.d0/2.d0/x(1)*(3.d0/2.d0*a(1)+1.d0) f1i= 0.d0 else if(abs(x(1)).lt.1.d-3) then f1= -3.d0/2.d0*z3-z2/2.d0+23.d0/16.d0+x(1)*(z3-z2/9.d0- # 25.d0/72.d0)+x(1)*x(1)/8.d0*(z2+25.d0/24.d0) f1i= 0.d0 else if(x(1).lt.1.d0.and.((1.d0-x(1)).lt.1.d-3)) then f1= -z3/2.d0-z2/12.d0+7.d0/8.d0-(1.d0-x(1))*(z3+z2+ # 13.d0/24.d0)+(1.d0-x(1))**2*(3.d0/8.d0*b(1)*b(1)- # (z2+9.d0/8.d0)*b(1)-3.d0/2.d0*z3-z2/3.d0+5.d0/24.d0) f1i= 0.d0 else if(x(1).gt.1.d0.and.((x(1)-1.d0).lt.1.d-3)) then f1= -z3/2.d0-z2/12.d0+7.d0/8.d0+(x(1)-1.d0)*(z3+z2+ # 13.d0/24.d0)+(x(1)-1.d0)**2*(3.d0/8.d0*bet(1)*bet(1)- # (z2+9.d0/8.d0)*bet(1)-3.d0/2.d0*z3-31.d0/12.d0*z2+ # 5.d0/24.d0) f1i= pi*(x(1)-1)**2*(-3.d0/4.d0*bet(1)+z2+9.d0/8.d0) else if(x(1).gt.1.d4) then f1= -x(1)*alp(1)/4.d0+3.d0/2.d0*(-alp(1)*alp(1)/4.d0- # alp(1)/4.d0+z2)+1.d0/2.d0/x(1)*(3.d0/2.d0*alp(1)+1.d0) f1i= pi*(x(1)/4.d0+3.d0/4.d0*(alp(1)+1.d0/2.d0)- # 3.d0/4.d0/x(1)) else f1= f1x f1i= f1xi endif * v1p= v1pr a1p= a1pr * return end * *-----cli3----------------------------------------------------------- * computes li_3(x) for complex x * subroutine zztocli3(xr,xi,cli3r,cli3i) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta * 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.d0,-0.75d0,0.236111111111111111111111111111111d0, # -3.472222222222222222222222222222222d-2, # 6.481481481481481481481481481481482d-4, # 4.861111111111111111111111111111111d-4, # -2.393550012597631645250692869740488d-5, # -1.062925170068027210884353741496599d-5, # 7.794784580498866213151927437641723d-7, # 2.526087595532039976484420928865373d-7, # -2.359163915200471237027273583310139d-8, # -6.168132746415574698402981231264060d-9, # 6.824456748981078267312315451125495d-10, # 1.524285616929084572552216019859487d-10, # -1.916909414174054295837274763110831d-11/ data z3/1.20205690315959428540d0/ data qpi/3.141592653589793238462643d0/ * qpis= qpi*qpi z2= qpis/6.d0 do n=0,14 bf(n)= b(n)/(n+1.d0) 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.d0 if(xtst.le.1.d-20) then y(1)= x(1) y(2)= x(2) addx(1)= 0.d0 addx(2)= 0.d0 else if(xm.gt.1.d-20) then y(1)= x(1)/xm2 y(2)= -x(2)/xm2 ox(1)= -x(1) ox(2)= -x(2) call zztocqlnx(ox,clnx) rlnx2= clnx(1)*clnx(1) ailnx2= clnx(2)*clnx(2) addx(1)= -clnx(1)*(z2+1.d0/6.d0*(rlnx2-3.d0*ailnx2)) addx(2)= -clnx(2)*(z2+1.d0/6.d0*(3.d0*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.d0.or.y2r.lt.0.d0) then ytst= y(1)-0.5d0 if(ytst.le.1.d-20) then * *-----li_3(y) is computed * omy(1)= 1.d0-y(1) omy(2)= -y(2) call zztocqlnomx(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.d-20) then ym2= y(1)*y(1)+y(2)*y(2) u1(1)= 1.d0-y(1)/ym2 u1(2)= y(2)/ym2 u2(1)= 1.d0-y(1) u2(2)= -y(2) call zztocqlnx(y,clny) omy(1)= 1.d0-y(1) omy(2)= -y(2) call zztocqlnomx(y,omy,clnomy) addy(1)= z3+z2*clny(1)+1.d0/6.d0*clny(1)* # (clny(1)*clny(1)-3.d0*clny(2)*clny(2))- # 0.5d0*clnomy(1)*(clny(1)*clny(1)-clny(2)* # clny(2))+clny(1)*clny(2)*clnomy(2) addy(2)= z2*clny(2)+1.d0/6.d0*clny(2)*(3.d0* # clny(1)*clny(1)-clny(2)*clny(2))-0.5d0* # 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.d0-u1(1) omu1(2)= -u1(2) call zztocqlnomx(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.d0-u2(1) omu2(2)= -u2(2) call zztocqlnomx(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.d0) then * *-----first t * t(1)= -y(1) t(2)= -y(2) if(t(1).le.0.5d0) then * *-----li_3(t) is computed * omt(1)= 1.d0-t(1) omt(2)= -t(2) call zztocqlnomx(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.5d0) then tm2= t(1)*t(1)+t(2)*t(2) u1(1)= 1.d0-t(1)/tm2 u1(2)= t(2)/tm2 u2(1)= 1.d0-t(1) u2(2)= -t(2) call zztocqlnx(t,clnt) omt(1)= 1.d0-t(1) omt(2)= -t(2) call zztocqlnomx(t,omt,clnomt) addt(1)= z3+z2*clnt(1)+1.d0/6.d0*clnt(1)* # (clnt(1)*clnt(1)-3.d0*clnt(2)*clnt(2))- # 0.5d0*clnomt(1)*(clnt(1)*clnt(1)-clnt(2)* # clnt(2))+clnt(1)*clnt(2)*clnomt(2) addt(2)= z2*clnt(2)+1.d0/6.d0*clnt(2)*(3.d0* # clnt(1)*clnt(1)-clnt(2)*clnt(2))-0.5d0* # 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.d0-u1(1) omu1(2)= -u1(2) call zztocqlnomx(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.d0-u2(1) omu2(2)= -u2(2) call zztocqlnomx(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.d0*y(1)*y(2) if(t(1).le.0.5d0) then * *-----li_3(t^2) is computed * omt(1)= 1.d0-t(1) omt(2)= -t(2) call zztocqlnomx(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.5d0) then tm2= t(1)*t(1)+t(2)*t(2) u1(1)= 1.d0-t(1)/tm2 u1(2)= t(2)/tm2 u2(1)= 1.d0-t(1) u2(2)= -t(2) call zztocqlnx(t,clnt) omt(1)= 1.d0-t(1) omt(2)= -t(2) call zztocqlnomx(t,omt,clnomt) addt2(1)= z3+z2*clnt(1)+1.d0/6.d0*clnt(1)* # (clnt(1)*clnt(1)-3.d0*clnt(2)*clnt(2))- # 0.5d0*clnomt(1)*(clnt(1)*clnt(1)-clnt(2)* # clnt(2))+clnt(1)*clnt(2)*clnomt(2) addt2(2)= z2*clnt(2)+1.d0/6.d0*clnt(2)*(3.d0* # clnt(1)*clnt(1)-clnt(2)*clnt(2))-0.5d0* # clnomt(2)*(clnt(1)*clnt(1)-clnt(2)*clnt(2))- # clnt(1)*clnomt(1)*clnt(2) * *-----li_3(1-1/t^2) is computed * omu1(1)= 1.d0-u1(1) omu1(2)= -u1(2) call zztocqlnomx(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 res3(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)))))))))))))))) res3(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-t^2) is computed * omu2(1)= 1.d0-u2(1) omu2(2)= -u2(2) call zztocqlnomx(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 res4(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)))))))))))))))) res4(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)))))))))))))))) resb(1)= -res3(1)-res4(1)+addt2(1) resb(2)= -res3(2)-res4(2)+addt2(2) endif cli3r= -resa(1)+0.25d0*resb(1)+addx(1) cli3i= -resa(2)+0.25d0*resb(2)+addx(2) return endif end * *--------------------------------------------------------------- * subroutine zztospence(xr,xi,omxr,cli2r,cli2i) implicit real*8(a-h,o-z) * dimension b(0:14),bf(0:14) dimension x(2),omx(2),y(2),oy(2),omy(2),z(2),omz(2),t(2),omt(2) dimension clnx(2),clnomx(2),clnoy(2),clnz(2),clnomz(2) dimension add1(2),add2(2),add3(2),par(2),res(2),ct(15),sn(15) * data qpi/3.141592653589793238462643d0/ * qpis= qpi*qpi x(1)= xr x(2)= xi omx(1)= omxr omx(2)= -xi if (xr.lt.0.d0) then y(1)= omxr y(2)= -xi sign1= -1.d0 call zztocqlnx(x,clnx) call zztocqlnomx(x,omx,clnomx) add1(1)= qpis/6.d0-clnx(1)*clnomx(1)+clnx(2)*clnomx(2) add1(2)= -clnx(1)*clnomx(2)-clnx(2)*clnomx(1) else y(1)= x(1) y(2)= x(2) sign1= 1.d0 add1(1)= 0.d0 add1(2)= 0.d0 endif omy(1)= 1.d0-y(1) omy(2)= -y(2) ym2= y(1)*y(1)+y(2)*y(2) ym= sqrt(ym2) if (ym.gt.1.d0) then z(1)= y(1)/ym2 z(2)= -y(2)/ym2 sign2= -1.d0 oy(1)= -y(1) oy(2)= -y(2) call zztocqlnx(oy,clnoy) add2(1)= -qpis/6.d0-0.5d0*((clnoy(1))**2-(clnoy(2))**2) add2(2)= -clnoy(1)*clnoy(2) else z(1)= y(1) z(2)= y(2) sign2= 1.d0 add2(1)= 0.d0 add2(2)= 0.d0 endif omz(1)= 1.d0-z(1) omz(2)= -z(2) zr= z(1) if (zr.gt.0.5d0) then t(1)= 1.d0-z(1) t(2)= -z(2) omt(1)= 1.d0-t(1) omt(2)= -t(2) sign3= -1.d0 call zztocqlnx(z,clnz) call zztocqlnomx(z,omz,clnomz) add3(1)= qpis/6.d0-clnz(1)*clnomz(1)+clnz(2)*clnomz(2) add3(2)= -clnz(1)*clnomz(2)-clnz(2)*clnomz(1) else t(1)= z(1) t(2)= z(2) omt(1)= 1.d0-t(1) omt(2)= -t(2) sign3= 1.d0 add3(1)= 0.d0 add3(2)= 0.d0 endif call zztocqlnomx(t,omt,par) b(0)= 1.d0 b(1)= -1.d0/2.d0 b(2)= 1.d0/6.d0 b(4)= -1.d0/30.d0 b(6)= 1.d0/42.d0 b(8)= -1.d0/30.d0 b(10)= 5.d0/66.d0 b(12)= -691.d0/2730.d0 b(14)= 7.d0/6.d0 fact= 1.d0 do n=0,14 bf(n)= b(n)/fact fact= fact*(n+2.d0) enddo parr= par(1) pari= par(2) parm2= parr*parr+pari*pari parm= sqrt(parm2) ct(1)= parr/parm sn(1)= pari/parm 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)= -((((((((bf(14)*ct(15)*parm2+bf(12)*ct(13))*parm2+ # bf(10)*ct(11))*parm2+bf(8)*ct(9))*parm2+ # bf(6)*ct(7))*parm2+bf(4)*ct(5))*parm2+ # bf(2)*ct(3))*(-parm)+bf(1)*ct(2))*(-parm)+ # bf(0)*ct(1))*parm res(2)= -((((((((bf(14)*sn(15)*parm2+bf(12)*sn(13))*parm2+ # bf(10)*sn(11))*parm2+bf(8)*sn(9))*parm2+ # bf(6)*sn(7))*parm2+bf(4)*sn(5))*parm2+ # bf(2)*sn(3))*(-parm)+bf(1)*sn(2))*(-parm)+ # bf(0)*sn(1))*parm cli2r= sign1*(sign2*(sign3*res(1)+add3(1))+add2(1))+add1(1) cli2i= sign1*(sign2*(sign3*res(2)+add3(2))+add2(2))+add1(2) * return end * *-------------------------------------------------- * subroutine zztocqlnx(arg,res) implicit real*8(a-h,o-z) * dimension arg(2),aarg(2),res(2) * data qpi/3.141592653589793238462643d0/ * qpis= qpi*qpi do i= 1,2 aarg(i)= abs(arg(i)) enddo zm2= (arg(1))**2+(arg(2))**2 zm= sqrt(zm2) res(1)= log(zm) if (arg(1).eq.0.d0) then if (arg(2).gt.0.d0) then teta= qpi/2.d0 else teta= -qpi/2.d0 endif res(2)= teta return else if (arg(2).eq.0.d0) then if (arg(1).gt.0.d0) then teta= 0.d0 else teta= qpi endif res(2)= teta return else tnteta= aarg(2)/aarg(1) teta= atan(tnteta) sr= arg(1)/aarg(1) si= arg(2)/aarg(2) if (sr.gt.0.d0) then res(2)= si*teta else res(2)= si*(qpi-teta) endif return endif end * *-------------------------------------- * subroutine zztocqlnomx(arg,omarg,res) implicit real*8(a-h,o-z) * dimension arg(2),omarg(2),res(2),ares(2),ct(10),sn(10) * data qpi/3.141592653589793238462643d0/ * qpis= qpi*qpi 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,10 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 ares(1)= ct(10)/10.d0 ares(2)= sn(10)/10.d0 do k=9,1,-1 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 zztocqlnx(omarg,ares) endif do i= 1,2 res(i)= ares(i) enddo return end