subroutine wtocalljs implicit real*8(a-h,o-z) * parameter (nout=99) parameter(nmxhep=2000) parameter (itmxmx=100000) * common/wtgnum/itmx common/wttopt/ios,iosf common/hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep), # jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep) * dimension aidhep(6) dimension aphep(5,6) dimension amps(itmxmx),amps1(itmxmx),amps2(itmxmx) * if(iosf.gt.0) then do i=1,6 isthep(i)= 1 enddo else do i=1,4 isthep(i)= 1 enddo endif open(nout,file='mom',status='old') * do i=1,itmx do i=1,4 nevhep= i read(nout,fmt=*) nhep read(nout,fmt=*) amps(i),amps1(i),amps2(i) if(iosf.gt.0) then read(nout,fmt=*) (aidhep(j),j=1,6) else read(nout,fmt=*) (aidhep(j),j=1,4) endif read(nout,fmt=*) (aphep(j,1),j=1,5) read(nout,fmt=*) (aphep(j,2),j=1,5) read(nout,fmt=*) (aphep(j,3),j=1,5) read(nout,fmt=*) (aphep(j,4),j=1,5) if(iosf.gt.0) then read(nout,fmt=*) (aphep(j,5),j=1,5) read(nout,fmt=*) (aphep(j,6),j=1,5) endif do j=1,4 idhep(j)= aidhep(j) enddo if(iosf.gt.0) then idhep(5)= aidhep(5) idhep(6)= aidhep(6) endif do j=1,5 do l=1,4 phep(j,l)= aphep(j,l) enddo if(iosf.gt.0) then phep(j,5)= aphep(j,5) phep(j,6)= aphep(j,6) endif enddo * istrat= 0 atotsq= amps(i) a1sq= amps1(i) a2sq= amps2(i) call lu4frm(atotsq,a1sq,a2sq,istrat,irad,itau,ierr) call lulist(1) enddo close(nout) * return end * *--------------------------------------------------------------------------- * *-- Author : * The following subroutine illustrates how to write a generic * interface between electroweak generators and QCD generators * for LEP 2 applications. * * Author of this routine: T.Sjostrand * * In this particular case, an electroweak generator is supposed to * have produced two fermions, two antifermions and an arbitrary * number of photons. These particles are stored in the HEPEVT * common block. The allowed order is specified by a standard. * In brief, the final fermions should appear in the order * fermion (1) - antifermion (2) - fermion (3) - antifermion (4). * The flavour pairs should be arranged so that, if possible, the * first two could come from a W+ and the second two from a W-; * else each pair should have flavours consistent with a Z0. * The subroutine LU4FRM is supposed to read the configuration, * and call JETSET to do parton showers and fragmentation. * * Since the colour flow need not be unique, three real numbers * should be given when LU4FRM is called: * ATOTSQ = total squared amplitude for the event, irrespective of * colour flow; * A1SQ = squared amplitude for the configuration with 1 + 2 and * 3 + 4 as the two colour singlets; and * A2SQ = squared amplitude for the configuration with 1 + 4 and * 3 + 2 as the two colour singlets. * * The choice of strategy is determined by an integer input: * ISTRAT = 0 : pick configurations according to A1SQ : A2SQ; * = 1 : assign interference to maximize 1 + 2 and 3 + 4; or * = 2 : assign interference to maximize 1 + 4 and 3 + 2. * * Final-state QED radiation may be allowed or inhibited: * IRAD = 0 : no final-state photon radiation. * = 1 : photon radiation inside each final fermion pair. * * tau lepton decay may be handled by QCD generator or not. * ITAU = 0 : taus are considered stable by QCD generator. * = 1 : taus are allowed to decay by QCD generator. * * IERR is an error flag, used both as input and output. * At input, 0 means leave routine in case of error, * nonzero means stop program execution. * At output, 0 means acceptable fermion configuration, * nonzero means treatment aborted for some reason. * It is up to the writer of the main program to pick error strategy, * i.e. to let the program crash or try to fix errors. subroutine lu4frm(atotsq,a1sq,a2sq,istrat,irad,itau,ierr) * real*8 atotsq,a1sq,a2sq * common/lujets/n,k(4000,5),p(4000,5),v(4000,5) common/ludat1/mstu(200),paru(200),mstj(200),parj(200) * dimension ijoin(2) * external rlu * * call luhepc to convert from hepevt to lujets common. * inerr= 0 mstu(28)= 0 call luhepc(2) if(mstu(28).eq.8) inerr= 1 * * loop through entries and pick up all final fermions/antifermions. * i1= 0 i2= 0 i3= 0 i4= 0 do 100 i=1,n if(k(i,1).le.0.or.k(i,1).gt.10) goto 100 kfa= iabs(k(i,2)) if((kfa.ge.1.and.kfa.le.6).or.(kfa.ge.11.and.kfa.le.16)) then if(k(i,2).gt.0) then if(i1.eq.0) then i1= i elseif(i3.eq.0) then i3= i else inerr= 2 call luerrm(6,'(lu4frm:) more than two fermions') endif else if(i2.eq.0) then i2= i elseif(i4.eq.0) then i4= i else inerr= 3 call luerrm(6,'(lu4frm:) more than two antifermions') endif endif endif 100 continue * * check that event is arranged according to conventions. * if(i1.eq.0.or.i2.eq.0.or.i3.eq.0.or.i4.eq.0) then inerr= 4 call luerrm(6,'(lu4frm:) event contains too few fermions') endif if(i2.lt.i1.or.i3.lt.i2.or.i4.lt.i3) then inerr= 5 call luerrm(6,'(lu4frm:) fermions arranged in wrong order') endif * * check which fermion pairs are quarks and which leptons. * if(iabs(k(i1,2)).lt.10.and.iabs(k(i2,2)).lt.10) then iql12= 1 elseif(iabs(k(i1,2)).gt.10.and.iabs(k(i2,2)).gt.10) then iql12= 2 else inerr= 6 call luerrm(6,'(lu4frm:) first fermion pair inconsistent') endif if(iabs(k(i3,2)).lt.10.and.iabs(k(i4,2)).lt.10) then iql34= 1 elseif(iabs(k(i3,2)).gt.10.and.iabs(k(i4,2)).gt.10) then iql34= 2 else inerr= 7 call luerrm(6,'(lu4frm:) second fermion pair inconsistent') endif * * return or stop program in case of problems. * if(inerr.eq.0) then ierr= 0 elseif(ierr.eq.0) then ierr= inerr return else write(6,*) ' error: listing of faulty event follows:' call lulist(2) write(6,*) ' fermions found in lines ',i1,i2,i3,i4 write(6,*) ' error type in event above is ',inerr write(6,*) ' program execution will be stopped now' write(6,*) ' since main program does not correct errors!' stop endif * * decide whether to allow or not photon radiation in showers. * mstj(41)= 2 if(irad.eq.0) mstj(41)= 1 * * decide on colour pairing. * if(iql12.eq.2.and.iql34.eq.2) then npair= 0 elseif(iql12.eq.1.and.iql34.eq.2) then npair= 1 ip1= i1 ip2= i2 elseif(iql12.eq.2.and.iql34.eq.1) then npair= 1 ip1= i3 ip2= i4 else npair= 2 ip1= i1 ip3= i3 r1sq= real(a1sq) r2sq= real(a2sq) delta= real(atotsq-a1sq-a2sq) if(istrat.eq.1) then if(delta.gt.0.) r1sq= r1sq+delta if(delta.lt.0.) r2sq= max(0.,r2sq+delta) elseif(istrat.eq.2) then if(delta.gt.0.) r2sq= r2sq+delta if(delta.lt.0.) r1sq= max(0.,r1sq+delta) endif if(r2sq.lt.rlu(0)*(r1sq+r2sq)) then ip2= i2 ip4= i4 else ip2= i4 ip4= i2 endif endif * * do colour joining and parton showers. * if(npair.ge.1) then ijoin(1)= ip1 ijoin(2)= ip2 call lujoin(2,ijoin) pm12s= (p(ip1,4)+p(ip2,4))**2-(p(ip1,1)+p(ip2,1))**2- # (p(ip1,2)+p(ip2,2))**2-(p(ip1,3)+p(ip2,3))**2 call lushow(ip1,ip2,sqrt(max(0.,pm12s))) endif if(npair.eq.2) then ijoin(1)= ip3 ijoin(2)= ip4 call lujoin(2,ijoin) pm34s= (p(ip3,4)+p(ip4,4))**2-(p(ip3,1)+p(ip4,1))**2- # (p(ip3,2)+p(ip4,2))**2-(p(ip3,3)+p(ip4,3))**2 call lushow(ip3,ip4,sqrt(max(0.,pm34s))) endif * * do fragmentation and decays. possibly except tau decay. * if(itau.eq.0) then if(iabs(k(i1,2)).eq.15) k(i1,1)= 11 if(iabs(k(i2,2)).eq.15) k(i2,1)= 11 if(iabs(k(i3,2)).eq.15) k(i3,1)= 11 if(iabs(k(i4,2)).eq.15) k(i4,1)= 11 endif call luexec if(itau.eq.0) then if(iabs(k(i1,2)).eq.15) k(i1,1)= 1 if(iabs(k(i2,2)).eq.15) k(i2,1)= 1 if(iabs(k(i3,2)).eq.15) k(i3,1)= 1 if(iabs(k(i4,2)).eq.15) k(i4,1)= 1 endif * 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 * *-------------------------------------------------------------------------- * subroutine wtocorrqcd(scal,als,bqm2,cqm2,rbm2,rcm2,rsm2) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * 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= wtoralphas(wm,scal1g,als,nf)/pi alsc4= wtoralphas(wm,cqm,als,nf)/pi alsc42= alsc4*alsc4 als1g2= als1g*als1g alsb= wtoralphas(wm,bqm,als,nf)/pi alsb2= alsb*alsb alsz= wtoralphas(wm,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= wtorrunm(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= wtorrunm(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 * *-----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/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy * dimension b0(5),b1(5),b2(5) * *-----limits for lambda_5 are 1 mev < lambda_5 < 10 gev * 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 * *-------------------------------------------------------------------- * real*8 function wtorals(qcdl,rs,nf) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy * 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) * wtorals= 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 * *-----wtorrunm-------------------------------------------------------- * computes the running quark mass at the pole mass * real*8 function wtorrunm(x1,x2,xacc,qm,als,rm1,rm2,fn) implicit real*8(a-h,o-z) * parameter (jmax=50) * fmid= wtoqcdmass(qm,als,rm1,rm2,fn,x2) f= wtoqcdmass(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 wtorrunm= x1 dx= x2-x1 else wtorrunm= x2 dx= x1-x2 endif do j=1,jmax dx= dx*0.5d0 xmid= wtorrunm+dx fmid= wtoqcdmass(qm,als,rm1,rm2,fn,xmid) if (fmid.le.0.d0) wtorrunm= xmid if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return enddo pause 'too many bisections' end * *-----wtoqcdmass------------------------------------------------------ * real*8 function wtoqcdmass(qm,als,rm1,rm2,fn,x) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * 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)) wtoqcdmass= qm-x*rhs * 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/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * 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 * subroutine wtohadr5(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 * *------------------------------------------------------------- * subroutine wtopself(p2x,pggf) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy * em2= em*em amm2= rmm*rmm tm2= tm*tm tqm2= tqm*tqm * call wtorbff(p2x,em2,em2,b0ee,b1ee,b21ee) call wtorbff(p2x,amm2,amm2,b0mm,b1mm,b21mm) call wtorbff(p2x,tm2,tm2,b0tau,b1tau,b21tau) call wtorbff(p2x,tqm2,tqm2,b0tt,b1tt,b21tt) call wtorbff0(em2,em2,b0ee0,b1ee0,b21ee0) call wtorbff0(amm2,amm2,b0mm0,b1mm0,b21mm0) call wtorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0) call wtorbff0(tqm2,tqm2,b0tt0,b1tt0,b21tt0) * bfe= 2.d0*b21ee-b0ee bfm= 2.d0*b21mm-b0mm bftau= 2.d0*b21tau-b0tau bft= 2.d0*b21tt-b0tt bfe0= 2.d0*b21ee0-b0ee0 bfm0= 2.d0*b21mm0-b0mm0 bftau0= 2.d0*b21tau0-b0tau0 bft0= 2.d0*b21tt0-b0tt0 * sllc= 1.d0+3.d0/4.d0/alphai/pi pggf= 4.d0*(bfe+bfm+bftau)*sllc+16.d0/3.d0*bft- # 4.d0*(bfe0+bfm0+bftau0)*sllc-16.d0/3.d0*bft0 * return end * *------------------------------------------------------------------ * subroutine wtopselfnp(op2x,opggnp) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy * em2= em*em amm2= rmm*rmm tm2= tm*tm tqm2= tqm*tqm sllc= 1.d0+3.d0/4.d0/alphai/pi * p2x= op2x p2z= -zm*zm qzm2= -p2z call wtorbff(p2z,em2,em2,b0eez,b1eez,b21eez) call wtorbff(p2z,amm2,amm2,b0mmz,b1mmz,b21mmz) call wtorbff(p2z,tm2,tm2,b0tauz,b1tauz,b21tauz) call wtorbff(p2z,tqm2,tqm2,b0ttz,b1ttz,b21ttz) call wtorbff0(em2,em2,b0ee0,b1ee0,b21ee0) call wtorbff0(amm2,amm2,b0mm0,b1mm0,b21mm0) call wtorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0) call wtorbff0(tqm2,tqm2,b0tt0,b1tt0,b21tt0) bfez= 2.d0*b21eez-b0eez bfmz= 2.d0*b21mmz-b0mmz bftauz= 2.d0*b21tauz-b0tauz bftz= 2.d0*b21ttz-b0ttz bfe0= 2.d0*b21ee0-b0ee0 bfm0= 2.d0*b21mm0-b0mm0 bftau0= 2.d0*b21tau0-b0tau0 bft0= 2.d0*b21tt0-b0tt0 pggl= 4.d0*(bfez+bfmz+bftauz)*sllc+16.d0/3.d0*bftz pggl0= 4.d0*(bfe0+bfm0+bftau0)*sllc+16.d0/3.d0*bft0 b= 0.00299d0 c= 1.d0 an= -b*log(1.d0+c*qzm2)-0.25d0/alphai/pi*(pggl-pggl0)+ # 1.d0-128.89d0/alphai da= an-0.00165d0 ap2x= abs(p2x) 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 * a= a+da opggnp= 4.d0*pi*alphai*(a+b*log(1.d0+c*ap2x)) * return end * *------------------------------------------------------------------- * subroutine wtorbff0(rm12,rm22,b0,b1,b21) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * 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= ddelta-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 wtorcg(n,yr,yi,omyr,fr) arg(1)= omyr arg(2)= -yi call wtoclnx(arg,cln) f1r= fr(1)-cln(1) f2r= 2.d0*fr(2)-cln(1) f3r= 3.d0*fr(3)-cln(1) fact= ddelta-log(arm22) b0= fact-f1r b1= 0.5d0*(-fact+f2r) b21= 1.d0/3.d0*(fact-f3r) return endif end * *---------------------------------------------------------- * subroutine wtorbff(p2r,rm12,rm22,rb0,rb1,rb21) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * dimension cmp2(2) dimension clnmp2(2) dimension gfpr(3),gfmr(3) * if(rm12.eq.0.d0.and.rm22.eq.0.d0) then cmp2(1)= p2r cmp2(2)= -eps call wtoclnx(cmp2,clnmp2) rb0= ddelta+2.d0-clnmp2(1) rb1= -0.5d0*rb0 rb21= 1.d0/18.d0+1.d0/3.d0*rb0 else cmp2(1)= -p2r cmp2(2)= -eps call wtoroots(p2r,rm12,rm22,rpr,rpi,rmr,rmi,omrpr,omrmr) call wtoclnx(cmp2,clnmp2) n= 3 call wtorcg(n,rpr,rpi,omrpr,gfpr) call wtorcg(n,rmr,rmi,omrmr,gfmr) auxdel= ddelta-clnmp2(1) rb0= auxdel-gfpr(1)-gfmr(1) rb1= -auxdel/2.d0+gfpr(2)+gfmr(2) rb21= auxdel/3.d0-gfpr(3)-gfmr(3) endif * return end * *---------------------------------------------------------- * subroutine wtoclnx(arg,res) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * 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 wtoclnomx(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 wtoclnx(omarg,ares) endif do i= 1,2 res(i)= ares(i) enddo * return end * *---------------------------------------------------------- * subroutine wtocqlnx(arg,res) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * 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.q0) then if (arg(2).gt.0.q0) then teta= qpi/2.q0 else teta= -qpi/2.q0 endif res(2)= teta return else if (arg(2).eq.0.q0) then if (arg(1).gt.0.q0) then teta= 0.q0 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.q0) then res(2)= si*teta else res(2)= si*(qpi-teta) endif return endif end * *------------------------------------------------------------- * subroutine wtocqlnomx(arg,omarg,res) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * 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.q-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.q0 ares(2)= sn(10)/10.q0 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.q0) then if(yi.gt.0.q0) then teta= qpi/2.q0 else teta= -qpi/2.q0 endif ares(2)= teta else if(yi.eq.0.q0) then if(yr.gt.0.q0) then teta= 0.q0 else teta= qpi 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.q0) then ares(2)= si*teta else ares(2)= si*(qpi-teta) endif endif endif do i= 1,2 res(i)= ares(i) enddo * return end * *-------------------------------------------------------------------- * subroutine wtorcg(n,zr,zi,omzr,gfr) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * 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.2d+0,0.666666666666667d-1, # -0.952380952380952d-2,-0.396825396825397d-3, # 0.317460317460317d-3,-0.132275132275132d-4, # -0.962000962000962d-5,0.105218855218855d-5, # 0.266488361726450d-6,-0.488745528428000d-7, # -0.675397500794000d-8,0.190720263471000d-8, # 0.153663007690000d-9,-0.679697905790000d-10, # -0.293683556000000d-11,0.228836696000000d-11/ * 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 wtoclnx(omz,clnomz) * |z| < 4 if (zmod.lt.4.d0) then oz(1)= -z(1) oz(2)= -z(2) call wtoclnx(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*zi)/zmod2 aux(2)= zi/zmod2 call wtoclnx(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 * * *-----rspence----------------------------------------------------- * gives the real part of li_2(z) in double precision. accuracy * is about 8 digits * real*8 function wtorspence(xr,xi) implicit real*8(a,b,d-h,o-z) implicit complex*16 (c) * data pi/3.141592653589793238462643d0/ * 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 * wtorspence= real(cli2) * return end * *---------------------------------------------------------- * real*8 function wtorfsr(x,r) implicit real*8(a,b,d-h,o-z) implicit complex*16(c) * x2= x*x r2= r*r * a1r= 1.d0 a1i= 1.d0/r a2r= 1.d0/(1.d0+r2*x2) a2i= r*x/(1.d0+r2*x2) a3r= 1.d0/(1.d0+r2) a3i= r/(1.d0+r2) a4r= (r2*x+1.d0)/(1.d0+r2) a4i= r*(1.d0-x)/(1.d0+r2) * ca1= cmplx(a1r,a1i) ca2= cmplx(a2r,a2i) * cpar= -log(ca1)*log(ca2) * wtorfsr= 2.d0*(real(cpar)+wtorspence(a3r,a3i)- # wtorspence(a4r,a4i)) * return end * *-------------------------------------------------------------- * subroutine funct1(n,xc,fc) implicit real*8(a-h,o-z) character*1 omssm character*4 otype parameter(npos=512) * common/wtkr/ireg common/wtmssmo/omssm common/wtochannel/otype common/wthx0/xshmx0(npos) common/wtcx0/xscmx0(npos) * dimension xc(n) * if(otype.eq.'cc11'.or.otype.eq.'cc20') then fc= -wtoxsc(n,xc)/xscmx0(ireg) else if(otype.eq.'cc12') then fc= -wtoxssc(n,xc)/xshmx0(ireg) else if(otype.eq.'nc33') then fc= -wtoxsnh32(n,xc)/xshmx0(ireg) else if(otype.eq.'nc32') then fc= -wtoxsn32(n,xc)/xshmx0(ireg) else if(otype.eq.'nc50') then fc= -wtoxsnh49(n,xc)/xshmx0(ireg) else if(otype.eq.'nc25') then fc= -wtoxsnh24(n,xc)/xshmx0(ireg) else if(otype.eq.'nc21') then fc= -wtoxsnh19(n,xc)/xshmx0(ireg) else if(otype.eq.'nc19') then fc= -wtoxsn(n,xc)/xshmx0(ireg) else if(otype.eq.'nc24') then fc= -wtoxsn(n,xc)/xshmx0(ireg) else if(otype.eq.'nc26') then fc= -wtoxsnh26(n,xc)/xshmx0(ireg) else if(otype.eq.'nc48') then fc= -wtoxsn(n,xc)/xshmx0(ireg) else if(otype.eq.'nc64') then fc= -wtoxsn64(n,xc)/xshmx0(ireg) else if(otype.eq.'nc68') then if(omssm.eq.'n') then fc= -wtoxsnh64(n,xc)/xshmx0(ireg) else if(omssm.eq.'y') then fc= -wtoxsnsh64(n,xc)/xshmx0(ireg) endif endif * return end * *------------------------------------------------------ * integer*4 function iwtopos(n,x) implicit real*8(a-h,o-z) character*1,osm * common/wtsmod/osm common/wtspl/prx0(9) common/wtspln/prx0n(9) * dimension x(n),i(9),apx(9) * do j=1,9 if(osm.eq.'n') then apx(j)= prx0(j) else if(osm.eq.'g') then apx(j)= prx0n(j) endif enddo * do j=1,n if(x(j).lt.apx(j)) then i(j)= 1 else i(j)= 2 endif enddo * iwtopos= 2*(2*(2*(2*(2*(2*(2*(2*(i(1)-1)+i(2)-1)+ # i(3)-1)+i(4)-1)+i(5)-1)+i(6)-1)+i(7)-1)+i(8)-1)+i(9) * return end * *------------------------------------------------------------------------- * subroutine wtomxregion(ipos,bl,bu) implicit real*8(a-h,o-z) * common/wtspl/prx0(9) * dimension bl(9),bu(9) * data eps/1.d-6/ * do i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 do i5=1,2 do i6=1,2 do i7=1,2 do i8=1,2 do i9=1,2 ip= 2*(2*(2*(2*(2*(2*(2*(2*(i1-1)+i2-1)+ # i3-1)+i4-1)+i5-1)+i6-1)+i7-1)+i8-1)+i9 if(ip.eq.ipos) then if(i1.eq.1) then bl(1)= eps bu(1)= prx0(1) else bl(1)= prx0(1) bu(1)= 1.d0-eps endif if(i2.eq.1) then bl(2)= eps bu(2)= prx0(2) else bl(2)= prx0(2) bu(2)= 1.d0-eps endif if(i3.eq.1) then bl(3)= eps bu(3)= prx0(3) else bl(3)= prx0(3) bu(3)= 1.d0-eps endif if(i4.eq.1) then bl(4)= eps bu(4)= prx0(4) else bl(4)= prx0(4) bu(4)= 1.d0-eps endif if(i5.eq.1) then bl(5)= eps bu(5)= prx0(5) else bl(5)= prx0(5) bu(5)= 1.d0-eps endif if(i6.eq.1) then bl(6)= eps bu(6)= prx0(6) else bl(6)= prx0(6) bu(6)= 1.d0-eps endif if(i7.eq.1) then bl(7)= eps bu(7)= prx0(7) else bl(7)= prx0(7) bu(7)= 1.d0-eps endif if(i8.eq.1) then bl(8)= eps bu(8)= prx0(8) else bl(8)= prx0(8) bu(8)= 1.d0-eps endif if(i9.eq.1) then bl(9)= eps bu(9)= prx0(9) else bl(9)= prx0(9) bu(9)= 1.d0-eps endif endif enddo enddo enddo enddo enddo enddo enddo enddo enddo * return end * *---------------------------------------------------------------- * subroutine wtovolume(npos,nint,xsmx,oxsmx,gbl,gbu,vol,volt) implicit real*8(a-h,o-z) character*1 ord parameter(nposm=512) * common/wtspln/prx0n(9) * dimension vol(npos,nint),volt(nint) dimension avol(nposm),ir(nposm) dimension bl(nposm,9),bu(nposm,9) dimension xsmx(npos,nint),oxsmx(npos,nint),gbl(npos,9,nint), # gbu(npos,9,nint) * do i=1,npos do i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 do i5=1,2 do i6=1,2 do i7=1,2 do i8=1,2 do i9=1,2 ip= 2*(2*(2*(2*(2*(2*(2*(2*(i1-1)+i2-1)+ # i3-1)+i4-1)+i5-1)+i6-1)+i7-1)+i8-1)+i9 if(ip.eq.i) then if(i1.eq.1) then bl(i,1)= 0.d0 bu(i,1)= prx0n(1) else bl(i,1)= prx0n(1) bu(i,1)= 1.d0 endif if(i2.eq.1) then bl(i,2)= 0.d0 bu(i,2)= prx0n(2) else bl(i,2)= prx0n(2) bu(i,2)= 1.d0 endif if(i3.eq.1) then bl(i,3)= 0.d0 bu(i,3)= prx0n(3) else bl(i,3)= prx0n(3) bu(i,3)= 1.d0 endif if(i4.eq.1) then bl(i,4)= 0.d0 bu(i,4)= prx0n(4) else bl(i,4)= prx0n(4) bu(i,4)= 1.d0 endif if(i5.eq.1) then bl(i,5)= 0.d0 bu(i,5)= prx0n(5) else bl(i,5)= prx0n(5) bu(i,5)= 1.d0 endif if(i6.eq.1) then bl(i,6)= 0.d0 bu(i,6)= prx0n(6) else bl(i,6)= prx0n(6) bu(i,6)= 1.d0 endif if(i7.eq.1) then bl(i,7)= 0.d0 bu(i,7)= prx0n(7) else bl(i,7)= prx0n(7) bu(i,7)= 1.d0 endif if(i8.eq.1) then bl(i,8)= 0.d0 bu(i,8)= prx0n(8) else bl(i,8)= prx0n(8) bu(i,8)= 1.d0 endif if(i9.eq.1) then bl(i,9)= 0.d0 bu(i,9)= prx0n(9) else bl(i,9)= prx0n(9) bu(i,9)= 1.d0 endif endif enddo enddo enddo enddo enddo enddo enddo enddo enddo enddo * do ij=1,nint do i=1,npos vol(i,ij)= 1.d0 do j=1,9 vol(i,ij)= vol(i,ij)*(bu(i,j)-bl(i,j)) enddo vol(i,ij)= vol(i,ij)*xsmx(i,ij) enddo enddo do ij=1,nint volt(ij)= 0.d0 do i=1,npos volt(ij)= volt(ij)+vol(i,ij) enddo enddo * m1= 1 m2= npos ord= 'd' do ij=1,nint do i=1,npos avol(i)= vol(i,ij) enddo ifail= 0 call m01daf(avol,m1,m2,ord,ir,ifail) do i=1,npos vol(ir(i),ij)= avol(i) oxsmx(ir(i),ij)= xsmx(i,ij) do j=1,9 gbl(ir(i),j,ij)= bl(i,j) gbu(ir(i),j,ij)= bu(i,j) enddo enddo enddo * return end * *------------------------------------------------------------------ * function ipack(shat,sm,sp,su,sd,sf) implicit real*8 (a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * zl= (zm-25.d0)*(zm-25.d0)/shat zu= (zm+25.d0)*(zm+25.d0)/shat bc= 2.5d3/shat sg= 1.d0-sm-sp-su-sd-sf * if(zu.gt.sm.and.sm.gt.zl.and.sp.gt.bc) then ipack= 1 else if(zu.gt.sp.and.sp.gt.zl.and.sm.gt.bc) then ipack= 1 else if(zu.gt.sf.and.sf.gt.zl.and.sg.gt.bc) then ipack= 1 else if(zu.gt.sg.and.sg.gt.zl.and.sf.gt.bc) then ipack= 1 else if(zu.gt.su.and.su.gt.zl.and.sd.gt.bc) then ipack= 1 else if(zu.gt.sd.and.sd.gt.zl.and.su.gt.bc) then ipack= 1 else ipack= 0 endif * return end * *------------------------------------------------------------------ * function ievs(a,b,s1,s2,s3,s4,s5,s6) implicit real*8 (a-h,o-z) * common/wthiggs/hm * c= hm-10.d0 d= hm+10.d0 * if((b.gt.s1.and.s1.gt.a).and.(s2.lt.c.or.s2.gt.d).and. # (s3.lt.c.or.s3.gt.d).and.(s4.lt.c.or.s4.gt.d).and. # (s5.lt.c.or.s5.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s2.and.s2.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s3.lt.c.or.s3.gt.d).and.(s4.lt.c.or.s4.gt.d).and. # (s5.lt.c.or.s5.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s3.and.s3.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s2.lt.c.or.s2.gt.d).and.(s4.lt.c.or.s4.gt.d).and. # (s5.lt.c.or.s5.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s4.and.s4.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s2.lt.c.or.s2.gt.d).and.(s3.lt.c.or.s3.gt.d).and. # (s5.lt.c.or.s5.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s5.and.s5.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s2.lt.c.or.s2.gt.d).and.(s3.lt.c.or.s3.gt.d).and. # (s4.lt.c.or.s4.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s6.and.s6.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s2.lt.c.or.s2.gt.d).and.(s3.lt.c.or.s3.gt.d).and. # (s4.lt.c.or.s4.gt.d).and.(s5.lt.c.or.s5.gt.d)) then ievs= 1 else ievs= 0 endif * return end * *------------------------------------------------------------------- * subroutine wtopole(p2r,p2i,p3q,fz,fw) implicit real*8(a-h,o-z) character*1 rio * common/wtai/rio common/wtparam/eps,ddelta common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * dimension bfl(2),bfh(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) * call wtosbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l) call wtosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h) call wtosbff(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) * if(rio.eq.'a') then do i=1,2 p3q(i)= 10.d0*bfl(i)+2.d0*bfh(i) enddo else if(rio.eq.'i') then p3q(1)= 0.d0 p3q(2)= 10.d0*bfl(2)+2.d0*bfh(2) endif * if(rio.eq.'a') then fz(1)= -0.5d0*(p2r*bfsr-p2i*bfsi)-3.d0/2.d0*tqm2*b0h(1) else if(rio.eq.'i') then fz(1)= 0.d0 endif 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) * if(rio.eq.'a') then fw(1)= 2.d0*(p2r*auxr-p2i*auxi)-3.d0*tqm2*(b1c(1)+b0c(1)) else if(rio.eq.'i') then fw(1)= 0.d0 endif fw(2)= 2.d0*(p2r*auxi+p2i*auxr)-3.d0*tqm2*(b1c(2)+b0c(2)) * return end * *------------------------------------------------------------------------ * subroutine wtobff(p2,rm12,rm22,b0,b1,b21) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * dimension cmp2(2) dimension clnmp2(2) dimension b0(2),b1(2),b21(2) dimension gfpr(3),gfmr(3),gfpi(3),gfmi(3) * cmp2(1)= -p2 cmp2(2)= -qeps * call wtoqroots(p2,rm12,rm22,rpr,rpi,rmr,rmi,omrpr,omrmr) * call wtocqlnx(cmp2,clnmp2) n= 3 call wtocg(n,rpr,rpi,omrpr,gfpr,gfpi) call wtocg(n,rmr,rmi,omrmr,gfmr,gfmi) * auxdel= qdelta-clnmp2(1) b0(1)= auxdel-gfpr(1)-gfmr(1) b0(2)= -clnmp2(2)-gfpi(1)-gfmi(1) b1(1)= -auxdel/2.q0+gfpr(2)+gfmr(2) b1(2)= +clnmp2(2)/2.q0+gfpi(2)+gfmi(2) b21(1)= auxdel/3.q0-gfpr(3)-gfmr(3) b21(2)= -clnmp2(2)/3.q0-gfpi(3)-gfmi(3) * return end * *------------------------------------------------------------------------- * subroutine wtosbff(p2r,p2i,rm12,rm22,b0,b1,b21) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * dimension cmp2(2) dimension clnmp2(2) dimension b0(2),b1(2),b21(2) * p2m= p2r*p2r+p2i*p2i * if(rm12.eq.0.d0.and.rm22.eq.0.d0) then if(p2i.eq.0) then cmp2(1)= p2r cmp2(2)= -eps call wtoclnx(cmp2,clnmp2) b0(1)= ddelta+2.d0-clnmp2(1) b0(2)= -clnmp2(2) else cmp2(1)= p2r cmp2(2)= p2i call wtoclnx(cmp2,clnmp2) b0(1)= ddelta+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*(-ddelta-1.d0+log(rm12)) if(p2i.eq.0) then cmp2(1)= 1.d0+p2r/rm12 cmp2(2)= -eps/rm12 call wtoclnx(cmp2,clnmp2) b0(1)= ddelta+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 wtoclnx(cmp2,clnmp2) xr= p2r/p2m xi= -p2i/p2m fctr= 1.d0+rm12*xr fcti= rm12*xi b0(1)= ddelta+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*(-ddelta-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 wtoclnx(cmp2,clnmp2) b0(1)= ddelta+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 wtoclnx(cmp2,clnmp2) b0(1)= ddelta+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*(-ddelta-1.d0+log(rm12)) a02= rm22*(-ddelta-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 wtoclnx(cmp2,clnmp2) crr= (br*clnmp2(1)-bi*clnmp2(2))/p2r cri= (br*clnmp2(2)+bi*clnmp2(1))/p2r b0(1)= ddelta-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 wtoclnx(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)= ddelta+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 wtosbffs(p2r,p2i,rm12,rm22,b0) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * dimension b0(2) dimension cmp2(2) dimension clnmp2(2) * 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 wtoclnx(cmp2,clnmp2) b0(1)= ddelta+2.d0-clnmp2(1) b0(2)= -clnmp2(2) else cmp2(1)= p2r cmp2(2)= p2i call wtoclnx(cmp2,clnmp2) b0(1)= ddelta+2.d0-clnmp2(1) b0(2)= -clnmp2(2)+2.d0*pi endif else if(rm22.eq.0.d0.and.rm12.ne.0.d0) then a0= rm12*(-ddelta-1.d0+log(rm12)) if(p2i.eq.0) then cmp2(1)= 1.d0+p2r/rm12 cmp2(2)= -eps/rm12 call wtoclnx(cmp2,clnmp2) b0(1)= ddelta+2.d0-log(rm12)-(1.d0+rm12/p2r)*clnmp2(1) b0(2)= -(1.d0+rm12/p2r)*clnmp2(2) else cmp2(1)= 1.d0+p2r/rm12 cmp2(2)= p2i/rm12 call wtoclnx(cmp2,clnmp2) xr= p2r/p2m xi= -p2i/p2m fctr= 1.d0+rm12*xr fcti= rm12*xi b0(1)= ddelta+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 endif else if(rm12.ne.0.d0.and.rm12.eq.rm22) then a0= rm12*(-ddelta-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 wtoclnx(cmp2,clnmp2) b0(1)= ddelta+2.d0-log(rm12)-(br*clnmp2(1)- # bi*clnmp2(2)) b0(2)= -(br*clnmp2(2)+bi*clnmp2(1)) 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 wtoclnx(cmp2,clnmp2) b0(1)= ddelta+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 endif else a01= rm12*(-ddelta-1.d0+log(rm12)) a02= rm22*(-ddelta-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 wtoclnx(cmp2,clnmp2) crr= (br*clnmp2(1)-bi*clnmp2(2))/p2r cri= (br*clnmp2(2)+bi*clnmp2(1))/p2r b0(1)= ddelta-crr-0.5d0*log(rm22*rm12)+0.5d0*(rm22-rm12)/ # p2r*log(rm22/rm12)+2.d0 b0(2)= -cri 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 wtoclnx(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)= ddelta+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 endif endif * return end * *--------------------------------------------------------------------- * subroutine wtocff(jflag,op12,op22,os,orm12,orm22,orm32, # oc0,oc1,oc2,oc3) implicit real*16(a-h,p-z) implicit real*8(o) character*1 rio * common/wtai/rio common/wtqparam/qpi,qpis,qeps,qdelta * dimension co(2) dimension tmp3_56(2,2),r3_56(2,2) dimension b0_12(2),b1_12(2),b21_12(2) dimension b0_13(2),b1_13(2),b21_13(2) dimension b0_23(2),b1_23(2),b21_23(2) dimension b22_12(2),b22_13(2),b22_23(2) dimension oc0(2),oc1(2,2),oc2(2,4),oc3(2,6) dimension r3_34(2,2),r3_13(2,2),r3_42(2,2) dimension tmp3_34(2,2),tmp3_13(2,2),tmp3_42(2,2) dimension r2_32(2,2),tmp2_32(2,2),tmp0(2),tmp24(2) dimension xmi(2,2),r1(2,2),tmp1(2,2),r2_13(2,2),tmp2_13(2,2) * qpis= qpi*qpi p12= op12*1.d15*1.q-15 p22= op22*1.d15*1.q-15 s= os*1.d15*1.q-15 rm12= orm12*1.d15*1.q-15 rm22= orm22*1.d15*1.q-15 rm32= orm32*1.d15*1.q-15 co(1)= 1.q0 co(2)= 0.q0 if(jflag.eq.0) then rsm= sqrt(p12/s) rsp= sqrt(p22/s) testc= rsm+rsp if(testc.gt.1.q0) then jflag= -1 return else jflag= 1 endif endif * p1dp2= 0.5q0*(s-p12-p22) den= -s*s/4.q0*(1.q0+((p12-p22)*(p12-p22)/s- # 2.q0*(p12+p22))/s) dm12= rm12-rm22 dm13= rm12-rm32 dm23= rm22-rm32 ra= -p22 rb= -p12 rc= -2.q0*p1dp2 rd= dm23+p22 re= dm12+p12+2.q0*p1dp2 raa= -p22+dm23 rbb= p12+dm12 rcc= -2.q0*p1dp2-p22+dm23 rdd= -p12+dm12 discp= +2.q0*sqrt(abs(den)) discm= -discp if(den.lt.0.q0) then arga= (p22-2.q0*(p12+s))*p22+(s-p12)*(s-p12) rta= sqrt(arga) s1a= abs(s/p12) s2a= abs(s/p22) if(s1a.lt.1.q-10.or.s2a.lt.1.q-10) then bal= p12-p22+s if(bal.gt.0.q0) then qal= -0.5q0*(bal+rta) bep= -s/qal bem= -qal/p12 else qal= -0.5q0*(bal-rta) bep= -qal/p12 bem= -s/qal endif ap= 1.q0-bem am= 1.q0-bep qomap= bem qomam= bep else bal= p12+p22-s if(bal.gt.0.q0) then qal= -0.5q0*(bal+rta) ap= -p22/qal am= -qal/p12 else qal= -0.5q0*(bal-rta) ap= -qal/p12 am= -p22/qal endif qomap= 1.q0-ap qomam= 1.q0-am endif if(jflag.eq.1) then alp= ap qomalp= qomap disc= discp else if(jflag.eq.2) then alp= am qomalp= qomam disc= discm else print*,'jflag > 2' stop endif alpi= 0.q0 y0r= -(rd+re*alp)/disc y0i= 0.q0 y1r= -(raa+rbb*alp)/disc y1i= y0i qomy1r= (rcc+rdd*alp)/disc y2r= y0r/qomalp y2i= y0i/qomalp qomy2r= qomy1r/qomalp y3r= -y0r/alp y3i= -y0i/alp qomy3r= y1r/alp else if(den.gt.0.q0) then alp= (p12+p22-s)/2.q0/p12 qomalp= (p12-p22+s)/2.q0/p12 calpi= sqrt(den)/p12 if(jflag.eq.1) then alpi= -calpi disc= discp else if(jflag.eq.2) then alpi= calpi disc= discm else print*,'jflag > 2' stop endif y0r= -re*alpi/disc y0i= (rd+re*alp)/disc y1r= -rbb*alpi/disc y1i= (raa+rbb*alp)/disc qomy1r= rdd*alpi/disc qomalpm2= qomalp*qomalp+alpi*alpi y2r= (y0r*qomalp-y0i*alpi)/qomalpm2 y2i= (y0r*alpi+y0i*qomalp)/qomalpm2 qomy2r= (qomy1r*qomalp+y1i*alpi)/qomalpm2 alpm2= alp*alp+alpi*alpi y3r= -(y0r*alp+y0i*alpi)/alpm2 y3i= -(-y0r*alpi+y0i*alp)/alpm2 qomy3r= (y1r*alp+y1i*alpi)/alpm2 endif * call wtoqroots(p12,rm22,rm12,rp1r,rp1i,rm1r,rm1i, # qomrp1r,qomrm1r) call wtoqroots(s,rm32,rm12,rp2r,rp2i,rm2r,rm2i, # qomrp2r,qomrm2r) call wtoqroots(p22,rm32,rm22,rp3r,rp3i,rm3r,rm3i, # qomrp3r,qomrm3r) * a1= -p12 a2= -s a3= -p22 eb1r= -1.q0-dm12/p12 eb1i= 0.q0 ec1i= qeps/p12 eb2r= -1.q0-dm13/s eb2i= 0.q0 ec2i= qeps/s eb3r= -1.q0-dm23/p22 eb3i= 0.q0 ec3i= qeps/p22 call wtos2(y1r,y1i,qomy1r,a1,eb1r,eb1i,ec1i,rp1r,rp1i, # rm1r,rm1i,qomrp1r,qomrm1r,s21r,s21i) call wtos2(y2r,y2i,qomy2r,a2,eb2r,eb2i,ec2i,rp2r,rp2i, # rm2r,rm2i,qomrp2r,qomrm2r,s22r,s22i) call wtos2(y3r,y3i,qomy3r,a3,eb3r,eb3i,ec3i,rp3r,rp3i, # rm3r,rm3i,qomrp3r,qomrm3r,s23r,s23i) * *-----scalar 3-point function c0 * if(den.lt.0.q0) then tmp0(1)= (s21r-s22r+s23r)/disc tmp0(2)= (s21i-s22i+s23i)/disc else if(den.gt.0.q0) then tmp0(1)= (s21i-s22i+s23i)/disc tmp0(2)= -(s21r-s22r+s23r)/disc endif * f1= dm12-p12 f2= dm23-s+p12 * call wtobff(p12,rm12,rm22,b0_12,b1_12,b21_12) call wtobff(s,rm12,rm32,b0_13,b1_13,b21_13) call wtobff(p22,rm22,rm32,b0_23,b1_23,b21_23) * do i=1,2 b22_12(i)= -p12*b21_12(i)+0.5q0*(rm12-rm22-p12)* # b1_12(i)+0.5q0*co(i)*rm22*(-qdelta-1.q0+ # log(rm22)) b22_13(i)= -s*b21_13(i)+0.5q0*(rm12-rm32-s)* # b1_13(i)+0.5q0*co(i)*rm32*(-qdelta-1.q0+ # log(rm32)) b22_23(i)= -p22*b21_23(i)+0.5q0*(rm22-rm32-p22)* # b1_23(i)+0.5q0*co(i)*rm32*(-qdelta-1.q0+ # log(rm32)) enddo * xmi(1,1)= p22/den xmi(1,2)= -p1dp2/den xmi(2,1)= xmi(1,2) xmi(2,2)= p12/den * *-----c1 form factors * do i= 1,2 r1(i,1)= 0.5q0*(f1*tmp0(i)+b0_13(i)-b0_23(i)) r1(i,2)= 0.5q0*(f2*tmp0(i)+b0_12(i)-b0_13(i)) enddo do j= 1,2 do i= 1,2 tmp1(j,i)= xmi(i,1)*r1(j,1)+xmi(i,2)*r1(j,2) enddo enddo * *-----c2 form factors * tmp24(1)= 0.25q0-0.5q0*rm12*tmp0(1)+ # 0.25q0*(b0_23(1)-f1*tmp1(1,1)-f2*tmp1(1,2)) tmp24(2)= -0.5q0*rm12*tmp0(2)+ # 0.25q0*(b0_23(2)-f1*tmp1(2,1)-f2*tmp1(2,2)) do i= 1,2 r2_13(i,1)= 0.5q0*(f1*tmp1(i,1)+b1_13(i)+b0_23(i))-tmp24(i) r2_13(i,2)= 0.5q0*(f2*tmp1(i,1)+b1_12(i)-b1_13(i)) enddo do j= 1,2 do i= 1,2 tmp2_13(j,i)= xmi(i,1)*r2_13(j,1)+xmi(i,2)*r2_13(j,2) enddo enddo do i= 1,2 r2_32(i,1)= 0.5q0*(f1*tmp1(i,2)+b1_13(i)-b1_23(i)) r2_32(i,2)= 0.5q0*(f2*tmp1(i,2)-b1_13(i))-tmp24(i) enddo do j= 1,2 do i= 1,2 tmp2_32(j,i)= xmi(i,1)*r2_32(j,1)+xmi(i,2)*r2_32(j,2) enddo enddo * *-----c3 form factors * do i=1,2 r3_56(i,1)= 0.5q0*(f1*tmp24(i)+b22_13(i)-b22_23(i)) r3_56(i,2)= 0.5q0*(f2*tmp24(i)+b22_12(i)-b22_13(i)) enddo do j= 1,2 do i= 1,2 tmp3_56(j,i)= xmi(i,1)*r3_56(j,1)+xmi(i,2)*r3_56(j,2) enddo enddo * do i=1,2 r3_34(i,1)= 0.5q0*(f1*tmp2_13(i,2)+b21_13(i)+b1_23(i))- # tmp3_56(i,2) r3_34(i,2)= 0.5q0*(f2*tmp2_13(i,2)-b21_13(i))-tmp3_56(i,1) enddo do j= 1,2 do i= 1,2 tmp3_34(j,i)= xmi(i,1)*r3_34(j,1)+xmi(i,2)*r3_34(j,2) enddo enddo * do i=1,2 r3_13(i,1)= 0.5q0*(f1*tmp2_13(i,1)+b21_13(i)-b0_23(i))- # 2.q0*tmp3_56(i,1) r3_13(i,2)= 0.5q0*(f2*tmp2_13(i,1)+b21_12(i)-b21_13(i)) enddo do j= 1,2 do i= 1,2 tmp3_13(j,i)= xmi(i,1)*r3_13(j,1)+xmi(i,2)*r3_13(j,2) enddo enddo * do i=1,2 r3_42(i,1)= 0.5q0*(f1*tmp2_32(i,2)+b21_13(i)-b21_23(i)) r3_42(i,2)= 0.5q0*(f2*tmp2_32(i,2)-b21_13(i))-2.q0* # tmp3_56(i,2) enddo do j= 1,2 do i= 1,2 tmp3_42(j,i)= xmi(i,1)*r3_42(j,1)+xmi(i,2)*r3_42(j,2) enddo enddo * if(rio.eq.'a') then do i= 1,2 oc0(i)= tmp0(i) oc1(i,1)= tmp1(i,1) oc1(i,2)= tmp1(i,2) oc2(i,1)= tmp2_13(i,1) oc2(i,2)= tmp2_32(i,2) oc2(i,3)= tmp2_13(i,2) oc2(i,4)= tmp24(i) oc3(i,5)= tmp3_56(i,1) oc3(i,6)= tmp3_56(i,2) oc3(i,3)= tmp3_34(i,1) oc3(i,4)= tmp3_34(i,2) oc3(i,1)= tmp3_13(i,1) oc3(i,2)= tmp3_42(i,2) enddo else if(rio.eq.'i') then oc0(1)= 0.d0 oc1(1,1)= 0.d0 oc1(1,2)= 0.d0 oc2(1,1)= 0.d0 oc2(1,2)= 0.d0 oc2(1,3)= 0.d0 oc2(1,4)= 0.d0 oc3(1,5)= 0.d0 oc3(1,6)= 0.d0 oc3(1,3)= 0.d0 oc3(1,4)= 0.d0 oc3(1,1)= 0.d0 oc3(1,2)= 0.d0 oc0(2)= tmp0(2) oc1(2,1)= tmp1(2,1) oc1(2,2)= tmp1(2,2) oc2(2,1)= tmp2_13(2,1) oc2(2,2)= tmp2_32(2,2) oc2(2,3)= tmp2_13(2,2) oc2(2,4)= tmp24(2) oc3(2,5)= tmp3_56(2,1) oc3(2,6)= tmp3_56(2,2) oc3(2,3)= tmp3_34(2,1) oc3(2,4)= tmp3_34(2,2) oc3(2,1)= tmp3_13(2,1) oc3(2,2)= tmp3_42(2,2) endif * return end * *------------------------------------------------------------------- * subroutine wtoroots(p2,rm12,rm22,rpr,rpi,rmr,rmi,omrpr,omrmr) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * if (rm12.eq.rm22) then disc2= (p2+4.d0*rm12)*p2 else disc2= (rm22+2.d0*(p2-rm12))*rm22+(p2+rm12)*(p2+rm12) endif if (disc2.gt.0.d0) then disc= sqrt(disc2) rpi= +eps/disc rmi= -eps/disc r1= abs(rm22/p2) r2= abs(rm22/rm12) * * one of the roots very near to 1 * if (r1.lt.1.d-10.or.r2.lt.1.d-10) then a= -p2 b= p2+rm12-rm22 c= rm22 if (b.gt.0.d0) then q= -0.5d0*(b+disc) rmr= 1.d0-c/q rpr= 1.d0-q/a omrpr= q/a omrmr= c/q else q= -0.5d0*(b-disc) rmr= 1.d0-q/a rpr= 1.d0-c/q omrpr= c/q omrmr= q/a endif else a= -p2 b= p2+rm22-rm12 c= rm12 if (b.gt.0.d0) then q= -0.5d0*(b+disc) rpr= c/q rmr= q/a omrpr= 1.d0-rpr omrmr= 1.d0-rmr else q= -0.5d0*(b-disc) rpr= q/a rmr= c/q omrpr= 1.d0-rpr omrmr= 1.d0-rmr endif endif else disc= sqrt(-disc2) a= -p2 b= p2+rm22-rm12 rpr= -b/2.d0/a rmr= -b/2.d0/a omrpr= 1.d0-rpr omrmr= 1.d0-rmr rpi= +disc/2.d0/a rmi= -disc/2.d0/a endif return end * *------------------------------------------------------------------- * subroutine wtoqroots(p2,rm12,rm22,rpr,rpi,rmr,rmi,omrpr,omrmr) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * if (rm12.eq.rm22) then disc2= (p2+4.q0*rm12)*p2 else disc2= (rm22+2.q0*(p2-rm12))*rm22+(p2+rm12)*(p2+rm12) endif if (disc2.gt.0.q0) then disc= sqrt(disc2) rpi= +qeps/disc rmi= -qeps/disc r1= abs(rm22/p2) r2= abs(rm22/rm12) * * one of the roots very near to 1 * if (r1.lt.1.q-10.or.r2.lt.1.q-10) then a= -p2 b= p2+rm12-rm22 c= rm22 if (b.gt.0.q0) then q= -0.5q0*(b+disc) rmr= 1.q0-c/q rpr= 1.q0-q/a omrpr= q/a omrmr= c/q else q= -0.5q0*(b-disc) rmr= 1.q0-q/a rpr= 1.q0-c/q omrpr= c/q omrmr= q/a endif else a= -p2 b= p2+rm22-rm12 c= rm12 if (b.gt.0.q0) then q= -0.5q0*(b+disc) rpr= c/q rmr= q/a omrpr= 1.q0-rpr omrmr= 1.q0-rmr else q= -0.5q0*(b-disc) rpr= q/a rmr= c/q omrpr= 1.q0-rpr omrmr= 1.q0-rmr endif endif else disc= sqrt(-disc2) a= -p2 b= p2+rm22-rm12 rpr= -b/2.q0/a rmr= -b/2.q0/a omrpr= 1.q0-rpr omrmr= 1.q0-rmr rpi= +disc/2.q0/a rmi= -disc/2.q0/a endif return end * *------------------------------------------------------------------- * subroutine wtos2(y0r,y0i,omy0r,a,ebr,ebi,eci,rpr,rpi,rmr,rmi, # omrpr,omrmr,s2r,s2i) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * dimension carg(2),comarg(2),cln(2) * * rp and rm are the roots of a*x^2+b*x+c = 0 * eb= b/a and ec= c/a are kept to compute im(rp*rm)= im(ec) * and im{(y0-rp)(y0-rm)}= im{y0^2+eb*y0+ec} * a1= a*eci a2= a*(2.q0*y0r*y0i+ebr*y0i+ebi*y0r+eci) if (a1.lt.0.q0) then del= 1.q-90 else del= -1.q-90 endif if (a2.lt.0.q0) then delp= 1.q-90 else delp= -1.q-90 endif * * computes the log which occurs in association with eta-functions y0m2= y0r*y0r+y0i*y0i carg(1)= y0r/y0m2 carg(2)= -y0i/y0m2 comarg(1)= -(omy0r*y0r-y0i*y0i)/y0m2 comarg(2)= y0i/y0m2 call wtocqlnomx(carg,comarg,cln) * * in calling eta only the sign of the arguments is relevant * a11= -rpi a12= -rmi a1p= eci a21= y0i-rpi a22= y0i-rmi a2p= 2.q0*y0r*y0i+ebr*y0i+ebi*y0r+eci a31= -del a32= delp a3p= a*(delp-del) call wtorlog(y0r,y0i,omy0r,rpr,rpi,omrpr,rfpr,rfpi) call wtorlog(y0r,y0i,omy0r,rmr,rmi,omrmr,rfmr,rfmi) eta1= wtoeta(a11,a12,a1p) eta2= wtoeta(a21,a22,a2p) eta3= wtoeta(a31,a32,a3p) s2r= rfpr+rfmr-cln(2)*(eta1-eta2-eta3) s2i= rfpi+rfmi+cln(1)*(eta1-eta2-eta3) return end * *----------------------------------------------------------- * subroutine wtorlog(z0r,z0i,omz0r,z1r,z1i,omz1r,rfunr,rfuni) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * dimension ca1(2),ca2(2),ca3(2),add(2) dimension cln1(2),cln2(2),cln3(2),ct(10),sn(10) * z1m2= z1r*z1r+z1i*z1i z1m= sqrt(z1m2) z0m2= z0r*z0r+z0i*z0i z0m= sqrt(z0m2) zrt= z1m/z0m * * |z1| << |z0| * if (zrt.lt.1.q-10) then ct(1)= (z1r*z0r+z1i*z0i)/z1m/z0m sn(1)= (z0r*z1i-z1r*z0i)/z1m/z0m do k=2,10 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 sumr= ct(10) sumi= sn(10) do j=9,1,-1 sumr= sumr*zrt+ct(j) sumi= sumi*zrt+sn(j) enddo sumr= sumr*zrt sumi= sumi*zrt a1r= 1.q0+sumr a1i= sumi a2r= (z0i*z0i-omz0r*z0r)/z0m2*a1r-z0i/z0m2*a1i a2i= (z0i*z0i-omz0r*z0r)/z0m2*a1i+z0i/z0m2*a1r oma1r= -sumr oma2r= 1.q0-a2r z01r= z0r-z1r z01i= z0i-z1i den= z01r*z01r+z01i*z01i add(1)= 0.q0 add(2)= 0.q0 sign= +1.q0 else * * if z0r and/or z1r are roots very near to 1 * z0r-z1r = (1-z1r)-(1-z0r) * aomz0r= abs(omz0r) aomz1r= abs(omz1r) az1i= abs(z1i) if (aomz0r.lt.1.q-10.or.aomz1r.lt.1.q-10) then z01r= omz1r-omz0r else z01r= z0r-z1r endif z01i= z0i-z1i * * if z0 and z1r are roots very near to 1 * and |z1i| << 1 , z0i = 0 * if (aomz0r.lt.1.q-10.and.aomz1r.lt.1.q-10. # and.az1i.lt.1.q-10.and.z0i.eq.0.q0) then den= z01r*z01r+z1i*z1i a1r= z01r/z0r a1i= -z1i/z0r oma1r= z1r/z0r ca3(1)= -a1r ca3(2)= -a1i call wtocqlnx(ca3,cln3) add(1)= -qpis/6.q0-0.5q0*(cln3(1)*cln3(1)-cln3(2)*cln3(2)) add(2)= -cln3(1)*cln3(2) sign= -1.q0 else den= z01r*z01r+z01i*z01i a1r= (z0r*z01r+z0i*z01i)/den a1i= (-z0r*z01i+z0i*z01r)/den oma1r= -(z01i*z1i+z1r*z01r)/den add(1)= 0.q0 add(2)= 0.q0 sign= +1.q0 endif a2r= -(omz0r*z01r-z0i*z01i)/den a2i= (omz0r*z01i+z0i*z01r)/den oma2r= (omz1r*z01r-z01i*z1i)/den endif * * in calling eta only the sign of the arguments is relevant * a1= -z1i a2= -z01i api= z1r*z0i-z0r*z1i apii= -z0i*omz1r+z1i*omz0r call wtospence(a1r,a1i,oma1r,sp1r,sp1i) call wtospence(a2r,a2i,oma2r,sp2r,sp2i) ca1(1)= a1r ca1(2)= a1i ca2(1)= a2r ca2(2)= a2i call wtocqlnx(ca1,cln1) call wtocqlnx(ca2,cln2) etai= wtoeta(a1,a2,api) etaii= wtoeta(a1,a2,apii) rfunr= sign*sp1r-sp2r-etai*cln1(2)+etaii*cln2(2)+add(1) rfuni= sign*sp1i-sp2i+etai*cln1(1)-etaii*cln2(1)+add(2) return end * *---------------------------------------------------------------------- * real*16 function wtoeta(z1i,z2i,zpi) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * * only the sign of the arguments is relevant * if (z1i.lt.0.q0.and.z2i.lt.0.q0.and.zpi.gt.0.q0) then wtoeta= 2.q0*qpi return else if (z1i.gt.0.q0.and.z2i.gt.0.q0.and.zpi.lt.0.q0) then wtoeta= -2.q0*qpi return else wtoeta= 0.q0 return endif end * *--------------------------------------------------------------------- * subroutine wtospence(xr,xi,omxr,cli2r,cli2i) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * 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) * x(1)= xr x(2)= xi omx(1)= omxr omx(2)= -xi if (xr.lt.0.q0) then y(1)= omxr y(2)= -xi sign1= -1.q0 call wtocqlnx(x,clnx) call wtocqlnomx(x,omx,clnomx) add1(1)= qpis/6.q0-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.q0 add1(1)= 0.q0 add1(2)= 0.q0 endif omy(1)= 1.q0-y(1) omy(2)= -y(2) ym2= y(1)*y(1)+y(2)*y(2) ym= sqrt(ym2) if (ym.gt.1.q0) then z(1)= y(1)/ym2 z(2)= -y(2)/ym2 sign2= -1.q0 oy(1)= -y(1) oy(2)= -y(2) call wtocqlnx(oy,clnoy) add2(1)= -qpis/6.q0-0.5q0*((clnoy(1))**2-(clnoy(2))**2) add2(2)= -clnoy(1)*clnoy(2) else z(1)= y(1) z(2)= y(2) sign2= 1.q0 add2(1)= 0.q0 add2(2)= 0.q0 endif omz(1)= 1.q0-z(1) omz(2)= -z(2) zr= z(1) if (zr.gt.0.5q0) then t(1)= 1.q0-z(1) t(2)= -z(2) omt(1)= 1.q0-t(1) omt(2)= -t(2) sign3= -1.q0 call wtocqlnx(z,clnz) call wtocqlnomx(z,omz,clnomz) add3(1)= qpis/6.q0-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.q0-t(1) omt(2)= -t(2) sign3= 1.q0 add3(1)= 0.q0 add3(2)= 0.q0 endif call wtocqlnomx(t,omt,par) b(0)= 1.q0 b(1)= -1.q0/2.q0 b(2)= 1.q0/6.q0 b(4)= -1.q0/30.q0 b(6)= 1.q0/42.q0 b(8)= -1.q0/30.q0 b(10)= 5.q0/66.q0 b(12)= -691.q0/2730.q0 b(14)= 7.q0/6.q0 fact= 1.q0 do n=0,14 bf(n)= b(n)/fact fact= fact*(n+2.q0) 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 wtopoleg(p2r,p2i,pggf,pgglq,pggnp) implicit real*8(a-h,o-z) character*1 rio * common/wtai/rio common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * dimension bfl(2),bfh(2) dimension pggf(2),pgglq(2) dimension b0l(2),b1l(2),b21l(2) dimension b0h(2),b1h(2),b21h(2) * call wtosbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l) call wtosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h) * do i=1,2 bfl(i)= 2.d0*b21l(i)-b0l(i) bfh(i)= 2.d0*b21h(i)-b0h(i) enddo * if(rio.eq.'a') then do i=1,2 pggf(i)= 12.d0*bfl(i)+16.d0/3.d0*bfh(i) pgglq(i)= 44.d0/3.d0*bfl(i) enddo else if(rio.eq.'i') then pggf(1)= 0.d0 pgglq(1)= 0.d0 pggf(2)= 12.d0*bfl(2)+16.d0/3.d0*bfh(2) pgglq(2)= 44.d0/3.d0*bfl(2) endif * if(rio.eq.'i') then pggnp= 0.d0 else if(rio.eq.'a') then zm2= zm*zm p2z= -zm2 call wtorbff(p2z,em2,em2,b0eez,b1eez,b21eez) call wtorbff(p2z,amm2,amm2,b0mmz,b1mmz,b21mmz) call wtorbff(p2z,tm2,tm2,b0tauz,b1tauz,b21tauz) call wtorbff0(em2,em2,b0ee0,b1ee0,b21ee0) call wtorbff0(amm2,amm2,b0mm0,b1mm0,b21mm0) call wtorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0) bfez= 2.d0*b21eez-b0eez bfmz= 2.d0*b21mmz-b0mmz bftauz= 2.d0*b21tauz-b0tauz bfe0= 2.d0*b21ee0-b0ee0 bfm0= 2.d0*b21mm0-b0mm0 bftau0= 2.d0*b21tau0-b0tau0 pggl= 4.d0*(bfez+bfmz+bftauz) pggl0= 4.d0*(bfe0+bfm0+bftau0) b= 0.00299d0 c= 1.d0 an= -b*log(1.d0+c*zm2)-1.d0/alphai/4.d0/pi*(pggl-pggl0)+ # 1.d0-128.89d0/alphai ap2x= abs(p2r) px= sqrt(ap2x) if(px.lt.0.3d0) then a= 0.d0 b= 0.00835d0 c= 1.d0 else if(px.lt.3.d0) then a= 0.d0 b= 0.00238d0 c= 3.927d0 else if(px.lt.100.d0) then a= 0.00165d0 b= 0.00299d0 c= 1.d0 else if(px.gt.100.d0) then a= 0.00221d0 b= 0.00293d0 c= 1.d0 endif * da= an-0.00165q0 a= a+da pggnp= 4.d0*pi*alphai*(a+b*log(1.d0+c*ap2x)) endif * return end * *------------------------------------------------------------------------- * subroutine wtopolez0(pggf0,fw0) implicit real*8(a-h,o-z) character*1 rio * common/wtai/rio common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * if(rio.eq.'i') then pggf0= 0.d0 fw0= 0.d0 return else if(rio.eq.'a') then call wtorbff0(em2,em2,b0ee0,b1ee0,b21ee0) call wtorbff0(amm2,amm2,b0mm0,b1mm0,b21mm0) call wtorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0) call wtorbff0(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 wtorbff0(tqm2,dmy2,b0tb0,b1tb0,b21tb0) * fw0= -3.d0*tqm2*(b0tb0+b1tb0) return endif * end * *------------------------------------------------------------- * *------------------------------------------------------------- * subroutine wtofst(nt,st,fvect,iflag) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * dimension zs(2),zl(2) dimension st(nt),fvect(nt) dimension pggfz(2),pgglqz(2),p3qz(2),fzz(2),p3qw(2),fw(2), # fzw(2),fwz(2) * tqm2= st(1) * *-----alpha(zm) * wm2= wm*wm zm2= zm*zm pzr= -zm2 pzi= 0.d0 call wtopoleg(pzr,pzi,pggfz,pgglqz,pggnpz) call wtopolez0(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 call wtopole(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 call wtopole(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 wtofsw(n,sw,fvecw,iflag) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * 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 * call wtopolez0(pggf0,fw0) p2r= -wm2 p2i= 0.d0 call wtopole(p2r,p2i,p3qw,fzw,fw) p2r= -swr p2i= -swi call wtopole(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 wtofsz(n,sz,fvecz,iflag) implicit real*8(a-h,o-z) * common/wtafsz/g2z common/wtiaz/alzr,alzi,alzir,alzii common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * 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) * zm2= zm*zm apis= 16.d0*pis call wtopolez0(pggf0,fw0) pzr= -zm2 pzi= 0.d0 call wtopoleg(pzr,pzi,pggfz,pgglqz,pggnpz) call wtopole(pzr,pzi,p3qz,fzz,fwz) p2r= -szr p2i= -szi call wtopoleg(p2r,p2i,pggf,pgglq,pggnp) call wtopole(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 wtocg(n,zr,zi,omzr,gfr,gfi) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * dimension gfr(n),gfi(n) dimension ca(2,10),aux(2),zp(2),ct(16),sn(16) dimension a(2,4),ra(16),z(2),omz(2),oz(2),clnomz(2),clnoz(2) * data ra/-0.2q+0,0.666666666666667q-1, # -0.952380952380952q-2,-0.396825396825397q-3, # 0.317460317460317q-3,-0.132275132275132q-4, # -0.962000962000962q-5,0.105218855218855q-5, # 0.266488361726450q-6,-0.488745528428000q-7, # -0.675397500794000q-8,0.190720263471000q-8, # 0.153663007690000q-9,-0.679697905790000q-10, # -0.293683556000000q-11,0.228836696000000q-11/ * z(1)= zr z(2)= zi omz(1)= omzr omz(2)= -zi if (zr.eq.0.q0.and.zi.eq.0.q0) then do k=1,n gfr(k)= -1.q0/k/k gfi(k)= 0.q0 enddo return else if (zr.eq.1.q0.and.zi.eq.0.q0) then a(1,1)= -1.q0 a(2,1)= qpi do j=2,4 a(1,j)= ((j-1.q0)*a(1,j-1)-1.q0/j)/j a(2,j)= (j-1.q0)/j*a(2,j-1) enddo do k=1,n gfr(k)= a(1,k) gfi(k)= a(2,k) enddo return else zmod2= zr*zr+zi*zi zmod= sqrt(zmod2) call wtocqlnx(omz,clnomz) * * |z| < 4 * if (zmod.lt.4.q0) then oz(1)= -z(1) oz(2)= -z(2) call wtocqlnx(oz,clnoz) ca(1,1)= omz(1)*clnomz(1)-omz(2)*clnomz(2)+z(1)*clnoz(1)- # z(2)*clnoz(2)-1.q0 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) gfi(1)= ca(2,1) return else do j= 2,n jm1= j-1 ca(1,j)= ((j-1.q0)*(z(1)*ca(1,jm1)-z(2)*ca(2,jm1))+ # omz(1)*clnomz(1)-omz(2)*clnomz(2)-1.q0/j)/j ca(2,j)= ((j-1.q0)*(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) gfi(k)= ca(2,k) enddo return endif * * |z| > 4 * else aux(1)= (-zr*omzr+zi**2)/zmod2 aux(2)= zi/zmod2 call wtocqlnx(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.q0/jp1)*z(1)+ca(2,jp1)*z(2))/zmod2 ca(2,j)= (ca(2,jp1)*z(1)-(ca(1,jp1)+1.q0/jp1)*z(2))/zmod2 enddo do k=1,n gfr(k)= (ca(1,k)+clnomz(1))/k gfi(k)= (ca(2,k)+clnomz(2))/k enddo return endif endif * end *