~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to vendor/CutTools/src/qcdloop/ffxc0i.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
*--#[ log:
 
2
*       $Id: ffxc0i.f,v 1.3 1996/06/03 12:11:43 gj Exp $
 
3
*       $Log: ffxc0i.f,v $
 
4
c Revision 1.3  1996/06/03  12:11:43  gj
 
5
c Added an error message for ffxc0j with zero masses, which is ill-defined.
 
6
c
 
7
c Revision 1.2  1995/12/01  15:04:40  gj
 
8
c Fixed a ridiculous bug: wrong sign for p4^2=0, m2<m1.
 
9
c
 
10
*--#] log:
 
11
*###[ ffxc0i:
 
12
        subroutine ffxc0i(cc0,xpi,dpipj,ier)
 
13
***#[*comment:***********************************************************
 
14
*                                                                       *
 
15
*       Calculates the infrared finite part of a infrared divergent     *
 
16
*       threepoint function with the factor ipi^2.  The cutoff          *
 
17
*       parameter is assumed to be in a common block /ffcut/. (ugly)    *
 
18
*                                                                       *
 
19
*       Input:  xpi(6)          (real)  pi.pi (B&D)                     *
 
20
*               dpipj(6,6)      (real)  xpi(i)-xpi(j)                   *
 
21
*               /ffcut/delta    (real)  cutoff (either foton mass**2 or *
 
22
*                                               radiation limit).       *
 
23
*       Output: cc0     (complex)       C0, the threepoint function.    *
 
24
*               ier     (integer)       usual error flag                *
 
25
*       Calls:                                                          *
 
26
*                                                                       *
 
27
***#]*comment:***********************************************************
 
28
*  #[ declarations:
 
29
        implicit none
 
30
*
 
31
*       arguments
 
32
*
 
33
        integer ier
 
34
        DOUBLE COMPLEX cc0
 
35
        DOUBLE PRECISION xpi(6),dpipj(6,6)
 
36
*
 
37
*       local variables
 
38
*
 
39
        integer init,ipi12,i,ilogi(3),irota,n
 
40
        integer j,inew(6,6)
 
41
        DOUBLE COMPLEX cs(15),csum,c,clogi(3)
 
42
        DOUBLE PRECISION xqi(6),dqiqj(6,6),qiDqj(6,6),sdel2,xmax,absc,
 
43
     +          dum66(6,6),del2
 
44
        save init,inew,ilogi
 
45
*
 
46
*       common blocks etc
 
47
*
 
48
        include 'ff.h'
 
49
        DOUBLE PRECISION delta
 
50
        common /ffcut/ delta
 
51
*
 
52
*       data
 
53
*
 
54
        data init /0/
 
55
        data inew /1,2,3,4,5,6,
 
56
     +             2,3,1,5,6,4,
 
57
     +             3,1,2,6,4,5,
 
58
     +             1,3,2,6,5,4,
 
59
     +             3,2,1,5,4,6,
 
60
     +             2,1,3,4,6,5/
 
61
        data ilogi /3*0/
 
62
*
 
63
*       statement function
 
64
*
 
65
        absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
66
*
 
67
*       initialisations
 
68
*
 
69
        do 1 i=1,15
 
70
            cs(i) = 0
 
71
    1   continue
 
72
        ipi12 = 0
 
73
*  #] declarations:
 
74
*  #[ check input:
 
75
        if ( init .eq. 0 .and. .not.lsmug ) then
 
76
            init = 1
 
77
            print *,'ffxc0i: infra-red divergent threepoint function, ',
 
78
     +          'working with a cutoff ',delta
 
79
        endif
 
80
        if ( .not.lsmug .and. delta .eq. 0 ) then
 
81
            call fferr(59,ier)
 
82
            return
 
83
        endif
 
84
        if ( lwrite ) then
 
85
*           print input
 
86
            print *,'ffxc0i: infrared divergent threepoint function'
 
87
            if ( .not.lsmug ) then
 
88
                print *,'  cutoff parameter:',delta
 
89
            endif
 
90
        endif
 
91
*  #] check input:
 
92
*  #[ groundwork:
 
93
*
 
94
*       rotate to xpi(3)=0, xpi(1)=xpi(6), xpi(2)=xpi(5)
 
95
*
 
96
        call ffrot3(irota,xqi,dqiqj,qiDqj,xpi,dpipj,dum66,6,3,3,ier)
 
97
*
 
98
*       get some dotproducts
 
99
*
 
100
        if ( ldot ) then
 
101
            call ffdot3(qiDqj,xqi,dqiqj,6,ier)
 
102
            do 5 i=1,6
 
103
                do 4 j=1,6
 
104
                    fpij3(j,i) = qiDqj(inew(i,irota),inew(j,irota))
 
105
    4           continue
 
106
    5       continue
 
107
        else
 
108
        if ( abs(xqi(4)) .lt. xqi(1) ) then
 
109
            qiDqj(4,1) = dqiqj(2,1) - xqi(4)
 
110
            xmax = abs(xqi(4))
 
111
        else
 
112
            qiDqj(4,1) = dqiqj(2,4) - xqi(1)
 
113
            xmax = xqi(1)
 
114
        endif
 
115
        if ( lwarn .and. abs(qiDqj(4,1)) .lt. xloss*xmax )
 
116
     +          call ffwarn(156,ier,qiDqj(4,1),xmax)
 
117
        qiDqj(4,1) = qiDqj(4,1)/2
 
118
        qiDqj(1,4) = qiDqj(4,1)
 
119
 
 
120
        if ( abs(xqi(4)) .lt. xqi(2) ) then
 
121
            qiDqj(4,2) = dqiqj(2,1) + xqi(4)
 
122
            xmax = abs(xqi(4))
 
123
        else
 
124
            qiDqj(4,2) = xqi(2) - dqiqj(1,4)
 
125
            xmax = xqi(2)
 
126
        endif
 
127
        if ( lwarn .and. abs(qiDqj(4,2)) .lt. xloss*xmax )
 
128
     +          call ffwarn(156,ier,qiDqj(4,2),xmax)
 
129
        qiDqj(4,2) = qiDqj(4,2)/2
 
130
        qiDqj(2,4) = qiDqj(4,2)
 
131
 
 
132
        if ( (xqi(1)) .lt. (xqi(2)) ) then
 
133
            qiDqj(1,2) = xqi(1) + dqiqj(2,4)
 
134
            xmax = xqi(1)
 
135
        else
 
136
            qiDqj(1,2) = xqi(2) + dqiqj(1,4)
 
137
            xmax = xqi(2)
 
138
        endif
 
139
        if ( lwarn .and. abs(qiDqj(1,2)) .lt. xloss*xmax )
 
140
     +          call ffwarn(156,ier,qiDqj(1,2),xmax)
 
141
        qiDqj(1,2) = qiDqj(1,2)/2
 
142
        qiDqj(2,1) = qiDqj(1,2)
 
143
 
 
144
        qiDqj(1,1) = xqi(1)
 
145
        qiDqj(2,2) = xqi(2)
 
146
        qiDqj(4,4) = xqi(4)
 
147
        endif
 
148
*  #] groundwork:
 
149
*  #[ calculations:
 
150
*
 
151
        call ffdel2(del2,qiDqj,6,1,2,4,1,ier)
 
152
        if ( ldot ) fdel2 = del2
 
153
*
 
154
*       the case del2=0 is hopeless - this is really a two-point function
 
155
*
 
156
        if ( del2 .eq. 0 ) then
 
157
            call fferr(58,ier)
 
158
            return
 
159
        endif
 
160
*
 
161
*       we cannot yet handle the complex case
 
162
*
 
163
        if ( del2 .gt. 0 ) then
 
164
            call fferr(83,ier)
 
165
            return
 
166
        endif
 
167
*
 
168
        sdel2 = isgnal*sqrt(-del2)
 
169
*
 
170
        call ffxc0j(cs,ipi12,sdel2,clogi,ilogi,xqi,dqiqj,qiDqj,
 
171
     +          delta,3,ier)
 
172
*  #] calculations:
 
173
*  #[ sum:
 
174
*
 
175
*       Sum
 
176
*
 
177
        xmax = 0
 
178
        csum = 0
 
179
        if ( .not.lsmug ) then
 
180
            n = 10
 
181
        else
 
182
            n = 15
 
183
        endif
 
184
        do 10 i=1,n
 
185
            csum = csum + cs(i)
 
186
            xmax = max(xmax,absc(csum))
 
187
   10   continue
 
188
        csum = csum + ipi12*DBLE(pi12)
 
189
        if ( lwarn .and. 2*absc(csum) .lt. xloss*xmax ) then
 
190
            call ffwarn(157,ier,absc(csum),xmax)
 
191
        endif
 
192
        cc0 = -csum*DBLE(1/(2*sdel2))
 
193
*  #] sum:
 
194
*  #[ debug:
 
195
  900   continue
 
196
        if (lwrite) then
 
197
            print '(a)','cs(i) = '
 
198
            print '(i3,2g20.10,1x)',(i,cs(i),i=1,n)
 
199
            print '(a3,2g20.10,1x)','pi ',ipi12*pi12
 
200
            print '(a)','+-----------'
 
201
            print '(a3,2g20.10,1x)','som  :',csum
 
202
            print '(a)',' '
 
203
            print *,'cc0  :',cc0,ier
 
204
        endif
 
205
*  #] debug:
 
206
*###] ffxc0i:
 
207
        end
 
208
*###[ ffxc0j:
 
209
        subroutine ffxc0j(cs,ipi12,sdel2i,clogi,ilogi,
 
210
     +                                  xpi,dpipj,piDpj,delta,npoin,ier)
 
211
***#[*comment:***********************************************************
 
212
*                                                                       *
 
213
*       Calculates the infrared finite part of a infrared divergent     *
 
214
*       threepoint function with the factor ipi^2.                      *
 
215
*                                                                       *
 
216
***#]*comment:***********************************************************
 
217
*  #[ declarations:
 
218
        implicit none
 
219
*
 
220
*       arguments
 
221
*
 
222
        integer ipi12,ilogi(3),npoin,ier
 
223
        DOUBLE COMPLEX cs(15),clogi(3)
 
224
        DOUBLE PRECISION xpi(6),dpipj(6,6),piDpj(6,6),delta,sdel2i
 
225
*
 
226
*       local variables
 
227
*
 
228
        integer i,ieps,ieps1,n,ier0
 
229
        DOUBLE COMPLEX clog1,clog2,cdum(2),cel3,cdyzm,cdyzp,cli,chulp,
 
230
     +          carg1,carg2,chulp1
 
231
        DOUBLE COMPLEX zfflog,zxfflg,cc
 
232
        DOUBLE PRECISION del2,zm,zp,zm1,zp1,sdel2,hulp,xheck,dum(3),
 
233
     +          dfflo1,dyzp,dyzm,wm,wp,absc,arg1,arg2,del3
 
234
*
 
235
*       common blocks etc
 
236
*
 
237
        include 'ff.h'
 
238
*
 
239
*       statement function
 
240
*
 
241
        absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
 
242
*
 
243
*  #] declarations:
 
244
*  #[ check input:
 
245
        if ( ltest ) then
 
246
            call ffxhck(xpi,dpipj,6,ier)
 
247
        endif
 
248
        if ( lwrite ) then
 
249
            print '(a)','  ##[ ffxc0j:'
 
250
            print *,'ffxc0j: input: '
 
251
            print *,'xpi   = ',xpi
 
252
            if ( .not.lsmug ) then
 
253
                print *,'delta = ',delta
 
254
            else
 
255
                print *,'cmipj(2,2) = ',cmipj(2,2)
 
256
                print *,'cmipj(1,3) = ',cmipj(1,3)
 
257
            endif
 
258
        endif
 
259
*  #] check input:
 
260
*  #[ get determinants, roots, ieps:
 
261
*
 
262
        if ( lsmug ) then
 
263
            del3 = (- DBLE(xpi(1))*DBLE(cmipj(2,2))**2
 
264
     +              - DBLE(xpi(2))*DBLE(cmipj(1,3))**2
 
265
     +              + 2*DBLE(piDpj(1,2))*DBLE(cmipj(2,2))*
 
266
     +                          DBLE(cmipj(1,3)) )/4
 
267
            if ( nschem .ge. 3 ) then
 
268
                cel3 = (- DBLE(xpi(1))*cmipj(2,2)**2
 
269
     +                  - DBLE(xpi(2))*cmipj(1,3)**2
 
270
     +                  + 2*DBLE(piDpj(1,2))*cmipj(2,2)*cmipj(1,3) )/4
 
271
            else
 
272
                cel3 = DBLE(del3)
 
273
            endif
 
274
            if ( lwrite ) print *,'cel3 = ',cel3
 
275
        endif
 
276
        del2 = -sdel2i**2
 
277
*
 
278
*       the routine as it stands can not handle sdel2<0.
 
279
*       the simplest solution seems to be to switch to sdel2>0 for
 
280
*       the time being - we calculate a complete 3point function so it
 
281
*       should not be a problem (just a sign).  Of course this spoils a
 
282
*       good check on the correctness.
 
283
*
 
284
        sdel2 = abs(sdel2i)
 
285
        if ( sdel2i .gt. 0 .and. lwrite ) print *,
 
286
     +          'ffxc0j: cannot handle sdel2>0, switched to sdel2<0'
 
287
*
 
288
        if ( xpi(4).eq.0 ) then
 
289
            zm = xpi(2)/dpipj(2,1)
 
290
            zm1 = -xpi(1)/dpipj(2,1)
 
291
        else
 
292
            call ffroot(zm,zp,xpi(4),piDpj(4,2),xpi(2),sdel2,ier)
 
293
            if ( dpipj(1,2) .ne. 0 ) then
 
294
                call ffroot(zp1,zm1,xpi(4),-piDpj(4,1),xpi(1),sdel2,ier)
 
295
            else
 
296
                zm1 = zp
 
297
                zp1 = zm
 
298
            endif
 
299
        endif
 
300
        if ( lwrite ) then
 
301
            print *,'ffxc0j: found roots:'
 
302
            print *,'  zm = ',zm,zm1
 
303
            if ( xpi(4).ne.0 ) print *,'  zp = ',zp,zp1
 
304
        endif
 
305
        if ( ltest ) then
 
306
            xheck = zm + zm1 - 1
 
307
            if ( xloss*abs(xheck) .gt. precx*max(x1,abs(zm)) ) print *,
 
308
     +          'ffxc0j: zm + zm1 <> 1: ',zm,zm1,xheck
 
309
            if ( xpi(4).ne.0 ) then
 
310
                xheck = zp + zp1 - 1
 
311
                if ( xloss*abs(xheck) .gt. precx*max(x1,abs(zp)) )
 
312
     +              print *,'ffxc0j: zp + zp1 <> 1: ',zp,zp1,xheck
 
313
            endif
 
314
        endif
 
315
 
 
316
*       imag sign ok 30-oct-1989
 
317
        ieps = -1
 
318
        if ( xpi(4).ne.0 ) dyzp = -2*sdel2/xpi(4)
 
319
*
 
320
*  #] get determinants, roots, ieps:
 
321
*  #[ the finite+divergent S1:
 
322
*
 
323
        if ( xpi(4).ne.0 ) then
 
324
            call ffcxr(cs(1),ipi12,zm,zm1,zp,zp1,dyzp,
 
325
     +          .FALSE.,x0,x0,x0,.FALSE.,dum,ieps,ier)
 
326
        endif
 
327
*
 
328
*       Next the divergent piece
 
329
*
 
330
        if ( .not.lsmug ) then
 
331
*
 
332
*           Here we dropped the term log(lam/delta)*log(-zm/zm1)
 
333
*
 
334
            if ( abs(zm1) .gt. 1/xloss ) then
 
335
                clog1 = dfflo1(1/zm1,ier)
 
336
            elseif ( zm.ne.0 ) then
 
337
                clog1 = zxfflg(-zm/zm1,-2,x0,ier)
 
338
            else
 
339
                call fferr(97,ier)
 
340
                return
 
341
            endif
 
342
            hulp = zm*zm1*4*del2/delta**2
 
343
*
 
344
*           14-jan-1994: do not count when this is small, this was 
 
345
*           meant to be so by the user carefully adjusting delta
 
346
*
 
347
            ier0 = ier
 
348
            if ( hulp.eq.0 ) call fferr(97,ier)
 
349
            clog2 = zxfflg(hulp,2,x0,ier0)
 
350
            cs(8) = -clog1*clog2/2
 
351
            if ( lwrite ) then
 
352
*               print *,'arg1 = ',-zm/zm1,1/zm1
 
353
                print *,'log1 = ',clog1
 
354
*               print *,'arg2 = ',hulp
 
355
                print *,'log2 = ',clog2
 
356
                print *,'cs(8)= ',cs(8)
 
357
            endif
 
358
        else
 
359
*
 
360
*           checked 4-aug-1992, but found Yet Another Bug 30-sep-1992
 
361
*
 
362
            cdyzm = cel3*DBLE(1/(-2*sdel2*del2))
 
363
             dyzm = del3/(-2*sdel2*del2)
 
364
            carg1 = +cdyzm*DBLE(1/zm)
 
365
             arg1 = +dyzm/zm
 
366
            clog1 = zfflog(-carg1,+ieps,DCMPLX(DBLE(zm),DBLE(0)),ier)
 
367
            if (DIMAG(cdyzm) .lt. 0 .and. arg1 .gt. 0 ) then
 
368
                clog1 = clog1 - c2ipi
 
369
                if ( lwrite ) then
 
370
                    print *,'added -2*pi*i to log1 S1'
 
371
                    print *,' arg1,zm = ',arg1,zm
 
372
                    print *,'carg1    = ',carg1
 
373
                endif
 
374
*               ier = ier + 50
 
375
            endif
 
376
            cs(8) = -clog1**2/2
 
377
            carg2 = -cdyzm*DBLE(1/zm1)
 
378
             arg2 = -dyzm/zm1
 
379
            clog2 = zfflog(-carg2,ieps,DCMPLX(DBLE(-zm1),DBLE(0)),ier)
 
380
            if ( DIMAG(cdyzm) .lt. 0 .and. arg2 .gt. 0 ) then
 
381
                clog2 = clog2 + c2ipi
 
382
                if ( lwrite ) then
 
383
                    print *,'added +2*pi*i to log2 S1'
 
384
                    print *,' arg2,zm = ',arg2,zm
 
385
                    print *,'carg2    = ',carg2
 
386
                endif
 
387
            endif
 
388
            cs(9) = +clog2**2/2
 
389
            if ( lwrite ) then
 
390
                print *,'y=zm = ',zm,zm1
 
391
                if ( xpi(4).ne.0 ) print *,'  zp = ',zp,zp1
 
392
                print *,'cdyzm= ',cdyzm
 
393
                print *,'arg1 = ',1/carg1
 
394
                print *,'log1 = ',clog1
 
395
                print *,'cs(8)= ',cs(8)
 
396
                print *,'arg2 = ',1/carg2
 
397
                print *,'log2 = ',clog2
 
398
                print *,'cs(9)= ',cs(9)
 
399
                print *,'ipi12= ',ipi12
 
400
                print *,'S1   = ',cs(1)+cs(2)+cs(3)+cs(4)+cs(5)+cs(6)+
 
401
     +                  cs(7)+cs(8)+cs(9)+ipi12*DBLE(pi12)
 
402
                print *,' '
 
403
            endif
 
404
        endif
 
405
*  #] the finite+divergent S1:
 
406
*  #[ log(1) for npoin=4:
 
407
        if ( npoin .eq. 4 ) then
 
408
            if ( ilogi(1) .eq. -999 ) then
 
409
                if ( .not.lsmug ) then
 
410
                    hulp = xpi(4)*delta/(4*del2)
 
411
                    ier0 = ier
 
412
                    if ( hulp.eq.0 ) call fferr(97,ier)
 
413
                    clogi(1) = -zxfflg(abs(hulp),0,x0,ier0)
 
414
                    if ( hulp .lt. 0 ) then
 
415
                        if ( xpi(4) .gt. 0 ) then
 
416
                            ilogi(1) = -1
 
417
                        else
 
418
                            ilogi(1) = +1
 
419
                        endif
 
420
                        if ( ltest ) then
 
421
                            print *,'ffxc0j: I am not 100% sure of the',
 
422
     +                          ' terms pi^2, please check against the',
 
423
     +                          ' limit lam->0 (id=',id,')'
 
424
                            ier = ier + 50
 
425
                        endif
 
426
                    else
 
427
                        ilogi(1) = 0
 
428
                    endif
 
429
                else
 
430
                    if ( xpi(4).eq.0 ) then
 
431
                        print *,'ffxc0i: cannot handle t=0 yet, sorry'
 
432
                        print *,'Please regularize with a small mass'
 
433
                        stop
 
434
                    endif
 
435
                    chulp = -cdyzm*DBLE(1/dyzp)
 
436
                    chulp1 = 1+chulp
 
437
                    if ( absc(chulp1) .lt. xloss )
 
438
     +                  call ffwarn(129,ier,absc(chulp1),x1)
 
439
                    call ffxclg(clogi(1),ilogi(1),chulp,chulp1,dyzp,
 
440
     +                  ier)
 
441
                endif
 
442
            endif
 
443
        endif
 
444
*  #] log(1) for npoin=4:
 
445
*  #[ the log(lam) Si:
 
446
        if ( .not.lsmug ) then
 
447
*
 
448
*           Next the divergent S_i (easy).
 
449
*           The term -2*log(lam/delta)*log(xpi(2)/xpi(1)) has been discarded
 
450
*           with lam the photon mass (regulator).
 
451
*           If delta = sqrt(xpi(1)*xpi(2)) the terms cancel as well
 
452
*
 
453
            if ( dpipj(1,2).ne.0 .and. xloss*abs(xpi(1)*xpi(2)-delta**2)
 
454
     +                  .gt.precx*delta**2 ) then
 
455
                if ( xpi(1) .ne. delta ) then
 
456
                    ier0 = ier
 
457
                    if ( xpi(1).eq.0 ) call fferr(97,ier)
 
458
                    cs(9) = -zxfflg(xpi(1)/delta,0,x0,ier0)**2 /4
 
459
                endif
 
460
                if ( xpi(2) .ne. delta ) then
 
461
                    ier0 = ier
 
462
                    if ( xpi(2).eq.0 ) call fferr(97,ier)
 
463
                    cs(10) = zxfflg(xpi(2)/delta,0,x0,ier0)**2 /4
 
464
                endif
 
465
            endif
 
466
            if ( lwrite ) then
 
467
                print *,'cs(9)= ',cs(9)
 
468
                print *,'cs(10)=',cs(10)
 
469
            endif
 
470
*  #] the log(lam) Si:
 
471
*  #[ the logs for A_i<0:
 
472
            if ( npoin.eq.4 ) then
 
473
                clogi(2) = 0
 
474
                ilogi(2) = 0
 
475
                clogi(3) = 0
 
476
                ilogi(3) = 0
 
477
            endif
 
478
*  #] the logs for A_i<0:
 
479
*  #[ the off-shell S3:
 
480
        else
 
481
*
 
482
*           the divergent terms in the offshell regulator scheme - not
 
483
*           quite as easy
 
484
*           wm = p3.p2/sqrtdel - 1 = -s1.s2/sqrtdel - 1
 
485
*           wp = p3.p2/sqrtdel + 1 = -s1.s2/sqrtdel + 1
 
486
*           Note that we took the choice sdel2<0 in S1 when
 
487
*           \delta^{p1 s2}_{p1 p2} < 0 by using yp=zm
 
488
*
 
489
            wm = -1 - piDpj(1,2)/sdel2
 
490
            wp = wm + 2
 
491
            if ( lwrite ) print *,'wm,wp = ',wm,wp
 
492
            if ( abs(wm) .lt. abs(wp) ) then
 
493
                wm = -xpi(5)*xpi(6)/(del2*wp)
 
494
                if ( lwrite ) print *,'wm+   = ',wm
 
495
            else
 
496
                wp = -xpi(5)*xpi(6)/(del2*wm)
 
497
                if ( lwrite ) print *,'wp+   = ',wp
 
498
            endif
 
499
*
 
500
*           the im sign
 
501
*
 
502
            if ( -DBLE(cmipj(1,3)) .gt. 0 ) then
 
503
                ieps = -1
 
504
            else
 
505
                ieps = +1
 
506
            endif
 
507
*
 
508
            if ( nschem .lt. 3 .or. DIMAG(cmipj(1,3)).eq.0 .and.
 
509
     +                  DIMAG(cmipj(2,2)).eq.0 ) then
 
510
*  #[           real case:
 
511
                if ( lwrite ) print *,'ffxc0i: Real S3'
 
512
*
 
513
*               first z-,z+
 
514
*
 
515
                dyzp = -DBLE(cmipj(1,3))*DBLE(wm)/(2*DBLE(xpi(6))) -
 
516
     +                  DBLE(cmipj(2,2))/(2*DBLE(sdel2))
 
517
                dyzm = -DBLE(cmipj(1,3))*DBLE(wp)/(2*DBLE(xpi(6))) -
 
518
     +                  DBLE(cmipj(2,2))/(2*DBLE(sdel2))
 
519
*
 
520
*               the (di)logs
 
521
*
 
522
                clog1 = zxfflg(-dyzp,-ieps,x1,ier)
 
523
                cs(10) = -clog1**2/2
 
524
                ipi12 = ipi12 - 4
 
525
                clog2 = zxfflg(-dyzm,+ieps,x1,ier)
 
526
                cs(11) = -clog2**2/2
 
527
                ipi12 = ipi12 - 2
 
528
                hulp = dyzp/dyzm
 
529
                if ( dyzp .lt. 0 ) then
 
530
                    ieps1 = -ieps
 
531
                else
 
532
                    ieps1 = +ieps
 
533
                endif
 
534
                call ffzxdl(cli,i,cdum(1),hulp,+ieps1,ier)
 
535
                cs(12) = -cli
 
536
                ipi12 = ipi12 - i
 
537
*
 
538
*               the log for npoin=4
 
539
*
 
540
                if ( npoin.eq.4 ) then
 
541
                    if ( ilogi(3) .eq. -999 ) then
 
542
                        if ( DBLE(cmipj(1,3)) .eq. 0 ) then
 
543
                            chulp = -1
 
544
                            chulp1 = 0
 
545
                        elseif ( dyzp .lt. dyzm ) then
 
546
                            chulp = -dyzm/dyzp
 
547
                            chulp1 = +DBLE(cmipj(1,3))/DBLE(xpi(6)*dyzp)
 
548
                        else
 
549
                            chulp = -dyzp/dyzm
 
550
                            chulp1 = -DBLE(cmipj(1,3))/DBLE(xpi(6)*dyzm)
 
551
                        endif
 
552
                        call ffxclg(clogi(3),ilogi(3),chulp,chulp1,dyzp,
 
553
     +                          ier)
 
554
                    endif
 
555
                endif
 
556
*
 
557
*               and some debug output:
 
558
*
 
559
                if ( lwrite ) then
 
560
                    print *,'z      = 1,0'
 
561
                    print *,'y-zm   = ',dyzm
 
562
                    print *,'y-zp   = ',dyzp
 
563
                    print *,'+Li2(y/(y-zp))    = ',cs(10)
 
564
                    print *,'+Li2(y/(y-zm))    = ',cs(11)
 
565
                    print *,'-Li2((y-1)/(y-zm))= ',cs(12)
 
566
                    print *,'ipi12   = ',ipi12
 
567
                endif
 
568
*  #]           real case:
 
569
            else
 
570
*  #[           complex case:
 
571
                if ( lwrite ) print *,'ffxc0i: Complex S3'
 
572
*
 
573
*               first z+
 
574
*
 
575
                cdyzp = -cmipj(1,3)*DBLE(wm)/(2*DBLE(xpi(6))) -
 
576
     +                  cmipj(2,2)/(2*DBLE(sdel2))
 
577
                clog1 = zfflog(-cdyzp,-ieps,c1,ier)
 
578
                if ( ieps*DIMAG(cdyzp).lt.0.and.DBLE(cdyzp).gt.0 ) then
 
579
                    if ( lwrite ) then
 
580
                        print *,'added ',-ieps,'*2*pi*i to log1 S3'
 
581
                        print *,'carg1    = ',-cdyzp
 
582
                        print *,'clog1 was  ',clog1
 
583
                        print *,'clog1 is   ',clog1 - ieps*c2ipi
 
584
                    endif
 
585
                    clog1 = clog1 - ieps*c2ipi
 
586
                else
 
587
                    if ( lwrite ) then
 
588
                        print *,'carg1    = ',-cdyzp
 
589
                        print *,'clog1 is   ',clog2
 
590
                    endif
 
591
                endif
 
592
                cs(10) = -clog1**2/2
 
593
                ipi12 = ipi12 - 4
 
594
*
 
595
*               now z-
 
596
*
 
597
                cdyzm = -cmipj(1,3)*DBLE(wp)/(2*DBLE(xpi(6))) -
 
598
     +                  cmipj(2,2)/(2*DBLE(sdel2))
 
599
                clog2 = zfflog(-cdyzm,+ieps,c1,ier)
 
600
                if ( ieps*DIMAG(cdyzm).gt.0.and.DBLE(cdyzm).gt.0 ) then
 
601
                    if ( lwrite ) then
 
602
                        print *,'added ',ieps,'*2*pi*i to log2 S3'
 
603
                        print *,'carg2    = ',-cdyzm
 
604
                        print *,'clog2 was  ',clog2
 
605
                        print *,'clog2 is   ',clog2 + ieps*c2ipi
 
606
                    endif
 
607
                    clog2 = clog2 + ieps*c2ipi
 
608
*                   ier = ier + 50
 
609
                else
 
610
                    if ( lwrite ) then
 
611
                        print *,'carg2    = ',-cdyzm
 
612
                        print *,'clog2 is   ',clog2
 
613
                    endif
 
614
                endif
 
615
                cs(11) = -clog2**2/2
 
616
                ipi12 = ipi12 - 2
 
617
*
 
618
*               the dilog
 
619
*
 
620
                chulp = cdyzp/cdyzm
 
621
                hulp = DBLE(cdyzp)/DBLE(cdyzm)
 
622
                if ( DBLE(cdyzp) .lt. 0 ) then
 
623
                    ieps1 = -ieps
 
624
                else
 
625
                    ieps1 = +ieps
 
626
                endif
 
627
                if ( DIMAG(chulp) .eq. 0 ) then
 
628
                    hulp = DBLE(chulp)
 
629
                    call ffzxdl(cli,i,cdum(1),hulp,+ieps1,ier)
 
630
                else
 
631
                    call ffzzdl(cli,i,cdum(1),chulp,ier)
 
632
                    if ( hulp.gt.1 .and. ieps1*DIMAG(chulp).lt.0 ) then
 
633
                        if ( lwrite ) then
 
634
                            print *,'addded 2ipi*log(z) to Li'
 
635
                            print *,'chulp = ',chulp
 
636
                            print *,'cli was ',cli
 
637
                            print *,'cli is  ',cli +
 
638
     +                          ieps1*c2ipi*zfflog(chulp,0,c0,ier)
 
639
                            call ffzxdl(cdum(2),i,cdum(1),hulp,+ieps1,
 
640
     +                          ier)
 
641
                            print *,'vgl   ',cdum(2)
 
642
                        endif
 
643
                        cli = cli + ieps1*c2ipi*zfflog(chulp,0,c0,ier)
 
644
                    endif
 
645
                endif
 
646
                cs(12) = -cli
 
647
                ipi12 = ipi12 - i
 
648
*
 
649
*               the log for npoin=4
 
650
*
 
651
                if ( npoin.eq.4 ) then
 
652
                    if ( ilogi(3) .eq. -999 ) then
 
653
                        if ( cmipj(1,3) .eq. 0 ) then
 
654
                            chulp = -1
 
655
                            chulp1 = 0
 
656
                        elseif ( DBLE(cdyzp) .lt. DBLE(cdyzm) ) then
 
657
                            chulp = -cdyzm/cdyzp
 
658
                            chulp1 = +cmipj(1,3)/cdyzp*DBLE(1/xpi(6))
 
659
                        else
 
660
                            chulp = -cdyzp/cdyzm
 
661
                            chulp1 = -cmipj(1,3)/cdyzm*DBLE(1/xpi(6))
 
662
                        endif
 
663
                        dyzp = DBLE(cdyzp)
 
664
                        call ffxclg(clogi(3),ilogi(3),chulp,chulp1,dyzp,
 
665
     +                          ier)
 
666
                    endif
 
667
                endif
 
668
*
 
669
*               and some debug output:
 
670
*
 
671
                if ( lwrite ) then
 
672
                    print *,'z      = 1,0'
 
673
                    print *,'y-zm   = ',cdyzm
 
674
                    print *,'y-zp   = ',cdyzp
 
675
                    print *,'+Li2(y/(y-zp))    = ',cs(10)
 
676
                    print *,'+Li2(y/(y-zm))    = ',cs(11)
 
677
                    print *,'-Li2((y-1)/(y-zm))= ',cs(12)
 
678
                    print *,'ipi12   = ',ipi12
 
679
                endif
 
680
*  #]           complex case:
 
681
            endif
 
682
*  #] the off-shell S3:
 
683
*  #[ the off-shell S2:
 
684
*
 
685
*           the im sign
 
686
*
 
687
            if ( -DBLE(cmipj(2,2)) .gt. 0 ) then
 
688
                ieps = -1
 
689
            else
 
690
                ieps = +1
 
691
            endif
 
692
*
 
693
            if ( nschem .lt. 3 ) then
 
694
*  #[           real case:
 
695
                if ( lwrite ) print *,'ffxc0i: Real S2'
 
696
*
 
697
*               first z-
 
698
*
 
699
                dyzm = -DBLE(cmipj(2,2))*DBLE(wp)/(2*DBLE(xpi(5))) -
 
700
     +                  DBLE(cmipj(1,3))/(2*DBLE(sdel2))
 
701
                clog1 = zxfflg(+dyzm,-ieps,x1,ier)
 
702
                cs(13) = +clog1**2/2
 
703
                ipi12 = ipi12 + 4
 
704
*
 
705
*               now z+
 
706
*
 
707
                dyzp = -DBLE(cmipj(2,2))*DBLE(wm)/(2*DBLE(xpi(5))) -
 
708
     +                  DBLE(cmipj(1,3))/(2*DBLE(sdel2))
 
709
                clog2 = zxfflg(+dyzp,+ieps,x1,ier)
 
710
                cs(14) = +clog2**2/2
 
711
                ipi12 = ipi12 + 2
 
712
                hulp = dyzm/dyzp
 
713
                if ( dyzm .lt. 0 ) then
 
714
                    ieps1 = -ieps
 
715
                else
 
716
                    ieps1 = +ieps
 
717
                endif
 
718
                call ffzxdl(cli,i,cdum(1),hulp,-ieps1,ier)
 
719
                cs(15) = +cli
 
720
                ipi12 = ipi12 + i
 
721
*
 
722
*               the log for npoin=4
 
723
*
 
724
                if ( npoin.eq.4 ) then
 
725
                    if ( ilogi(2) .eq. -999 ) then
 
726
                        if ( DBLE(cmipj(2,2)) .eq. 0 ) then
 
727
                            chulp = -1
 
728
                            chulp1 = 0
 
729
                        elseif ( dyzp .lt. dyzm ) then
 
730
                            chulp = -dyzm/dyzp
 
731
                            chulp1 = +DBLE(cmipj(2,2))/DBLE(xpi(5)*dyzp)
 
732
                        elseif ( dyzp .gt. dyzm ) then
 
733
                            chulp = -dyzp/dyzm
 
734
                            chulp1 = -DBLE(cmipj(2,2))/DBLE(xpi(5)*dyzm)
 
735
                        endif
 
736
                        call ffxclg(clogi(2),ilogi(2),chulp,chulp1,dyzp,
 
737
     +                          ier)
 
738
                    endif
 
739
                endif
 
740
*
 
741
*               and some debug output:
 
742
*
 
743
                if ( lwrite ) then
 
744
                    print *,'z      = 0,1'
 
745
                    print *,'y-zm   = ',dyzm
 
746
                    print *,'y-zp   = ',dyzp
 
747
                    print *,'-Li2((y-1)/(y-zm))= ',cs(13)
 
748
                    print *,'-Li2((y-1)/(y-zp))= ',cs(14)
 
749
                    print *,'+Li2(y/(y-zp))    = ',cs(15)
 
750
                    print *,'ipi12   = ',ipi12
 
751
                endif
 
752
*  #]           real case:
 
753
            else
 
754
*  #[           complex case:
 
755
                if ( lwrite ) print *,'ffxc0i: Complex S2'
 
756
*
 
757
*               first z-
 
758
*
 
759
                cdyzm = -cmipj(2,2)*DBLE(wp)/(2*DBLE(xpi(5))) -
 
760
     +                  cmipj(1,3)/(2*DBLE(sdel2))
 
761
                clog1 = zfflog(+cdyzm,-ieps,c1,ier)
 
762
                if ( DBLE(cdyzm).lt.0.and.ieps*DIMAG(cdyzm).gt.0 ) then
 
763
                    if ( lwrite ) print *,'added 2*i*pi to log1'
 
764
                    clog1 = clog1 - ieps*c2ipi
 
765
                endif
 
766
                cs(13) = +clog1**2/2
 
767
                ipi12 = ipi12 + 4
 
768
*
 
769
*               now z+
 
770
*
 
771
                cdyzp = -cmipj(2,2)*DBLE(wm)/(2*DBLE(xpi(5))) -
 
772
     +                  cmipj(1,3)/(2*DBLE(sdel2))
 
773
                clog2 = zfflog(+cdyzp,+ieps,c1,ier)
 
774
                if ( DBLE(cdyzp).lt.0.and.ieps*DIMAG(cdyzp).lt.0 ) then
 
775
                    if ( lwrite ) then
 
776
                        print *,'added ',ieps,'*2*pi*i to log2 S2'
 
777
                        print *,'carg1    = ',+cdyzp
 
778
                    endif
 
779
                    clog2 = clog2 + ieps*c2ipi
 
780
                endif
 
781
                cs(14) = +clog2**2/2
 
782
                ipi12 = ipi12 + 2
 
783
*               
 
784
*               and ghe dilog
 
785
*               
 
786
                chulp = cdyzm/cdyzp
 
787
                hulp = DBLE(dyzm)/DBLE(dyzp)
 
788
                if ( DBLE(cdyzm) .lt. 0 ) then
 
789
                    ieps1 = -ieps
 
790
                else
 
791
                    ieps1 = +ieps
 
792
                endif
 
793
                if ( DIMAG(chulp ) .eq. 0 ) then
 
794
                    hulp = DBLE(chulp)
 
795
                    call ffzxdl(cli,i,cdum(1),hulp,-ieps1,ier)
 
796
                else
 
797
                    call ffzzdl(cli,i,cdum(1),chulp,ier)
 
798
                    if ( hulp.gt.1 .and. ieps1*DIMAG(chulp).gt.0 ) then
 
799
                        if ( lwrite ) then
 
800
                            print *,'addded 2ipi*log(z) to Li'
 
801
                            print *,'chulp = ',chulp
 
802
                            print *,'cli was ',cli
 
803
                            print *,'cli is  ',cli -
 
804
     +                          ieps1*c2ipi*zfflog(chulp,0,c0,ier)
 
805
                            call ffzxdl(cdum(2),i,cdum(1),hulp,-ieps1,
 
806
     +                          ier)
 
807
                            print *,'vgl   ',cdum(2)
 
808
                        endif
 
809
                        cli = cli - ieps1*c2ipi*zfflog(chulp,0,c0,ier)
 
810
                    endif
 
811
                endif
 
812
                cs(15) = +cli
 
813
                ipi12 = ipi12 + i
 
814
*
 
815
*               the log for npoin=4
 
816
*
 
817
                if ( npoin.eq.4 ) then
 
818
                    if ( ilogi(2) .eq. -999 ) then
 
819
                        if ( cmipj(2,2) .eq. 0 ) then
 
820
                            chulp = -1
 
821
                            chulp1 = 0
 
822
                        elseif ( DBLE(cdyzp) .lt. DBLE(cdyzm) ) then
 
823
                            chulp = -cdyzm/cdyzp
 
824
                            chulp1 = +cmipj(2,2)/cdyzp*DBLE(1/xpi(5))
 
825
                        elseif ( DBLE(cdyzp) .gt. DBLE(cdyzm) ) then
 
826
                            chulp = -cdyzp/cdyzm
 
827
                            chulp1 = -cmipj(2,2)/cdyzm*DBLE(1/xpi(5))
 
828
                        endif
 
829
                        dyzp = DBLE(cdyzp)
 
830
                        call ffxclg(clogi(2),ilogi(2),chulp,chulp1,dyzp,
 
831
     +                          ier)
 
832
                    endif
 
833
                endif
 
834
*
 
835
*               and some debug output:
 
836
*
 
837
                if ( lwrite ) then
 
838
                    print *,'z      = 0,1'
 
839
                    print *,'y-zm   = ',cdyzm
 
840
                    print *,'y-zp   = ',cdyzp
 
841
                    print *,'-Li2((y-1)/(y-zm))= ',cs(13)
 
842
                    print *,'-Li2((y-1)/(y-zp))= ',cs(14)
 
843
                    print *,'+Li2(y/(y-zp))    = ',cs(15)
 
844
                    print *,'ipi12   = ',ipi12
 
845
                endif
 
846
*  #]           complex case:
 
847
            endif
 
848
        endif
 
849
*  #] the off-shell S2:
 
850
*  #[ sdel2<0!:
 
851
        if ( sdel2i.gt.0 .neqv. xpi(4).eq.0.and.xpi(1).gt.xpi(2) ) then
 
852
            if ( .not.lsmug ) then
 
853
                n = 10
 
854
            else
 
855
                n = 15
 
856
            endif
 
857
            do 10 i=1,n
 
858
                cs(i) = -cs(i)
 
859
   10       continue
 
860
            ipi12 = -ipi12
 
861
            if ( npoin.eq.4 ) then
 
862
                do 20 i=1,3
 
863
                    ilogi(i) = -ilogi(i)
 
864
                    clogi(i) = -clogi(i)
 
865
   20           continue
 
866
            endif
 
867
        endif
 
868
        if ( lwrite ) print '(a)','  ##] ffxc0j:'
 
869
*  #] sdel2<0!:
 
870
*###] ffxc0j:
 
871
        end
 
872
*###[ ffxclg:
 
873
        subroutine ffxclg(clg,ilg,chulp,chulp1,dyzp,ier)
 
874
***#[*comment:***********************************************************
 
875
*                                                                       *
 
876
*       compute the extra logs for npoin=4 given chulp=-cdyzm/cdyzp     *
 
877
*       all flagchecking has already been done.                         *
 
878
*                                                                       *
 
879
*       Input:  chulp   (complex)       see above                       *
 
880
*               chulp1  (complex)       1+chulp (in case chulp ~ -1)    *
 
881
*               dyzp    (real)          (real part of) y-z+ for im part *
 
882
*       Output: clg     (complex)       the log                         *
 
883
*               ilg     (integer)       factor i*pi split off clg       *
 
884
*                                                                       *
 
885
***#]*comment:***********************************************************
 
886
*  #[ declarations:
 
887
        implicit none
 
888
*
 
889
*       arguments
 
890
*
 
891
        integer ilg,ier
 
892
        DOUBLE PRECISION dyzp
 
893
        DOUBLE COMPLEX clg,chulp,chulp1
 
894
*
 
895
*       local variables
 
896
*
 
897
        DOUBLE PRECISION hulp,hulp1,dfflo1
 
898
        DOUBLE COMPLEX zxfflg,zfflog,zfflo1,check
 
899
*
 
900
*       common blocks
 
901
*
 
902
        include 'ff.h'
 
903
*
 
904
*  #] declarations:
 
905
*  #[ check:
 
906
        if ( ltest ) then
 
907
            check = c1 + chulp - chulp1
 
908
            if ( xloss*abs(check) .gt. precc*max(abs(c1),abs(chulp)) )
 
909
     +          print *,'ffxclg: error: chulp1 != 1+chulp: ',chulp1,
 
910
     +          c1+chulp,check
 
911
        endif
 
912
*  #] check:
 
913
*  #[ work:
 
914
*
 
915
        if ( DIMAG(chulp) .eq. 0 ) then
 
916
            hulp = DBLE(chulp)
 
917
            hulp1 = DBLE(chulp1)
 
918
            if ( abs(hulp1) .lt. xloss ) then
 
919
                clg = DBLE(dfflo1(hulp1,ier))
 
920
            else
 
921
                clg = zxfflg(abs(hulp),0,x0,ier)
 
922
            endif
 
923
            if ( hulp .lt. 0 ) then
 
924
                if ( dyzp.lt.0 ) then
 
925
                    ilg = +1
 
926
                else
 
927
                    ilg = -1
 
928
                endif
 
929
            else
 
930
                ilg = 0
 
931
            endif
 
932
            if ( lwrite ) print *,'clg(real) = ',clg+c2ipi*ilg/2
 
933
        else
 
934
*
 
935
*           may have to be improved
 
936
*
 
937
            if ( abs(DBLE(chulp1))+abs(DIMAG(chulp1)) .lt. xloss ) then
 
938
                clg = zfflo1(chulp1,ier)
 
939
            else
 
940
                clg = zfflog(chulp,0,c0,ier)
 
941
            endif
 
942
            ilg = 0
 
943
            if ( DBLE(chulp) .lt. 0 ) then
 
944
                if ( dyzp.lt.0 .and. DIMAG(clg).lt.0 ) then
 
945
                    if ( lwrite ) print *,'ffxclg: added -2*pi to log'
 
946
                    ilg = +2
 
947
                elseif ( dyzp.gt.0 .and. DIMAG(clg).gt.0 ) then
 
948
                    if ( lwrite ) print *,'ffxclg: added +2*pi to log'
 
949
                    ilg = -2
 
950
                endif
 
951
            endif
 
952
            if ( lwrite ) print *,'clg(cmplx)= ',clg+c2ipi*ilg/2
 
953
        endif
 
954
*  #] work:
 
955
*###] ffxclg:
 
956
        end