subroutine ww(iosf,inpts,inrand,wmi,umi,dmi,rsi, # s0i,tc,alszi,sig,esig) implicit real*8(a-h,o-z) character*1 oqcd,oral * parameter (nin=98,nout=99,ndimmx=6) * common/wtcut/tcs,s0 common/wpro/ipro,isf common/wtifz/ifz(7),ifzzf,ik common/wtapar/s,rwm,rwg,g2,sth2 common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis, # pih,tfact,eps,em,em2,um,dm * dimension vk(ndimmx) * external d01gcf,wtoregion,wwcc20 * do i=1,7 ifz(i)= 0 enddo ifzzf= 0 ik= 0 iqcd= 1 oqcd= 'n' oral= 'r' ios= 2 ipro= 2 * eps= 1.d-37 cfct= 0.38937966d9 alphai= 137.0359895d0 * em= 0.51099906d-3 em= 0.5116d-3 um= umi dm= dmi pi= 3.141592653589793238462643d0 pis= pi*pi pih= 0.5d0*pi ge= 0.5772156649d0 gf= 8.2476172696d-6 zm= 91.189d0 rs= rsi s= rs*rs 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) * 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 * 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/2048.d0/pi**5*fcol*cfct * if(iosf.eq.1) then ndim= 6 else if(iosf.eq.0) then ndim= 4 endif if(inpts.le.6) then jnpts= inpts else if(inpts.eq.7) then jnpts= 2011*101 if(ndim.eq.4) then vk(1)= 0.1d+01 vk(2)= 0.12294d+05 vk(3)= 0.27852d+05 vk(4)= 0.170453d+06 else if(ndim.eq.6) then vk(1)= 0.1d+01 vk(2)= 0.9463d+05 vk(3)= 0.79132d+05 vk(4)= 0.167923d+06 vk(5)= 0.164405d+06 vk(6)= 0.154994d+06 endif else if(inpts.eq.8) then jnpts= 10193*101 if(ndim.eq.4) then vk(1)= 0.1d+01 vk(2)= 0.9601d+06 vk(3)= 0.449688d+06 vk(4)= 0.792432d+06 endif endif * 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,7 ifzt= ifzt+ifz(i) enddo print 6,ifail print 1,ik print 2,1.d2*(1.d0-1.d0*ifzt/ik) print 3,ifzzf do i=1,7 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) 6 format(/' on exit IFAIL= ',i1) return end * *------------------------------------------------------------------ * real*8 function wwcc20(ndim,x) implicit real*8(a-h,o-z) * common/wtcut/tcs,s0 common/wpro/ipro,isf common/wtifz/ifz(7),ifzzf,ik common/wtapar/s,rwm,rwg,g2,sth2 common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis, # pih,tfact,eps,em,em2,um,dm * dimension x(ndim),ww(5) * external s07aaf * *-----order of integration is: uv,vv; ym; yp; zv; xv * do ii=1,5 ww(ii)= 0.d0 enddo ik= ik+1 em4= em2*em2 em6= em4*em2 s2= s*s s3= s2*s rem= em/sqrt(s) rwms= rwm*rwm rwmq= rwms*rwms rwgs= rwg*rwg rwgm= rwg*rwm rwgms= rwgs*rwms rwmpgs= rwms+rwgs * if(ndim.eq.4) then ymx= x(1) ypx= x(2) zx= x(3) xx= x(4) else if(ndim.eq.6) then uvx= x(1) vvx= x(2) ymx= x(3) ypx= x(4) zx= x(5) xx= x(6) endif * iz= 1 * if(ndim.eq.4) 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= 0.25d0*(rem+sqrt(rem*rem+4.d0*s0))* # (rem+sqrt(rem*rem+4.d0*s0)) 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 shs= sh*sh ssh= sqrt(sh) vvs= vv*vv vvi= 1.d0/vv vvis= vvi*vvi vvic= vvis*vvi vviq= vvic*vvi * 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= (4.d0*s0+tcs)/(4.d0+tcs)/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 ymvvs= ymvv*ymvv ymvvc= ymvvs*ymvv ymvvq= ymvvc*ymvv omym= (1.d0-s0/vv)*(1.d0-ymx)+em/ssh*ymx omyms= omym*omym omymc= omyms*omym * *-----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 ypiv= ypiq*ypi ypivi= ypiv*ypi omyp= 1.d0-yp omyps= omyp*omyp omypc= omyps*omyp * qh0s= em2*(xm*ym+1.d0-xm)*ym/omym qhcs= qh0s+0.25d0*xm*xm*omym*tcs*s * z0= qh0s/ym/sh zc0= qhcs/ym/sh zc1= 1.d0-s0/ymvv zc2= yp-s0/ymvv zc= dmin1(zc0,zc1,zc2) * if(z0.gt.zc) then iz= 0 ifz(5)= ifz(5)+1 go to 1 endif * *-----Loop for integration starts here * z is initialized * do ii=1,5 * do ii=1,3 * if(ii.eq.1) then zjac= (zc-z0)/zc/z0/ym/sh zv= z0*zc/((zc-z0)*zx+z0) else if(ii.eq.2.or.ii.eq.3) then zjac= log(zc/z0) zv= exp(zjac*zx+log(z0)) else if(ii.eq.4) then zjac= (zc-z0)/zc/z0/(ym*sh)**2 zv= z0*zc/((zc-z0)*zx+z0) else if(ii.eq.5) then zjac= log(zc/z0)/ym/sh zv= exp(zjac*zx+log(z0)) endif * *-----x is initialized * xl1= s0/ymvv+zv xl2= (1.d0-0.5d0/zv)*yp xl= dmax1(xl1,xl2) xu= yp if(xl.gt.xu) then iz= 0 ifz(6)= ifz(6)+1 go to 1 endif * atma= ((xl-zv)*ymvv-rwms)/rwgm atmb= ((xu-zv)*ymvv-rwms)/rwgm * if(atma.gt.1.d0.and.atmb.gt.1.d0) then atma= 1.d0/atma atma= atan(atma) zmat= (pih-atma)/rwgm atmb= 1.d0/atmb atmb= atan(atmb) zmbt= (pih-atmb)/rwgm xjac0= (-atmb+atma)/rwgm else if(atma.gt.1.d0.and.atmb.lt.-1.d0) then atma= 1.d0/atma atma= atan(atma) zmat= (pih-atma)/rwgm atmb= -1.d0/atmb atmb= atan(atmb) zmbt= (-pih+atmb)/rwgm xjac0= (-pi+atmb+atma)/rwgm else if(atma.gt.1.d0.and.abs(atmb).lt.1.d0) then atma= 1.d0/atma atma= atan(atma) zmat= (pih-atma)/rwgm atmb= atan(atmb) zmbt= atmb/rwgm xjac0= (-pih+atmb+atma)/rwgm else if(atma.lt.-1.d0.and.atmb.gt.1.d0) then atma= -1.d0/atma atma= atan(atma) zmat= (-pih+atma)/rwgm atmb= 1.d0/atmb atmb= atan(atmb) zmbt= (pih-atmb)/rwgm xjac0= (pi-atmb-atma)/rwgm else if(atma.lt.-1.d0.and.atmb.lt.-1.d0) then atma= -1.d0/atma atma= atan(atma) zmat= (-pih+atma)/rwgm atmb= -1.d0/atmb atmb= atan(atmb) zmbt= (-pih+atmb)/rwgm xjac0= (atmb-atma)/rwgm else if(atma.lt.-1.d0.and.abs(atmb).lt.1.d0) then atma= -1.d0/atma atma= atan(atma) zmat= (-pih+atma)/rwgm atmb= atan(atmb) zmbt= atmb/rwgm xjac0= (pih+atmb-atma)/rwgm else if(abs(atma).lt.1.d0.and.atmb.gt.1.d0) then atma= atan(atma) zmat= atma/rwgm atmb= 1.d0/atmb atmb= atan(atmb) zmbt= (pih-atmb)/rwgm xjac0= (pih-atmb-atma)/rwgm else if(abs(atma).lt.1.d0.and.atmb.lt.-1.d0) then atma= atan(atma) zmat= atma/rwgm atmb= -1.d0/atmb atmb= atan(atmb) zmbt= (-pih+atmb)/rwgm xjac0= (-pih+atmb-atma)/rwgm else if(abs(atma).lt.1.d0.and.abs(atmb).lt.1.d0) then atma= atan(atma) zmat= atma/rwgm atmb= atan(atmb) zmbt= atmb/rwgm xjac0= (atmb-atma)/rwgm endif * xjac= xjac0/ymvv argx= rwgm*(xjac0*xx+zmat) iftn= 1 s07= s07aaf(argx,iftn) if(iftn.ne.0) then print*,' call to s07aaf, IFAIL = ',iftn endif xv= rwms/ymvv+zv+rwgm/ymvv*s07 xvs= xv*xv xvc= xvs*xv xvq= xvc*xv xvv= xvq*xv xvvi= xvv*xv * pwtm2= ((yp-zv)*ymvv-rwgm*s07)**2+rwms*rwgs * ums= um*um/s dms= dm*dm/s qms= zv*ymvv qmss= qms*qms qps= (yp-xv)*ymvv qpm= -0.5d0*yp*ymvv qpms= qpm*qpm qra= qps/qpm qras= qra*qra ophqra= 1.d0+0.5d0*qra opqra= 1.d0+qra tpqra= 2.d0+qra ophqras= ophqra*ophqra bqs0= qps+2.d0*qpm bqs= bqs0*(1.d0+qms/bqs0) qex= qms/qpm qexs= qex*qex qexx= qms*qex qt= qps*qms/qpms grm2= qpms*(1.d0-qt) grm= sqrt(grm2) opgrm= 1.d0+grm/qpm omgrm= 1.d0-grm/qpm ophqrae= ophqra+0.25d0*qras * umds= ums-dms dmus= dms-ums upds= ums+dms * if(ums.ne.0.d0.and.dms.ne.0.d0) then if(qt.lt.1.d-4) then denu= vvs*(bqs0+qms-umds)/(4.d0*ums*qpms- # bqs0*qps*qms-qps*qmss-2.d0*qps*qms* # upds+2.d0*qpm*qms*umds) dend= vvs*(bqs0+qms+umds)/(4.d0*dms*qpms- # bqs0*qps*qms-qps*qmss-2.d0*qps*qms* # upds-2.d0*qpm*qms*umds) denu= (1.d0+2.d0*bqs*ums/(bqs+umds)**2)*denu dend= (1.d0+2.d0*bqs*dms/(bqs-umds)**2)*dend * unum= qpm*tpqra+2.d0*dms-ums+ # qms*(1.d0-0.5d0*qra*ophqra)-0.25d0*qra* # qexx*ophqrae-0.5d0*ums* # qex+0.5d0*dms*(1.d0-qra)*qex uden= -ums+0.5d0*qms*qra*ophqra+0.25d0*qra* # qexx*ophqrae-0.5d0*ums* # qex+0.5d0*dms*opqra*qex rlncu= log(unum/uden) dnum= qpm*tpqra+2.d0*ums-dms+ # qms*(1.d0-0.5d0*qra*ophqra)-0.25d0*qra* # qexx*ophqrae-0.5d0*dms* # qex+0.5d0*ums*(1.d0-qra)*qex dden= -dms+0.5d0*qms*qra*ophqra+0.25d0*qra* # qexx*ophqrae-0.5d0*dms* # qex+0.5d0*ums*opqra*qex rlncd= log(dnum/dden) else denu= vvs*(bqs-umds)/(4.d0*ums*qpms-qms*qpm*( # qra*tpqra*qpm+2.d0*opqra*dmus+qra*qms)) dend= vvs*(bqs+umds)/(4.d0*dms*qpms-qms*qpm*( # qra*tpqra*qpm+2.d0*opqra*umds+qra*qms)) unum= ophqra*qpm+grm*ophqra+0.5d0*qms*(opgrm+ # dmus/qpm)+dms*opgrm-ums uden= ophqra*(qpm-grm)+0.5d0*qms*(omgrm+ # dmus/qpm)+dms*omgrm-ums dnum= ophqra*qpm+grm*ophqra+0.5d0*qms*(opgrm+ # umds/qpm)+ums*opgrm-dms dden= ophqra*(qpm-grm)+0.5d0*qms*(omgrm+ # umds/qpm)+ums*omgrm-dms rlncu= -log(unum/uden) rlncd= -log(dnum/dden) endif else if(qt.lt.1.d-4) then denu= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms)) dend= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms)) unum= qpm*tpqra+qms*(1.d0-0.5d0*qra*ophqra)- # 0.25d0*qexx*qra*ophqras uden= 0.5d0*qms*qra*(ophqra+0.5d0*qex*ophqras+ # qra/8.d0*qexs*(opqra*ophqra)) if((unum/uden).le.0.d0) then iz= 0 ifz(6)= ifz(6)+1 go to 1 endif rlncu= log(unum/uden) rlncd= rlncd else denu= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms)) dend= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms)) unum= ophqra*qpm+grm*ophqra+0.5d0*qms*opgrm uden= ophqra*(qpm-grm)+0.5d0*qms*omgrm dnum= ophqra*qpm+grm*ophqra+0.5d0*qms*opgrm dden= ophqra*(qpm-grm)+0.5d0*qms*omgrm if((unum/uden).le.0.d0) then iz= 0 ifz(6)= ifz(6)+1 go to 1 endif rlncu= -log(unum/uden) rlncd= -log(dnum/dden) endif endif * bd= -xp*xm+(xp+xm)*em2/s qyjc= 1.d0-em4/s2*(1.d0+(xp+xm*omxm)/bd)+em6/s3/bd*(xp+xm) * if(ii.eq.1) then psf= em2*(2.d0+omxm*(2.d0/ym-1.d0)) else if(ii.eq.2) then psf= -1.d0+2.d0/ym*(1.d0-1.d0/ym) else if(ii.ge.3) then psf= 1.d0 endif * tjac0= ujc*vjc*ymjc*ypjc*zjac*xjac/qyjc tjac= tjac0*stf*psf/sh * if(ii.le.2) then * wcomu1= ypiq*rwms*(-8./9.*xvq) wcomu1= wcomu1+ypic*rwms*(16./9.*xvc) wcomu1= wcomu1+ypis*rwms*(-4./3.*xvs) wcomu1= wcomu1+ypi*rwms*(4./9.*xv) wcomu1= wcomu1+ymvv*ypiq*(4./9.*xvv) wcomu1= wcomu1+ymvv*ypic*(-8./9.*xvq) wcomu1= wcomu1+ymvv*ypis*(2./3.*xvc) wcomu1= wcomu1+ymvv*ypi*(-2./9.*xvs) wcomu1= ymvv*wcomu1 * wcomu2= ypic*(4./9.*xvc) wcomu2= wcomu2+ypis*(-8./9.*xvs) wcomu2= wcomu2+ypi*(2./3.*xv) wcomu2= wcomu2-2./9. wcomu2= ypi*rwms*rwmpgs*wcomu2 * wcomu= wcomu1+wcomu2 * wcomd1= ypiq*rwms*(-2./9.*xvq) wcomd1= wcomd1+ypic*rwms*(4./9.*xvc+4./9.*xvq) wcomd1= wcomd1+ypis*rwms*(-1./3.*xvs-8./ # 9.*xvc-2./9.*xvq) wcomd1= wcomd1+ypi*rwms*(1./9.*xv+2./3.* # xvs+4./9.*xvc) wcomd1= wcomd1+rwms*omyp*(-1./9.*xv) wcomd1= wcomd1+rwms*(-1./9.*xv-1./3.*xvs) wcomd1= wcomd1+ymvv*ypiq*(1./9.*xvv) wcomd1= wcomd1+ymvv*ypic*(-2./9.*xvq-2./9.*xvv) wcomd1= wcomd1+ymvv*ypis*(1./6.*xvc+4./9.*xvq+1./9.*xvv) wcomd1= wcomd1+ymvv*ypi*(-1./18.*xvs-1./3.*xvc-2./9.*xvq) wcomd1= wcomd1+ymvv*omyp*(1./18.*xvs) wcomd1= wcomd1+ymvv*(1./18.*xvs+1./6.*xvc) wcomd1= ymvv*wcomd1 * wcomd2= ypiq*(1./9.*xvc) wcomd2= wcomd2+ypic*(-2./9.*xvs-2./9.*xvc) wcomd2= wcomd2+ypis*(1./6.*xv+4./9.* # xvs+1./9.*xvc) wcomd2= wcomd2+ypi*(-1./18.-1./3.*xv-2./9.*xvs) wcomd2= wcomd2+omyp*(1./18.) wcomd2= wcomd2+(1./18.+1./6.*xv) wcomd2= rwms*rwmpgs*wcomd2 * wcomd= wcomd1+wcomd2 * wcom01= ypiq*rwms*(14./3.*xvq) wcom01= wcom01+ypic*rwms*(-28./3.*xvc-10./3.*xvq) wcom01= wcom01+ypis*rwms*(62./9.*xvs+20./3.* # xvc+11./9.*xvq) wcom01= wcom01+ypi*rwms*(-20./9.*xv-50./9.* # xvs-22./9.*xvc-2./3.*xvq) wcom01= wcom01+rwms*omyp*(11./18.*xv+xvs) wcom01= wcom01+rwms*omyps*(1./3.*xv) wcom01= wcom01+rwms*(23./18.*xv+3./2.*xvs+4./3.*xvc) wcom01= wcom01+ymvv*ypiq*(-7./3.*xvv) wcom01= wcom01+ymvv*ypic*(5*xvq+5./3.*xvv) wcom01= wcom01+ymvv*ypis*(-71./18.*xvc-11./3.*xvq-11./18.*xvv) wcom01= wcom01+ymvv*ypi*(29./18.*xvs+59./18.*xvc # +14./9.*xvq+1./3.*xvv) wcom01= wcom01+ymvv*omyp*(-1./6.*xv-5./36.*xvs-7./6.*xvc) wcom01= wcom01+ymvv*omyps*(-2./3.*xvs) wcom01= wcom01+ymvv*omypc*(-1./6.*xv) wcom01= wcom01+ymvv*(-41./36.*xvs-3./4.*xvc-xvq) wcom01= ymvv*wcom01 * wcom02= ypiq*(-7./3.*xvc) wcom02= wcom02+ypic*(13./3.*xvs+5./3.*xvc) wcom02= wcom02+ypis*(-53./18.*xv-3*xvs-11./18.*xvc) wcom02= wcom02+ypi*(11./18.+41./18.*xv+8./9.*xvs+1./3.*xvc) wcom02= wcom02+omyp*(-5./36.-1./6.*xv) wcom02= wcom02+(-5./36.-3./4.*xv-1./3.*xvs) wcom02= wcom02*rwms*rwmpgs * wcom0= wcom01+wcom02 * denum= denu*ums/vv dendm= dend*dms/vv * wcomum= ypis/vv*rwms*rwmpgs*(4./9.*xvs) wcomum= wcomum+ypi/vv*rwms*rwmpgs*(-8./9.*xv) wcomum= wcomum+1./vv*rwms*rwmpgs*(4./9.) wcomum= wcomum+ymvv*ypis/vv*rwms*(-8./9.*xvc) wcomum= wcomum+ymvv*ypi/vv*rwms*(16./9.*xvs) wcomum= wcomum+ymvv/vv*rwms*(-8./9.*xv) wcomum= wcomum+ymvvs*ypis/vv*(4./9.*xvq) wcomum= wcomum+ymvvs*ypi/vv*(-8./9.*xvc) wcomum= wcomum+ymvvs/vv*(4./9.*xvs) * wcomdm= ypis/vv*rwms*rwmpgs*(1./9.*xvs) wcomdm= wcomdm+ypi/vv*rwms*rwmpgs # *(-2./9.*xv-2./9.*xvs) wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyp*(2./9.*xv) wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyps*(1./9.) wcomdm= wcomdm+1./vv*rwms*rwmpgs*(2./9.*xv+1./9.*xvs) wcomdm= wcomdm+ymvv*ypis/vv*rwms*(-2./9.*xvc) wcomdm= wcomdm+ymvv*ypi/vv*rwms*(4./9.*xvs+4./9.*xvc) wcomdm= wcomdm+ymvv/vv*rwms*omyp*(-4./9.*xvs) wcomdm= wcomdm+ymvv/vv*rwms*omyps*(-2./9.*xv) wcomdm= wcomdm+ymvv/vv*rwms*(-4./9.*xvs-2./9.*xvc) wcomdm= wcomdm+ymvvs*ypis/vv*(1./9.*xvq) wcomdm= wcomdm+ymvvs*ypi/vv*(-2./9.*xvc-2./9.*xvq) wcomdm= wcomdm+ymvvs/vv*omyp*(2./9.*xvc) wcomdm= wcomdm+ymvvs/vv*omyps*(1./9.*xvs) wcomdm= wcomdm+ymvvs/vv*(2./9.*xvc+1./9.*xvq) * wcom= wcomu*rlncu+wcomd*rlncd+wcom0+ # ymvv*(wcomum*denum+wcomdm*dendm) * else if(ii.ge.3) then * dw1= ymvv*ypivi*rlncu*vvis*rwms*(56./3.*xvv) dw1= dw1+ymvv*ypivi*rlncu*rwms*(64./3.*xvv) dw1= dw1+ymvv*ypivi*rlncd*vvis*rwms*(16./3.*xvv) dw1= dw1+ymvv*ypivi*rlncd*rwms*(8*xvv) dw1= dw1+ymvv*ypivi*vvis*rwms*(-500./9.*xvv) dw1= dw1+ymvv*ypivi*rwms*(-232./3.*xvv) dw1= dw1+ymvv*ypiv*rlncu*vvis*rwms*(-48*xvq-136./9.*xvv) dw1= dw1+ymvv*ypiv*rlncu*rwms*(-64*xvq-40./3.*xvv) dw1= dw1+ymvv*ypiv*rlncd*vvis*rwms*(-40./3.* # xvq-104./9.*xvv) dw1= dw1+ymvv*ypiv*rlncd*rwms*(-24.*xvq-16.*xvv) dw1= dw1+ymvv*ypiv*vvis*rwms*(400./3.*xvq+680./9.*xvv) dw1= dw1+ymvv*ypiv*rwms*(232.*xvq+260./3.*xvv) dw1= dw1+ymvv*ypiq*rlncu*vvis*rwms*(362./9.* # xvc+364./9.*xvq) dw1= dw1+ymvv*ypiq*rlncu*rwms*(226./3.*xvc+116./3.*xvq) dw1= dw1+ymvv*ypiq*rlncd*vvis*rwms*(95./9.* # xvc+278./9.*xvq+8*xvv) dw1= dw1+ymvv*ypiq*rlncd*rwms*(85./3.*xvc+142./ # 3.*xvq+28./3.*xvv) dw1= dw1+ymvv*ypiq*vvis*rwms*(-101*xvc- # 1754./9.*xvq-28*xvv) dw1= dw1+ymvv*ypiq*rwms*(-2411./9.*xvc-758./ # 3.*xvq-254./9.*xvv) dw1= dw1+ymvv*ypic*rlncu*vvis*rwms*(-100./9. # *xvs-328./9.*xvc) dw1= dw1+ymvv*ypic*rlncu*rwms*(-376./9.*xvs-380./9.*xvc) dw1= dw1+ymvv*ypic*rlncd*vvis*rwms*(-22./9.* # xvs-250./9.*xvc-212./9.*xvq-16./9.*xvv) dw1= dw1+ymvv*ypic*rlncd*rwms*(-142./9.*xvs- # 496./9.*xvc-80./3.*xvq-4./3.*xvv) dw1= dw1+ymvv*ypic*vvis*rwms*(199./9.*xvs+ # 505./3.*xvc+242./3.*xvq+40./9.*xvv) dw1= dw1+ymvv*ypic*rwms*(1285./9.*xvs+2519./9. # *xvc+236./3.*xvq+10./3.*xvv) dw1= dw1+ymvv*ypis*rlncu*vvis*rwms*(2./9.*xv+12*xvs) dw1= dw1+ymvv*ypis*rlncu*rwms*(86./9.*xv+184./9.*xvs) dw1= dw1+ymvv*ypis*rlncd*vvis*rwms*(-1./9.* # xv+82./9.*xvs+223./9.*xvc+6*xvq) dw1= dw1+ymvv*ypis*rlncd*rwms*(32./9.*xv+274./9. # *xvs+265./9.*xvc+10./3.*xvq) dw1= dw1+ymvv*ypis*vvis*rwms*(29./18.*xv-941. # /18.*xvs-497./6.*xvc-15*xvq-2./3.*xvv) dw1= dw1+ymvv*ypis*rwms*(-593./18.*xv-1270./9. # *xvs-169./2.*xvc-25./3.*xvq) dw1= dw1+ymvv*ypi*rlncu*vvis*rwms*(-8./9.*xv) dw1= dw1+ymvv*ypi*rlncu*rwms*(-4./9.-32./9.*xv) dw1= dw1+ymvv*ypi*rlncd*vvis*rwms*(-2./3.* # xv-98./9.*xvs-68./9.*xvc) dw1= dw1+ymvv*ypi*rlncd*rwms*(-1./9.-62./9.* # xv-46./3.*xvs-8./3.*xvc) dw1= dw1+ymvv*ypi*vvis*rwms*(-1./3.+25./9. # *xv+319./9.*xvs+173./9.*xvc+8./3.*xvq) dw1= dw1+ymvv*ypi*rwms*(23./9.+277./9.*xv+383. # /9.*xvs+29./3.*xvc-2./3.*xvq) dw1= dw1+ymvv*rlncd*vvis*rwms*omyp*(8./9.*xv) dw1= dw1+ymvv*rlncd*vvis*rwms*(7./9.*xv+38./9.*xvs) dw1= dw1+ymvv*rlncd*rwms*omyp*(1./9.) dw1= dw1+ymvv*rlncd*rwms*(1./9.+10./3.*xv+2./3.*xvs) dw1= dw1+ymvv*vvis*rwms*omyp*(-11./9.*xv-8./3.*xvs) dw1= dw1+ymvv*vvis*rwms*omyps*(-2./3.*xv) dw1= dw1+ymvv*vvis*rwms*(1./3.-61./18.*xv-77./ # 9.*xvs-4*xvc) dw1= dw1+ymvv*rwms*omyp*(-11./18.-7./3.*xv+2./3.*xvs) dw1= dw1+ymvv*rwms*omyps*(-1./3.+1./3.*xv) dw1= dw1+ymvv*rwms*(-13./9.-73./9.*xv-20./3.*xvs+xvc) dw1= dw1+ymvvs*ypivi*rlncu*vvis*(-28./3.*xvvi) dw1= dw1+ymvvs*ypivi*rlncu*(-32./3.*xvvi) dw1= dw1+ymvvs*ypivi*rlncd*vvis*(-8./3.*xvvi) dw1= dw1+ymvvs*ypivi*rlncd*(-4*xvvi) dw1= dw1+ymvvs*ypivi*vvis*(250./9.*xvvi) dw1= dw1+ymvvs*ypivi*(116./3.*xvvi) dw1= dw1+ymvvs*ypiv*rlncu*vvis*(24.*xvv+68./9.*xvvi) dw1= dw1+ymvvs*ypiv*rlncu*(104./3.*xvv+20./3.*xvvi) dw1= dw1+ymvvs*ypiv*rlncd*vvis*(20./3.*xvv+52./9.*xvvi) dw1= dw1+ymvvs*ypiv*rlncd*(40./3.*xvv+8*xvvi) dw1= dw1+ymvvs*ypiv*vvis*(-200./3.*xvv-340./9.*xvvi) dw1= dw1+ymvvs*ypiv*(-386./3.*xvv-130./3.*xvvi) dw1= dw1+ymvvs*ypiq*rlncu*vvis*(-178./9.*xvq-188./9.*xvv) dw1= dw1+ymvvs*ypiq*rlncu*(-400./9.*xvq-64./3.*xvv) dw1= dw1+ymvvs*ypiq*rlncd*vvis*(-46./9.* # xvq-142./9.*xvv-4*xvvi) dw1= dw1+ymvvs*ypiq*rlncd*(-157./9.*xvq-80./ # 3.*xvv-14./3.*xvvi) dw1= dw1+ymvvs*ypiq*vvis*(146./3.*xvq+910./9.*xvv+14*xvvi) dw1= dw1+ymvvs*ypiq*(1495./9.*xvq+428./3.* # xvv+127./9.*xvvi) dw1= dw1+ymvvs*ypic*rlncu*vvis*(44./9.*xvc+176./9.*xvq) dw1= dw1+ymvvs*ypic*rlncu*(244./9.*xvc+232./9.*xvq) dw1= dw1+ymvvs*ypic*rlncd*vvis*(8./9.*xvc+ # 131./9.*xvq+112./9.*xvv+8./9.*xvvi) dw1= dw1+ymvvs*ypic*rlncd*(97./9.*xvc+313./9.* # xvq+46./3.*xvv+2./3.*xvvi) dw1= dw1+ymvvs*ypic*vvis*(-68./9.*xvc-91 # *xvq-130./3.*xvv-20./9.*xvvi) dw1= dw1+ymvvs*ypic*(-908./9.*xvc-1612./9.* # xvq-137./3.*xvv-5./3.*xvvi) dw1= dw1+ymvvs*ypis*rlncu*vvis*(2./9.*xvs-20./3.*xvc) dw1= dw1+ymvvs*ypis*rlncu*(-64./9.*xvs-128./9.*xvc) dw1= dw1+ymvvs*ypis*rlncd*vvis*(2./9.*xvs- # 44./9.*xvc-128./9.*xvq-10./3.*xvv) dw1= dw1+ymvvs*ypis*rlncd*(-25./9.*xvs-64./3. # *xvc-58./3.*xvq-2*xvv) dw1= dw1+ymvvs*ypis*vvis*(-23./9.*xvs+ # 523./18.*xvc+49*xvq+25./3.*xvv+1./3.*xvvi) dw1= dw1+ymvvs*ypis*(28*xvs+935./9.*xvc+511./9.*xvq+5*xvv) dw1= dw1+ymvvs*ypi*rlncu*vvis*(4./9.*xvs) dw1= dw1+ymvvs*ypi*rlncu*(4./9.*xv+28./9.*xvs) dw1= dw1+ymvvs*ypi*rlncd*vvis*(1./3.*xvs+ # 64./9.*xvc+43./9.*xvq) dw1= dw1+ymvvs*ypi*rlncd*(1./9.*xv+49./9.*xvs # +101./9.*xvc+2*xvq) dw1= dw1+ymvvs*ypi*vvis*(1./3.*xv-19./18.* # xvs-215./9.*xvc-221./18.*xvq-5./3.*xvv) dw1= dw1+ymvvs*ypi*(-32./9.*xv-163./6.*xvs # -203./6.*xvc-20./3.*xvq+1./3.*xvv) dw1= dw1+ymvvs*rlncd*vvis*omyp*(-7./9.*xvs) dw1= dw1+ymvvs*rlncd*vvis*(-5./9.*xvs-28./9.*xvc) dw1= dw1+ymvvs*rlncd*omyp*(-1./9.*xv) dw1= dw1+ymvvs*rlncd*(-1./9.*xv-8./3.*xvs-2./3.*xvc) dw1= dw1+ymvvs*vvis*omyp*(-1./24.+5./12.*xv- # 7./18.*xvs+10./3.*xvc) dw1= dw1+ymvvs*vvis*omyps*(-2./3.*xv+5./3.*xvs) dw1= dw1+ymvvs*vvis*omypc*(1./3.*xv) dw1= dw1+ymvvs*vvis*(-5./12.*xv+109./36.*xvs # +49./9.*xvc+10./3.*xvq) dw1= dw1+ymvvs*omyp*(5./24.+1./36.*xv+25./6.*xvs-5./6.*xvc) dw1= dw1+ymvvs*omyps*(3./2.*xv-1./2.*xvs) dw1= dw1+ymvvs*omypc*(1./6.-1./6.*xv) dw1= dw1+ymvvs*(85./36.*xv+251./36.*xvs+13./2.*xvc-5./6.*xvq) dw1= dw1+ypivi*rlncu*vvis*rwms*rwmpgs*(-28./3.*xvq) dw1= dw1+ypivi*rlncu*rwms*rwmpgs*(-32./3.*xvq) dw1= dw1+ypivi*rlncd*vvis*rwms*rwmpgs*(-8./3.*xvq) dw1= dw1+ypivi*rlncd*rwms*rwmpgs*(-4*xvq) dw1= dw1+ypivi*vvis*rwms*rwmpgs*(250./9.*xvq) dw1= dw1+ypivi*rwms*rwmpgs*(116./3.*xvq) dw1= dw1+ypiv*rlncu*vvis*rwms*rwmpgs*(24.*xvc+68./9.*xvq) dw1= dw1+ypiv*rlncu*rwms*rwmpgs*(88./3.*xvc+20./3.*xvq) dw1= dw1+ypiv*rlncd*vvis*rwms*rwmpgs*(20./3.*xvc+52./9.*xvq) dw1= dw1+ypiv*rlncd*rwms*rwmpgs*(32./3.*xvc+8*xvq) dw1= dw1+ypiv*vvis*rwms*rwmpgs*(-200./3.* # xvc-340./9.*xvq) dw1= dw1+ypiv*rwms*rwmpgs*(-310./3.*xvc-130./3.*xvq) dw1= dw1+ypiq*rlncu*vvis*rwms*rwmpgs*(-184./ # 9.*xvs-176./9.*xvc) dw1= dw1+ypiq*rlncu*rwms*rwmpgs*(-278./9.*xvs-52./3.*xvc) dw1= dw1+ypiq*rlncd*vvis*rwms*rwmpgs*(-49./9. # *xvs-136./9.*xvc-4*xvq) dw1= dw1+ypiq*rlncd*rwms*rwmpgs*(-98./9.*xvs # -62./3.*xvc-14./3.*xvq) dw1= dw1+ypiq*vvis*rwms*rwmpgs*(157./3.*xvs # +844./9.*xvc+14*xvq) dw1= dw1+ypiq*rwms*rwmpgs*(916./9.*xvs+110*xvc+127./9.*xvq) dw1= dw1+ypic*rlncu*vvis*rwms*rwmpgs*(56./9.* # xv+152./9.*xvs) dw1= dw1+ypic*rlncu*rwms*rwmpgs*(44./3.*xv+148./9.*xvs) dw1= dw1+ypic*rlncd*vvis*rwms*rwmpgs*(14./9.* # xv+119./9.*xvs+100./9.*xvc+8./9.*xvq) dw1= dw1+ypic*rlncd*rwms*rwmpgs*(5*xv+61./3.* # xvs+34./3.*xvc+2./3.*xvq) dw1= dw1+ypic*vvis*rwms*rwmpgs*(-131./9.*xv # -232./3.*xvs-112./3.*xvc-20./9.*xvq) dw1= dw1+ypic*rwms*rwmpgs*(-383./9.*xv-907./9. # *xvs-33*xvc-5./3.*xvq) dw1= dw1+ypis*rlncu*vvis*rwms*rwmpgs*(-4./9.-16./3.*xv) dw1= dw1+ypis*rlncu*rwms*rwmpgs*(-22./9.-56./9.*xv) dw1= dw1+ypis*rlncd*vvis*rwms*rwmpgs*(-1./9. # -38./9.*xv-95./9.*xvs-8./3.*xvc) dw1= dw1+ypis*rlncd*rwms*rwmpgs*(-7./9.-82./9. # *xv-91./9.*xvs-4./3.*xvc) dw1= dw1+ypis*vvis*rwms*rwmpgs*(17./18.+209./ # 9.*xv+203./6.*xvs+20./3.*xvc+1./3.*xvq) dw1= dw1+ypis*rwms*rwmpgs*(107./18.+341./9.*xv # +499./18.*xvs+10./3.*xvc) dw1= dw1+ypi*rlncu*vvis*rwms*rwmpgs*(4./9.) dw1= dw1+ypi*rlncu*rwms*rwmpgs*(4./9.) dw1= dw1+ypi*rlncd*vvis*rwms*rwmpgs*(1./3.+ # 34./9.*xv+25./9.*xvs) dw1= dw1+ypi*rlncd*rwms*rwmpgs*(13./9.+37./9.*xv+2./3.*xvs) dw1= dw1+ypi*vvis*rwms*rwmpgs*(-31./18.- # 104./9.*xv-125./18.*xvs-xvc) dw1= dw1+ypi*rwms*rwmpgs*(-83./18.-169./18.* # xv-3*xvs+1./3.*xvc) dw1= dw1+rlncd*vvis*rwms*rwmpgs*omyp*(-1./9.) dw1= dw1+rlncd*vvis*rwms*rwmpgs*(-2./9.-10./9.*xv) dw1= dw1+rlncd*rwms*rwmpgs*(-2./3.) dw1= dw1+vvis*rwms*rwmpgs*omyp*(5./18.+1./3.*xv) dw1= dw1+vvis*rwms*rwmpgs*(7./9.+22./9.*xv+xvs) dw1= dw1+rwms*rwmpgs*omyp*(1./6.-1./6.*xv) dw1= dw1+rwms*rwmpgs*(11./9.+5./6.*xv-1./6.*xvs) dw1= dw1+denu*ymvvs*ypiq*vviq*rwms*rwmpgs*(-4./9.*xvq) dw1= dw1+denu*ymvvs*ypic*vviq*rwms*rwmpgs*(4./3.*xvc) dw1= dw1+denu*ymvvs*ypis*vviq*rwms*rwmpgs*(-4./3.*xvs) dw1= dw1+denu*ymvvs*ypis*vvis*rwms*rwmpgs*(-2./9.*xvs) dw1= dw1+denu*ymvvs*ypi*vviq*rwms*rwmpgs*(4./9.*xv) dw1= dw1+denu*ymvvs*ypi*vvis*rwms*rwmpgs*(4./9.*xv) dw1= dw1+denu*ymvvs*vvis*rwms*rwmpgs*(-2./9.) dw1= dw1+denu*ymvvc*ypiq*vviq*rwms*(8./9.*xvv) dw1= dw1+denu*ymvvc*ypic*vviq*rwms*(-8./3.*xvq) dw1= dw1+denu*ymvvc*ypis*vviq*rwms*(8./3.*xvc) dw1= dw1+denu*ymvvc*ypis*vvis*rwms*(4./9.*xvc) dw1= dw1+denu*ymvvc*ypi*vviq*rwms*(-8./9.*xvs) dw1= dw1+denu*ymvvc*ypi*vvis*rwms*(-8./9.*xvs) dw1= dw1+denu*ymvvc*vvis*rwms*(4./9.*xv) dw1= dw1+denu*ymvvq*ypiq*vviq*(-4./9.*xvvi) dw1= dw1+denu*ymvvq*ypic*vviq*(4./3.*xvv) dw1= dw1+denu*ymvvq*ypis*vviq*(-4./3.*xvq) dw1= dw1+denu*ymvvq*ypis*vvis*(-2./9.*xvq) dw1= dw1+denu*ymvvq*ypi*vviq*(4./9.*xvc) dw1= dw1+denu*ymvvq*ypi*vvis*(4./9.*xvc) dw1= dw1+denu*ymvvq*vvis*(-2./9.*xvs) dw1= dw1+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-1./9.*xvq) dw1= dw1+dend*ymvvs*ypic*vviq*rwms*rwmpgs # *(1./3.*xvc+2./9.*xvq) dw1= dw1+dend*ymvvs*ypis*vviq*rwms*rwmpgs # *(-1./3.*xvs-2./3.*xvc-1./9.*xvq) dw1= dw1+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./18.*xvs) dw1= dw1+dend*ymvvs*ypi*vviq*rwms*rwmpgs # *(1./9.*xv+2./3.*xvs+1./3.*xvc) dw1= dw1+dend*ymvvs*ypi*vvis*rwms*rwmpgs # *(1./9.*xv+1./9.*xvs) dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-1./9.*xv) dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*(-1./9. # *xv-1./3.*xvs) dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-1./9.*xv) dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./18.) dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*(-1./9.*xv-1./18.*xvs) dw1= dw1+dend*ymvvc*ypiq*vviq*rwms*(2./9.*xvv) dw1= dw1+dend*ymvvc*ypic*vviq*rwms*(-2./3.*xvq-4./9.*xvv) dw1= dw1+dend*ymvvc*ypis*vviq*rwms*(2./3.* # xvc+4./3.*xvq+2./9.*xvv) dw1= dw1+dend*ymvvc*ypis*vvis*rwms*(1./9.*xvc) dw1= dw1+dend*ymvvc*ypi*vviq*rwms*(-2./ # 9.*xvs-4./3.*xvc-2./3.*xvq) dw1= dw1+dend*ymvvc*ypi*vvis*rwms*(-2./9.*xvs-2./9.*xvc) dw1= dw1+dend*ymvvc*vviq*rwms*omyp*(2./9.*xvs) dw1= dw1+dend*ymvvc*vviq*rwms*(2./9.*xvs+2./3.*xvc) dw1= dw1+dend*ymvvc*vvis*rwms*omyp*(2./9.*xvs) dw1= dw1+dend*ymvvc*vvis*rwms*omyps*(1./9.*xv) dw1= dw1+dend*ymvvc*vvis*rwms*(2./9.*xvs+1./9.*xvc) dw1= dw1+dend*ymvvq*ypiq*vviq*(-1./9.*xvvi) dw1= dw1+dend*ymvvq*ypic*vviq*(1./3.*xvv+2./9.*xvvi) dw1= dw1+dend*ymvvq*ypis*vviq*(-1./3.* # xvq-2./3.*xvv-1./9.*xvvi) dw1= dw1+dend*ymvvq*ypis*vvis*(-1./18.*xvq) dw1= dw1+dend*ymvvq*ypi*vviq*(1./9.*xvc # +2./3.*xvq+1./3.*xvv) dw1= dw1+dend*ymvvq*ypi*vvis*(1./9.*xvc+1./9.*xvq) dw1= dw1+dend*ymvvq*vviq*omyp*(-1./9.*xvc) dw1= dw1+dend*ymvvq*vviq*(-1./9.*xvc-1./3.*xvq) dw1= dw1+dend*ymvvq*vvis*omyp*(-1./9.*xvc) dw1= dw1+dend*ymvvq*vvis*omyps*(-1./18.*xvs) dw1= dw1+dend*ymvvq*vvis*(-1./9.*xvc-1./18.*xvq) * dw2= ymvv*ypivi*rlncu*vvis*rwms*(56*xvv) dw2= dw2+ymvv*ypivi*rlncu*rwms*(64./3.*xvv) dw2= dw2+ymvv*ypivi*rlncd*vvis*rwms*(16.*xvv) dw2= dw2+ymvv*ypivi*rlncd*rwms*(8*xvv) dw2= dw2+ymvv*ypivi*vvis*rwms*(-500./3.*xvv) dw2= dw2+ymvv*ypivi*rwms*(-232./3.*xvv) dw2= dw2+ymvv*ypiv*rlncu*vvis*rwms*(-144*xvq-136./3.*xvv) dw2= dw2+ymvv*ypiv*rlncu*rwms*(-64*xvq-40./3.*xvv) dw2= dw2+ymvv*ypiv*rlncd*vvis*rwms*(-40*xvq-104./3.*xvv) dw2= dw2+ymvv*ypiv*rlncd*rwms*(-24.*xvq-16.*xvv) dw2= dw2+ymvv*ypiv*vvis*rwms*(400.*xvq+680./3.*xvv) dw2= dw2+ymvv*ypiv*rwms*(232.*xvq+260./3.*xvv) dw2= dw2+ymvv*ypiq*rlncu*vvis*rwms*(362./3.* # xvc+364./3.*xvq) dw2= dw2+ymvv*ypiq*rlncu*rwms*(226./3.*xvc+116./3.*xvq) dw2= dw2+ymvv*ypiq*rlncd*vvis*rwms*(95./3.* # xvc+278./3.*xvq+24.*xvv) dw2= dw2+ymvv*ypiq*rlncd*rwms*(85./3.*xvc+142./ # 3.*xvq+28./3.*xvv) dw2= dw2+ymvv*ypiq*vvis*rwms*(-303*xvc- # 1754./3.*xvq-84*xvv) dw2= dw2+ymvv*ypiq*rwms*(-2411./9.*xvc-758./ # 3.*xvq-254./9.*xvv) dw2= dw2+ymvv*ypic*rlncu*vvis*rwms*(-100./3. # *xvs-328./3.*xvc) dw2= dw2+ymvv*ypic*rlncu*rwms*(-376./9.*xvs-380./9.*xvc) dw2= dw2+ymvv*ypic*rlncd*vvis*rwms*(-22./3.* # xvs-250./3.*xvc-212./3.*xvq-16./3.*xvv) dw2= dw2+ymvv*ypic*rlncd*rwms*(-142./9.*xvs- # 496./9.*xvc-80./3.*xvq-4./3.*xvv) dw2= dw2+ymvv*ypic*vvis*rwms*(199./3.*xvs+ # 505*xvc+242.*xvq+40./3.*xvv) dw2= dw2+ymvv*ypic*rwms*(1285./9.*xvs+2519./9. # *xvc+236./3.*xvq+10./3.*xvv) dw2= dw2+ymvv*ypis*rlncu*vvis*rwms*(2./3.*xv+36*xvs) dw2= dw2+ymvv*ypis*rlncu*rwms*(86./9.*xv+184./9.*xvs) dw2= dw2+ymvv*ypis*rlncd*vvis*rwms*(-1./3.* # xv+82./3.*xvs+223./3.*xvc+18*xvq) dw2= dw2+ymvv*ypis*rlncd*rwms*(32./9.*xv+274./9. # *xvs+265./9.*xvc+10./3.*xvq) dw2= dw2+ymvv*ypis*vvis*rwms*(29./6.*xv-941./ # 6.*xvs-497./2.*xvc-45*xvq-2*xvv) dw2= dw2+ymvv*ypis*rwms*(-593./18.*xv-1270./9. # *xvs-169./2.*xvc-25./3.*xvq) dw2= dw2+ymvv*ypi*rlncu*vvis*rwms*(-8./3.*xv) dw2= dw2+ymvv*ypi*rlncu*rwms*(-4./9.-32./9.*xv) dw2= dw2+ymvv*ypi*rlncd*vvis*rwms*(-2*xv- # 98./3.*xvs-68./3.*xvc) dw2= dw2+ymvv*ypi*rlncd*rwms*(-1./9.-62./9.* # xv-46./3.*xvs-8./3.*xvc) dw2= dw2+ymvv*ypi*vvis*rwms*(-1+25./3.*xv # +319./3.*xvs+173./3.*xvc+8*xvq) dw2= dw2+ymvv*ypi*rwms*(23./9.+277./9.*xv+383. # /9.*xvs+29./3.*xvc-2./3.*xvq) dw2= dw2+ymvv*rlncd*vvis*rwms*omyp*(8./3.*xv) dw2= dw2+ymvv*rlncd*vvis*rwms*(7./3.*xv+38./3.*xvs) dw2= dw2+ymvv*rlncd*rwms*omyp*(1./9.) dw2= dw2+ymvv*rlncd*rwms*(1./9.+10./3.*xv+2./3.*xvs) dw2= dw2+ymvv*vvis*rwms*omyp*(-11./3.*xv-8*xvs) dw2= dw2+ymvv*vvis*rwms*omyps*(-2*xv) dw2= dw2+ymvv*vvis*rwms*(1.-61./6.*xv-77./3.*xvs-12*xvc) dw2= dw2+ymvv*rwms*omyp*(-11./18.-7./3.*xv+2./3.*xvs) dw2= dw2+ymvv*rwms*omyps*(-1./3.+1./3.*xv) dw2= dw2+ymvv*rwms*(-13./9.-73./9.*xv-20./3.*xvs+xvc) dw2= dw2+ymvvs*ypivi*rlncu*vvis*(-28*xvvi) dw2= dw2+ymvvs*ypivi*rlncu*(-32./3.*xvvi) dw2= dw2+ymvvs*ypivi*rlncd*vvis*(-8*xvvi) dw2= dw2+ymvvs*ypivi*rlncd*(-4*xvvi) dw2= dw2+ymvvs*ypivi*vvis*(250./3.*xvvi) dw2= dw2+ymvvs*ypivi*(116./3.*xvvi) dw2= dw2+ymvvs*ypiv*rlncu*vvis*(72*xvv+68./3.*xvvi) dw2= dw2+ymvvs*ypiv*rlncu*(104./3.*xvv+20./3.*xvvi) dw2= dw2+ymvvs*ypiv*rlncd*vvis*(20*xvv+52./3.*xvvi) dw2= dw2+ymvvs*ypiv*rlncd*(40./3.*xvv+8*xvvi) dw2= dw2+ymvvs*ypiv*vvis*(-200*xvv-340./3.*xvvi) dw2= dw2+ymvvs*ypiv*(-386./3.*xvv-130./3.*xvvi) dw2= dw2+ymvvs*ypiq*rlncu*vvis*(-178./3.*xvq-188./3.*xvv) dw2= dw2+ymvvs*ypiq*rlncu*(-400./9.*xvq-64./3.*xvv) dw2= dw2+ymvvs*ypiq*rlncd*vvis*(-46./3.* # xvq-142./3.*xvv-12*xvvi) dw2= dw2+ymvvs*ypiq*rlncd*(-157./9.*xvq-80./ # 3.*xvv-14./3.*xvvi) dw2= dw2+ymvvs*ypiq*vvis*(146*xvq+910./3.*xvv+42.*xvvi) dw2= dw2+ymvvs*ypiq*(1495./9.*xvq+428./3.*xvv+127./9.*xvvi) dw2= dw2+ymvvs*ypic*rlncu*vvis*(44./3.*xvc+176./3.*xvq) dw2= dw2+ymvvs*ypic*rlncu*(244./9.*xvc+232./9.*xvq) dw2= dw2+ymvvs*ypic*rlncd*vvis*(8./3.*xvc+ # 131./3.*xvq+112./3.*xvv+8./3.*xvvi) dw2= dw2+ymvvs*ypic*rlncd*(97./9.*xvc+313./9.* # xvq+46./3.*xvv+2./3.*xvvi) dw2= dw2+ymvvs*ypic*vvis*(-68./3.*xvc- # 273*xvq-130*xvv-20./3.*xvvi) dw2= dw2+ymvvs*ypic*(-908./9.*xvc-1612./9.* # xvq-137./3.*xvv-5./3.*xvvi) dw2= dw2+ymvvs*ypis*rlncu*vvis*(2./3.*xvs-20*xvc) dw2= dw2+ymvvs*ypis*rlncu*(-64./9.*xvs-128./9.*xvc) dw2= dw2+ymvvs*ypis*rlncd*vvis*(2./3.*xvs- # 44./3.*xvc-128./3.*xvq-10*xvv) dw2= dw2+ymvvs*ypis*rlncd*(-25./9.*xvs-64./3. # *xvc-58./3.*xvq-2*xvv) dw2= dw2+ymvvs*ypis*vvis*(-23./3.*xvs+ # 523./6.*xvc+147*xvq+25*xvv+xvvi) dw2= dw2+ymvvs*ypis*(28*xvs+935./9.*xvc+511./9.*xvq+5*xvv) dw2= dw2+ymvvs*ypi*rlncu*vvis*(4./3.*xvs) dw2= dw2+ymvvs*ypi*rlncu*(4./9.*xv+28./9.*xvs) dw2= dw2+ymvvs*ypi*rlncd*vvis*(xvs+64./3.*xvc+43./3.*xvq) dw2= dw2+ymvvs*ypi*rlncd*(1./9.*xv+49./9.*xvs # +101./9.*xvc+2*xvq) dw2= dw2+ymvvs*ypi*vvis*(xv-19./6.*xvs- # 215./3.*xvc-221./6.*xvq-5*xvv) dw2= dw2+ymvvs*ypi*(-32./9.*xv-163./6.*xvs # -203./6.*xvc-20./3.*xvq+1./3.*xvv) dw2= dw2+ymvvs*rlncd*vvis*omyp*(-7./3.*xvs) dw2= dw2+ymvvs*rlncd*vvis*(-5./3.*xvs-28./3.*xvc) dw2= dw2+ymvvs*rlncd*omyp*(-1./9.*xv) dw2= dw2+ymvvs*rlncd*(-1./9.*xv-8./3.*xvs-2./3.*xvc) dw2= dw2+ymvvs*vvis*omyp*(-1./8.+5./4.*xv-7./6.*xvs+10*xvc) dw2= dw2+ymvvs*vvis*omyps*(-2*xv+5*xvs) dw2= dw2+ymvvs*vvis*omypc*(xv) dw2= dw2+ymvvs*vvis*(-5./4.*xv+109./12.*xvs # +49./3.*xvc+10*xvq) dw2= dw2+ymvvs*omyp*(5./24.+1./36.*xv+25./6.*xvs-5./6.*xvc) dw2= dw2+ymvvs*omyps*(3./2.*xv-1./2.*xvs) dw2= dw2+ymvvs*omypc*(1./6.-1./6.*xv) dw2= dw2+ymvvs*(85./36.*xv+251./36.*xvs+13./2.*xvc-5./6.*xvq) dw2= dw2+ypivi*rlncu*vvis*rwms*rwmpgs*(-28*xvq) dw2= dw2+ypivi*rlncu*rwms*rwmpgs*(-32./3.*xvq) dw2= dw2+ypivi*rlncd*vvis*rwms*rwmpgs*(-8*xvq) dw2= dw2+ypivi*rlncd*rwms*rwmpgs*(-4.*xvq) dw2= dw2+ypivi*vvis*rwms*rwmpgs*(250./3.*xvq) dw2= dw2+ypivi*rwms*rwmpgs*(116./3.*xvq) dw2= dw2+ypiv*rlncu*vvis*rwms*rwmpgs*(72*xvc+68./3.*xvq) dw2= dw2+ypiv*rlncu*rwms*rwmpgs*(88./3.*xvc+20./3.*xvq) dw2= dw2+ypiv*rlncd*vvis*rwms*rwmpgs*(20*xvc+52./3.*xvq) dw2= dw2+ypiv*rlncd*rwms*rwmpgs*(32./3.*xvc+8*xvq) dw2= dw2+ypiv*vvis*rwms*rwmpgs*(-200*xvc-340./3.*xvq) dw2= dw2+ypiv*rwms*rwmpgs*(-310./3.*xvc-130./3.*xvq) dw2= dw2+ypiq*rlncu*vvis*rwms*rwmpgs*(-184./ # 3.*xvs-176./3.*xvc) dw2= dw2+ypiq*rlncu*rwms*rwmpgs*(-278./9.*xvs-52./3.*xvc) dw2= dw2+ypiq*rlncd*vvis*rwms*rwmpgs*(-49./3. # *xvs-136./3.*xvc-12*xvq) dw2= dw2+ypiq*rlncd*rwms*rwmpgs*(-98./9.*xvs # -62./3.*xvc-14./3.*xvq) dw2= dw2+ypiq*vvis*rwms*rwmpgs*(157*xvs+844./3.*xvc+42.*xvq) dw2= dw2+ypiq*rwms*rwmpgs*(916./9.*xvs+110*xvc+127./9.*xvq) dw2= dw2+ypic*rlncu*vvis*rwms*rwmpgs*(56./3.*xv+152./3.*xvs) dw2= dw2+ypic*rlncu*rwms*rwmpgs*(44./3.*xv+148./9.*xvs) dw2= dw2+ypic*rlncd*vvis*rwms*rwmpgs*(14./3.* # xv+119./3.*xvs+100./3.*xvc+8./3.*xvq) dw2= dw2+ypic*rlncd*rwms*rwmpgs*(5*xv+61./3.* # xvs+34./3.*xvc+2./3.*xvq) dw2= dw2+ypic*vvis*rwms*rwmpgs*(-131./3.*xv # -232.*xvs-112*xvc-20./3.*xvq) dw2= dw2+ypic*rwms*rwmpgs*(-383./9.*xv-907./9. # *xvs-33*xvc-5./3.*xvq) dw2= dw2+ypis*rlncu*vvis*rwms*rwmpgs*(-4./3.-16.*xv) dw2= dw2+ypis*rlncu*rwms*rwmpgs*(-22./9.-56./9.*xv) dw2= dw2+ypis*rlncd*vvis*rwms*rwmpgs*(-1./3. # -38./3.*xv-95./3.*xvs-8*xvc) dw2= dw2+ypis*rlncd*rwms*rwmpgs*(-7./9.-82./9. # *xv-91./9.*xvs-4./3.*xvc) dw2= dw2+ypis*vvis*rwms*rwmpgs*(17./6.+209./ # 3.*xv+203./2.*xvs+20*xvc+xvq) dw2= dw2+ypis*rwms*rwmpgs*(107./18.+341./9.*xv # +499./18.*xvs+10./3.*xvc) dw2= dw2+ypi*rlncu*vvis*rwms*rwmpgs*(4./3.) dw2= dw2+ypi*rlncu*rwms*rwmpgs*(4./9.) dw2= dw2+ypi*rlncd*vvis*rwms*rwmpgs*(1+34./3.*xv+25./3.*xvs) dw2= dw2+ypi*rlncd*rwms*rwmpgs*(13./9.+37./9.*xv+2./3.*xvs) dw2= dw2+ypi*vvis*rwms*rwmpgs*(-31./6.- # 104./3.*xv-125./6.*xvs-3*xvc) dw2= dw2+ypi*rwms*rwmpgs*(-83./18.-169./18.*xv-3*xvs+1./3.*xvc) dw2= dw2+rlncd*vvis*rwms*rwmpgs*omyp*(-1./3.) dw2= dw2+rlncd*vvis*rwms*rwmpgs*(-2./3.-10./3.*xv) dw2= dw2+rlncd*rwms*rwmpgs*(-2./3.) dw2= dw2+vvis*rwms*rwmpgs*omyp*(5./6.+xv) dw2= dw2+vvis*rwms*rwmpgs*(7./3.+22./3.*xv+3*xvs) dw2= dw2+rwms*rwmpgs*omyp*(1./6.-1./6.*xv) dw2= dw2+rwms*rwmpgs*(11./9.+5./6.*xv-1./6.*xvs) dw2= dw2+denu*ymvvs*ypiq*vviq*rwms*rwmpgs*(-4./3.*xvq) dw2= dw2+denu*ymvvs*ypic*vviq*rwms*rwmpgs*(4.*xvc) dw2= dw2+denu*ymvvs*ypis*vviq*rwms*rwmpgs*(-4.*xvs) dw2= dw2+denu*ymvvs*ypis*vvis*rwms*rwmpgs*(-2./9.*xvs) dw2= dw2+denu*ymvvs*ypi*vviq*rwms*rwmpgs*(4./3.*xv) dw2= dw2+denu*ymvvs*ypi*vvis*rwms*rwmpgs*(4./9.*xv) dw2= dw2+denu*ymvvs*vvis*rwms*rwmpgs*(-2./9.) dw2= dw2+denu*ymvvc*ypiq*vviq*rwms*(8./3.*xvv) dw2= dw2+denu*ymvvc*ypic*vviq*rwms*(-8*xvq) dw2= dw2+denu*ymvvc*ypis*vviq*rwms*(8*xvc) dw2= dw2+denu*ymvvc*ypis*vvis*rwms*(4./9.*xvc) dw2= dw2+denu*ymvvc*ypi*vviq*rwms*(-8./3.*xvs) dw2= dw2+denu*ymvvc*ypi*vvis*rwms*(-8./9.*xvs) dw2= dw2+denu*ymvvc*vvis*rwms*(4./9.*xv) dw2= dw2+denu*ymvvq*ypiq*vviq*(-4./3.*xvvi) dw2= dw2+denu*ymvvq*ypic*vviq*(4.*xvv) dw2= dw2+denu*ymvvq*ypis*vviq*(-4.*xvq) dw2= dw2+denu*ymvvq*ypis*vvis*(-2./9.*xvq) dw2= dw2+denu*ymvvq*ypi*vviq*(4./3.*xvc) dw2= dw2+denu*ymvvq*ypi*vvis*(4./9.*xvc) dw2= dw2+denu*ymvvq*vvis*(-2./9.*xvs) dw2= dw2+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-1./3.*xvq) dw2= dw2+dend*ymvvs*ypic*vviq*rwms*rwmpgs*(xvc+2./3.*xvq) dw2= dw2+dend*ymvvs*ypis*vviq*rwms*rwmpgs*(-xvs-2*xvc-1./3.*xvq) dw2= dw2+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./18.*xvs) dw2= dw2+dend*ymvvs*ypi*vviq*rwms*rwmpgs*(1./3.*xv+2*xvs+xvc) dw2= dw2+dend*ymvvs*ypi*vvis*rwms*rwmpgs*(1./9.*xv+1./9.*xvs) dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-1./3.*xv) dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*(-1./3.*xv-xvs) dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-1./9.*xv) dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./18.) dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*(-1./9.*xv-1./18.*xvs) dw2= dw2+dend*ymvvc*ypiq*vviq*rwms*(2./3.*xvv) dw2= dw2+dend*ymvvc*ypic*vviq*rwms*(-2*xvq-4./3.*xvv) dw2= dw2+dend*ymvvc*ypis*vviq*rwms*(2*xvc+4.*xvq+2./3.*xvv) dw2= dw2+dend*ymvvc*ypis*vvis*rwms*(1./9.*xvc) dw2= dw2+dend*ymvvc*ypi*vviq*rwms*(-2./3.*xvs-4.*xvc-2*xvq) dw2= dw2+dend*ymvvc*ypi*vvis*rwms*(-2./9.*xvs-2./9.*xvc) dw2= dw2+dend*ymvvc*vviq*rwms*omyp*(2./3.*xvs) dw2= dw2+dend*ymvvc*vviq*rwms*(2./3.*xvs+2*xvc) dw2= dw2+dend*ymvvc*vvis*rwms*omyp*(2./9.*xvs) dw2= dw2+dend*ymvvc*vvis*rwms*omyps*(1./9.*xv) dw2= dw2+dend*ymvvc*vvis*rwms*(2./9.*xvs+1./9.*xvc) dw2= dw2+dend*ymvvq*ypiq*vviq*(-1./3.*xvvi) dw2= dw2+dend*ymvvq*ypic*vviq*(xvv+2./3.*xvvi) dw2= dw2+dend*ymvvq*ypis*vviq*(-xvq-2*xvv-1./3.*xvvi) dw2= dw2+dend*ymvvq*ypis*vvis*(-1./18.*xvq) dw2= dw2+dend*ymvvq*ypi*vviq*(1./3.*xvc+2*xvq+xvv) dw2= dw2+dend*ymvvq*ypi*vvis*(1./9.*xvc+1./9.*xvq) dw2= dw2+dend*ymvvq*vviq*omyp*(-1./3.*xvc) dw2= dw2+dend*ymvvq*vviq*(-1./3.*xvc-xvq) dw2= dw2+dend*ymvvq*vvis*omyp*(-1./9.*xvc) dw2= dw2+dend*ymvvq*vvis*omyps*(-1./18.*xvs) dw2= dw2+dend*ymvvq*vvis*(-1./9.*xvc-1./18.*xvq) * if(ii.eq.3) then wcom= -dw1+2.d0/ym*(1.d0-1.d0/ym)*dw2 wcom= zv*wcom else if(ii.eq.4) then wcom= 0.5d0*em4*omxm*omxm*(dw1-dw2) else if(ii.eq.5) then wcom= em2*((2.d0+0.5d0*omxm)*dw1+omxm*(2.d0/ym- # 3.d0/2.d0)*dw2) endif * endif * www= wcom*ymvvs/pwtm2 * 1 if(iz.eq.0) then ww(ii)= 0.d0 else ww(ii)= www*tjac*tfact endif * enddo * wwt= ww(1)+ww(2)+ww(3)+ww(4)+ww(5) * if(wwt.lt.0.d0) then iz= 0 ifz(7)= ifz(7)+1 ifzzf= ifzzf+1 wwcc20= 0.d0 else wwcc20= wwt 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,em,em2,um,dm * 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,em,em2,um,dm * 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,em,em2,um,dm * 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,em,em2,um,dm * 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,em,em2,um,dm * 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