2
subroutine ffxdb0(cdb0,cdb0p,xp,xma,xmb,ier)
3
***#[*comment:***********************************************************
5
* Calculates the the derivative of the two-point function with *
6
* respect to p2 and the same times p2 (one is always well-defined)*
8
* Input: xp (real) k2, in B&D metric *
12
* Output: cdb0 (complex) dB0/dxp *
13
* cdb0p (complex) xp*dB0/dxp *
14
* ier (integer) # of digits lost, if >=100: error *
18
***#]*comment:***********************************************************
25
DOUBLE COMPLEX cdb0,cdb0p
26
DOUBLE PRECISION xp,xma,xmb
31
DOUBLE PRECISION dmamb,dmap,dmbp
40
print *,'ffxdb0: input:'
41
print *,'xma,xmb,xp,ier = ',xma,xmb,xp,ier
44
if ( xma .lt. 0 .or. xmb .lt. 0 ) then
45
print *,'ffxdb0: error: xma,b < 0: ',xma,xmb
56
if ( abs(dmamb) .lt. xloss*abs(xma) .and. xma .ne. xmb )
57
+ call ffwarn(97,ier0,dmamb,xma)
58
if ( abs(dmap) .lt. xloss*abs(xp) .and. xp .ne. xma )
59
+ call ffwarn(98,ier0,dmap,xp)
60
if ( abs(dmbp) .lt. xloss*abs(xp) .and. xp .ne. xmb )
61
+ call ffwarn(99,ier0,dmbp,xp)
65
call ffxdbp(cdb0,cdb0p,xp,xma,xmb,dmap,dmbp,dmamb,ier)
66
if ( lwrite ) print *,'B0'' = ',cdb0,cdb0p,ier
71
subroutine ffxdbp(cdb0,cdb0p,xp,xma,xmb,dmap,dmbp,dmamb,ier)
72
***#[*comment:***********************************************************
74
* calculates the derivatives of the two-point function *
75
* Veltman) for all possible cases: masses equal, unequal, *
78
* Input: xp (real) p.p, in B&D metric *
81
* dm[ab]p (real) xm[ab] - xp *
82
* dmamb (real) xma - xmb *
84
* Output: cdb0 (complex) B0' = dB0/dxp *
85
* cdb0p (complex) xp*dB0/dxp *
86
* ier (integer) 0=ok,>0=numerical problems,>100=error *
90
***#]*comment:***********************************************************
97
DOUBLE COMPLEX cdb0,cdb0p
98
DOUBLE PRECISION xp,xma,xmb,dmap,dmbp,dmamb
102
integer i,initeq,jsign,initir
103
DOUBLE PRECISION ax,ffbnd,
104
+ xprceq,bdeq01,bdeq05,bdeq11,bdeq17,bdeq25,
105
+ xprcn3,bdn301,bdn305,bdn310,bdn315,
106
+ xprcn5,bdn501,bdn505,bdn510,bdn515,
107
+ xprec0,bdn001,bdn005,bdn010,bdn015,bdn020
108
DOUBLE PRECISION xcheck,xm,dmp,xm1,xm2,dm1m2,dm1p,
109
+ dm2p,s,s1,s1a,s1b,s1p,s2,s2a,s2b,s2p,x,y,som,
110
+ xlam,slam,xlogmm,alpha,alph1,xnoe,xpneq(30),
111
+ xx,dfflo1,dfflo3,d1,d2,diff,h,a,b,c,d,beta,
112
+ betm2n,xmax,s1c,s1d,s1e,s1f,s3
113
DOUBLE COMPLEX cc,zxfflg
114
save initeq,xpneq,initir,
115
+ xprceq,bdeq01,bdeq05,bdeq11,bdeq17,bdeq25,
116
+ xprcn3,bdn301,bdn305,bdn310,bdn315,
117
+ xprcn5,bdn501,bdn505,bdn510,bdn515,
118
+ xprec0,bdn001,bdn005,bdn010,bdn015,bdn020
123
DOUBLE PRECISION delta
137
xcheck = xma - xmb - dmamb
138
if ( abs(xcheck) .gt. precx*max(abs(xma),abs(xmb),abs(
139
+ dmamb))/xloss ) then
140
print *,'ffxdbp: input not OK, dmamb <> xma-xmb',xcheck
142
xcheck = -xp + xma - dmap
143
if ( abs(xcheck) .gt. precx*max(abs(xp),abs(xma),abs(
144
+ dmap))/xloss ) then
145
print *,'ffxdbp: input not OK, dmap <> xma - xp',xcheck
147
xcheck = -xp + xmb - dmbp
148
if ( abs(xcheck) .gt. precx*max(abs(xp),abs(xmb),abs(
149
+ dmbp))/xloss ) then
150
print *,'ffxdbp: input not OK, dmbp <> xmb - xp',xcheck
156
* sort according to the type of masscombination encountered:
157
* 100: both masses zero, 200: one equal to zero, 300: both equal
160
if ( xma .eq. 0 ) then
161
if ( xmb .eq. 0 ) then
168
if ( xmb .eq. 0 ) then
172
elseif ( dmamb .eq. 0 ) then
176
elseif ( xma .gt. xmb ) then
191
* #[ both masses equal to zero:
193
if ( xp.ne.0 ) cdb0 = -1/xp
196
* #] both masses equal to zero:
197
* #[ one mass equal to zero:
200
* special case xp = 0
202
if ( xp .eq. 0 ) then
207
* special case xp = xm
209
elseif ( dmp.eq.0 ) then
211
if ( DBLE(cmipj(1,3)).lt.DBLE(cmipj(2,3)) ) then
212
cdb0p = -1 - log(cmipj(1,3)*DBLE(1/xm))
214
cdb0p = -1 - log(cmipj(2,3)*DBLE(1/xm))
217
if ( initir.eq.0 ) then
219
print *,'ffxdb0: IR divergent B0'', using cutoff ',
222
if ( delta.eq.0 ) then
226
cdb0p = -1 + log(xm/delta)/2
229
cdb0 = cdb0p*(1/DBLE(xp))
237
if ( ax .lt. xloss ) then
238
* #[ Taylor expansion:
239
if ( xprec0 .ne. precx ) then
241
bdn001 = ffbnd(2,1,xninv)
242
bdn005 = ffbnd(2,5,xninv)
243
bdn010 = ffbnd(2,10,xninv)
244
bdn015 = ffbnd(2,15,xninv)
245
bdn020 = ffbnd(2,20,xninv)
247
if ( lwarn .and. ax .gt. bdn020 ) then
248
call ffwarn(15,ier,precx,xninv(21)*ax**20)
250
if ( ax .gt. bdn015 ) then
251
som = x*(xninv(17) + x*(xninv(18) + x*(xninv(19) +
252
+ x*(xninv(20) + x*(xninv(21) )))))
256
if ( ax .gt. bdn010 ) then
257
som = x*(xninv(12) + x*(xninv(13) + x*(xninv(14) +
258
+ x*(xninv(15) + x*(xninv(16) + som )))))
260
if ( ax .gt. bdn005 ) then
261
som = x*(xninv(7) + x*(xninv(8) + x*(xninv(9) +
262
+ x*(xninv(10) + x*(xninv(11) + som )))))
264
if ( ax .gt. bdn001 ) then
265
som = x*(xninv(3) + x*(xninv(4) + x*(xninv(5) +
266
+ x*(xninv(6) + som ))))
268
cdb0p = x*(xninv(2) + som)
270
print *,'cdb0p = ',cdb0p
271
print *,'verg ',-1 - xm/xp*dfflo1(x,ier),1
273
* #] Taylor expansion:
277
cdb0p = -(1 + s*xm/xp)
278
if ( xp.gt.xm ) cdb0p = cdb0p+DCMPLX(DBLE(0),DBLE(xm/xp*pi))
281
cdb0 = cdb0p*(1/DBLE(xp))
283
* #] one mass equal to zero:
284
* #[ both masses equal:
287
* Both masses are equal. Not only this speeds up things, some
288
* cancellations have to be avoided as well.
290
* first a special case
292
if ( abs(xp) .lt. 8*xloss*xm ) then
293
* -#[ taylor expansion:
295
* a Taylor expansion seems appropriate as the result will go
296
* as k^2 but seems to go as 1/k !!
298
*--#[ data and bounds:
299
if ( initeq .eq. 0 ) then
303
xpneq(i) = - xpneq(i-1)*DBLE(i)/DBLE(2*(2*i+1))
306
if (xprceq .ne. precx ) then
308
* calculate the boundaries for the number of terms to be
309
* included in the taylorexpansion
312
bdeq01 = ffbnd(1,1,xpneq)
313
bdeq05 = ffbnd(1,5,xpneq)
314
bdeq11 = ffbnd(1,11,xpneq)
315
bdeq17 = ffbnd(1,17,xpneq)
316
bdeq25 = ffbnd(1,25,xpneq)
318
*--#] data and bounds:
321
if ( lwarn .and. ax .gt. bdeq25 ) then
322
call ffwarn(15,ier,precx,abs(xpneq(25))*ax**25)
324
if ( ax .gt. bdeq17 ) then
325
som = x*(xpneq(18) + x*(xpneq(19) + x*(xpneq(20) +
326
+ x*(xpneq(21) + x*(xpneq(22) + x*(xpneq(23) +
327
+ x*(xpneq(24) + x*(xpneq(25) ))))))))
331
if ( ax .gt. bdeq11 ) then
332
som = x*(xpneq(12) + x*(xpneq(13) + x*(xpneq(14) +
333
+ x*(xpneq(15) + x*(xpneq(16) + x*(xpneq(17) + som ))))
336
if ( ax .gt. bdeq05 ) then
337
som = x*(xpneq(6) + x*(xpneq(7) + x*(xpneq(8) + x*(
338
+ xpneq(9) + x*(xpneq(10) + x*(xpneq(11) + som ))))))
340
if ( ax .gt. bdeq01 ) then
341
som = x*(xpneq(2) + x*(xpneq(3) + x*(xpneq(4) + x*(
342
+ xpneq(5) + som ))))
344
cdb0p = -x*(xpneq(1)+som)
346
print *,'ffxdbp: m1 = m2, Taylor expansion in ',x
347
print *,'cdb0p = ',cdb0p
350
cdb0 = cdb0p*(1/DBLE(xp))
355
* -#] taylor expansion:
361
call ffxlmb(xlam,-xp,-xm,-xm,dmp,dmp,x0,ier)
362
if ( xlam .eq. 0 ) then
365
elseif ( xlam .gt. 0 ) then
370
if ( abs(s2) .gt. xloss*slam ) then
378
if ( ax .lt. xalogm ) then
379
if ( lwarn ) call ffwarn(16,ier,ax,xalogm)
381
elseif( ax-1 .lt. .1 .and. s2 .gt. 0 ) then
382
* In this case a quicker and more accurate way is to
383
* calculate log(1-x).
385
* the following line is superfluous.
386
if ( lwarn .and. abs(s2) .lt. xloss*slam )
387
+ call ffwarn(17,ier,s2,slam)
388
s = 2*xm/slam*dfflo1(s2/(2*xm),ier)
390
* finally the normal case
391
s = 2*xm/slam*log(ax)
392
if ( jsign .eq. -1 ) s = -s
394
if ( xp .gt. 2*xm ) then
395
* in this case ( xlam>0, so xp>(2*m)^2) ) there also
396
* is an imaginary part
402
* the root is complex (k^2 between 0 and (2*m1)^2)
404
s = 4*xm/slam*atan2(xp,slam)
407
if (lwrite) print *,'s = ',s
409
if ( lwarn .and. abs(xx).lt.xloss ) call ffwarn(18,ier,xx,x1)
410
cdb0p = DCMPLX(DBLE(xx),DBLE(y))
411
cdb0 = cdb0p*(1/DBLE(xp))
415
* #] both masses equal:
416
* #[ unequal nonzero masses:
417
* -#[ get log(xm2/xm1):
420
if ( 1 .lt. xalogm*x ) then
423
elseif ( abs(x-1) .lt. xloss ) then
424
xlogmm = dfflo1(dm1m2/xm1,ier)
428
* -#] get log(xm2/xm1):
431
* first a special case
433
if ( xp .eq. 0 ) then
435
* repaired 19-nov-1993, see b2.frm
437
s1 = xm1*xm2*xlogmm/dm1m2**3
438
s2 = (xm1+xm2)/(2*dm1m2**2)
440
if ( abs(s) .lt. xloss**2*s2 ) then
444
h = dfflo3(dm1m2/xm1,ier)
447
s3 = xm1**2*h/dm1m2**3
449
if ( abs(s) .lt. xloss*max(abs(s2),abs(s3)) ) then
450
call ffwarn(228,ier,s,s2)
460
* proceeding with the normal case
462
call ffxlmb(xlam,-xp,-xm2,-xm1,dm2p,dm1p,dm1m2,ier)
463
diff = xlam + xp*(dm2p+xm1)
464
if ( lwrite ) print *,'diff = ',diff
465
if ( abs(diff) .lt. xloss*xlam ) then
466
h = dm1m2**2 - xp*(xm1+xm2)
467
if ( lwrite ) print *,'diff+= ',h
468
if ( abs(h) .lt. xloss*dm1m2**2 ) then
469
if ( dm1m2**2 .lt. abs(xlam) ) diff = h
471
call ffwarn(221,ier,diff,min(dm1m2**2,abs(xlam)))
475
if ( xlam .eq. 0 ) then
478
elseif ( xlam .gt. 0 ) then
479
* cases k^2 < -(m2+m1)^2 or k^2 > -(m2-m1)^2:
481
* first try the normal way
485
if ( abs(s2) .gt. xloss*slam ) then
492
s2 = s2**2/(4*xm1*xm2)
493
if ( abs(s2) .lt. xalogm ) then
496
elseif ( abs(s2-1) .lt. xloss ) then
497
if ( jsign.eq.1 ) then
498
if (lwrite) print *,'s2 ',-diff/(2*slam*xp)*log(s2)
499
s2 = -slam*(s2a+slam)/(2*xm1*xm2)
500
s2 = -diff/(2*slam*xp)*dfflo1(s2,ier)
503
print *,'ffxdb0: untested: s2 better in first try'
504
if (lwrite) print *,'s2 ',+diff/(2*slam*xp)*log(s2)
505
s2 = +slam*(s2a-slam)/(2*xm1*xm2)
506
s2 = +diff/(2*slam*xp)*dfflo1(s2,ier)
508
if ( lwrite ) print *,'s2+ ',s2,jsign
510
s2 = -diff/(2*slam*xp)*log(s2)
511
if ( jsign .eq. -1 ) s2 = -s2
513
s1 = -dm1m2*xlogmm/(2*xp)
516
print *,'ffxdbp: lam>0, first try, xx = ',xx,s1,s2,-1
519
if ( abs(xx) .lt. xloss**2*max(abs(s1),abs(s2)) ) then
521
* this is unacceptable, try a better solution
522
s1a = diff + slam*dm1m2
523
if (lwrite) print *,'s1 = ',-s1a/(2*xp*slam),diff/
525
if ( abs(s1a) .gt. xloss*diff ) then
527
s1 = -s1a/(2*xp*slam)
529
* by division a more accurate form can be found
530
s1 = -2*xm1*xm2*xp/(slam*(diff - slam*dm1m2))
531
if (lwrite) print *,'s1+= ',s1
535
if ( abs(xp) .lt. xm2 ) then
541
if (lwrite) print *,'s2 = ',s2/(2*xm2),slam/(2*xm2)
542
if ( abs(s2) .gt. xloss*slam ) then
543
* at least reasonable
547
s2 = (2*xp) / (s2a+slam)
548
if (lwrite) print *,'s2+= ',s2
550
if ( abs(s2) .lt. .1 ) then
551
* choose a quick way to get the logarithm
553
elseif ( s2.eq.1 ) then
554
print *,'ffxdbp: error: arg log would be 0!'
555
print *,' xp,xma,xmb = ',xp,xma,xmb
559
s2 = zxfflg(h,0,c0,ier)
561
s2 = -diff/(slam*xp)*s2
564
print *,'ffxdbp: lam>0, 2nd try, xx = ',xx,s1,s2,-1
567
if ( abs(xx) .lt. xloss**2*max(abs(s1),abs(s2)) ) then
569
* (we accept two times xloss because that's the same
571
* A Taylor expansion might work. We expand
572
* inside the logs. Only do the necessary work.
579
d = sqrt((2/xnoe)**2 + 1/xp**2)
580
call ffroot(d1,d2,a,b,c,d,ier)
586
alpha = beta*diff/slam
588
if ( alph1 .lt. xloss ) then
589
s1a = 4*xp**2*xm1*xm2/(slam*dm1m2*(diff-slam*
591
s1b = -diff/slam*4*xm1*xp/(dm1m2*xnoe*(2*xp-
595
call ffroot(d1,d2,a,b,c,d,ier)
601
d1 = s1a + s1b - diff/slam*betm2n
603
print *,'alph1 = ',d1,s1a,s1b,-diff/slam*
605
print *,'verg ',1-alpha
607
xmax = max(abs(s1a),abs(s1b))
608
if ( xmax .lt. 1 ) then
613
if ( lwarn .and. abs(alph1).lt.xloss*xmax ) then
614
call ffwarn(222,ier,alph1,xmax)
615
if ( lwrite ) print *,'d1,s1a,s2b,... = ',
616
+ d1,s1a,s1b,diff/slam*betm2n
619
betm2n = beta - 2/xnoe
622
print *,' s1 - alph1 = ',s1-alph1
623
print *,' s2 - alpha = ',s2-alpha
632
if ( abs(s2p) .lt. xloss*abs(s2) ) then
634
* determine the boundaries for 1,5,10,15 terms
635
if ( xprcn5 .ne. precx ) then
637
bdn501 = ffbnd(3,1,xinfac)
638
bdn505 = ffbnd(3,5,xinfac)
639
bdn510 = ffbnd(3,10,xinfac)
640
bdn515 = ffbnd(3,15,xinfac)
645
if ( lwarn .and. ax .gt. bdn515 ) then
646
* do not do the Taylor expansion
647
call ffwarn(23,ier,s2p,s2)
650
if ( ax .gt. bdn510 ) then
651
s2a = x*(xinfac(13) + x*(xinfac(14) + x*(
652
+ xinfac(15) + x*(xinfac(16) + x*(
657
if ( ax .gt. bdn505 ) then
658
s2a = x*(xinfac(8) + x*(xinfac(9) + x*(
659
+ xinfac(10) + x*(xinfac(11) + x*(
660
+ xinfac(12) + s2a)))))
662
if ( ax .gt. bdn501 ) then
663
s2a = x*(xinfac(4) + x*(xinfac(5) + x*(
664
+ xinfac(6) + x*(xinfac(7) + s2a))))
666
s2a = x**3*(xinfac(3)+s2a)
667
s2b = 2*xp/xnoe*(s2a + x**2/2)
669
if ( lwarn .and. abs(s2p).lt.xloss*abs(s2a) )
670
+ call ffwarn(24,ier,s2p,s2a)
671
s2p = -diff/(xp*slam)*dfflo1(s2p,ier)
673
print *,'ffxdbp: Taylor expansion of s2-a'
675
print *,' gives s2p = ',s2p
685
if ( abs(s1p) .lt. xloss*abs(s1) ) then
687
* determine the boundaries for 1,5,10,15 terms
688
if ( xprcn3 .ne. precx ) then
690
bdn301 = ffbnd(3,1,xinfac)
691
bdn305 = ffbnd(3,5,xinfac)
692
bdn310 = ffbnd(3,10,xinfac)
693
bdn315 = ffbnd(3,15,xinfac)
697
x = slam*(diff-slam*dm1m2)*alph1/(2*xp*xm1*xm2)
698
h = (2*xp*(xm1+xm2) - xp**2)/(slam-dm1m2)
700
if ( lwarn .and. ax .gt. bdn315 ) then
701
* do not do the Taylor expansion
702
call ffwarn(21,ier,s1p,s1)
706
* see form job gets1.frm
708
s1b = diff*(diff-slam*dm1m2)*betm2n/(2*xp*xm1*
710
s1c = 1/(xm1*xnoe*(2*xp-xnoe))*(
711
+ xp*( 4*xp*xm2 + 2*dm1m2**2/xm2*(xp-h) +
712
+ 2*dm1m2*(3*xp-h) - 8*dm1m2**2 )
713
+ - 2*dm1m2**3/xm2*(3*xp-h)
717
print *,'s1c was ',-2*xp/dm1m2 + 2*diff*
718
+ (diff-slam*dm1m2)/(xm2*dm1m2*xnoe*(2*xp-
720
print *,' en is ',s1c
721
print *,'s1b+s1c was ',dm1m2/xm1-x
722
print *,' en is ',s1b+s1c
726
if ( ax .gt. bdn310 ) then
727
s1a = x*(xinfac(13) + x*(xinfac(14) + x*(
728
+ xinfac(15) + x*(xinfac(16) + x*(
733
if ( ax .gt. bdn305 ) then
734
s1a = x*(xinfac(8) + x*(xinfac(9) + x*(
735
+ xinfac(10) + x*(xinfac(11) + x*(
736
+ xinfac(12) + s1a)))))
738
if ( ax .gt. bdn301 ) then
739
s1a = x*(xinfac(4) + x*(xinfac(5) + x*(
740
+ xinfac(6) + x*(xinfac(7) + s1a))))
742
s1a = -x**3 *(xinfac(3) + s1a)
743
s1f = dm1m2/xm1*(x**2/2 - s1a)
744
s1p = s1e + s1d + s1c + s1b + s1a + s1f
745
xmax = max(abs(s1a),abs(s1b),abs(s1c),abs(s1d),
747
if ( lwarn .and. abs(s1p).lt.xloss*xmax ) then
748
call ffwarn(223,ier,s1p,xmax)
750
+ print *,'s1p,s1e,s1d,s1c,s1b,s1a,s1f = '
751
+ ,s1p,s1e,s1d,s1c,s1b,s1a,s1f
753
s1p = s*dfflo1(s1p,ier)
762
print *,'ffxdbp: Taylor exp. of s1-(1-a)'
764
print *,' gives s1p = ',s1p
765
print *,' verg ',s*log(xm2/xm1
775
if ( lwarn .and. abs(xx) .lt. xloss*abs(s1p) ) then
776
call ffwarn(25,ier,xx,s1p)
782
if ( xp .gt. xm1+xm2 ) then
783
*--#[ imaginary part:
784
* in this case ( xlam>0, so xp>(m1+m2)^2) ) there also
785
* is an imaginary part
786
y = -pi*diff/(slam*xp)
789
*--#] imaginary part:
792
* the root is complex (k^2 between -(m1+m2)^2 and -(m2-m1)^2)
796
s1 = -(dm1m2/(2*xp))*xlogmm
797
s2 = -diff/(slam*xp)*atan2(slam,xnoe)
800
print *,'ffxdbp: lam<0, first try, xx = ',xx,s1,s2,-1
801
* alpha = -xlam/(2*xp*xnoe)
802
* alph1 = -(xp**2-dm1m2**2)/(2*xp*xnoe)
803
* print *,' alpha = ',alpha
804
* print *,' s1 = ',s1,' - 2alph1 = ',s1-2*alph1
805
* print *,' s2 = ',s2,' - 2alpha = ',s2-2*alpha
808
if ( lwarn .and. abs(xx).lt.xloss**2*max(abs(s1),abs(s2)) )
810
call ffwarn(224,ier,xx,max(abs(s1),abs(s2)))
815
cdb0p = DCMPLX(DBLE(xx),DBLE(y))
816
cdb0 = cdb0p*(1/DBLE(xp))
819
* #] unequal nonzero masses:
823
print *,'cdb0 = ',cdb0,cdb0p
added
removed
vendor/CutTools/src/qcdloop/ffxdb0.f
