subroutine ww(inpts,inrand,scalei,wmi,rsi,s0i,tc,alszi,sig,esig) implicit real*8(a-h,o-z) character*1 oqcd,oral * parameter (nin=98,nout=99,ndimmx=4,liw=ndimmx+2, # lw=ndimmx*(ndimmx-1)/2+12*ndimmx) * common/wpro/ipro common/wtcut/tcs,s0 common/wtifz/ifz(5),ifzzf,ik common/wtapar/s,rwm,rwg,g2,sth2,isf common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis, # pih,tfact,eps,scale * dimension vk(ndimmx) * external d01gcf,wtoregion,wwcc20 * do i=1,5 ifz(i)= 0 enddo ifzzf= 0 ik= 0 iqcd= 1 * oqcd= 'y' oqcd= 'n' oral= 'r' iosf= 1 ios= 2 ipro= 2 * if(ipro.eq.1) then scale= 0.10565839D0 else if(ipro.eq.2) then scale= scalei endif eps= 1.d-37 cfct= 0.38937966d9 alphai= 137.0359895d0 em= 0.51099906d-3 pi= 3.141592653589793238462643d0 pis= pi*pi pih= 0.5d0*pi ge= 0.5772156649d0 gf= 8.2476172696d-6 zm= 91.1866d0 rs= rsi s= rs*rs scale= scale*scale/s tcs= tc*tc s0= s0i * wm= wmi alsz= alszi * nf= 5 if(iqcd.eq.2) then x1= 0.001d0 x2= 10.0d0 xacc= 1.d-12 nf= 5 qcdl= wtoqcdlam(nf,als,wm,x1,x2,xacc) else qcdl= 0.d0 endif als= wtoralphas(zm,wm,alsz,nf) * qgfopi= 0.25d0*gf/pi wm2= wm*wm wmc= wm2*wm zm2= zm*zm wg0= 1.5d0*gf*wmc/pi if(oqcd.eq.'n') then wg= wg0 else wg= wg0*(1.d0+2.d0/3.d0*als/pi) endif wg2= wg*wg * rwm= wm/rs rwm2= wm2/s rwg= wg/rs * alwi= 128.07d0 * *-----e.m. quantities are initialized * em2= em*em alpha= 1.d0/alphai alw= 1.d0/alwi aexp= alpha/pi rln= log(s/em2) beta= 2.d0*aexp*(rln-1.d0) isf= iosf if(iosf.eq.1) then eta= beta eob= 1.d0 else if(iosf.gt.1) then eta= 2.d0*aexp*rln eob= rln/(rln-1.d0) endif feta= eta/16.d0 hbet= beta/2.d0 heta= eta/2.d0 hbeti= 2.d0/beta omhb= 1.d0-hbet ophb= 1.d0+hbet if(iosf.lt.3) then arge= hbet*(3.d0/4.d0-ge) else if(iosf.eq.3) then arge= -hbet*ge+3.d0/4.d0*heta endif ifg= 1 gamb= s14aaf(ophb,ifg) if(ifg.ne.0) then stop else d0gl= exp(arge)/gamb endif * if(ios.eq.1) then sth2= 0.5d0*pi*alw/gf/wm2 g2= 4.d0*pi*alw/sth2 else if(ios.eq.2) then sth2= 1.d0-wm2/zm2 g2= 8.d0*gf*wm2 endif g4= g2*g2 g8= g4*g4 sth4= sth2*sth2 if(ipro.eq.1) then fcol= 1.d0 else if(ipro.eq.2) then fcol= 3.d0 endif tfact= g8*sth4/2.d0/(2.d0*pi)**8*pi/4.d0/8.d0* # pi**2/4.d0*fcol*cfct * if(inpts.le.6) then jnpts= inpts else if(inpts.eq.7) then jnpts= 2011*101 vk(1)= 0.1d+1 vk(2)= 0.12294d+5 vk(3)= 0.27852d+5 vk(4)= 0.170453d+6 else if(inpts.eq.8) then jnpts= 10193*101 vk(1)= 0.1d+1 vk(2)= 0.9601d+6 vk(3)= 0.449688d+6 vk(4)= 0.792432d+6 else if(inpts.eq.9) then jnpts= 22807*151 vk(1)= 0.1d+01 vk(2)= 0.2670673d+7 vk(3)= 0.124894d+7 vk(4)= 0.521697d+6 else if(inpts.eq.10) then jnpts= 32771*181 vk(1)= 0.1d+01 vk(2)= 0.416383d+7 vk(3)= 0.1562225d+7 vk(4)= 0.39176d+7 endif * ndim= 4 npts= jnpts itrans= 0 nrand= inrand ifail= 0 call d01gcf(ndim,wwcc20,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) * ifzt= 0 do i=1,5 ifzt= ifzt+ifz(i) enddo print 1,ik print 2,1.d2*(1.d0-1.d0*ifzt/ik) print 3,ifzzf do i=1,5 print 5,ifz(i) enddo print 4 1 format(/' # of calls = ',i10) 2 format(/' 1-nc/int (%) = ',f10.3) 3 format(/' # of neg M^2 = ',i10) 4 format(/' End of Call DUMP -----------------------------',/) 5 format(1x,i10) return end * *------------------------------------------------------------------ * real*8 function wwcc20(ndim,x) implicit real*8(a-h,o-z) * common/wpro/ipro common/wtcut/tcs,s0 common/wtifz/ifz(5),ifzzf,ik common/wtapar/s,rwm,rwg,g2,sth2,isf common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis, # pih,tfact,eps,scale * dimension x(ndim),zrt1(2),zrt2(2),arg(2) * external s07aaf * *-----order of integration is: uv,vv; ym; yp * ik= ik+1 em= 0.51099906d-3 em2= em*em rem= em/sqrt(s) rwms= rwm*rwm rwmc= rwms*rwm rwmq= rwmc*rwm rwmv= rwmq*rwm rwmvi= rwmv*rwm rwmvii= rwmvi*rwm rwgs= rwg*rwg rwgm= rwg*rwm rwgms= rwgs*rwms s0s= s0*s0 * if(ndim.eq.2) then ymx= x(1) ypx= x(2) else if(ndim.eq.4) then uvx= x(1) vvx= x(2) ymx= x(3) ypx= x(4) endif * iz= 1 * if(ndim.eq.2) then ueps= 0.d0 uv= 1.d0 uvs= uv*uv ujc= 1.d0 veps= 0.d0 vv= 1.d0 vjc= 1.d0 else if(uvx.gt.1.d0) then iz= 0 ifz(1)= ifz(1)+1 go to 1 endif vvl= (sqrt(s0)+0.5d0*rem)*(sqrt(s0)+0.5d0*rem) omul= 1.d0-vvl ueps= omul*(1.d0-uvx)**hbeti uv= 1.d0-ueps ujc= omul**hbet uvs= uv*uv * uvl= uv-vvl veps= uvl*(1.d0-vvx)**hbeti vv= uv-veps vjc0= 1.d0-vvl/uv if(vjc0.le.0.d0) then iz= 0 ifz(2)= ifz(2)+1 go to 1 else vjc= vjc0**hbet endif endif * if(vv.lt.0.d0) then iz= 0 ifz(2)= ifz(2)+1 go to 1 endif svv= sqrt(vv) xm= uv xp= vv/uv xmop= xm/xp sh= vv*s ssh= sqrt(sh) * opxp= 1.d0+xp opxm= 1.d0+xm omxp= veps/uv omxm= ueps if(isf.eq.0) then stfp= 1.d0 stfm= 1.d0 else if(isf.gt.0) then if(omxp.eq.0) then stfp= d0gl else rcpx= 0.25d0*opxp*opxp rcpy= xp iflp= 1 rclp= s21baf(rcpx,rcpy,iflp) stfp= d0gl+eob*omxp**omhb*(-0.5d0*opxp+ # feta*(-4.d0*opxp*log(omxp)+ # 3.d0*opxp*log(xp)+4.d0*rclp-5.d0-xp)) endif if(omxm.eq.0) then stfm= d0gl else rcmx= 0.25d0*opxm*opxm rcmy= xm iflm= 1 rclm= s21baf(rcmx,rcmy,iflm) stfm= d0gl+eob*omxm**omhb*(-0.5d0*opxm+ # feta*(-4.d0*opxm*log(omxm)+ # 3.d0*opxm*log(xm)+4.d0*rclm-5.d0-xm)) endif endif * stf= stfp*stfm * *-----ym is initialized * yml= s0/vv ymu= 1.d0-em/ssh if(yml.gt.ymu) then iz= 0 ifz(3)= ifz(3)+1 go to 1 endif ymjc= ymu-yml ym= ymjc*ymx+yml ymvv= ym*vv ymvvi= 1.d0/ymvv omym= (1.d0-s0/vv)*(1.d0-ymx)+em/ssh*ymx omyms= omym*omym omymc= omyms*omym ymi= 1.d0/ym * *-----yp is initialized * ypl= s0/ymvv ypu= 1.d0 if(ypl.gt.ypu) then iz= 0 ifz(4)= ifz(4)+1 go to 1 endif ypjc= ypu-ypl yp= ypjc*ypx+ypl ypi= 1.d0/yp ypis= ypi*ypi ypic= ypis*ypi ypiq= ypic*ypi ypmwi= 1.d0/(yp**2*ymvv**2+4.d0*rwms*rwgs) * qyjc= 1.d0-2.d0*em2/s*(xp+xm)/vv * qh0s= em2*ym/(1.d0-ym)*(xm*ym+1.d0-xm) qhcs= qh0s+0.25d0*xm*xm*(1.d0-ym)*tcs*s by= 1.d0-xm*(1.d0-ym) bk= 0.25d0*xm*xm*tcs psf= -omymc/ym/by*bk*s/(omyms*bk*s+ym*by*em2)* # (1.d0+xm+2.d0*(1.d0-xm)/ym)+ # (1.d0-2.d0/ym+2.d0/(ym*ym))*log(1.d0+omyms*bk*s/ # (em2*ym*by)) * tjac0= ujc*vjc*ymjc*ypjc/qyjc tjac= tjac0*stf*psf/sh * ypmvv= yp*ymvv aux1= rwms+rwgs aux2= s0-ypmvv-rwms aux3= ypmvv-rwms aux4= ypmvv+rwms aux5= s0-rwms aux6= rwms+ypmvv-s0 auxg2= aux2*aux2+rwgms auxg3= aux3*aux3+rwgms auxg4= aux4*aux4+rwgms auxg5= aux5*aux5+rwgms * arg_lns= auxg3/auxg5 rlns= log(arg_lns) arg_lnt= rwms*aux1/auxg2 rlnt= log(arg_lnt) * arg_atns1= aux3/rwgm arg_atns2= aux5/rwgm * if(arg_atns1.gt.1.d0) then arg_atns1= 1.d0/arg_atns1 atns= pih-atan(arg_atns1) else if(arg_atns1.lt.-1.d0) then arg_atns1= -1.d0/arg_atns1 atns= -pih+atan(arg_atns1) else if(abs(arg_atns1).lt.1.d0) then atns= atan(arg_atns1) endif * if(arg_atns2.gt.1.d0) then arg_atns2= 1.d0/arg_atns2 atns= atns-pih+atan(arg_atns2) else if(arg_atns2.lt.-1.d0) then arg_atns2= -1.d0/arg_atns2 atns= atns+pih-atan(arg_atns2) else if(abs(arg_atns2).lt.1.d0) then atns= atns-atan(arg_atns2) endif * arg_atnt1= -rwm/rwg arg_atnt2= aux2/rwgm * if(arg_atnt1.gt.1.d0) then arg_atnt1= 1.d0/arg_atnt1 atnt= pih-atan(arg_atnt1) else if(arg_atnt1.lt.-1.d0) then arg_atnt1= -1.d0/arg_atnt1 atnt= -pih+atan(arg_atnt1) else if(abs(arg_atnt1).lt.1.d0) then atnt= atan(arg_atnt1) endif * if(arg_atnt2.gt.1.d0) then arg_atnt2= 1.d0/arg_atnt2 atnt= atnt-pih+atan(arg_atnt2) else if(arg_atnt2.lt.-1.d0) then arg_atnt2= -1.d0/arg_atnt2 atnt= atnt+pih-atan(arg_atnt2) else if(abs(arg_atnt2).lt.1.d0) then atnt= atnt-atan(arg_atnt2) endif * aln= log(s0) bln= log(ypmvv) * zsb= ypmvv-rwms arga= rwgm argbs= rwms fh0s= fh0(zsb,arga,argbs) fh1s= fh1(zsb,arga,argbs) * zsa= s0-rwms fh0s= fh0s-fh0(zsa,arga,argbs) fh1s= fh1s-fh1(zsa,arga,argbs) * ztb= -rwms argbt= ypmvv+rwms fh0t= fh0(ztb,arga,argbt) fh1t= fh1(ztb,arga,argbt) * zta= s0-ypmvv-rwms fh0t= fh0t-fh0(zta,arga,argbt) fh1t= fh1t-fh1(zta,arga,argbt) * rlnc0= log(scale) * if(ipro.eq.1) then * ww_m1= # +ymvvi*rwm**5*atns*ymi*(-2./3.*ypi+2./3.*ypi**2) ww_m1= ww_m1+ymvvi*rwm**5*atnt*ymi*rlnc0*(2*ypi**2-4* # ypi**3+2*ypi**4) ww_m1= ww_m1+ymvvi*rwm**5*atnt*ymi*(11./3.*ypi**2-10* # ypi**3+6*ypi**4) ww_m1= ww_m1+ymvvi*rwm**5*ymi*fh0t*(-2*ypi**2+4*ypi**3 # -2*ypi**4) ww_m1= ww_m1+vv*rwm*atns*(1-1./3.*yp-4./3.*ypi+2./3.* # ypi**2) ww_m1= ww_m1+vv*rwm*atnt*rlnc0*(1-2*ypi+ypi**2) ww_m1= ww_m1+vv*rwm*atnt*(23./6.-23./3.*ypi+11./3.* # ypi**2) ww_m1= ww_m1+vv*rwm*fh0t*(-1+2*ypi-ypi**2) ww_m1= ww_m1+rwm**3*atns*ymi*(2./3.-2./3.*ypi) ww_m1= ww_m1+rwm**3*atnt*ymi*rlnc0*(2*ypi-4*ypi**2+2* # ypi**3) ww_m1= ww_m1+rwm**3*atnt*ymi*(17./3.*ypi-14*ypi**2+8* # ypi**3) ww_m1= ww_m1+rwm**3*ymi*fh0t*(-2*ypi+4*ypi**2-2* # ypi**3) * ww_0= # +ymvv*vv*rwm**2*rlns*ypmwi*(-1./2.+1./3.*yp+1./6.*ypi # ) ww_0= ww_0+ymvv*vv*rwm**2*rlnt*ypmwi*(1./2.-1./3.*yp-1./ # 6.*ypi) ww_0= ww_0+ymvvi*rwm**2*ymi*rlnc0*(-4*s0*ypi**2+8*s0* # ypi**3-4*s0*ypi**4) ww_0= ww_0+ymvvi*rwm**2*ymi*aln*(4*s0*ypi**2-8*s0*ypi**3 # +4*s0*ypi**4) ww_0= ww_0+ymvvi*rwm**2*ymi*(4./3.*s0*ypi-14*s0*ypi**2+ # 28*s0*ypi**3-44./3.*s0*ypi**4) ww_0= ww_0+ymvvi*rwm**4*rlns*ymi*(-ypi+3*ypi**2-3* # ypi**3+ypi**4) ww_0= ww_0+ymvvi*rwm**4*rlnt*ymi*rlnc0*(3*ypi**2-6*ypi**3 # +3*ypi**4) ww_0= ww_0+ymvvi*rwm**4*rlnt*ymi*(15./2.*ypi**2-18*ypi**3 # +10*ypi**4) ww_0= ww_0+ymvvi*rwm**4*ymi*fh1t*(-6*ypi**2+12*ypi**3 # -6*ypi**4) ww_0= ww_0+ymvvi*ymi*rlnc0*(-s0**2*ypi**2+2*s0**2* # ypi**3-s0**2*ypi**4) ww_0= ww_0+ymvvi*ymi*aln*(s0**2*ypi**2-2*s0**2*ypi**3+ # s0**2*ypi**4) ww_0= ww_0+ymvvi*ymi*(1./3.*s0**2*ypi-7./2.*s0**2*ypi**2 # +7*s0**2*ypi**3-11./3.*s0**2*ypi**4) ww_0= ww_0+vv*rlns*(1./2.-1./6.*yp-2./3.*ypi+1./3.* # ypi**2) ww_0= ww_0+vv*rlnt*rlnc0*(1./2.-ypi+1./2.*ypi**2) ww_0= ww_0+vv*rlnt*(23./12.-23./6.*ypi+11./6.*ypi**2) ww_0= ww_0+vv*rlnc0*(1-2*ypi+ypi**2) ww_0= ww_0+vv*fh1t*(-1+2*ypi-ypi**2) ww_0= ww_0+vv*bln*(-1+2*ypi-ypi**2) ww_0= ww_0+vv*(7./2.+1./3.*yp-9*ypi+5*ypi**2) ww_0= ww_0+rwm**2*rlns*ymi*(2./3.-5./3.*ypi+3./2.* # ypi**2-1./2.*ypi**3) ww_0= ww_0+rwm**2*rlnt*ymi*rlnc0*(2*ypi-4*ypi**2+2* # ypi**3) ww_0= ww_0+rwm**2*rlnt*ymi*(20./3.*ypi-31./2.*ypi**2+17. # /2.*ypi**3) ww_0= ww_0+rwm**2*ymi*rlnc0*(4*ypi-8*ypi**2+4*ypi**3) ww_0= ww_0+rwm**2*ymi*fh1t*(-4*ypi+8*ypi**2-4*ypi**3 # ) ww_0= ww_0+rwm**2*ymi*bln*(-4*ypi+8*ypi**2-4*ypi**3) ww_0= ww_0+rwm**2*ymi*(-4./3.+14*ypi-28*ypi**2+44./ # 3.*ypi**3) ww_0= ww_0+rwm**6*rlns*ypmwi*ymi*(-4./3.*ypi+2*ypi**2 # -2./3.*ypi**3) ww_0= ww_0+rwm**6*rlnt*ypmwi*ymi*(4./3.*ypi-2*ypi**2+2./ # 3.*ypi**3) ww_0= ww_0+ymi*(2*s0*ypi**2-4./3.*s0*ypi**3-2./3.*s0) * ww_1= # +ymvvi*rwm**3*atns*ymi*(2*ypi-6*ypi**2+6*ypi**3-2* # ypi**4) ww_1= ww_1+ymvvi*rwm**3*atnt*ymi*rlnc0*(-6*ypi**2+12* # ypi**3-6*ypi**4) ww_1= ww_1+ymvvi*rwm**3*atnt*ymi*(-15*ypi**2+36*ypi**3 # -20*ypi**4) ww_1= ww_1+ymvvi*rwm**3*ymi*fh0t*(6*ypi**2-12*ypi**3+6* # ypi**4) ww_1= ww_1+ymvvi*rwm**7*atns*ypmwi*ymi*(16./3.*ypi**2-8* # ypi**3+8./3.*ypi**4) ww_1= ww_1+ymvvi*rwm**7*atnt*ypmwi*ymi*(16./3.*ypi**2-8* # ypi**3+8./3.*ypi**4) ww_1= ww_1+vv*rwm**3*atns*ypmwi*(-4./3.+2*ypi-2./3.* # ypi**2) ww_1= ww_1+vv*rwm**3*atnt*ypmwi*(-4./3.+2*ypi-2./3.* # ypi**2) ww_1= ww_1+rwm*atns*ymi*(-2./3.+2*ypi-2*ypi**2+2./3. # *ypi**3) ww_1= ww_1+rwm*atnt*ymi*rlnc0*(-2*ypi+4*ypi**2-2* # ypi**3) ww_1= ww_1+rwm*atnt*ymi*(-7*ypi+16*ypi**2-26./3.* # ypi**3) ww_1= ww_1+rwm*ymi*fh0t*(2*ypi-4*ypi**2+2*ypi**3) * ww_2= # +ymvvi*rwm**2*rlns*ymi*(1./3.*ypi-ypi**2+ypi**3-1./3. # *ypi**4) ww_2= ww_2+ymvvi*rwm**2*rlnt*ymi*rlnc0*(-ypi**2+2* # ypi**3-ypi**4) ww_2= ww_2+ymvvi*rwm**2*rlnt*ymi*(-5./2.*ypi**2+6* # ypi**3-10./3.*ypi**4) ww_2= ww_2+ymvvi*rwm**2*ymi*fh1t*(2*ypi**2-4*ypi**3+2* # ypi**4) * else if(ipro.eq.2) then ** * ww_m1= * # +ymvvi*rwm**5*atns*ymi*(-2./3.*ypi+2./3.*ypi**2) * ww_m1= ww_m1+ymvvi*rwm**5*atnt*ymi*rlnc0*(2./9.*ypi**2-4./ * # 9.*ypi**3+10./9.*ypi**4) * ww_m1= ww_m1+ymvvi*rwm**5*atnt*ymi*(5./9.*ypi**2-2*ypi**3 * # +10./3.*ypi**4) * ww_m1= ww_m1+ymvvi*rwm**5*ymi*fh0t*(-2./9.*ypi**2+4./9. * # *ypi**3-10./9.*ypi**4) * ww_m1= ww_m1+vv*rwm*atns*(1-1./3.*yp-4./3.*ypi+2./3.* * # ypi**2) * ww_m1= ww_m1+vv*rwm*atnt*rlnc0*(1./9.-2./9.*ypi+5./9.* * # ypi**2) * ww_m1= ww_m1+vv*rwm*atnt*(1./2.-ypi+7./3.*ypi**2) * ww_m1= ww_m1+vv*rwm*fh0t*(-1./9.+2./9.*ypi-5./9.* * # ypi**2) * ww_m1= ww_m1+rwm**3*atns*ymi*(2./3.-2./3.*ypi) * ww_m1= ww_m1+rwm**3*atnt*ymi*rlnc0*(2./9.*ypi-4./9.* * # ypi**2+10./9.*ypi**3) * ww_m1= ww_m1+rwm**3*atnt*ymi*(7./9.*ypi-22./9.*ypi**2+ * # 40./9.*ypi**3) * ww_m1= ww_m1+rwm**3*ymi*fh0t*(-2./9.*ypi+4./9.*ypi**2 * # -10./9.*ypi**3) ** * ww_0= * # +ymvv*vv*rwm**2*rlns*ypmwi*(-1./6.+1./6.*ypi) * ww_0= ww_0+ymvv*vv*rwm**2*rlnt*ypmwi*(1./6.-1./6.*ypi) * ww_0= ww_0+ymvvi*rwm**2*ymi*rlnc0*(-4./9.*s0*ypi**2+8./ * # 9.*s0*ypi**3-20./9.*s0*ypi**4) * ww_0= ww_0+ymvvi*rwm**2*ymi*aln*(4./9.*s0*ypi**2-8./9.*s0 * # *ypi**3+20./9.*s0*ypi**4) * ww_0= ww_0+ymvvi*rwm**2*ymi*(4./3.*s0*ypi-22./9.*s0* * # ypi**2+20./3.*s0*ypi**3-28./3.*s0*ypi**4) * ww_0= ww_0+ymvvi*rwm**4*rlns*ymi*(-ypi+ypi**2-ypi**3 * # +ypi**4) * ww_0= ww_0+ymvvi*rwm**4*rlnt*ymi*rlnc0*(1./3.*ypi**2-2./3. * # *ypi**3+5./3.*ypi**4) * ww_0= ww_0+ymvvi*rwm**4*rlnt*ymi*(5./6.*ypi**2-4*ypi**3 * # +6*ypi**4) * ww_0= ww_0+ymvvi*rwm**4*ymi*fh1t*(-2./3.*ypi**2+4./3.* * # ypi**3-10./3.*ypi**4) * ww_0= ww_0+ymvvi*ymi*rlnc0*(-1./9.*s0**2*ypi**2+2./9.* * # s0**2*ypi**3-5./9.*s0**2*ypi**4) * ww_0= ww_0+ymvvi*ymi*aln*(1./9.*s0**2*ypi**2-2./9.*s0**2* * # ypi**3+5./9.*s0**2*ypi**4) * ww_0= ww_0+ymvvi*ymi*(1./3.*s0**2*ypi-11./18.*s0**2* * # ypi**2+5./3.*s0**2*ypi**3-7./3.*s0**2*ypi**4) * ww_0= ww_0+vv*rlns*(1./2.-1./6.*yp-2./3.*ypi+1./3.* * # ypi**2) * ww_0= ww_0+vv*rlnt*rlnc0*(1./18.-1./9.*ypi+5./18.* * # ypi**2) * ww_0= ww_0+vv*rlnt*(1./4.-1./2.*ypi+7./6.*ypi**2) * ww_0= ww_0+vv*rlnc0*(1./9.-2./9.*ypi+5./9.*ypi**2) * ww_0= ww_0+vv*fh1t*(-1./9.+2./9.*ypi-5./9.*ypi**2) * ww_0= ww_0+vv*bln*(-1./9.+2./9.*ypi-5./9.*ypi**2) * ww_0= ww_0+vv*(1./6.+1./3.*yp-13./9.*ypi+25./9.* * # ypi**2) * ww_0= ww_0+rwm**2*rlns*ymi*(2./3.-2./3.*ypi+1./2.* * # ypi**2-1./2.*ypi**3) * ww_0= ww_0+rwm**2*rlnt*ymi*rlnc0*(2./9.*ypi-4./9.*ypi**2 * # +10./9.*ypi**3) * ww_0= ww_0+rwm**2*rlnt*ymi*(7./9.*ypi-53./18.*ypi**2+89. * # /18.*ypi**3) * ww_0= ww_0+rwm**2*ymi*rlnc0*(4./9.*ypi-8./9.*ypi**2+20./ * # 9.*ypi**3) * ww_0= ww_0+rwm**2*ymi*fh1t*(-4./9.*ypi+8./9.*ypi**2- * # 20./9.*ypi**3) * ww_0= ww_0+rwm**2*ymi*bln*(-4./9.*ypi+8./9.*ypi**2-20. * # /9.*ypi**3) * ww_0= ww_0+rwm**2*ymi*(-4./3.+22./9.*ypi-20./3.* * # ypi**2+28./3.*ypi**3) * ww_0= ww_0+rwm**6*rlns*ypmwi*ymi*(2./3.*ypi**2-2./3.* * # ypi**3) * ww_0= ww_0+rwm**6*rlnt*ypmwi*ymi*(-2./3.*ypi**2+2./3.* * # ypi**3) * ww_0= ww_0+ymi*(4./9.*s0*ypi-2./9.*s0*ypi**2-4./9.*s0* * # ypi**3-2./3.*s0) ** * ww_1= * # +ymvvi*rwm**3*atns*ymi*(2*ypi-2*ypi**2+2*ypi**3-2* * # ypi**4) * ww_1= ww_1+ymvvi*rwm**3*atnt*ymi*rlnc0*(-2./3.*ypi**2+4. * # /3.*ypi**3-10./3.*ypi**4) * ww_1= ww_1+ymvvi*rwm**3*atnt*ymi*(-5./3.*ypi**2+8* * # ypi**3-12*ypi**4) * ww_1= ww_1+ymvvi*rwm**3*ymi*fh0t*(2./3.*ypi**2-4./3.* * # ypi**3+10./3.*ypi**4) * ww_1= ww_1+ymvvi*rwm**7*atns*ypmwi*ymi*(-8./3.*ypi**3+8. * # /3.*ypi**4) * ww_1= ww_1+ymvvi*rwm**7*atnt*ypmwi*ymi*(-8./3.*ypi**3+8. * # /3.*ypi**4) * ww_1= ww_1+vv*rwm**3*atns*ypmwi*(2./3.*ypi-2./3.*ypi**2) * ww_1= ww_1+vv*rwm**3*atnt*ypmwi*(2./3.*ypi-2./3.*ypi**2) * ww_1= ww_1+rwm*atns*ymi*(-2./3.+2./3.*ypi-2./3.* * # ypi**2+2./3.*ypi**3) * ww_1= ww_1+rwm*atnt*ymi*rlnc0*(-2./9.*ypi+4./9.*ypi**2 * # -10./9.*ypi**3) * ww_1= ww_1+rwm*atnt*ymi*(-7./9.*ypi+28./9.*ypi**2-46./ * # 9.*ypi**3) * ww_1= ww_1+rwm*ymi*fh0t*(2./9.*ypi-4./9.*ypi**2+10./9.* * # ypi**3) ** * ww_2= * # +ymvvi*rwm**2*rlns*ymi*(1./3.*ypi-1./3.*ypi**2+1./3.* * # ypi**3-1./3.*ypi**4) * ww_2= ww_2+ymvvi*rwm**2*rlnt*ymi*rlnc0*(-1./9.*ypi**2+2. * # /9.*ypi**3-5./9.*ypi**4) * ww_2= ww_2+ymvvi*rwm**2*rlnt*ymi*(-5./18.*ypi**2+4./3.* * # ypi**3-2*ypi**4) * ww_2= ww_2+ymvvi*rwm**2*ymi*fh1t*(2./9.*ypi**2-4./9.* * # ypi**3+10./9.*ypi**4) ** ww_m1= # +ymvv*vv*atnt/rwm*(-1+ypi) ww_m1= ww_m1+ymvvi*rwm**5*atnt*ymi*rlnc0*(-4./9.*ypi**2 # +8./9.*ypi**3-8./9.*ypi**4) ww_m1= ww_m1+ymvvi*rwm**5*atnt*ymi*(-4./9.*ypi**2+8./3. # *ypi**3-8./3.*ypi**4) ww_m1= ww_m1+ymvvi*rwm**5*ymi*fh0t*(4./9.*ypi**2-8./9.* # ypi**3+8./9.*ypi**4) ww_m1= ww_m1+vv*rwm*atnt*rlnc0*(1./9.-2./9.*ypi+5./9.* # ypi**2) ww_m1= ww_m1+vv*rwm*atnt*(1./3.-5./3.*ypi+8./3.*ypi**2 # ) ww_m1= ww_m1+vv*rwm*fh0t*(-1./9.+2./9.*ypi-5./9.* # ypi**2) ww_m1= ww_m1+rwm**3*atnt*ymi*rlnc0*(-4./9.*ypi+8./9.* # ypi**2-8./9.*ypi**3) ww_m1= ww_m1+rwm**3*atnt*ymi*(-2./9.*ypi+20./9.*ypi**2 # -14./9.*ypi**3) ww_m1= ww_m1+rwm**3*ymi*fh0t*(4./9.*ypi-8./9.*ypi**2+8./ # 9.*ypi**3) * ww_0= # +ymvvi*rwm**2*ymi*rlnc0*(8./9.*s0*ypi**2-16./9.*s0*ypi**3 # +16./9.*s0*ypi**4) ww_0= ww_0+ymvvi*rwm**2*ymi*aln*(-8./9.*s0*ypi**2+16./9. # *s0*ypi**3-16./9.*s0*ypi**4) ww_0= ww_0+ymvvi*rwm**2*ymi*(8./9.*s0*ypi**2-16./3.*s0* # ypi**3+16./3.*s0*ypi**4) ww_0= ww_0+ymvvi*rwm**4*rlnt*ymi*rlnc0*(-2./3.*ypi**2+4. # /3.*ypi**3-4./3.*ypi**4) ww_0= ww_0+ymvvi*rwm**4*rlnt*ymi*(-2./3.*ypi**2+4* # ypi**3-4*ypi**4) ww_0= ww_0+ymvvi*rwm**4*ymi*fh1t*(4./3.*ypi**2-8./3.* # ypi**3+8./3.*ypi**4) ww_0= ww_0+ymvvi*ymi*rlnc0*(2./9.*s0**2*ypi**2-4./9.* # s0**2*ypi**3+4./9.*s0**2*ypi**4) ww_0= ww_0+ymvvi*ymi*aln*(-2./9.*s0**2*ypi**2+4./9.* # s0**2*ypi**3-4./9.*s0**2*ypi**4) ww_0= ww_0+ymvvi*ymi*(2./9.*s0**2*ypi**2-4./3.*s0**2* # ypi**3+4./3.*s0**2*ypi**4) ww_0= ww_0+vv*rlnt*rlnc0*(1./18.-1./9.*ypi+5./18.* # ypi**2) ww_0= ww_0+vv*rlnt*(1./6.-5./6.*ypi+4./3.*ypi**2) ww_0= ww_0+vv*rlnc0*(-2./9.+4./9.*ypi-4./9.*ypi**2) ww_0= ww_0+vv*fh1t*(-1./9.+2./9.*ypi-5./9.*ypi**2) ww_0= ww_0+vv*bln*(2./9.-4./9.*ypi+4./9.*ypi**2) ww_0= ww_0+vv*(8./9.*ypi-2./9.*ypi**2) ww_0= ww_0+rwm**2*rlnt*ymi*rlnc0*(-4./9.*ypi+8./9.* # ypi**2-8./9.*ypi**3) ww_0= ww_0+rwm**2*rlnt*ymi*(-2./9.*ypi+20./9.*ypi**2- # 14./9.*ypi**3) ww_0= ww_0+rwm**2*ymi*rlnc0*(-8./9.*ypi+16./9.*ypi**2 # -16./9.*ypi**3) ww_0= ww_0+rwm**2*ymi*fh1t*(8./9.*ypi-16./9.*ypi**2+16./ # 9.*ypi**3) ww_0= ww_0+rwm**2*ymi*bln*(8./9.*ypi-16./9.*ypi**2+16./ # 9.*ypi**3) ww_0= ww_0+rwm**2*ymi*(-8./9.*ypi+16./3.*ypi**2-16./3. # *ypi**3) ww_0= ww_0+ymi*(-2./9.*s0*ypi+4./9.*s0*ypi**2-10./9.* # s0*ypi**3) * ww_1= # +ymvvi*rwm**3*atnt*ymi*rlnc0*(4./3.*ypi**2-8./3.*ypi**3 # +8./3.*ypi**4) ww_1= ww_1+ymvvi*rwm**3*atnt*ymi*(4./3.*ypi**2-8*ypi**3 # +8*ypi**4) ww_1= ww_1+ymvvi*rwm**3*ymi*fh0t*(-4./3.*ypi**2+8./3.* # ypi**3-8./3.*ypi**4) ww_1= ww_1+rwm*atnt*ymi*rlnc0*(4./9.*ypi-8./9.*ypi**2+8. # /9.*ypi**3) ww_1= ww_1+rwm*atnt*ymi*(2./9.*ypi-20./9.*ypi**2+14./9. # *ypi**3) ww_1= ww_1+rwm*ymi*fh0t*(-4./9.*ypi+8./9.*ypi**2-8./9. # *ypi**3) * ww_2= # +ymvvi*rwm**2*rlnt*ymi*rlnc0*(2./9.*ypi**2-4./9.*ypi**3 # +4./9.*ypi**4) ww_2= ww_2+ymvvi*rwm**2*rlnt*ymi*(2./9.*ypi**2-4./3.* # ypi**3+4./3.*ypi**4) ww_2= ww_2+ymvvi*rwm**2*ymi*fh1t*(-4./9.*ypi**2+8./9.* # ypi**3-8./9.*ypi**4) * endif * www= (ww_m1+rwg*(ww_0+rwg*(ww_1+rwg*ww_2)))/rwg * if((www*tjac*tfact).lt.0.d0) then iz= 0 ifz(5)= ifz(5)+1 ifzzf= ifzzf+1 endif * 1 if(iz.eq.0) then wwcc20= 0.d0 else wwcc20= www*tjac*tfact endif * return end * *-----wtoregion----------------------------------------------------- * subroutine wtoregion(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 * *-----wtoqcdlam----------------------------------------------------------------- *-----computes lamba^nf_msbar from alpha_s(rs0) * real*8 function wtoqcdlam(nf,als,rs,x1,x2,xacc) implicit real*8(a-h,o-z) * parameter (jmax=50,nout=21) * fmid= wtoqcdscale(nf,als,rs,x2) f= wtoqcdscale(nf,als,rs,x1) if (f*fmid.ge.0.d0) then write(nout,fmt=*) 'root must be bracketed for bisection' print*,'root must be bracketed for bisection' write(nout,fmt=1) als print 1,als 1 format(/' error detected by wtoqcdlam ',/ # ' current value of alpha_s = ',e20.5) stop endif if (f.lt.0.d0) then wtoqcdlam= x1 dx= x2-x1 else wtoqcdlam= x2 dx= x1-x2 endif do j=1,jmax dx= dx*0.5d0 xmid= wtoqcdlam+dx fmid= wtoqcdscale(nf,als,rs,xmid) if (fmid.le.0.d0) wtoqcdlam= xmid if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return enddo pause 'too many bisections' end * *-----wtoqcdscale------------------------------------------------------------- *-----computes lamba^nf_msbar from alpha_s(rs0) * real*8 function wtoqcdscale(nf,als,rs,x) implicit real*8(a-h,o-z) * common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis, # pih,tfact,eps,scale * 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) wtoqcdscale= 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 * *-----wtoralphas-------------------------------------------------------------- *-----computes alpha_s(rs) for nf flavors given als= alpha_s(wm) * real*8 function wtoralphas(rs0,rs,als,nf) implicit real*8(a-h,o-z) * common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis, # pih,tfact,eps,scale * dimension b0(5),b1(5),b2(5) * *-----limits for lambda_5 are 1 mev < lambda_5 < 10 gev * bqm= 4.7d0 cqm= 1.55d0 * if(als.eq.0.d0) then wtoralphas= 0.d0 else x1= 0.001d0 x2= 10.0d0 xacc= 1.d-12 qcdl= wtoqcdlam(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) * wtoralphas= 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 * *---------------------------------------------------------------------- * subroutine wclnx(arg,res) implicit real*8(a-h,o-z) * common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis, # pih,tfact,eps,scale * dimension arg(2),aarg(2),res(2) * 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= 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 wclnomx(arg,omarg,res) implicit real*8(a-h,o-z) * common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis, # pih,tfact,eps,scale * dimension arg(2),aarg(2),omarg(2),res(2),ares(2), # ct(10),sn(10) * 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 yr= omarg(1) yi= omarg(2) ym2= yr*yr+yi*yi ym= sqrt(ym2) do i= 1,2 aarg(i)= abs(omarg(i)) enddo res(1)= log(ym) if(yr.eq.0.d0) then if(yi.gt.0.d0) then teta= pi/2.d0 else teta= -pi/2.d0 endif ares(2)= teta else if(yi.eq.0.d0) then if(yr.gt.0.d0) then teta= 0.d0 else teta= pi endif ares(2)= teta else tnteta= aarg(2)/aarg(1) teta= atan(tnteta) sr= yr/aarg(1) si= yi/aarg(2) if(sr.gt.0.d0) then ares(2)= si*teta else ares(2)= si*(pi-teta) endif endif endif do i= 1,2 res(i)= ares(i) enddo * return end * *------------------------------------------------------------------------ * subroutine wcli2(arg,fres) implicit real*8(a-h,o-z) * common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis, # pih,tfact,eps,scale * dimension arg(2),fres(2) dimension b(0:14),bf(0:14) 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) dimension x(2),omx(2),y(2),oy(2),omy(2),z(2),omz(2),t(2),omt(2) * xr= arg(1) xi= arg(2) omxr= 1.d0-arg(1) 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 wclnx(x,clnx) call wclnomx(x,omx,clnomx) add1(1)= pis/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 wclnx(oy,clnoy) add2(1)= -pis/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 wclnx(z,clnz) call wclnomx(z,omz,clnomz) add3(1)= pis/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 wclnomx(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) * fres(1)= cli2r fres(2)= cli2i * return end * *------------------------------------------------------------------------ * real*8 function fh0(z,a,b) implicit real*8(a-h,o-z) * dimension arg(2),rln1(2),rln2(2),sp1(2),sp2(2) * a2b2= a*a+b*b * if((z+b).lt.0.d0) then stop else arg(1)= 1.d0-(z+b)/a2b2*b arg(2)= (z+b)/a2b2*a call wclnx(arg,rln1) * arg(1)= 1.d0-b/a2b2*b arg(2)= b/a2b2*a call wclnx(arg,rln2) * arg(1)= (z+b)/a2b2*b arg(2)= -(z+b)/a2b2*a call wcli2(arg,sp1) * arg(1)= b*b/a2b2 arg(2)= -b*a/a2b2 call wcli2(arg,sp2) * fh0= log(z+b)*rln1(2)+sp1(2)-log(b)*rln2(2)-sp2(2) endif * return end * *------------------------------------------------------------------------ * real*8 function fh1(z,a,b) implicit real*8(a-h,o-z) * dimension arg(2),rln1(2),sp1(2),sp2(2) * a2b2= a*a+b*b if((z+b).lt.0.d0) then stop else arg(1)= -b arg(2)= a call wclnx(arg,rln1) * arg(1)= (z+b)/a2b2*b arg(2)= (z+b)/a2b2*a call wcli2(arg,sp1) * arg(1)= b*b/a2b2 arg(2)= b*a/a2b2 call wcli2(arg,sp2) * fh1= 0.5d0*(log(z+b)*log(z*z+a*a)-log(b)*log(a*a))- # (rln1(1)*(log(z+b)-log(b))-sp1(1)+sp2(1)) endif * return end