2
subroutine ffxb0(cb0,d0,xmu,xp,xma,xmb,ier)
3
***#[*comment:***********************************************************
5
* Calculates the the two-point function (cf 't Hooft and Veltman) *
6
* we include an overall factor 1/(i*pi^2) relative to FormF *
8
* Input: d0 (real) infinity arising from renormalization *
9
* xmu (real) renormalization mass *
10
* xp (real) k2, in B&D metric *
14
* Output: cb0 (complex) B0, the two-point function, *
15
* ier (integer) # of digits lost, if >=100: error *
19
***#]*comment:***********************************************************
27
DOUBLE PRECISION d0,xmu,xp,xma,xmb
33
DOUBLE PRECISION dmamb,dmap,dmbp,xm,absc
41
absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
46
print *,'ffxb0: nevent,id = ',nevent,id,' input:'
47
print *,'xma,xmb,xp,ier = ',xma,xmb,xp,ier
50
if ( xma .lt. 0 .or. xmb .lt. 0 ) then
51
print *,'ffxb0: error: xma,b < 0: ',xma,xmb
62
if ( abs(dmamb) .lt. xloss*abs(xma) .and. xma .ne. xmb )
63
+ call ffwarn(97,ier0,dmamb,xma)
64
if ( abs(dmap) .lt. xloss*abs(xp) .and. xp .ne. xma )
65
+ call ffwarn(98,ier0,dmap,xp)
66
if ( abs(dmbp) .lt. xloss*abs(xp) .and. xp .ne. xmb )
67
+ call ffwarn(99,ier0,dmbp,xp)
71
call ffxb0p(cb0p,xp,xma,xmb,dmap,dmbp,dmamb,ier)
72
if ( xma .eq. 0 ) then
73
if ( xmb .eq. 0 ) then
78
elseif ( xmb .eq. 0 ) then
83
if ( xmu .ne. 0 ) xm = xm/xmu**2
84
if ( abs(xm) .gt. xalogm ) then
85
cb0 = DBLE(d0 - log(xm)/2) - cb0p
86
if ( lwarn .and. absc(cb0).lt.xloss*max(abs(d0),absc(cb0p)))
87
+ call ffwarn(150,ier,absc(cb0),max(abs(d0),absc(cb0p)))
92
if ( lwrite ) print *,'B0 = ',cb0,ier
97
subroutine ffxb0p(cb0p,xp,xma,xmb,dmap,dmbp,dmamb,ier)
98
***#[*comment:***********************************************************
100
* calculates the two-point function (see 't Hooft and *
101
* Veltman) for all possible cases: masses equal, unequal, *
104
* Input: xp (real) p.p, in B&D metric *
105
* xma (real) mass2, *
106
* xmb (real) mass2, *
107
* dm[ab]p (real) xm[ab] - xp *
108
* dmamb (real) xma - xmb *
110
* Output: cb0p (complex) B0, the two-point function, minus *
111
* log(xm1*xm2)/2, delta and ipi^2 *
112
* ier (integer) 0=ok, 1=numerical problems, 2=error *
116
***#]*comment:***********************************************************
124
DOUBLE PRECISION xp,xma,xmb,dmap,dmbp,dmamb
128
integer i,initeq,initn1,iflag,jsign,init
129
DOUBLE PRECISION ax,ay,ffbnd,
130
+ xprceq,bdeq01,bdeq05,bdeq11,bdeq17,bdeq25,
131
+ xprcn1,bdn101,bdn105,bdn110,bdn115,bdn120,
132
+ xprnn2,bdn201,bdn205,bdn210,bdn215,bdn220,
133
+ xprcn3,bdn301,bdn305,bdn310,bdn315,
134
+ xprcn5,bdn501,bdn505,bdn510,bdn515,
136
DOUBLE PRECISION xcheck,xm,dmp,xm1,xm2,dm1m2,dm1p,
137
+ dm2p,s,s1,s1a,s1b,s1p,s2,s2a,s2b,s2p,x,y,som,
138
+ xlam,slam,xlogmm,alpha,alph1,xnoe,xpneq(30),
139
+ xpnn1(30),xx,xtel,dfflo1
140
DOUBLE COMPLEX cs2a,cs2b,cs2p,c,cx
141
save initeq,initn1,init,xpneq,xpnn1,
142
+ xprceq,bdeq01,bdeq05,bdeq11,bdeq17,bdeq25,
143
+ xprcn1,bdn101,bdn105,bdn110,bdn115,bdn120,
144
+ xprnn2,bdn201,bdn205,bdn210,bdn215,bdn220,
145
+ xprcn3,bdn301,bdn305,bdn310,bdn315,
146
+ xprcn5,bdn501,bdn505,bdn510,bdn515
165
absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
169
xcheck = xma - xmb - dmamb
170
if ( abs(xcheck) .gt. precx*max(abs(xma),abs(xmb),abs(
171
+ dmamb))/xloss ) then
172
print *,'ffxb0q: input not OK, dmamb <> xma-xmb',xcheck
174
xcheck = -xp + xma - dmap
175
if ( abs(xcheck) .gt. precx*max(abs(xp),abs(xma),abs(
176
+ dmap))/xloss ) then
177
print *,'ffxb0q: input not OK, dmap <> xma - xp',xcheck
179
xcheck = -xp + xmb - dmbp
180
if ( abs(xcheck) .gt. precx*max(abs(xp),abs(xmb),abs(
181
+ dmbp))/xloss ) then
182
print *,'ffxb0q: input not OK, dmbp <> xmb - xp',xcheck
186
* #[ fill some dotproducts:
188
call ffdot2(fpij2,xp,xma,xmb,dmap,dmbp,dmamb,ier)
190
* #] fill some dotproducts:
193
* sort according to the type of masscombination encountered:
194
* 100: both masses zero, 200: one equal to zero, 300: both equal
197
if ( xma .eq. 0 ) then
198
if ( xmb .eq. 0 ) then
205
if ( xmb .eq. 0 ) then
209
elseif ( dmamb .eq. 0 ) then
213
elseif ( xma .gt. xmb ) then
228
* #[ both masses equal to zero:
230
if ( xp .lt. -xalogm ) then
232
elseif ( xp .gt. xalogm ) then
233
cb0p = DCMPLX( DBLE(log(xp) - 2), DBLE(-pi) )
239
* #] both masses equal to zero:
240
* #[ one mass equal to zero:
243
* special case xp = 0
245
if ( xp .eq. 0 ) then
249
* special case xp = xm
251
elseif ( dmp.eq.0 ) then
259
if ( abs(s1) .lt. xloss ) then
267
+ cb0p = cb0p - DCMPLX(DBLE(x0),DBLE(-(dmp/xp)*pi))
268
if ( lwarn .and. absc(cb0p) .lt. xloss*x2 )
269
+ call ffwarn(14,ier,absc(cb0p),x2)
271
* #] one mass equal to zero:
272
* #[ both masses equal:
275
* Both masses are equal. Not only this speeds up things, some
276
* cancellations have to be avoided as well.
278
* first a special case
280
if ( abs(xp) .lt. 8*xloss*xm ) then
281
* -#[ taylor expansion:
283
* a Taylor expansion seems appropriate as the result will go
284
* as k^2 but seems to go as 1/k !!
286
*--#[ data and bounds:
287
if ( initeq .eq. 0 ) then
291
xpneq(i) = - xpneq(i-1)*DBLE(i-1)/DBLE(2*(2*i+1))
294
if (xprceq .ne. precx ) then
296
* calculate the boundaries for the number of terms to be
297
* included in the taylorexpansion
300
bdeq01 = ffbnd(1,1,xpneq)
301
bdeq05 = ffbnd(1,5,xpneq)
302
bdeq11 = ffbnd(1,11,xpneq)
303
bdeq17 = ffbnd(1,17,xpneq)
304
bdeq25 = ffbnd(1,25,xpneq)
306
*--#] data and bounds:
309
if ( lwarn .and. ax .gt. bdeq25 ) then
310
call ffwarn(15,ier,precx,abs(xpneq(25))*ax**25)
312
if ( ax .gt. bdeq17 ) then
313
som = x*(xpneq(18) + x*(xpneq(19) + x*(xpneq(20) +
314
+ x*(xpneq(21) + x*(xpneq(22) + x*(xpneq(23) +
315
+ x*(xpneq(24) + x*(xpneq(25) ))))))))
319
if ( ax .gt. bdeq11 ) then
320
som = x*(xpneq(12) + x*(xpneq(13) + x*(xpneq(14) +
321
+ x*(xpneq(15) + x*(xpneq(16) + x*(xpneq(17) + som ))))
324
if ( ax .gt. bdeq05 ) then
325
som = x*(xpneq(6) + x*(xpneq(7) + x*(xpneq(8) + x*(
326
+ xpneq(9) + x*(xpneq(10) + x*(xpneq(11) + som ))))))
328
if ( ax .gt. bdeq01 ) then
329
som = x*(xpneq(2) + x*(xpneq(3) + x*(xpneq(4) + x*(
330
+ xpneq(5) + som ))))
332
cb0p = x*(xpneq(1)+som)
334
print *,'ffxb0q: m1 = m2, Taylor expansion in ',x
337
* -#] taylor expansion:
343
call ffxlmb(xlam,-xp,-xm,-xm,dmp,dmp,x0,ier)
344
if ( xlam .ge. 0 ) then
349
if ( abs(s2) .gt. xloss*slam ) then
357
if ( ax .lt. xalogm ) then
358
if ( lwarn ) call ffwarn(16,ier,ax,xalogm)
360
elseif( ax-1 .lt. .1 .and. s2 .gt. 0 ) then
361
* In this case a quicker and more accurate way is to
362
* calculate log(1-x).
364
* the following line is superfluous.
365
if ( lwarn .and. abs(s2) .lt. xloss*slam )
366
+ call ffwarn(17,ier,s2,slam)
367
s = -slam/xp*dfflo1(s2/(2*xm),ier)
369
* finally the normal case
371
if ( jsign .eq. -1 ) s = -s
373
if ( xp .gt. 2*xm ) then
374
* in this case ( xlam>0, so xp>(2*m)^2) ) there also
375
* is an imaginary part
381
* the root is complex (k^2 between 0 and (2*m1)^2)
383
s = 2*slam/xp*atan2(xp,slam)
386
if (lwrite) print *,'s = ',s
388
if ( lwarn .and. abs(xx).lt.xloss*2 ) call ffwarn(18,ier,xx,x2)
389
cb0p = DCMPLX(DBLE(xx),DBLE(y))
393
* #] both masses equal:
394
* #[ unequal nonzero masses:
395
* -#[ get log(xm2/xm1):
398
if ( 1 .lt. xalogm*x ) then
401
elseif ( abs(x-1) .lt. xloss ) then
402
xlogmm = dfflo1(dm1m2/xm1,ier)
406
* -#] get log(xm2/xm1):
409
* first a special case
411
if ( xp .eq. 0 ) then
412
s2 = ((xm2+xm1) / dm1m2)*xlogmm
414
* save the factor 1/2 for the end
415
if (lwrite) print *,'s = ',s/2
416
* save the factor 1/2 for the end
417
if ( abs(s) .lt. xloss*2 ) then
418
* Taylor expansions: choose which one
421
if ( ax .lt. .15 .or. precx .gt. 1.E-8 .and. ax
424
* This is the simple Taylor expansion 'n1'
426
*--#[ data and bounds:
427
* get the coefficients of the taylor expansion
428
if ( initn1 .eq. 0 ) then
431
410 xpnn1(i) = DBLE(i)/DBLE((i+1)*(i+2))
433
* determine the boundaries for 1,5,10,15 terms
434
if ( xprcn1 .ne. precx ) then
436
bdn101 = ffbnd(1,1,xpnn1)
437
bdn105 = ffbnd(1,5,xpnn1)
438
bdn110 = ffbnd(1,10,xpnn1)
439
bdn115 = ffbnd(1,15,xpnn1)
440
bdn120 = ffbnd(1,20,xpnn1)
442
*--#] data and bounds:
444
if ( lwarn .and. ax .gt. bdn120 )
445
+ call ffwarn(19,ier,precx,abs(xpnn1(20))*ax**19)
446
if ( ax .gt. bdn115 ) then
447
s = x*(xpnn1(16) + x*(xpnn1(17) + x*(xpnn1(18) +
448
+ x*(xpnn1(19) + x*(xpnn1(20)) ))))
452
if ( ax .gt. bdn110 ) then
453
s = x*(xpnn1(11) + x*(xpnn1(12) + x*(xpnn1(13) +
454
+ x*(xpnn1(14) + x*(xpnn1(15)) + s))))
456
if ( ax .gt. bdn105 ) then
457
s = x*(xpnn1(6) + x*(xpnn1(7) + x*(xpnn1(8) + x*
458
+ (xpnn1(9) + x*(xpnn1(10) + s)))))
460
if ( ax .gt. bdn101 ) then
461
s = x*(xpnn1(2) + x*(xpnn1(3) + x*(xpnn1(4) + x*
464
s = x*x*(xpnn1(1) + s)
466
print *,'ffxb0q: xp = 0, simple Taylor exp'
468
print *,' gives s ',s/2
472
* This is the more complicated Taylor expansion 'fc'
475
* determine the boundaries for 1,5,10,15 terms for
476
* the exponential taylor expansion, assuming it
479
if ( xprnn2 .ne. precx ) then
481
bdn201 = ffbnd(4,1,xinfac)
482
bdn205 = ffbnd(4,5,xinfac)
483
bdn210 = ffbnd(4,10,xinfac)
484
bdn215 = ffbnd(4,15,xinfac)
485
bdn220 = ffbnd(4,20,xinfac)
491
if ( lwarn .and. ay .gt. bdn220 )
492
+ call ffwarn(20,ier,precx,xinfac(23)*ay**23)
493
if ( ay .gt. bdn220 ) then
494
s = y*(xinfac(19) + y*(xinfac(20) + y*(xinfac(
495
+ 21) + y*(xinfac(22) + y*(xinfac(
500
if ( ay .gt. bdn215 ) then
501
s = y*(xinfac(14) + y*(xinfac(15) + y*(xinfac(
502
+ 16) + y*(xinfac(17) + y*(xinfac(
505
if ( ay .gt. bdn210 ) then
506
s = y*(xinfac(9) + y*(xinfac(10) + y*(xinfac(11)
507
+ + y*(xinfac(12) + y*(xinfac(13) + s)))))
509
if ( ay .gt. bdn205 ) then
510
s = y*(xinfac(5) + y*(xinfac(6) + y*(xinfac(7) +
511
+ y*(xinfac(8) + s))))
513
s = (1-x)*y**4*(xinfac(4)+s)
514
s = x*y**2*(1+y)/12 - s
515
s = - 2*dfflo1(s,ier)/y
517
print *,'ffxb0q: xp = 0, other Taylor expansion'
529
* proceeding with the normal case
531
call ffxlmb(xlam,-xp,-xm2,-xm1,dm2p,dm1p,dm1m2,ier)
532
if ( xlam .gt. 0 ) then
533
* cases k^2 < -(m2+m1)^2 or k^2 > -(m2-m1)^2:
535
* first try the normal way
540
if ( abs(s2) .gt. xloss*slam ) then
547
s2 = s2**2/(4*xm1*xm2)
548
if ( abs(s2) .lt. xalogm ) then
551
elseif ( abs(s2-1) .lt. xloss ) then
552
if ( jsign.eq.1 ) then
553
if ( lwrite ) print *,'s2 was ',-slam/(2*xp)*log(s2)
554
s2 = -slam*(s2a+slam)/(2*xm1*xm2)
555
s2 = -slam/(2*xp)*dfflo1(s2,ier)
557
if ( lwrite ) print *,'s2 was ',+slam/(2*xp)*log(s2)
558
s2 = +slam*(s2a-slam)/(2*xm1*xm2)
559
s2 = +slam/(2*xp)*dfflo1(s2,ier)
561
if ( lwrite ) print *,'s2 is ',s2,jsign
563
s2 = -slam/(2*xp)*log(s2)
564
if ( jsign .eq. -1 ) s2 = -s2
566
s1 = -dm1m2*xlogmm/(2*xp)
569
print *,'ffxb0q: lam>0, first try, xx = ',xx,s1,s2,-2
572
if ( abs(xx) .lt. xloss*max(abs(s1),abs(s2)) ) then
574
* this is unacceptable, try a better solution
576
if (lwrite) print *,'s1 = ',-s1a/(2*xp),slam/(2*xp)
577
if ( abs(s1a) .gt. xloss*slam ) then
578
* (strangely) this works
581
* by division a more accurate form can be found
582
s1 = ( -xp/2 + xm1 + xm2 ) / ( slam - dm1m2 )
583
if (lwrite) print *,'s1+= ',s1
586
if ( abs(xp) .lt. xm2 ) then
592
if (lwrite) print *,'s2 = ',s2/(2*xm2),slam/(2*xm2)
593
if ( abs(s2) .gt. xloss*slam ) then
594
* at least reasonable
598
s2 = (2*xp) / (s2a+slam)
599
if (lwrite) print *,'s2+= ',s2
601
if ( abs(s2) .lt. .1 ) then
602
* choose a quick way to get the logarithm
611
print *,'ffxb0q: lam>0, second try, xx = ',xx
612
alpha = slam/(slam-dm1m2)
613
alph1 = -dm1m2/(slam-dm1m2)
614
print *,' alpha = ',alpha
615
print *,' s1 = ',s1,',- 2alph1 = ',s1-2*alph1
616
print *,' s2 = ',s2,',- 2alpha = ',s2-2*alpha
619
if ( abs(xx) .lt. xloss**2*max(abs(s1),abs(s2)) ) then
621
* (we accept two times xloss because that's the same
623
* A Taylor expansion might work. We expand
624
* inside the logs. Only do the necessary work.
626
alpha = slam/(slam-dm1m2)
627
alph1 = -dm1m2/(slam-dm1m2)
632
if ( abs(s1p) .lt. xloss*abs(s1) ) then
634
* determine the boundaries for 1,5,10,15 terms
635
if ( xprcn3 .ne. precx ) then
637
bdn301 = ffbnd(3,1,xinfac)
638
bdn305 = ffbnd(3,5,xinfac)
639
bdn310 = ffbnd(3,10,xinfac)
640
bdn315 = ffbnd(3,15,xinfac)
643
xnoe = -xp + 2*xm1 + 2*xm2
646
if ( lwarn .and. ax .gt. bdn315 ) then
647
call ffwarn(21,ier,precx,xinfac(17)*ax**14)
649
if ( ax .gt. bdn310 ) then
650
s1a = x*(xinfac(13) + x*(xinfac(14) + x*(
651
+ xinfac(15) + x*(xinfac(16) + x*(
656
if ( ax .gt. bdn305 ) then
657
s1a = x*(xinfac(8) + x*(xinfac(9) + x*(
658
+ xinfac(10) + x*(xinfac(11) + x*(
659
+ xinfac(12) + s1a)))))
661
if ( ax .gt. bdn301 ) then
662
s1a = x*(xinfac(4) + x*(xinfac(5) + x*(
663
+ xinfac(6) + x*(xinfac(7) + s1a))))
665
s1a = x**3 *(xinfac(3) + s1a) *xm2/xm1
666
s1b = dm1m2*(4*dm1m2**2 - xp*(4*xm1-xp))/
669
if ( lwarn .and. abs(s1p).lt.xloss*abs(s1a) )
670
+ call ffwarn(22,ier,s1p,s1a)
671
s1p = xnoe*dfflo1(s1p,ier)/(slam - dm1m2)/2
673
print *,'ffxb0q: Taylor exp. of s1-2(1-a)'
675
print *,' gives s1p = ',s1p
681
490 s2p = s2 - 2*alpha
682
if ( abs(s2p) .lt. xloss*abs(s2) ) then
684
* determine the boundaries for 1,5,10,15 terms
685
if ( xprcn5 .ne. precx ) then
687
bdn501 = ffbnd(4,1,xinfac)
688
bdn505 = ffbnd(4,5,xinfac)
689
bdn510 = ffbnd(4,10,xinfac)
690
bdn515 = ffbnd(4,15,xinfac)
696
if ( ax .gt. bdn515 ) then
697
* do not do the Taylor expansion
698
if ( lwarn ) call ffwarn(23,ier,s2p,s2)
701
if ( ax .gt. bdn510 ) then
702
s2a = x*(xinfac(14) + x*(xinfac(15) + x*(
703
+ xinfac(16) + x*(xinfac(17) + x*(
708
if ( ax .gt. bdn505 ) then
709
s2a = x*(xinfac(9) + x*(xinfac(10) + x*(
710
+ xinfac(11) + x*(xinfac(12) + x*(
711
+ xinfac(13) + s2a)))))
713
if ( ax .gt. bdn501 ) then
714
s2a = x*(xinfac(5) + x*(xinfac(6) + x*(
715
+ xinfac(7) + x*(xinfac(8) + s2a))))
717
s2a = x**4*(xinfac(4)+s2a)*(1-2*xp/(xnoe+xp))
718
s2b = -2*xp**3 *(-2*xp - xnoe)/(3*(xnoe+xp)*
721
if ( lwarn .and. abs(s2p).lt.xloss*abs(s2a) )
722
+ call ffwarn(24,ier,s2p,s2a)
723
s2p = -slam/xp*dfflo1(s2p,ier)
725
print *,'ffxb0q: Taylor expansion of s2-2a'
727
print *,' gives s2p = ',s2p
734
if ( lwarn .and. abs(xx) .lt. xloss*abs(s1p) ) then
735
call ffwarn(25,ier,xx,s1p)
740
if ( xp .gt. xm1+xm2 ) then
741
*--#[ imaginary part:
742
* in this case ( xlam>0, so xp>(m1+m2)^2) ) there also
743
* is an imaginary part
747
*--#] imaginary part:
750
* the root is complex (k^2 between -(m1+m2)^2 and -(m2-m1)^2)
754
s1 = -(dm1m2/(2*xp))*xlogmm
755
s2 = (slam/xp)*atan2(slam,xnoe)
758
print *,'ffxb0q: lam<0, first try, xx = ',xx
759
alpha = -xlam/(2*xp*xnoe)
760
alph1 = -(xp**2-dm1m2**2)/(2*xp*xnoe)
761
print *,' alpha = ',alpha
762
print *,' s1 = ',s1,' - 2alph1 = ',s1-2*alph1
763
print *,' s2 = ',s2,' - 2alpha = ',s2-2*alpha
766
if ( abs(xx) .lt. xloss**2*max(abs(s1),abs(s2)) ) then
768
* Again two times xloss as we'll accept that in the next
771
xtel = dm1m2**2 - xp**2
772
alpha = -xlam/(2*xp*xnoe)
773
alph1 = xtel/(2*xp*xnoe)
775
* try a taylor expansion on the terms. First s1:
778
if ( abs(s1p) .lt. xloss*abs(s1) ) then
780
* determine the boundaries for 1,5,10,15 terms
781
if ( xprcn3 .ne. precx ) then
783
bdn301 = ffbnd(3,1,xinfac)
784
bdn305 = ffbnd(3,5,xinfac)
785
bdn310 = ffbnd(3,10,xinfac)
786
bdn315 = ffbnd(3,15,xinfac)
789
x = 2*xtel/(dm1m2*xnoe)
791
if ( ax .gt. bdn315 ) then
792
* do not do the Taylor expansion
793
if ( lwarn ) call ffwarn(21,ier,s1p,s1)
796
if ( ax .gt. bdn310 ) then
797
s1a = x*(xinfac(13) + x*(xinfac(14) + x*(
798
+ xinfac(15) + x*(xinfac(16) + x*(
803
if ( ax .gt. bdn305 ) then
804
s1a = x*(xinfac(8) + x*(xinfac(9) + x*(
805
+ xinfac(10) + x*(xinfac(11) + x*(
806
+ xinfac(12) + s1a)))))
808
if ( ax .gt. bdn301 ) then
809
s1a = x*(xinfac(4) + x*(xinfac(5) + x*(
810
+ xinfac(6) + x*(xinfac(7) + s1a))))
812
s1a = x**3 *(xinfac(3) + s1a) *xm2/xm1
813
s1b = (dm1m2**3*(dm1m2**2-2*xp*xm1) + xp**2*(4*
814
+ dm1m2*xm1**2-dm1m2**2*(dm1m2+2*xm1))-2*xm2*
815
+ xp**3*(dm1m2+xp))/(xm1*dm1m2**2*xnoe**2)
817
if ( lwarn .and. abs(s1p) .lt. xloss*abs(s1a) )
818
+ call ffwarn(22,ier,s1p,s1a)
819
s1p = -dm1m2*dfflo1(s1p,ier)/(2*xp)
821
print *,'ffxb0q: Taylor expansion of s1-2(1-a)'
823
print *,' gives s1p = ',s1p
831
if ( abs(s2p) .lt. xloss*abs(s2) ) then
833
* determine the boundaries for 1,5,10,15 terms
834
if ( xprcn3 .ne. precx ) then
836
bdn301 = ffbnd(3,1,xinfac)
837
bdn305 = ffbnd(3,5,xinfac)
838
bdn310 = ffbnd(3,10,xinfac)
839
bdn315 = ffbnd(3,15,xinfac)
842
cx = DCMPLX(DBLE(x0),DBLE(-slam/xnoe))
844
if ( ax .gt. bdn315 ) then
845
if ( lwarn ) call ffwarn(23,ier,s2p,s2)
848
if ( ax .gt. bdn310 ) then
849
cs2a = cx*(DBLE(xinfac(13)) + cx*(DBLE(xinfac(14
850
+ )) + cx*(DBLE(xinfac(15)) + cx*(DBLE(xinfac(16
851
+ )) + cx*(DBLE(xinfac(17)))))))
855
if ( ax .gt. bdn305 ) then
856
cs2a = cx*(DBLE(xinfac(8)) + cx*(DBLE(xinfac(9))
857
+ + cx*(DBLE(xinfac(10)) + cx*(DBLE(xinfac(11))
858
+ + cx*(DBLE(xinfac(12)) + cs2a)))))
860
if ( ax .gt. bdn301 ) then
861
cs2a = cx*(DBLE(xinfac(4)) + cx*(DBLE(xinfac(5))
862
+ + cx*(DBLE(xinfac(6)) + cx*(DBLE(xinfac(7))
865
cs2a = cx**3*(DBLE(xinfac(3))+cs2a)*
866
+ DCMPLX(DBLE(xnoe),DBLE(slam))
867
cs2b = DCMPLX(DBLE(xnoe-xlam/xnoe/2),
868
+ -DBLE(slam**3/xnoe**2/2))
870
if ( lwarn .and. absc(cs2p) .lt. xloss*absc(cs2a) )
871
+ call ffwarn(24,ier,absc(cs2p),absc(cs2b))
872
s2p = slam*atan2(DIMAG(cs2p),DBLE(cs2p))/xp
874
print *,'ffxb0q: Taylor expansion of s2-2a'
875
print *,' in x = ',cx
876
print *,' gives s2p = ',s2p
881
if ( lwarn .and. abs(xx) .lt. xloss*abs(s1p) ) then
882
call ffwarn(25,ier,xx,s1p)
888
cb0p = DCMPLX(DBLE(xx),DBLE(y))
891
* #] unequal nonzero masses:
895
print *,'cb0p = ',cb0p,ier
901
subroutine ffxlmb(xlambd,a1,a2,a3,a12,a13,a23,ier)
902
***#[*comment:***********************************************************
903
* calculate in a numerically stable way *
904
* lambda(a1,a2,a3) = *
905
* a1**2 + a2**2 + a3**2 - 2*a2*a3 - 2*a3*a1 - 2*a1*a2 *
906
* aij = ai - aj are required for greater accuracy at times *
907
* ier is the usual error flag. *
908
***#]*comment:***********************************************************
915
DOUBLE PRECISION xlambd,a1,a2,a3,a12,a13,a23
919
DOUBLE PRECISION aa1,aa2,aa3,aa12,aa13,aa23,
935
xcheck = a1 - a2 - a12
936
if ( xloss*abs(xcheck) .gt. precx*max(aa1,aa2,aa12) )
937
+ print *,'ffxlmb: input not OK, a12 /= a1 - a2',a12,a1,
939
xcheck = a1 - a3 - a13
940
if ( xloss*abs(xcheck) .gt. precx*max(aa1,aa3,aa13) )
941
+ print *,'ffxlmb: input not OK, a13 /= a1 - a3',a13,a3,
943
xcheck = a2 - a3 - a23
944
if ( xloss*abs(xcheck) .gt. precx*max(aa2,aa3,aa23) )
945
+ print *,'ffxlmb: input not OK, a23 /= a2 - a3',a23,a2,
949
* first see if there are input parameters with opposite sign:
951
if ( a1 .lt. 0 .and. a2 .gt. 0 .or.
952
+ a1 .gt. 0 .and. a2 .lt. 0 ) then
954
elseif ( a1 .lt. 0 .and. a3 .gt. 0 .or.
955
+ a1 .gt. 0 .and. a3 .lt. 0 ) then
958
* all have the same sign, choose the smallest 4*ai*aj term
960
elseif ( aa1 .gt. aa2 .and. aa1 .gt. aa3 ) then
962
elseif ( aa2 .gt. aa3 ) then
968
if ( aa1 .gt. aa2 ) then
976
if ( aa1 .gt. aa3 ) then
984
if ( aa2 .gt. aa3 ) then
993
if ( lwarn .and. abs(xlambd) .lt. xloss*asq )
994
+ call ffwarn(69,ier,xlambd,asq)
999
subroutine ffclmb(clambd,cc1,cc2,cc3,cc12,cc13,cc23,ier)
1000
***#[*comment:***********************************************************
1001
* calculate in cc numerically stable way *
1002
* lambda(cc1,cc2,cc3) = *
1003
* cc1**2 + cc2**2 + cc3**2 - 2*cc2*cc3 - 2*cc3*cc1 - 2*cc1*cc2 *
1004
* cij = ci - cj are required for greater accuracy at times *
1005
* ier is the usual error flag. *
1006
***#]*comment:***********************************************************
1013
DOUBLE COMPLEX clambd,cc1,cc2,cc3,cc12,cc13,cc23
1017
DOUBLE PRECISION aa1,aa2,aa3,aa12,aa13,aa23,absc
1018
DOUBLE COMPLEX check,cc,cff,csq,c
1024
* statement function
1026
absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
1028
* #[ calculations (rather old style ...):
1037
check = cc1 - cc2 - cc12
1038
if ( xloss*absc(check) .gt. precc*max(aa1,aa2,aa12) ) then
1039
print *,'ffclmb: input not OK, cc12 /= cc1 - cc2',check
1041
check = cc1 - cc3 - cc13
1042
if ( xloss*absc(check) .gt. precc*max(aa1,aa3,aa13) ) then
1043
print *,'ffclmb: input not OK, cc13 /= cc1 - cc3',check
1045
check = cc2 - cc3 - cc23
1046
if ( xloss*absc(check) .gt. precc*max(aa2,aa3,aa23) ) then
1047
print *,'ffclmb: input not OK, cc23 /= cc2 - cc3',check
1051
* first see if there are input parameters with opposite sign:
1053
if ( DBLE(cc1) .lt. 0 .and. DBLE(cc2) .gt. 0 .or.
1054
+ DBLE(cc1) .gt. 0 .and. DBLE(cc2) .lt. 0 ) then
1056
elseif ( DBLE(cc1) .lt. 0 .and. DBLE(cc3) .gt. 0 .or.
1057
+ DBLE(cc1) .gt. 0 .and. DBLE(cc3) .lt. 0 ) then
1060
* all have the same sign, choose the smallest 4*ci*cj term
1062
elseif ( aa1 .gt. aa2 .and. aa1 .gt. aa3 ) then
1064
elseif ( aa2 .gt. aa3 ) then
1070
if ( aa1 .gt. aa2 ) then
1078
if ( aa1 .gt. aa3 ) then
1086
if ( aa2 .gt. aa3 ) then
1095
if ( lwarn .and. absc(clambd) .lt. xloss*absc(csq) )
1096
+ call ffwarn(68,ier,absc(clambd),absc(csq))
1097
* #] calculations (rather old style ...):
1101
subroutine ffdot2(piDpj,xp,xma,xmb,dmap,dmbp,dmamb,ier)
1102
***#[*comment:***********************************************************
1104
* Store the 3 dotproducts in the common block ffdot. *
1106
* Input: see ffxb0p *
1108
***#]*comment:***********************************************************
1115
DOUBLE PRECISION piDpj(3,3),xp,xma,xmb,dmap,dmbp,dmamb
1125
* statement function
1127
* absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
1134
if ( abs(dmap) .lt. abs(dmbp) ) then
1135
piDpj(1,2) = (dmap + xmb)/2
1137
piDpj(1,2) = (dmbp + xma)/2
1139
piDpj(2,1) = piDpj(1,2)
1140
if ( lwarn .and. abs(piDpj(1,2)) .lt. xloss*min(xma,xmb)/2
1142
call ffwarn(24,ier1,piDpj(1,2),min(xma,xmb)/2)
1144
if ( abs(dmamb) .lt. abs(dmbp) ) then
1145
piDpj(1,3) = (-dmamb - xp)/2
1147
piDpj(1,3) = (dmbp - xma)/2
1149
piDpj(3,1) = piDpj(1,3)
1150
if ( lwarn .and. abs(piDpj(1,3)) .lt. xloss*min(xma,
1153
call ffwarn(25,ier0,piDpj(1,3),min(xma,abs(xp))/2)
1154
ier1 = max(ier0,ier1)
1156
if ( abs(dmamb) .lt. abs(dmap) ) then
1157
piDpj(2,3) = (-dmamb + xp)/2
1159
piDpj(2,3) = (-dmap + xmb)/2
1161
piDpj(3,2) = piDpj(2,3)
1162
if ( lwarn .and. abs(piDpj(2,3)) .lt. xloss*min(xmb,
1165
call ffwarn(25,ier0,piDpj(2,3),min(xmb,abs(xp))/2)
1166
ier1 = max(ier0,ier1)
added
removed
vendor/CutTools/src/qcdloop/ffxb0.f
