bub=4.d0*tauu/TC
alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10)))
z=pi*4.d0*tauu*tauu*alf*(1.+nf/6.d0)
+ if (z < 0.) write(6,*) "Warning in PLJETR: z < 0.", z
bubs=dsqrt(abs(z))/TC
alfs=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bubs,1.d-10)))
phmin2=z
* if quark production outside the QGP then
* arg=(((dsin(uu)*cosh(uu))**2)+((dcos(uu)*sinh(uu))**2))/(2.d0*uu*uu);
* here quark production inside the QGP
- arg=((dcos(uu)*cosh(uu))**2)+((dsin(uu)*sinh(uu))**2)
- if (z < 0.) write(6,*) "Warning in PLFUN1: z < 0.", z
+c AM: better numerical stability
+ arg = ((dcos(uu)*cosh(uu))**2)+((dsin(uu)*sinh(uu))**2)
+ if (uu .lt. 5.) then
+ xlogarg = dlog(max(arg,1.d-20))
+ else
+ xlogarg = log(0.25) + 2. * uu
+ endif
+
+
gl1=(ygl/(cfac*abs(z)))**0.3333333d0
gl2=(um/epa)**1.333333d0
- dc=1.d0/((1.d0+((gl1*gl2*or)**1.5d0))**2) ! massive parton
+ dc=1.d0/((1.d0+((gl1*gl2*or)**1.5d0))**2) ! massive parton
c dc=1.d0 !massless parton
- plfun1=dc*3.d0*alfs*ygl*dlog(max(arg,1.d-20))*spinf/(pi*al*or)
+c plfun1 = dc*3.d0*alfs*ygl*dlog(max(arg,1.d-20))*spinf/(pi*al*or)
+ plfun1 = dc*3.d0*alfs*ygl*xlogarg*spinf/(pi*al*or)
return
end