← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

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

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/CutTools/src/qcdloop/ffcxs3.f
 
1
*###[ ffcxs3:
 
2
        subroutine ffcxs3(cs3,ipi12,y,z,dyz,d2yzz,dy2z,xpi,piDpj,ii,ns,
 
3
     +                                  isoort,ier)
 
4
***#[*comment:***********************************************************
 
5
*                                                                       *
 
6
*       calculates the s3 as defined in appendix b.                     *
 
7
*               (ip = ii+3, is1 = ii, is2 = ii+1)                       *
 
8
*                                                                       *
 
9
*                 log( xk*y^2 + (-xk+xm1-xm2)*y + xm2 - i*eps )         *
 
10
*            /1                                   - log( ... ) |y=yi    *
 
11
*       s3 = \ dy --------------------------------------------------    *
 
12
*            /0                         y - yi                          *
 
13
*                                                                       *
 
14
*           = r(yi,y-,+) + r(yi,y+,-)                                   *
 
15
*                                                                       *
 
16
*       with y+- the roots of the argument of the logarithm.            *
 
17
*       the sign of the argument to the logarithms in r is passed       *
 
18
*       in ieps                                                         *
 
19
*                                                                       *
 
20
*       input:  y(4),z(4)       (real)  roots in form (z-,z+,1-z-,1-z+) *
 
21
*               dyz(2,2),d2yzz, (real)  y() - z(), y+ - z- - z+         *
 
22
*               dy2z(4)         (real)  y() - 2z()                      *
 
23
*               xpi     (real(ns))      p(i).p(i) (B&D metric)  i=1,3   *
 
24
*                                       m(i)^2 = si.si          i=4,6   *
 
25
*               ii      (integer)       xk = xpi(ii+3) etc              *
 
26
*               ns      (integer)       size of arrays                  *
 
27
*               isoort  (integer)       returns kind of action taken    *
 
28
*               cs3     (complex)(20)   assumed zero.                   *
 
29
*               ccy     (complex)(3)    if i0 != 0: complex y           *
 
30
*                                                                       *
 
31
*       output: cs3     (complex)       mod factors pi^2/12, in array   *
 
32
*               ipi12   (integer)       these factors                   *
 
33
*               ier     (integer)       0=ok 1=inaccurate 2=error       *
 
34
*                                                                       *
 
35
*       calls:  ffcrr,ffcxr,real/dble,DCMPLX,log,ffadd1,ffadd2,ffadd3   *
 
36
*                                                                       *
 
37
***#]*comment:***********************************************************
 
38
*  #[ declarations:
 
39
        implicit none
 
40
*
 
41
*       arguments:
 
42
*
 
43
        integer ipi12(2),ii,ns,isoort(2),ier
 
44
        DOUBLE COMPLEX cs3(20)
 
45
        DOUBLE PRECISION y(4),z(4),dyz(2,2),d2yzz,dy2z(4),
 
46
     +          xpi(ns),piDpj(ns,ns)
 
47
*
 
48
*       local variables:
 
49
*
 
50
        integer i,ip,ieps(2),ipi12p(2),ier0,i2,i3
 
51
        DOUBLE COMPLEX c,csum,cs3p(14)
 
52
        DOUBLE PRECISION yy,yy1,zz,zz1,dyyzz,xdilog,xlog,x00(3)
 
53
        DOUBLE PRECISION absc,xmax
 
54
        logical ld2yzz
 
55
*
 
56
*       common blocks
 
57
*
 
58
        include 'ff.h'
 
59
        absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
60
*
 
61
*  #] declarations:
 
62
*  #[ get counters:
 
63
        if ( ltest .and. ns .ne. 6 )
 
64
     +      print *,'ffcxs3: error: only for ns=6, not ',ns
 
65
        ip = ii+3
 
66
        if ( isoort(2) .ne. 0 ) then
 
67
            if ( (z(2).gt.z(1) .or.  z(1).eq.z(2) .and. z(4).lt.z(3) )
 
68
     +                  .eqv. xpi(ip) .gt. 0 ) then
 
69
                ieps(1) = +1
 
70
                ieps(2) = -1
 
71
            else
 
72
                ieps(1) = -1
 
73
                ieps(2) = +1
 
74
            endif
 
75
        else
 
76
            if ( piDpj(ip,ii) .gt. 0 ) then
 
77
                ieps(1) = +1
 
78
            else
 
79
                ieps(1) = -1
 
80
            endif
 
81
        endif
 
82
        i2 = mod(ii,3) + 1
 
83
        i3 = mod(i2,3) + 1
 
84
*  #] get counters:
 
85
*  #[ special case |z| >> |y|:
 
86
        if ( xpi(ip).lt.0 .and. max(abs(y(2)),abs(y(4))) .lt.
 
87
     +          xloss*min(abs(z(1)), abs(z(2)))/2 ) then
 
88
*
 
89
*           we will obtain cancellations of the type Li_2(x) + Li_2(-x)
 
90
*           with x small.
 
91
*
 
92
            if ( lwrite ) then
 
93
                print *,'ffcxs3: special case |z| >> |y|'
 
94
                print *,'  y,y1  = ',y(2),y(4)
 
95
                print *,'  z,z1- = ',z(1),z(3)
 
96
                print *,'  z,z1+ = ',z(2),z(4)
 
97
            endif
 
98
            yy = dyz(2,1)/d2yzz
 
99
            yy1 = dyz(2,2)/d2yzz
 
100
            if ( y(2) .eq. 0 ) goto 10
 
101
            zz = z(2)*yy/y(2)
 
102
            zz1 = 1-zz
 
103
            if ( lwarn .and. abs(zz) .lt. xloss )
 
104
     +          call ffwarn(44,ier,abs(zz),x1)
 
105
            dyyzz = dyz(2,2)*yy/y(2)
 
106
            call ffcxr(cs3(1),ipi12(1),yy,yy1,zz,zz1,dyyzz,.FALSE.,x0,
 
107
     +                  x0,x0,.FALSE.,x00,0,ier)
 
108
   10       continue
 
109
            if ( y(4) .eq. 0 ) goto 30
 
110
            zz = yy*z(4)/y(4)
 
111
            zz1 = 1-zz
 
112
            if ( lwarn .and. abs(zz) .lt. xloss )
 
113
     +          call ffwarn(44,ier,abs(zz),x1)
 
114
            dyyzz = -yy*dyz(2,2)/y(4)
 
115
            call ffcxr(cs3(8),ipi12(2),yy,yy1,zz,zz1,dyyzz,.FALSE.,x0,
 
116
     +                  x0,x0,.FALSE.,x00,0,ier)
 
117
            do 20 i=8,14
 
118
   20           cs3(i) = -cs3(i)
 
119
   30       continue
 
120
*           And now the remaining Li_2(x^2) terms
 
121
            call ffxli2(xdilog,xlog,(y(2)/dyz(2,1))**2,ier)
 
122
            cs3(15) = +xdilog/2
 
123
            call ffxli2(xdilog,xlog,(y(4)/dyz(2,1))**2,ier)
 
124
            cs3(16) = -xdilog/2
 
125
            if ( lwrite ) then
 
126
                lwrite = .FALSE.
 
127
                ipi12p(1) = 0
 
128
                ipi12p(2) = 0
 
129
                ier0 = 0
 
130
                do 40 i=1,14
 
131
   40               cs3p(i) = 0
 
132
                call ffcxr(cs3p(1),ipi12p(1),y(2),y(4),z(1),z(3),
 
133
     +              dyz(2,1),.FALSE.,x0,x0,x0,.FALSE.,x00,ieps(1),ier0)
 
134
                call ffcxr(cs3p(8),ipi12p(2),y(2),y(4),z(2),z(4),
 
135
     +              dyz(2,2),.FALSE.,x0,x0,x0,.FALSE.,x00,ieps(2),ier0)
 
136
                csum = 0
 
137
                xmax = 0
 
138
                do 50 i=1,14
 
139
                    csum = csum + cs3p(i)
 
140
                    xmax = max(xmax,absc(csum))
 
141
   50           continue
 
142
                csum = csum + (ipi12p(1)+ipi12(2))*DBLE(pi12)
 
143
                print '(a,3g20.10,3i3)','cmp',csum,xmax,ipi12p,ier0
 
144
                lwrite = .TRUE.
 
145
            endif
 
146
            goto 900
 
147
        endif
 
148
*  #] special case |z| >> |y|:
 
149
*  #[ normal:
 
150
        if ( xpi(ip) .eq. 0 ) then
 
151
            ld2yzz = .FALSE.
 
152
        else
 
153
            ld2yzz = .TRUE.
 
154
        endif
 
155
        if ( lwrite ) print *,  'ieps = ',ieps
 
156
        if ( isoort(1) .ne. 0 ) call ffcxr(cs3(1),ipi12(1),y(2),y(4),
 
157
     +      z(1),z(3),dyz(2,1),ld2yzz,d2yzz,z(2),z(4),.TRUE.,dy2z(1),
 
158
     +      ieps(1),ier)
 
159
        if ( isoort(2) .ne. 0 ) then
 
160
            if ( mod(isoort(2),10) .eq. 2 ) then
 
161
*               both roots are equal: multiply by 2
 
162
                if ( lwrite ) print *,'ffcxs3: skipped next R as it ',
 
163
     +                  'is the conjugate'
 
164
                do 60 i=1,7
 
165
                    cs3(i) = 2*DBLE(cs3(i))
 
166
   60           continue
 
167
                ipi12(1) = 2*ipi12(1)
 
168
            else
 
169
                call ffcxr(cs3(8),ipi12(2),y(2),y(4),z(2),z(4),dyz(2,2),
 
170
     +              ld2yzz,d2yzz,z(1),z(3),.TRUE.,dy2z(2),ieps(2),ier)
 
171
            endif
 
172
        endif
 
173
*
 
174
*  #] normal:
 
175
*  #[ print output:
 
176
  900   if (lwrite) then
 
177
            print *,'  cs3 ='
 
178
            do 905 i=1,20
 
179
                if ( cs3(i).ne.0 ) print '(i3,2g20.10,1x)',i,cs3(i)
 
180
  905       continue
 
181
            print '(a3,2g20.10,1x)','pi ',(ipi12(1)+ipi12(2))*pi12
 
182
            print *,'+-----------'
 
183
            csum = 0
 
184
            do 910 i=1,20
 
185
  910       csum = csum + cs3(i)
 
186
            csum = csum+(ipi12(1)+ipi12(2))*DBLE(pi12)
 
187
            print '(a,2g20.10)','Si ',csum
 
188
            print *,'  ipi12,ier=',ipi12,ier
 
189
            print *,' '
 
190
        endif
 
191
*  #] print output:
 
192
*###] ffcxs3:
 
193
        end
 
194
*###[ ffcs3:
 
195
        subroutine ffcs3(cs3,ipi12,cy,cz,cdyz,cd2yzz,cpi,cpiDpj,ii,ns,
 
196
     +          isoort,ier)
 
197
***#[*comment:***********************************************************
 
198
*                                                                       *
 
199
*       calculates the s3 as defined in appendix b.                     *
 
200
*                                                                       *
 
201
*                 log( cpi(ii+3)*y^2 + (cpi(ii+3)+cpi(ii)-cpi(ii+1))*y  *
 
202
*            /1                      +  cpi(ii+1))  - log( ... ) |y=cyi *
 
203
*       s3 = \ dy ----------------------------------------------------  *
 
204
*            /0                         y - cyi                         *
 
205
*                                                                       *
 
206
*          = r(cyi,cy+) + r(cyi,cy-) +  ( eta(-cy-,-cy+) -              *
 
207
*               eta(1-cy-,1-cy+) - eta(...) )*log(1-1/cyi)              *
 
208
*                                                                       *
 
209
*       with y+- the roots of the argument of the logarithm.            *
 
210
*                                                                       *
 
211
*       input:  cy(4)    (complex)  cy(1)=y^-,cy(2)=y^+,cy(i+2)=1-cy(1) *
 
212
*               cz(4)    (complex)  cz(1)=z^-,cz(2)=z^+,cz(i+2)=1-cz(1) *
 
213
*               cpi(6)   (complex)  masses & momenta (B&D)              *
 
214
*               ii       (integer)  position of cp,cma,cmb in cpi       *
 
215
*               ns       (integer)  size of arrays                      *
 
216
*               isoort(2)(integer)  returns the kind of action taken    *
 
217
*               cs3      (complex)(14)  assumed zero.                   *
 
218
*                                                                       *
 
219
*       output: cs3      (complex)  mod factors ipi12                   *
 
220
*               ipi12(2) (integer)  these factors                       *
 
221
*               ier      (integer)  0=ok, 1=numerical problems, 2=error *
 
222
*                                                                       *
 
223
*       calls:  ffcrr,DIMAG,DBLE,zfflog                                 *
 
224
*                                                                       *
 
225
***#]*comment:***********************************************************
 
226
*  #[ declarations:
 
227
        implicit none
 
228
*
 
229
*       arguments:
 
230
*
 
231
        integer ipi12(2),ii,ns,isoort(2),ier
 
232
        DOUBLE COMPLEX cs3(20),cpi(ns),cpiDpj(ns,ns)
 
233
        DOUBLE COMPLEX cy(4),cz(4),cdyz(2,2),cd2yzz
 
234
*
 
235
*       local variables:
 
236
*
 
237
        integer i,ip,ieps(2),ieps0,ni(4),ipi12p(2),ier0,ntot,i2,i3
 
238
        logical ld2yzz
 
239
        DOUBLE COMPLEX c,csum,zdilog,zlog,cyy,cyy1,czz,czz1,cdyyzz
 
240
     +          ,cs3p(14)
 
241
        DOUBLE PRECISION absc,xmax,y,y1,z,z1,dyz,d2yzz,zz,zz1,
 
242
     +          x00(3),sprec
 
243
*
 
244
*       common blocks
 
245
*
 
246
        include 'ff.h'
 
247
*
 
248
*       statement function
 
249
*
 
250
        absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
251
*  #] declarations:
 
252
*  #[ get ieps:
 
253
        if ( ltest ) then
 
254
            if ( ns .ne. 6 ) then
 
255
                print *,'ffcs3: error: only for ns=6, not ',ns
 
256
                stop
 
257
            endif
 
258
        endif
 
259
        ip = ii+3
 
260
        call ffieps(ieps,cz(1),cpi(ip),cpiDpj(ip,ii),isoort)
 
261
        i2 = mod(ii,3) + 1
 
262
        i3 = mod(i2,3) + 1
 
263
*  #] get ieps:
 
264
*  #[ special case |cz| >> |cy|:
 
265
        if ( isoort(2) .ne. 0 .and. max(absc(cy(2)),absc(cy(4))) .lt.
 
266
     +          xloss*min(absc(cz(1)),absc(cz(2)))/2 ) then
 
267
*
 
268
*           we will obtain cancellations of the type Li_2(x) + Li_2(-x)
 
269
*           with x small.
 
270
*
 
271
            if ( lwrite ) print *,'Special case |cz| >> |cy|'
 
272
            cyy = cdyz(2,1)/cd2yzz
 
273
            cyy1 = cdyz(2,2)/cd2yzz
 
274
            if ( absc(cy(2)) .lt. xclogm ) then
 
275
                if ( DIMAG(cy(2)) .eq. 0 .and. abs(DBLE(cy(2))) .gt.
 
276
     +                  xalogm ) then
 
277
                    czz = cz(2)*cyy*DCMPLX(1/DBLE(cy(2)))
 
278
                    cdyyzz = cyy*cdyz(2,2)*DCMPLX(1/DBLE(cy(2)))
 
279
                elseif ( cy(2) .eq. 0 .and. cz(2) .ne. 0 .and. cyy
 
280
     +                  .ne. 0 ) then
 
281
*                   the answer IS zero
 
282
                    goto 30
 
283
                else
 
284
*                   the answer is rounded off to zero
 
285
                    if (lwarn) call ffwarn(42,ier,absc(cy(2)),xclogm)
 
286
                endif
 
287
            else
 
288
                czz = cz(2)*cyy/cy(2)
 
289
                cdyyzz = cyy*cdyz(2,2)/cy(2)
 
290
            endif
 
291
            czz1 = 1-czz
 
292
            if ( lwarn .and. absc(czz) .lt. xloss )
 
293
     +          call ffwarn(43,ier,absc(czz),x1)
 
294
            if ( isoort(1) .eq. -10 ) then
 
295
*               no eta terms.
 
296
                ieps0 = 99
 
297
            else
 
298
*               do not know the im part
 
299
                ieps0 = 0
 
300
            endif
 
301
            call ffcrr(cs3(1),ipi12(1),cyy,cyy1,czz,czz1,cdyyzz,.FALSE.,
 
302
     +          c0,c0,c0,-1,ieps0,ier)
 
303
   30       continue
 
304
            if ( absc(cy(4)) .lt. xclogm ) then
 
305
                if ( DIMAG(cy(4)) .eq. 0 .and. abs(DBLE(cy(4))) .gt.
 
306
     +                  xalogm ) then
 
307
                    czz = cz(4)*cyy*DCMPLX(1/DBLE(cy(4)))
 
308
                    cdyyzz = -cyy*cdyz(2,2)*DCMPLX(1/DBLE(cy(4)))
 
309
                elseif ( cy(4) .eq. 0 .and. cz(4) .ne. 0 .and. cyy
 
310
     +                  .ne. 0 ) then
 
311
*                   the answer IS zero
 
312
                    goto 50
 
313
                else
 
314
*                   the answer is rounded off to zero
 
315
                    if (lwarn) call ffwarn(42,ier,absc(cy(4)),xclogm)
 
316
                endif
 
317
            else
 
318
                czz = cz(4)*cyy/cy(4)
 
319
                cdyyzz = -cyy*cdyz(2,2)/cy(4)
 
320
            endif
 
321
            czz1 = 1-czz
 
322
            if ( lwarn .and. absc(czz) .lt. xloss )
 
323
     +          call ffwarn(43,ier,absc(czz),x1)
 
324
            call ffcrr(cs3(8),ipi12(2),cyy,cyy1,czz,czz1,cdyyzz,.FALSE.,
 
325
     +          c0,c0,c0,-1,ieps0,ier)
 
326
            do 40 i=8,14
 
327
                cs3(i) = -cs3(i)
 
328
   40       continue
 
329
   50       continue
 
330
*
 
331
*           And now the remaining Li_2(x^2) terms
 
332
*           stupid Gould NP1
 
333
*
 
334
            c = cy(2)*cy(2)/(cdyz(2,1)*cdyz(2,1))
 
335
            call ffzli2(zdilog,zlog,c,.FALSE.,ier)
 
336
            cs3(15) = +zdilog/2
 
337
*           stupid Gould NP1
 
338
            c = cy(4)*cy(4)/(cdyz(2,1)*cdyz(2,1))
 
339
            call ffzli2(zdilog,zlog,c,.FALSE.,ier)
 
340
            cs3(16) = -zdilog/2
 
341
            if ( lwrite ) then
 
342
                lwrite = .FALSE.
 
343
                ipi12p(1) = 0
 
344
                ipi12p(2) = 0
 
345
                ier0 = 0
 
346
                do 60 i=1,14
 
347
                    cs3p(i) = 0
 
348
   60           continue
 
349
                call ffcrr(cs3p(1),ipi12p(1),cy(2),cy(4),cz(1),
 
350
     +                  cz(3),cdyz(2,1),.TRUE.,cd2yzz,cz(2),
 
351
     +                  cz(4),isoort(1),ieps(1),ier0)
 
352
                call ffcrr(cs3p(8),ipi12p(2),cy(2),cy(4),cz(2),
 
353
     +                  cz(4),cdyz(2,2),.TRUE.,cd2yzz,cz(1),
 
354
     +                  cz(3),isoort(2),ieps(2),ier0)
 
355
                csum = 0
 
356
                xmax = 0
 
357
                do 70 i=1,14
 
358
                    csum = csum + cs3p(i)
 
359
                    xmax = max(xmax,absc(csum))
 
360
   70           continue
 
361
                csum = csum + (ipi12p(1)+ipi12(2))*DBLE(pi12)
 
362
                print '(a,3g20.10,3i3)','cmp',csum,xmax,ipi12p,ier0
 
363
                lwrite = .TRUE.
 
364
            endif
 
365
            goto 900
 
366
        endif
 
367
*  #] special case |cz| >> |cy|:
 
368
*  #[ normal:
 
369
        if ( isoort(2) .eq. 0 ) then
 
370
            ld2yzz = .FALSE.
 
371
        else
 
372
            ld2yzz = .TRUE.
 
373
        endif
 
374
        if ( isoort(1) .eq. 0 ) then
 
375
*           do nothing
 
376
        elseif ( mod(isoort(1),10).eq.0 .or. mod(isoort(1),10).eq.-1
 
377
     +          .or. mod(isoort(1),10).eq.-3 ) then
 
378
            call ffcrr(cs3(1),ipi12(1),cy(2),cy(4),cz(1),cz(3),
 
379
     +          cdyz(2,1),ld2yzz,cd2yzz,cz(2),cz(4),isoort(1),
 
380
     +          ieps(1),ier)
 
381
        elseif ( mod(isoort(1),10) .eq. -5 .or. mod(isoort(1),10) .eq.
 
382
     +          -6 ) then
 
383
            y = DBLE(cy(2))
 
384
            y1 = DBLE(cy(4))
 
385
            z = DBLE(cz(1))
 
386
            z1 = DBLE(cz(3))
 
387
            dyz = DBLE(cdyz(2,1))
 
388
            d2yzz = DBLE(cd2yzz)
 
389
            zz = DBLE(cz(2))
 
390
            zz1 = DBLE(cz(4))
 
391
            sprec = precx
 
392
            precx = precc
 
393
            call ffcxr(cs3(1),ipi12(1),y,y1,z,z1,dyz,ld2yzz,d2yzz,zz,zz1
 
394
     +                          ,.FALSE.,x00,ieps(1),ier)
 
395
            precx = sprec
 
396
        else
 
397
            call fferr(12,ier)
 
398
        endif
 
399
        if ( isoort(2) .eq. 0 ) then
 
400
*           do nothing
 
401
        elseif ( mod(isoort(2),5) .eq. 0 ) then
 
402
            if ( lwrite ) print *,'ffcs3: skipped next R as it is the ',
 
403
     +          'conjugate'
 
404
            do 100 i=1,7
 
405
  100           cs3(i) = 2*DBLE(cs3(i))
 
406
            ipi12(1) = 2*ipi12(1)
 
407
        elseif ( mod(isoort(2),10).eq.-1 .or. mod(isoort(1),10).eq.-3 )
 
408
     +          then
 
409
            call ffcrr(cs3(8),ipi12(2),cy(2),cy(4),cz(2),cz(4),
 
410
     +          cdyz(2,2),ld2yzz,cd2yzz,cz(1),cz(3),isoort(2),
 
411
     +          ieps(2),ier)
 
412
        elseif ( mod(isoort(2),10) .eq. -6 ) then
 
413
            y = DBLE(cy(2))
 
414
            y1 = DBLE(cy(4))
 
415
            z = DBLE(cz(2))
 
416
            z1 = DBLE(cz(4))
 
417
            dyz = DBLE(cdyz(2,2))
 
418
            d2yzz = DBLE(cd2yzz)
 
419
            zz = DBLE(cz(1))
 
420
            zz1 = DBLE(cz(3))
 
421
            sprec = precx
 
422
            precx = precc
 
423
            call ffcxr(cs3(8),ipi12(2),y,y1,z,z1,dyz,ld2yzz,d2yzz,zz,zz1
 
424
     +                          ,.FALSE.,x00,ieps(2),ier)
 
425
            precx = sprec
 
426
        else
 
427
            call fferr(13,ier)
 
428
        endif
 
429
*  #] normal:
 
430
*  #[ eta's:
 
431
        if ( mod(isoort(1),10).eq.-5 .or. mod(isoort(1),10).eq.-6 )
 
432
     +          then
 
433
            if ( mod(isoort(2),10).ne.-5 .and. mod(isoort(1),10).ne.-6
 
434
     +                  ) then
 
435
                print *,'ffcxs3: error: I assumed both would be real!'
 
436
                ier = ier + 50
 
437
            endif
 
438
*           we called ffcxr - no eta's
 
439
        elseif ( DIMAG(cpi(ip)).eq.0 ) then
 
440
            call ffgeta(ni,cz(1),cdyz(1,1),cd2yzz,
 
441
     +          cpi(ip),cpiDpj(ii,ip),ieps,isoort,ier)
 
442
            if ( lwrite ) print *,'ffcs3: eta''s are ',ni
 
443
            ntot = ni(1) + ni(2) + ni(3) + ni(4)
 
444
            if ( ntot .ne. 0 ) call ffclgy(cs3(15),ipi12(2),ntot,
 
445
     +          cy(1),cz(1),cd2yzz,ier)
 
446
        else
 
447
*
 
448
*           cpi(ip) is really complex (occurs in transformed
 
449
*           4pointfunction)
 
450
*
 
451
            print *,'THIS PART IS NOT READY ',
 
452
     +                  'and should not be reached'
 
453
            stop
 
454
        endif
 
455
*  #] eta's:
 
456
*  #[ print output:
 
457
  900   if (lwrite) then
 
458
            print *,'  cs3 ='
 
459
            do 905 i=1,20
 
460
                if ( cs3(i).ne.0 ) print '(i3,2g20.10,1x)',i,cs3(i)
 
461
  905       continue
 
462
            print '(a3,2g20.10,1x)','pi ',(ipi12(1)+ipi12(2))*pi12
 
463
            print *,'+-----------'
 
464
            csum = 0
 
465
            do 910 i=1,20
 
466
  910       csum = csum + cs3(i)
 
467
            csum = csum+(ipi12(1)+ipi12(2))*DBLE(pi12)
 
468
            print '(a,2g20.10)','Si ',csum
 
469
            print *,'  ipi12,ier=',ipi12,ier
 
470
        endif
 
471
*  #] print output:
 
472
*###] ffcs3:
 
473
        end
 
474
*###[ ffclgy:
 
475
        subroutine ffclgy(cs3,ipi12,ntot,cy,cz,cd2yzz,ier)
 
476
***#[*comment:***********************************************************
 
477
*                                                                       *
 
478
*       calculates the the difference of two S's with cy(3,4),cz(3,4),  *
 
479
*       cy(4)cz(3)-cy(3)cz(4) given.  Note the difference with ffdcs4,  *
 
480
*       in which the cy's are the same and only the cz's different.     *
 
481
*       Here both can be different.     Also we skip an intermediat     *
 
482
*       level.                                                          *
 
483
*                                                                       *
 
484
*       input:  cy(4)        (complex)  cy,1-cy in S with s3,s4         *
 
485
*               cz(4)        (complex)  cz,1-cz in S with s3,s4         *
 
486
*               cdyz(2,2)    (complex)  cy - cz                         *
 
487
*               cd2yzz       (complex)  2*cy - cz+ - cz-                *
 
488
*               cdyzzy(4)    (complex)  cy(i,4)*cz(i,4)-cy(i,3)*cz(i,4) *
 
489
*               cpiDpj(6,6)  (complex)  usual                           *
 
490
*               cs3          (complex)  assumed zero.                   *
 
491
*                                                                       *
 
492
*       output: cs3          (complex)  mod factors pi^2/12, in array   *
 
493
*               ipi12        (integer)  these factors                   *
 
494
*               isoort       (integer)  returns kind of action taken    *
 
495
*               ier          (integer)  number of digits lost           *
 
496
*                                                                       *
 
497
***#]*comment:***********************************************************
 
498
*  #[ declarations:
 
499
        implicit none
 
500
*
 
501
*       arguments
 
502
*
 
503
        DOUBLE COMPLEX cs3
 
504
        DOUBLE COMPLEX cy(4),cz(4),cd2yzz
 
505
        integer ipi12,ntot,ier
 
506
*
 
507
*       local variables
 
508
*
 
509
        integer ipi
 
510
        DOUBLE COMPLEX c,cc,clogy,c2y1,zfflog,zfflo1,csum
 
511
        DOUBLE PRECISION absc
 
512
*
 
513
*       common blocks
 
514
*
 
515
        include 'ff.h'
 
516
*
 
517
*       statement function
 
518
*
 
519
        absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
520
*  #] declarations:
 
521
*  #[ calculations:
 
522
        ipi = 0
 
523
        if ( 1 .lt. xloss*absc(cy(2)) ) then
 
524
            clogy = zfflo1(1/cy(2),ier)
 
525
        else
 
526
            if ( absc(cy(2)) .lt. xclogm .or. absc(cy(4)) .lt. xclogm )
 
527
     +                  then
 
528
                if ( ntot .ne. 0 ) call fferr(15,ier)
 
529
                clogy = 0
 
530
            else
 
531
                c = -cy(4)/cy(2)
 
532
                if ( DBLE(c) .gt. -abs(DIMAG(c)) ) then
 
533
                    clogy = zfflog(c,0,c0,ier)
 
534
                else
 
535
*                   take out the factor 2*pi^2
 
536
                    cc = c+1
 
537
                    if ( absc(cc) .lt. xloss ) then
 
538
                        c2y1 = -cd2yzz - cz(1) + cz(4)
 
539
                        if ( absc(c2y1) .lt. xloss*max(absc(cz(1)),
 
540
     +                          absc(cz(4))) ) then
 
541
                            c2y1 = -cd2yzz - cz(2) + cz(3)
 
542
                            if ( lwarn .and. absc(c2y1) .lt. xloss*max(
 
543
     +                          absc(cz(2)),absc(cz(3))) ) call ffwarn(
 
544
     +                          56,ier,absc(c2y1),absc(cy(2)))
 
545
                        endif
 
546
                        csum = -c2y1/cy(2)
 
547
                        clogy = zfflo1(csum,ier)
 
548
                        if ( lwrite ) then
 
549
                            print *,'c  = ',c
 
550
                            print *,'c+ = ',-1+csum
 
551
                        endif
 
552
                    else
 
553
                        csum = 0
 
554
                        clogy = zfflog(-c,0,c0,ier)
 
555
                    endif
 
556
                    if ( DIMAG(c) .lt. -precc*absc(c) .or.
 
557
     +                   DIMAG(csum) .lt. -precc*absc(csum) ) then
 
558
                        ipi = -1
 
559
                    elseif ( DIMAG(c) .gt. precc*absc(c) .or.
 
560
     +                       DIMAG(csum) .gt. precc*absc(csum) ) then
 
561
                        ipi = +1
 
562
                    else
 
563
                        call fferr(51,ier)
 
564
                        ipi = 0
 
565
                    endif
 
566
                endif
 
567
            endif
 
568
        endif
 
569
        if ( ltest .and. cs3 .ne. 0 ) then
 
570
            print *,'ffclgy: error: cs3 al bezet! ',cs3
 
571
        endif
 
572
        cs3 = cs3 + ntot*c2ipi*clogy
 
573
        if ( ipi .ne. 0 ) then
 
574
            ipi12 = ipi12 - 24*ntot*ipi
 
575
        endif
 
576
*  #] calculations:
 
577
*###] ffclgy:
 
578
        end
 
579
*###[ ffieps:
 
580
        subroutine ffieps(ieps,cz,cp,cpDs,isoort)
 
581
***#[*comment:***********************************************************
 
582
*                                                                       *
 
583
*       Get the ieps prescription in such a way that it is compatible   *
 
584
*       with the imaginary part of cz if non-zero, compatible with the  *
 
585
*       real case if zero.                                              *
 
586
*                                                                       *
 
587
*       Input:  cz      complex(4)      the roots z-,z+,1-z-,1-z+       *
 
588
*               cp      complex         p^2                             *
 
589
*               cpDs    complex         p.s                             *
 
590
*               isoort  integer(2)      which type of Ri                *
 
591
*                                                                       *
 
592
*       Output: ieps    integer(2)      z -> z-ieps*i*epsilon           *
 
593
*                                       will give correct im part       *
 
594
*                                                                       *
 
595
***#]*comment:***********************************************************
 
596
*  #[ declarations:
 
597
        implicit none
 
598
*
 
599
*       arguments:
 
600
*
 
601
        integer ieps(2),isoort(2)
 
602
        DOUBLE COMPLEX cp,cpDs,cz(4)
 
603
*
 
604
*  #] declarations:
 
605
*  #[ work:
 
606
        if ( DIMAG(cp) .ne. 0 ) then
 
607
*           do not calculate ANY eta terms, we'll do them ourselves.
 
608
            ieps(1) = 99
 
609
            ieps(2) = 99
 
610
        elseif ( isoort(2) .ne. 0 ) then
 
611
            if ( DIMAG(cz(1)) .lt. 0 ) then
 
612
                ieps(1) = +1
 
613
                if ( DIMAG(cz(2)) .lt. 0 ) then
 
614
                    ieps(2) = +1
 
615
                else
 
616
                    ieps(2) = -1
 
617
                endif
 
618
            elseif ( DIMAG(cz(1)) .gt. 0 ) then
 
619
                ieps(1) = -1
 
620
                if ( DIMAG(cz(2)) .le. 0 ) then
 
621
                    ieps(2) = +1
 
622
                else
 
623
                    ieps(2) = -1
 
624
                endif
 
625
            else
 
626
                if ( DIMAG(cz(2)) .lt. 0 ) then
 
627
                    ieps(1) = -1
 
628
                    ieps(2) = +1
 
629
                elseif ( DIMAG(cz(2)) .gt. 0 ) then
 
630
                    ieps(1) = +1
 
631
                    ieps(2) = -1
 
632
                else
 
633
                    if ( (DBLE(cz(2)).gt.DBLE(cz(1))
 
634
     +                  .or.  (DBLE(cz(1)).eq.DBLE(cz(2))
 
635
     +                          .and. DBLE(cz(4)).lt.DBLE(cz(3)))
 
636
     +                  ) .eqv. DBLE(cp).gt.0 ) then
 
637
                        ieps(1) = +1
 
638
                        ieps(2) = -1
 
639
                    else
 
640
                        ieps(1) = -1
 
641
                        ieps(2) = +1
 
642
                    endif
 
643
                endif
 
644
            endif
 
645
        else
 
646
            if ( DIMAG(cz(1)) .lt. 0 ) then
 
647
                ieps(1) = +1
 
648
            elseif ( DIMAG(cz(1)) .gt. 0 ) then
 
649
                ieps(1) = -1
 
650
            elseif ( DBLE(cpDs) .gt. 0 ) then
 
651
                ieps(1) = +1
 
652
            else
 
653
                ieps(1) = -1
 
654
            endif
 
655
            ieps(2) = -9999
 
656
        endif
 
657
*  #] work:
 
658
*###] ffieps:
 
659
        end
 
660
*###[ ffgeta:
 
661
        subroutine ffgeta(ni,cz,cdyz,cd2yzz,cp,cpDs,ieps,isoort,ier)
 
662
***#[*comment:***********************************************************
 
663
*                                                                       *
 
664
*       Get the eta terms which arise from splitting up                 *
 
665
*       log(p2(x-z-)(x-z+)) - log(p2(y-z-)(y-z+))                       *
 
666
*                                                                       *
 
667
*       Input:  cz      complex(4)      the roots z-,z+,1-z-,1-z+       *
 
668
*               cdyz    complex(2,2)    y-z                             *
 
669
*               cd2yzz  complex(2)      2y-(z-)-(z+)                    *
 
670
*               cp      complex         p^2                             *
 
671
*               cpDs    complex         p.s                             *
 
672
*               ieps    integer(2)      the assumed im part if Im(z)=0  *
 
673
*               isoort  integer(2)      which type of Ri                *
 
674
*                                                                       *
 
675
*       Output: ni      integer(4)      eta()/(2*pi*i)                  *
 
676
*                                                                       *
 
677
***#]*comment:***********************************************************
 
678
*  #[ declarations:
 
679
        implicit none
 
680
*
 
681
*       arguments:
 
682
*
 
683
        integer ni(4),ieps(2),isoort(2),ier
 
684
        DOUBLE COMPLEX cp,cpDs,cz(4),cdyz(2,2),cd2yzz
 
685
*
 
686
*       local variables
 
687
*
 
688
        integer i,nffeta,nffet1
 
689
        DOUBLE COMPLEX cmip
 
690
*
 
691
*       common
 
692
*
 
693
        include 'ff.h'
 
694
*
 
695
*  #] declarations:
 
696
*  #[ complex  masses or imaginary roots:
 
697
*
 
698
*       only complex because of complex roots in y or z
 
699
*       [checked and in agreement with ieps definition 23-sep-1991]
 
700
*
 
701
        if ( lwrite ) print *,'ffgeta: isoort = ',isoort
 
702
*
 
703
*       isoort = +1:        y is real, z is real
 
704
*       isoort = -1-n*10:   y is complex, possibly z as well
 
705
*       isoort = -3-n*10:   y,z complex, (y-z-)*(y-z+) real
 
706
*       isoort = 0:         y is complex, one z root only
 
707
*       isoort = -10-n*10:  y is real, z is complex
 
708
*       isoort = -5,6-n*10: y,z real
 
709
*
 
710
        if ( isoort(1) .gt. 0 ) then
 
711
*
 
712
*           really a real case
 
713
*
 
714
            ni(1) = 0
 
715
            ni(2) = 0
 
716
            ni(3) = 0
 
717
            ni(4) = 0
 
718
        elseif ( mod(isoort(1),10) .ne. 0 .and. isoort(2) .ne. 0 ) then
 
719
            cmip = DCMPLX(DBLE(x0),-DBLE(cp))
 
720
*
 
721
*           ni(1) = eta(p2,(x-z-)(x-z+)) = 0 by definition (see ni(3))
 
722
*           ni(2) = eta(x-z-,x-z+)
 
723
*
 
724
            ni(1) = 0
 
725
            if ( ieps(1) .gt. 0 .neqv. ieps(2) .gt. 0 ) then
 
726
                ni(2) = 0
 
727
            else
 
728
                ni(2) = nffet1(-cz(1),-cz(2),cmip,ier)
 
729
                if ( cz(3).ne.0 .and. cz(4).ne.0 ) then
 
730
                    i = nffet1(cz(3),cz(4),cmip,ier)
 
731
                    if ( i .ne. ni(2) ) call fferr(53,ier)
 
732
                endif
 
733
            endif
 
734
*
 
735
*           ni(3) compensates for whatever convention we chose in ni(1)
 
736
*           ni(4) = -eta(y-z-,y-z+)
 
737
*
 
738
***         if ( DBLE(cd2yzz).eq.0 .and. ( DIMAG(cz(1)).eq.0 .and.
 
739
***     +                DIMAG(cz(2)).eq.0 .or. DBLE(cdyz(2,1)).eq.0 .and.
 
740
***     +                DBLE(cdyz(2,2)) .eq. 0 ) ) then
 
741
            if ( mod(isoort(1),10).eq.-3 ) then
 
742
*               follow the i*epsilon prescription as (y-z-)(y-z+) real
 
743
                ni(3) = 0
 
744
                if ( ltest ) then
 
745
                    if ( DIMAG(cdyz(2,1)).eq.0 .or. DIMAG(cdyz(2,2))
 
746
     +              .eq.0 ) print *,'ffgeta: error: calling nffet1',
 
747
     +              ' with im(y-z-)=im(y-z+)=0: ',cdyz(2,1),cdyz(2,2)
 
748
                endif
 
749
                ni(4) = -nffet1(cdyz(2,1),cdyz(2,2),cmip,ier)
 
750
            else
 
751
                if ( DBLE(cp) .lt. 0 .and. DIMAG(cdyz(2,1)*
 
752
     +                  cdyz(2,2)) .lt. 0 ) then
 
753
                    ni(3) = -1
 
754
                else
 
755
                    ni(3) = 0
 
756
                endif
 
757
                ni(4) = -nffeta(cdyz(2,1),cdyz(2,2),ier)
 
758
            endif
 
759
        elseif ( (mod(isoort(1),10).eq.-1 .or. mod(isoort(1),10).eq.-3)
 
760
     +          .and. isoort(2) .eq. 0 ) then
 
761
            ni(1) = 0
 
762
            if ( DIMAG(cz(1)) .ne. 0 ) then
 
763
                ni(2) = nffet1(-cpDs,-cz(1),DCMPLX(DBLE(0),
 
764
     +                  DBLE(-1)),ier)
 
765
            else
 
766
                ni(2) = nffet1(-cpDs,DCMPLX(DBLE(0),DBLE(1)),
 
767
     +                  DCMPLX(DBLE(0),DBLE(-1)),ier)
 
768
            endif
 
769
            ni(3) = 0
 
770
            ni(4) = -nffeta(-cpDs,cdyz(2,1),ier)
 
771
        else
 
772
            ni(1) = 0
 
773
            ni(2) = 0
 
774
            ni(3) = 0
 
775
            ni(4) = 0
 
776
        endif
 
777
*  #] complex  masses or imaginary roots:
 
778
*###] ffgeta:
 
779
        end