← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.2

« back to all changes in this revision

Viewing changes to vendor/IREGI/src/qcdloop/ff/ffxdi.f

  • Committer: Marco Zaro
  • Date: 2014-01-27 16:54:10 UTC
  • mfrom: (78.124.55 MG5_aMC_2.1)
  • Revision ID: marco.zaro@gmail.com-20140127165410-5lma8c2hzbzm426j
merged with lp:~maddevelopers/madgraph5/MG5_aMC_2.1 r 267

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/IREGI/src/qcdloop/ff/ffxdi.f
 
1
*###[ ffxdi:
 
2
        subroutine ffxdi(cd4pppp,cd4ppdel,cd4deldel, cd3ppp,cd3pdel,
 
3
     +          cd2pp,cd2del, cd1p, dl2pij, cd0,cc0i,cb0ij,ca0i,
 
4
     +          del4s,del3p,del2pi, xpi,piDpj, d0,xmu, degree, ier)
 
5
***#[*comment:***********************************************************
 
6
*                                                                       *
 
7
*       Compute the tensor functions D1-D(degree) in the determinant    *
 
8
*       scheme, i.e. with basis p1-p3 and (instead of d_(mu,nu))        *
 
9
*       \delta_{p1 p2 p3 mu}^{p1 p2 p3 nu}.                             *
 
10
*                                                                       *
 
11
*       Input:  cd0        (complex)    D0                              *
 
12
*               cc0i(4)    (complex)    C0 with Ni=(Q+..)^2-mi^2 missing*
 
13
*               cb0ij(4,4) (complex)    B0 _with_ Ni,Nj (only for       *
 
14
*                                       degree>1)                       *
 
15
*               ca0i(4)    (complex)    A0 with Ni (only for degree>2)  *
 
16
*               del4s      (real)       delta(s1,s2,s3,s4)(s1,s2,s3,s4) *
 
17
*                                       (only needed when degree>1)     *
 
18
*               del3p      (real)       delta(p1,p2,p3,p1,p2,p3)        *
 
19
*               del2pi(4)  (real)       delta(pipj)(pi,pj) belonging to *
 
20
*                                       cc0i(i)                         *
 
21
*               xpi(13)    (real)       1-4: mi^2, 5-10: p(i-4)^2       *
 
22
*               piDpj(10,10) (re)       pi.pj                           *
 
23
*               d0         (real)       \ renormalization constants     *
 
24
*               xmu        (real)       / used in B0, A0                *
 
25
*               degree     (integer)    1-4                             *
 
26
*               ier        (integer)    number of unreliable digits in  *
 
27
*                                       input                           *
 
28
*                                                                       *
 
29
*       Output: ier                     number of digits lost in the    *
 
30
*                                       least stable result             *
 
31
*               dl2pij(6,6)(real)       determinants delta(pi,pj,pk,pl) *
 
32
*               cd1p(3)    (complex)    coeffs of p1,p2,p3              *
 
33
*          only when degree>1:                                          *
 
34
*               cd2pp(3,3) (complex)    coeffs of p1p1,(p1p2+p2p1),...  *
 
35
*               cd2del     (complex)    coeff of delta(p1,p2,p3,mu,..)  *
 
36
*          only when degree>2:                                          *
 
37
*               cd3ppp(3,3,3)(compl)    coeffs of p1p1p1,p1(p1p2+p2p1), *
 
38
*                                       (p1p2p3+p1p3p2+p2p1p3+p2p3p1+..)*
 
39
*               cd3pdel(3) (complex)    coeffs of pidel (symmetrized)   *
 
40
*          only when degree>3:                                          *
 
41
*               cd4pppp(3,3,3,3)(co)    you guessed it!                 *
 
42
*               cd4ppdel(3,3)(compl)                                    *
 
43
*               cd4deldel  (complex)                                    *
 
44
*                                                                       *
 
45
*       Note: at this moment (28-feb-1993) only D1 and D2 are coded.    *
 
46
*                                                                       *
 
47
***#]*comment:***********************************************************
 
48
*  #[ declarations:
 
49
        implicit none
 
50
*
 
51
*       arguments
 
52
*
 
53
        integer degree,ier
 
54
        DOUBLE PRECISION dl2pij(6,6),del4s,del3p,del2pi(4),xpi(13),
 
55
     +          piDpj(10,10),d0,xmu
 
56
        DOUBLE COMPLEX cd4pppp(3,3,3,3),cd4ppdel(3,3),cd4deldel,
 
57
     +          cd3ppp(3,3,3),cd3pdel(3),cd2pp(3,3),cd2del,
 
58
     +          cd1p(3),cd0,cc0i(4),cb0ij(4,4),ca0i(4)
 
59
*
 
60
*       local variables
 
61
*
 
62
        integer i,j,k,ier0,ier1,ier2,inx43(6,4),sgn43(6,4),i2p(5:8,5:8),
 
63
     +          isgnsa,ii4(6)
 
64
        logical lsave1,lsave2
 
65
        DOUBLE PRECISION a,xpi3(6),xlosn,dl3qi(7),xmax,vgl,xnul
 
66
        DOUBLE COMPLEX cc,cs(25),cnul
 
67
        save inx43,sgn43,i2p,ii4
 
68
*
 
69
*       common blocks
 
70
*
 
71
        include 'ff.h'
 
72
*
 
73
*       data
 
74
*
 
75
        data inx43 /2,3,4,6,7,10,
 
76
     +              1,3,4,9,7,8,
 
77
     +              1,2,4,5,10,8,
 
78
     +              1,2,3,5,6,9/
 
79
        data sgn43 /+1,+1,+1,+1,+1,-1,
 
80
     +              +1,+1,+1,-1,+1,+1,
 
81
     +              +1,+1,+1,+1,+1,+1,
 
82
     +              +1,+1,+1,+1,+1,+1/
 
83
        data i2p /0,0,0,0,
 
84
     +            1,0,0,0,
 
85
     +            2,4,0,0,
 
86
     +            3,5,6,0/
 
87
        data ii4 /5,6,7,8,9,10/
 
88
*
 
89
*  #] declarations:
 
90
*  #[ check input:
 
91
        if ( lwrite ) then
 
92
            print *,'ffxdi: input:'
 
93
            print *,'  degree ',degree
 
94
            print *,'  xpi  = ',xpi
 
95
            print *,'  ier  = ',ier
 
96
        endif
 
97
        if ( degree .gt. 2 ) then
 
98
            print *,'ffxdi: degree > 2 not yet supported: ',degree
 
99
            stop
 
100
        endif
 
101
        if ( del2pi(1).eq.0 .or. del2pi(2).eq.0 .or. del2pi(3).eq.0
 
102
     +          .or. del2pi(4).eq.0 ) then
 
103
            call fferr(87,ier)
 
104
            return
 
105
        endif
 
106
        if ( ltest ) then
 
107
*
 
108
*           the D0
 
109
*
 
110
            ier0 = ier
 
111
            lsave1 = ldot
 
112
            lsave2 = lwrite
 
113
            ldot = .TRUE.
 
114
            lwrite = .FALSE.
 
115
            isgnsa = isgnal
 
116
            call ffxd0(cc,xpi,ier0)
 
117
            isgnal = isgnsa
 
118
            ldot = lsave1
 
119
            lwrite = lsave2
 
120
            xlosn = xloss*DBLE(10)**(-mod(ier0,50))
 
121
            if ( xlosn*abs(cc-cd0) .gt. precc*abs(cd0) ) print *,
 
122
     +          'ffxdi: error: input D0 disagrees with recomputed: ',
 
123
     +          cd0,cc,cd0-cc,ier,ier0
 
124
            if ( xlosn*abs(del3p-fdel3) .gt. precx*abs(del3p) ) print *,
 
125
     +          'ffxdi: error: input del3p disagrees with recomputed: ',
 
126
     +          del3p,fdel3,del3p-fdel3,ier,ier0
 
127
            if ( xlosn*abs(del4s-fdel4s) .gt. precx*abs(del4s) ) print*,
 
128
     +          'ffxdi: error: input del4s disagrees with recomputed: ',
 
129
     +          del4s,fdel4s,del4s-fdel4s,ier,ier0
 
130
            do 20 i=1,10
 
131
                do 10 j=1,10
 
132
                    if ( xlosn*abs(piDpj(j,i)-fpij4(j,i)) .gt. precx*
 
133
     +                  abs(piDpj(j,i)) ) print *,'ffxdi: error: input '
 
134
     +                  ,'piDpj(',j,i,') disagrees with recomputed: ',
 
135
     +                  piDpj(j,i),fpij4(j,i),piDpj(j,i)-fpij4(j,i)
 
136
   10           continue
 
137
   20       continue
 
138
*
 
139
*           the C0s
 
140
*
 
141
            do 40 i=1,4
 
142
                do 30 j=1,6
 
143
                    xpi3(j) = xpi(inx43(j,i))
 
144
   30           continue
 
145
                if ( idot.gt.0 ) then
 
146
                    do 36 j=1,6
 
147
*                       distribute dotproducts
 
148
                        do 35 k=1,6
 
149
                            fpij3(k,j) = fpij4(inx43(k,i),inx43(j,i))*
 
150
     +                          sgn43(k,i)*sgn43(j,i)
 
151
   35                   continue
 
152
   36               continue
 
153
                endif
 
154
                ier0 = ier
 
155
                lsave1 = ldot
 
156
                lsave2 = lwrite
 
157
                ldot = .TRUE.
 
158
                lwrite = .FALSE.
 
159
                call ffxc0(cc,xpi3,ier0)
 
160
                isgnal = isgnsa
 
161
                ldot = lsave1
 
162
                lwrite = lsave2
 
163
                xlosn = xloss*DBLE(10)**(-mod(ier0,50))
 
164
                if ( xlosn*abs(cc-cc0i(i)) .gt. precc*abs(cc0i(i)) )
 
165
     +              print *,'ffxdi: error: input C0(',i,') disagrees ',
 
166
     +              'with recomputed: ',cc0i(i),cc,cc0i(i)-cc,ier,ier0
 
167
                if ( xlosn*abs(del2pi(i)-fdel2) .gt. precx*abs(del2pi(i)
 
168
     +              ) ) print *,'ffxdi: error: input del2pi(',i,
 
169
     +              ') disagrees with recomputed: ',del2pi(i),fdel2,
 
170
     +              del2pi(i)-fdel2
 
171
   40       continue
 
172
*
 
173
*           the B0s
 
174
*
 
175
            if ( degree .lt. 2 ) goto 80
 
176
            do 60 i=1,3
 
177
                do 50 j=i+1,4
 
178
                    ier0 = ier
 
179
                    lsave2 = lwrite
 
180
                    lwrite = .FALSE.
 
181
                    call ffxb0(cc,d0,xmu,xpi(inx(i,j)),xpi(i),xpi(j),
 
182
     +                  ier0)
 
183
                    lwrite = lsave2
 
184
                    xlosn = xloss*DBLE(10)**(-mod(ier0,50))
 
185
                    if ( cb0ij(i,j) .ne. cb0ij(j,i) ) print *,
 
186
     +                  'ffxdi: error: cb0ij(',i,j,') != cb0ij(',j,i,
 
187
     +                  ') : ',cb0ij(i,j),cb0ij(j,i)
 
188
                    if ( xlosn*abs(cc-cb0ij(i,j)) .gt. precc*abs(cb0ij(i
 
189
     +                  ,j)) ) print *,'ffxdi: error: input B0(',i,j,
 
190
     +                  ') disagrees with recomputed: ',cb0ij(i,j),cc,
 
191
     +                  cb0ij(i,j)-cc,ier,ier0
 
192
   50           continue
 
193
   60       continue
 
194
*
 
195
*           the A0s
 
196
*
 
197
            if ( degree .lt. 3 ) goto 80
 
198
            do 70 i=1,4
 
199
                ier0 = ier
 
200
                lsave2 = lwrite
 
201
                lwrite = .FALSE.
 
202
                call ffxa0(cc,d0,xmu,xpi(i),ier0)
 
203
                lwrite = lsave2
 
204
                xlosn = xloss*DBLE(10)**(-mod(ier0,50))
 
205
                if ( xlosn*abs(cc-ca0i(i)) .gt. precc*abs(ca0i(i)) )
 
206
     +              print *,'ffxdi: error: input A0(',i,') disagrees ',
 
207
     +              'with recomputed: ',ca0i(i),cc,ca0i(i)-cc,ier,ier0
 
208
   70       continue
 
209
   80       continue
 
210
        endif
 
211
        if ( .not.ltest ) then
 
212
*           to check when called from ffzfi, ffzei
 
213
            do i=1,10
 
214
                xnul = piDpj(i,5) + piDpj(i,6) + piDpj(i,9)
 
215
                xmax = max(abs(piDpj(i,6)),abs(piDpj(i,9)))
 
216
                if ( xloss*abs(xnul).gt.precx*xmax ) then
 
217
                    print *,'ffxdi: error: i569 does not add up to 0: ',
 
218
     +                  i,piDpj(i,5),piDpj(i,6),piDpj(i,9),xnul,ier
 
219
                endif
 
220
                xnul = piDpj(i,6) + piDpj(i,7) - piDpj(i,10)
 
221
                xmax = max(abs(piDpj(i,7)),abs(piDpj(i,10)))
 
222
                if ( xloss*abs(xnul).gt.precx*xmax ) then
 
223
                    print *,'ffxdi: error: i670 does not add up to 0: ',
 
224
     +                  i,piDpj(i,6),piDpj(i,7),piDpj(i,10),xnul,ier
 
225
                endif
 
226
                xnul = piDpj(i,7) + piDpj(i,8) - piDpj(i,9)
 
227
                xmax = max(abs(piDpj(i,8)),abs(piDpj(i,9)))
 
228
                if ( xloss*abs(xnul).gt.precx*xmax ) then
 
229
                    print *,'ffxdi: error: i789 does not add up to 0: ',
 
230
     +                  i,piDpj(i,7),piDpj(i,8),piDpj(i,9),xnul,ier
 
231
                endif
 
232
                xnul = piDpj(i,8) + piDpj(i,5) + piDpj(i,10)
 
233
                xmax = max(abs(piDpj(i,5)),abs(piDpj(i,10)))
 
234
                if ( xloss*abs(xnul).gt.precx*xmax ) then
 
235
                    print *,'ffxdi: error: i850 does not add up to 0: ',
 
236
     +                  i,piDpj(i,8),piDpj(i,5),piDpj(i,10),xnul,ier
 
237
                endif
 
238
            enddo
 
239
            ier0 = ier
 
240
            call ffdl3p(vgl,piDpj,10,ii4,ii4,ier0)
 
241
            xlosn = xloss*DBLE(10)**(-mod(ier0,50))
 
242
            if ( xlosn*abs(del3p-vgl).gt.precx*abs(vgl) ) then
 
243
                print *,'ffxdi: error: input del3p disagrees with '//
 
244
     +          'recomputed: ',del3p,vgl,del3p-vgl,ier,ier0
 
245
            endif
 
246
            do i=1,4
 
247
                ier0 = ier
 
248
                call ffdel2(vgl,piDpj,10,inx43(4,i),inx43(5,i),
 
249
     +                  inx43(6,i),0,ier0)
 
250
                xlosn = xloss*DBLE(10)**(-mod(ier0,50))
 
251
                if ( xlosn*abs(del2pi(i)-vgl).gt.precx*abs(vgl) ) then
 
252
                    print *,'ffxdi: error: input del2pi(',i,
 
253
     +                  ') disagrees with recomputed: ',del2pi(i),vgl,
 
254
     +                  del2pi(i)-vgl,ier,ier0
 
255
                endif
 
256
            enddo
 
257
        endif
 
258
        if ( degree .le. 0 ) then
 
259
            if ( ltest ) print *,'ffxdi: rather useless call to ffxdi'
 
260
            return
 
261
        endif
 
262
*  #] check input:
 
263
*  #[ preliminaries:
 
264
*       not needed?  security first!
 
265
        if ( lwrite ) then
 
266
            print *,'i2p(5,6) = ',i2p(5,6)
 
267
            print *,'i2p(6,7) = ',i2p(6,7)
 
268
            print *,'i2p(7,8) = ',i2p(7,8)
 
269
            print *,'i2p(5,8) = ',i2p(5,8)
 
270
        endif
 
271
        dl2pij(i2p(5,6),i2p(5,6)) = del2pi(4)
 
272
        dl2pij(i2p(6,7),i2p(6,7)) = del2pi(1)
 
273
        dl2pij(i2p(7,8),i2p(7,8)) = del2pi(2)
 
274
        dl2pij(i2p(5,8),i2p(5,8)) = del2pi(3)
 
275
*  #] preliminaries:
 
276
*  #[ get determinants:
 
277
*
 
278
        ier1 = ier
 
279
        call ffdl2i(dl2pij(i2p(6,7),i2p(7,8)),piDpj,10,
 
280
     +          6,7,10,+1,7,8,9,+1,ier1)
 
281
        dl2pij(i2p(7,8),i2p(6,7)) = dl2pij(i2p(6,7),i2p(7,8))
 
282
*
 
283
        ier0 = ier
 
284
        call ffdl2i(dl2pij(i2p(5,8),i2p(6,7)),piDpj,10,
 
285
     +          6,7,10,+1,5,8,10,-1,ier0)
 
286
        ier1 = max(ier1,ier0)
 
287
        dl2pij(i2p(6,7),i2p(5,8)) = dl2pij(i2p(5,8),i2p(6,7))
 
288
*
 
289
        ier0 = ier
 
290
        call ffdl2i(dl2pij(i2p(5,6),i2p(6,7)),piDpj,10,
 
291
     +          6,7,10,+1,5,6,9,-1,ier0)
 
292
        ier1 = max(ier1,ier0)
 
293
        dl2pij(i2p(6,7),i2p(5,6)) = dl2pij(i2p(5,6),i2p(6,7))
 
294
*
 
295
        ier0 = ier
 
296
        call ffdl2t(dl2pij(i2p(5,7),i2p(6,7)),piDpj,5,7,
 
297
     +          6,7,10,-1,-1, 10,ier0)
 
298
        ier1 = max(ier1,ier0)
 
299
        dl2pij(i2p(6,7),i2p(5,7)) = dl2pij(i2p(5,7),i2p(6,7))
 
300
*
 
301
        ier0 = ier
 
302
        call ffdl2t(dl2pij(i2p(5,7),i2p(7,8)),piDpj,5,7,
 
303
     +          7,8,9,-1,-1, 10,ier0)
 
304
        ier1 = max(ier1,ier0)
 
305
        dl2pij(i2p(7,8),i2p(5,7)) = dl2pij(i2p(5,7),i2p(7,8))
 
306
*
 
307
        ier0 = ier
 
308
        call ffdl2t(dl2pij(i2p(5,7),i2p(5,8)),piDpj,5,7,
 
309
     +          5,8,10,+1,-1, 10,ier0)
 
310
        ier1 = max(ier1,ier0)
 
311
        dl2pij(i2p(5,8),i2p(5,7)) = dl2pij(i2p(5,7),i2p(5,8))
 
312
*
 
313
        ier0 = ier
 
314
        call ffdl2t(dl2pij(i2p(5,6),i2p(5,7)),piDpj,5,7,
 
315
     +          5,6,9,+1,-1, 10,ier0)
 
316
        ier1 = max(ier1,ier0)
 
317
        dl2pij(i2p(5,7),i2p(5,6)) = dl2pij(i2p(5,6),i2p(5,7))
 
318
*
 
319
        ier0 = ier
 
320
        call ffdl2i(dl2pij(i2p(5,6),i2p(7,8)),piDpj,10,
 
321
     +          5,6,9,-1,7,8,9,+1,ier0)
 
322
        ier1 = max(ier1,ier0)
 
323
        dl2pij(i2p(7,8),i2p(5,6)) = dl2pij(i2p(5,6),i2p(7,8))
 
324
*
 
325
        ier0 = ier
 
326
        call ffdl2i(dl2pij(i2p(5,6),i2p(5,8)),piDpj,10,
 
327
     +          5,6,9,-1,5,8,10,-1,ier0)
 
328
        ier1 = max(ier1,ier0)
 
329
        dl2pij(i2p(5,8),i2p(5,6)) = dl2pij(i2p(5,6),i2p(5,8))
 
330
*
 
331
        ier0 = ier
 
332
        call ffdl3q(dl3qi(i2p(6,7)),piDpj, 1,6,7, 0,10,0, 0,-1,0,
 
333
     +          0,+1,0, ier0)
 
334
        ier1 = max(ier1,ier0)
 
335
*
 
336
        ier0 = ier
 
337
        call ffdl3q(dl3qi(i2p(5,7)),piDpj, 1,5,7, 2,0,0, -1,0,0,
 
338
     +          +1,0,0, ier0)
 
339
        ier1 = max(ier1,ier0)
 
340
*
 
341
        ier0 = ier
 
342
        call ffdl3q(dl3qi(i2p(5,6)),piDpj, 1,2,3, 5,6,9, +1,+1,+1,
 
343
     +          -1,-1,-1, ier0)
 
344
        ier1 = max(ier1,ier0)
 
345
*
 
346
        if ( degree.gt.1 ) then
 
347
*
 
348
        ier0 = ier
 
349
        call ffdl3q(dl3qi(i2p(5,8)),piDpj, 1,5,8, 2,10,4, -1,-1,+1,
 
350
     +          +1,-1,+1, ier0)
 
351
        ier1 = max(ier1,ier0)
 
352
*
 
353
        ier0 = ier
 
354
        call ffdl3q(dl3qi(i2p(7,8)),piDpj, 3,4,1, 7,8,10, +1,+1,+1,
 
355
     +          -1,-1,-1, ier0)
 
356
        ier1 = max(ier1,ier0)
 
357
*
 
358
        ier0 = ier
 
359
        call ffdl3q(dl3qi(7),piDpj, 2,3,4, 6,7,10, +1,+1,+1,
 
360
     +          -1,-1,+1, ier0)
 
361
        ier1 = max(ier1,ier0)
 
362
*
 
363
        endif
 
364
        ier = ier1
 
365
        if ( lwrite ) print *,'ier after determinants = ',ier
 
366
*
 
367
*  #] get determinants:
 
368
*  #[ D1:
 
369
*-      #[ D11:
 
370
*
 
371
*       see the Form job D1.frm
 
372
*
 
373
        if ( lwrite ) print *,'ffxdi: D11'
 
374
        cs(1) = - cc0i(1)*DBLE(del2pi(1))
 
375
        cs(2) = + cc0i(2)*DBLE(dl2pij(i2p(6,7),i2p(7,8)))
 
376
        cs(3) = + cc0i(3)*DBLE(dl2pij(i2p(5,8),i2p(6,7)))
 
377
        cs(4) = + cc0i(4)*DBLE(dl2pij(i2p(5,6),i2p(6,7)))
 
378
        cs(5) = + 2*cd0*DBLE(dl3qi(i2p(6,7)))
 
379
*
 
380
        cd1p(1) = 0
 
381
        xmax = 0
 
382
        do 110 i=1,5
 
383
            cd1p(1) = cd1p(1) + cs(i)
 
384
            a = abs(cs(i))
 
385
            xmax = max(xmax,a)
 
386
  110   continue
 
387
        if ( lwarn .and. abs(cd1p(1)) .lt. xloss*xmax ) then
 
388
            a = abs(cd1p(1))
 
389
            call ffwarn(164,ier1,a,xmax)
 
390
            if ( lwrite ) print *,'cs,cd1p(1) = ',(cs(i),i=1,5),cd1p(1)
 
391
        endif
 
392
        cd1p(1) = cd1p(1)*(1/DBLE(2*del3p))
 
393
*
 
394
*-      #] D11:
 
395
*-      #[ D12:
 
396
*
 
397
        if ( lwrite ) print *,'ffxdi: D12'
 
398
        cs(1) = + cc0i(1)*DBLE(dl2pij(i2p(5,7),i2p(6,7)))
 
399
        cs(2) = - cc0i(2)*DBLE(dl2pij(i2p(5,7),i2p(7,8)))
 
400
        cs(3) = - cc0i(3)*DBLE(dl2pij(i2p(5,7),i2p(5,8)))
 
401
        cs(4) = - cc0i(4)*DBLE(dl2pij(i2p(5,6),i2p(5,7)))
 
402
        cs(5) = - 2*cd0*DBLE(dl3qi(i2p(5,7)))
 
403
*
 
404
        cd1p(2) = 0
 
405
        xmax = 0
 
406
        do 120 i=1,5
 
407
            cd1p(2) = cd1p(2) + cs(i)
 
408
            a = abs(cs(i))
 
409
            xmax = max(xmax,a)
 
410
  120   continue
 
411
        if ( lwarn .and. abs(cd1p(2)) .lt. xloss*xmax ) then
 
412
            a = abs(cd1p(2))
 
413
            ier0 = ier
 
414
            call ffwarn(164,ier0,a,xmax)
 
415
            ier1 = max(ier1,ier0)
 
416
            if ( lwrite ) print *,'cs,cd1p(2) = ',(cs(i),i=1,5),cd1p(2)
 
417
        endif
 
418
        cd1p(2) = cd1p(2)*(1/DBLE(2*del3p))
 
419
*
 
420
*-      #] D12:
 
421
*-      #[ D13:
 
422
*
 
423
        if ( lwrite ) print *,'ffxdi: D13'
 
424
        cs(1) = - cc0i(1)*DBLE(dl2pij(i2p(5,6),i2p(6,7)))
 
425
        cs(2) = + cc0i(2)*DBLE(dl2pij(i2p(5,6),i2p(7,8)))
 
426
        cs(3) = + cc0i(3)*DBLE(dl2pij(i2p(5,6),i2p(5,8)))
 
427
        cs(4) = + cc0i(4)*DBLE(del2pi(4))
 
428
        cs(5) = + 2*cd0*DBLE(dl3qi(i2p(5,6)))
 
429
*
 
430
        cd1p(3) = 0
 
431
        xmax = 0
 
432
        do 130 i=1,5
 
433
            cd1p(3) = cd1p(3) + cs(i)
 
434
            a = abs(cs(i))
 
435
            xmax = max(xmax,a)
 
436
  130   continue
 
437
        if ( lwarn .and. abs(cd1p(3)) .lt. xloss*xmax ) then
 
438
            a = abs(cd1p(3))
 
439
            ier0 = ier
 
440
            call ffwarn(164,ier0,a,xmax)
 
441
            ier1 = max(ier1,ier0)
 
442
            if ( lwrite ) print *,'cs,cd1p(3) = ',(cs(i),i=1,5),cd1p(3)
 
443
        endif
 
444
        cd1p(3) = cd1p(3)*(1/DBLE(2*del3p))
 
445
*
 
446
*-      #] D13:
 
447
*-      #[ print output:
 
448
        if ( lwrite ) then
 
449
            print *,'ffxdi: D1:'
 
450
            print *,'cd1p = '
 
451
            print '(6e20.13)',cd1p
 
452
            print *,'ier  = ',ier1
 
453
        endif
 
454
*-      #] print output:
 
455
        if ( degree .eq. 1 ) then
 
456
            ier = ier1
 
457
            return
 
458
        endif
 
459
*  #] D1:
 
460
*  #[ D2:
 
461
*
 
462
*       see the form job d2.frm
 
463
*
 
464
*-      #[ D2del:
 
465
*
 
466
        if ( lwrite ) print *,'ffxdi: D2del'
 
467
        cs(1) = -2*DBLE(del4s)*cd0
 
468
        cs(2) = +DBLE(dl3qi(i2p(5,6)))*cc0i(4)
 
469
        cs(3) = +DBLE(dl3qi(i2p(5,8)))*cc0i(3)
 
470
        cs(4) = +DBLE(dl3qi(i2p(7,8)))*cc0i(2)
 
471
        cs(5) = -DBLE(dl3qi(7))*cc0i(1)
 
472
*
 
473
        cd2del = 0
 
474
        xmax = 0
 
475
        do 210 i=1,5
 
476
            cd2del = cd2del + cs(i)
 
477
            a = abs(cs(i))
 
478
            xmax = max(xmax,a)
 
479
  210   continue
 
480
        if ( lwarn .and. abs(cd2del) .lt. xloss*xmax ) then
 
481
            a = abs(cd2del)
 
482
            ier0 = ier
 
483
            call ffwarn(189,ier0,a,xmax)
 
484
            ier1 = max(ier1,ier0)
 
485
            if ( lwrite ) print *,'cs,cd2del = ',(cs(i),i=1,5),cd2del
 
486
        endif
 
487
        cd2del = cd2del*DBLE(1/(-2*Del3p**2))
 
488
*
 
489
*-      #] D2del:
 
490
*-      #[ D2pp(1,1):
 
491
*
 
492
        if ( lwrite ) print *,'D2pp(1,1)'
 
493
        cs(1) = -cb0ij(1,2)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*piDpj(5,6)*
 
494
     +          del3p/del2pi(4))
 
495
        cs(2) = -cb0ij(1,2)*DBLE(dl2pij(i2p(5,8),i2p(6,7))*piDpj(5,10)*
 
496
     +          del3p/del2pi(3))
 
497
        cs(3) = -cb0ij(1,3)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*piDpj(6,9)*
 
498
     +          del3p/del2pi(4))
 
499
        cs(4) = +cb0ij(1,3)*DBLE(dl2pij(i2p(6,7),i2p(7,8))*piDpj(7,9)*
 
500
     +          del3p/del2pi(2))
 
501
        cs(5) = -cb0ij(1,4)*DBLE(dl2pij(i2p(5,8),i2p(6,7))*piDpj(8,10)*
 
502
     +          del3p/del2pi(3))
 
503
        cs(6) = -cb0ij(1,4)*DBLE(dl2pij(i2p(6,7),i2p(7,8))*piDpj(7,8)*
 
504
     +          del3p/del2pi(2))
 
505
        cs(7) = -cb0ij(2,3)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*piDpj(6,6)*
 
506
     +          del3p/del2pi(4))
 
507
        cs(8) = -cb0ij(2,4)*DBLE(dl2pij(i2p(5,8),i2p(6,7))*piDpj(10,10)*
 
508
     +          del3p/del2pi(3))
 
509
        cs(9) = -cb0ij(3,4)*DBLE(dl2pij(i2p(6,7),i2p(7,8))*piDpj(7,7)*
 
510
     +          del3p/del2pi(2))
 
511
        cs(10) = -4*cc0i(1)*DBLE(dl3qi(i2p(6,7))*del2pi(1))
 
512
        cs(11) = +2*cc0i(1)*DBLE(dl3qi(7)*del2pi(1))
 
513
        cs(12) = -2*cc0i(2)*DBLE(dl2pij(i2p(6,7),i2p(7,8))*
 
514
     +          dl2pij(i2p(6,7),i2p(7,8))*dl3qi(i2p(7,8))/del2pi(2))
 
515
        cs(13) = +4*cc0i(2)*DBLE(dl2pij(i2p(6,7),i2p(7,8))*
 
516
     +          dl3qi(i2p(6,7)))
 
517
        cs(14) = -2*cc0i(3)*DBLE(dl2pij(i2p(5,8),i2p(6,7))*
 
518
     +          dl2pij(i2p(5,8),i2p(6,7))*dl3qi(i2p(5,8))/del2pi(3))
 
519
        cs(15) = +4*cc0i(3)*DBLE(dl2pij(i2p(5,8),i2p(6,7))*
 
520
     +          dl3qi(i2p(6,7)))
 
521
        cs(16) = -2*cc0i(4)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*
 
522
     +          dl2pij(i2p(5,6),i2p(6,7))*dl3qi(i2p(5,6))/del2pi(4))
 
523
        cs(17) = +4*cc0i(4)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*
 
524
     +          dl3qi(i2p(6,7)))
 
525
        cs(18) = +4*cd0*DBLE(dl3qi(i2p(6,7))*dl3qi(i2p(6,7)))
 
526
*
 
527
        cd2pp(1,1) = 0
 
528
        xmax = 0
 
529
        do 220 i=1,18
 
530
            cd2pp(1,1) = cd2pp(1,1) + cs(i)
 
531
            a = abs(cs(i))
 
532
            xmax = max(xmax,a)
 
533
  220   continue
 
534
        if ( lwarn .and. abs(cd2pp(1,1)) .lt. xloss*xmax ) then
 
535
            a = abs(cd2pp(1,1))
 
536
            ier0 = ier
 
537
            call ffwarn(190,ier0,a,xmax)
 
538
            ier1 = max(ier1,ier0)
 
539
            if ( lwrite ) print *,'cs,cd2pp(1,1) = ',(cs(i),i=1,18),
 
540
     +          cd2pp(1,1)
 
541
        endif
 
542
        cd2pp(1,1) = cd2pp(1,1)*DBLE(1/(4*Del3p**2))
 
543
*
 
544
*-      #] D2pp(1,1):
 
545
*-      #[ D2pp(1,2):
 
546
*
 
547
        if ( lwrite ) print *,'D2pp(1,2)'
 
548
        cs(1)=+cb0ij(1,2)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*piDpj(5,
 
549
     +  6)*del3p/del2pi(4))
 
550
        cs(2)=+cb0ij(1,2)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*piDpj(5,
 
551
     +  10)*del3p/del2pi(3))
 
552
        cs(3)=+cb0ij(1,3)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*piDpj(6,
 
553
     +  9)*del3p/del2pi(4))
 
554
        cs(4)=-cb0ij(1,3)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*piDpj(7,
 
555
     +  9)*del3p/del2pi(2))
 
556
        cs(5)=+cb0ij(1,4)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*piDpj(8,
 
557
     +  10)*del3p/del2pi(3))
 
558
        cs(6)=+cb0ij(1,4)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*piDpj(7,
 
559
     +  8)*del3p/del2pi(2))
 
560
        cs(7)=+cb0ij(2,3)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*piDpj(6,
 
561
     +  6)*del3p/del2pi(4))
 
562
        cs(8)=+cb0ij(2,4)*DBLE(dl2pij(i2p(5,8),i2p(6,7))*piDpj(5,
 
563
     +  10)*del3p/del2pi(3))
 
564
        cs(9)=-cb0ij(2,4)*DBLE(piDpj(7,10)*del3p)
 
565
        cs(10)=+cb0ij(3,4)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*piDpj(7,
 
566
     +  7)*del3p/del2pi(2))
 
567
        cs(11)=-2*cc0i(1)*DBLE(dl2pij(i2p(5,7),i2p(6,7))*del3p)
 
568
        cs(12)=+2*cc0i(1)*DBLE(dl3qi(i2p(5,7))*del2pi(1))
 
569
        cs(13)=+2*cc0i(2)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*dl2pij(i2p(6,
 
570
     +  7),i2p(7,8))*dl3qi(i2p(7,8))/del2pi(2))
 
571
        cs(14)=-2*cc0i(2)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*dl3qi(i2p(6,
 
572
     +  7)))
 
573
        cs(15)=-2*cc0i(2)*DBLE(dl2pij(i2p(6,7),i2p(7,8))*dl3qi(i2p(5,
 
574
     +  7)))
 
575
        cs(16)=+2*cc0i(3)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*dl2pij(i2p(5,
 
576
     +  8),i2p(6,7))*dl3qi(i2p(5,8))/del2pi(3))
 
577
        cs(17)=-2*cc0i(3)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*dl3qi(i2p(6,
 
578
     +  7)))
 
579
        cs(18)=-2*cc0i(3)*DBLE(dl2pij(i2p(5,8),i2p(6,7))*dl3qi(i2p(5,
 
580
     +  7)))
 
581
        cs(19)=+2*cc0i(4)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*dl2pij(i2p(5,
 
582
     +  6),i2p(6,7))*dl3qi(i2p(5,6))/del2pi(4))
 
583
        cs(20)=-2*cc0i(4)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*dl3qi(i2p(6,
 
584
     +  7)))
 
585
        cs(21)=-2*cc0i(4)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*dl3qi(i2p(5,
 
586
     +  7)))
 
587
        cs(22)=-4*cd0*DBLE(dl3qi(i2p(5,7))*dl3qi(i2p(6,7)))
 
588
*
 
589
        cd2pp(1,2) = 0
 
590
        xmax = 0
 
591
        do 230 i=1,22
 
592
            cd2pp(1,2) = cd2pp(1,2) + cs(i)
 
593
            a = abs(cs(i))
 
594
            xmax = max(xmax,a)
 
595
  230   continue
 
596
        if ( lwarn .and. abs(cd2pp(1,2)) .lt. xloss*xmax ) then
 
597
            a = abs(cd2pp(1,2))
 
598
            ier0 = ier
 
599
            call ffwarn(190,ier0,a,xmax)
 
600
            ier1 = max(ier1,ier0)
 
601
            if ( lwrite ) print *,'cs,cd2pp(1,2) = ',(cs(i),i=1,22),
 
602
     +          cd2pp(1,2)
 
603
        endif
 
604
        cd2pp(1,2) = cd2pp(1,2)*DBLE(1/(4*Del3p**2))
 
605
        cd2pp(2,1) = cd2pp(1,2)
 
606
*
 
607
*-      #] D2pp(1,2):
 
608
*-      #[ D2pp(1,3):
 
609
*
 
610
        if ( lwrite ) print *,'D2pp(1,3)'
 
611
        cs(1)=-cb0ij(1,2)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*piDpj(5,
 
612
     +  10)*del3p/del2pi(3))
 
613
        cs(2)=-cb0ij(1,2)*DBLE(piDpj(5,6)*del3p)
 
614
        cs(3)=+cb0ij(1,3)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*piDpj(7,
 
615
     +  9)*del3p/del2pi(2))
 
616
        cs(4)=-cb0ij(1,3)*DBLE(piDpj(6,9)*del3p)
 
617
        cs(5)=-cb0ij(1,4)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*piDpj(8,
 
618
     +  10)*del3p/del2pi(3))
 
619
        cs(6)=-cb0ij(1,4)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*piDpj(7,
 
620
     +  8)*del3p/del2pi(2))
 
621
        cs(7)=-cb0ij(2,3)*DBLE(piDpj(6,6)*del3p)
 
622
        cs(8)=-cb0ij(2,4)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*piDpj(10,
 
623
     +  10)*del3p/del2pi(3))
 
624
        cs(9)=-cb0ij(3,4)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*piDpj(7,
 
625
     +  7)*del3p/del2pi(2))
 
626
        cs(10)=+2*cc0i(1)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*del3p)
 
627
        cs(11)=-2*cc0i(1)*DBLE(dl3qi(i2p(5,6))*del2pi(1))
 
628
        cs(12)=-2*cc0i(2)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*dl2pij(i2p(6,
 
629
     +  7),i2p(7,8))*dl3qi(i2p(7,8))/del2pi(2))
 
630
        cs(13)=+2*cc0i(2)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*dl3qi(i2p(6,
 
631
     +  7)))
 
632
        cs(14)=+2*cc0i(2)*DBLE(dl2pij(i2p(6,7),i2p(7,8))*dl3qi(i2p(5,
 
633
     +  6)))
 
634
        cs(15)=-2*cc0i(3)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*dl2pij(i2p(5,
 
635
     +  8),i2p(6,7))*dl3qi(i2p(5,8))/del2pi(3))
 
636
        cs(16)=+2*cc0i(3)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*dl3qi(i2p(6,
 
637
     +  7)))
 
638
        cs(17)=+2*cc0i(3)*DBLE(dl2pij(i2p(5,8),i2p(6,7))*dl3qi(i2p(5,
 
639
     +  6)))
 
640
        cs(18)=+2*cc0i(4)*DBLE(dl3qi(i2p(6,7))*del2pi(4))
 
641
        cs(19)=+4*cd0*DBLE(dl3qi(i2p(5,6))*dl3qi(i2p(6,7)))
 
642
*
 
643
        cd2pp(1,3) = 0
 
644
        xmax = 0
 
645
        do 240 i=1,19
 
646
            cd2pp(1,3) = cd2pp(1,3) + cs(i)
 
647
            a = abs(cs(i))
 
648
            xmax = max(xmax,a)
 
649
  240   continue
 
650
        if ( lwarn .and. abs(cd2pp(1,3)) .lt. xloss*xmax ) then
 
651
            a = abs(cd2pp(1,3))
 
652
            ier0 = ier
 
653
            call ffwarn(190,ier0,a,xmax)
 
654
            ier1 = max(ier1,ier0)
 
655
            if ( lwrite ) print *,'cs,cd2pp(1,3) = ',(cs(i),i=1,19),
 
656
     +          cd2pp(1,3)
 
657
        endif
 
658
        cd2pp(1,3) = cd2pp(1,3)*DBLE(1/(4*Del3p**2))
 
659
        cd2pp(3,1) = cd2pp(1,3)
 
660
*
 
661
*-      #] D2pp(1,3):
 
662
*-      #[ D2pp(2,2):
 
663
*
 
664
        if ( lwrite ) print *,'D2pp(2,2)'
 
665
        cs(1)=-cb0ij(1,2)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*piDpj(5,
 
666
     +  5)*del3p/del2pi(4))
 
667
        cs(2)=-cb0ij(1,2)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*piDpj(5,
 
668
     +  5)*del3p/del2pi(3))
 
669
        cs(3)=-cb0ij(1,3)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*piDpj(5,
 
670
     +  9)*del3p/del2pi(4))
 
671
        cs(4)=-cb0ij(1,3)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*piDpj(7,
 
672
     +  9)*del3p/del2pi(2))
 
673
        cs(5)=-cb0ij(1,4)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*piDpj(5,
 
674
     +  8)*del3p/del2pi(3))
 
675
        cs(6)=+cb0ij(1,4)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*piDpj(7,
 
676
     +  8)*del3p/del2pi(2))
 
677
        cs(7)=-cb0ij(2,3)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*piDpj(5,
 
678
     +  6)*del3p/del2pi(4))
 
679
        cs(8)=-cb0ij(2,3)*DBLE(dl2pij(i2p(5,7),i2p(6,7))*piDpj(6,
 
680
     +  7)*del3p/del2pi(1))
 
681
        cs(9)=-cb0ij(2,4)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*piDpj(5,
 
682
     +  10)*del3p/del2pi(3))
 
683
        cs(10)=+cb0ij(2,4)*DBLE(dl2pij(i2p(5,7),i2p(6,7))*piDpj(7,
 
684
     +  10)*del3p/del2pi(1))
 
685
        cs(11)=-cb0ij(3,4)*DBLE(dl2pij(i2p(5,7),i2p(6,7))*piDpj(7,
 
686
     +  7)*del3p/del2pi(1))
 
687
        cs(12)=+cb0ij(3,4)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*piDpj(7,
 
688
     +  7)*del3p/del2pi(2))
 
689
        cs(13)=+2*cc0i(1)*DBLE(dl2pij(i2p(5,7),i2p(6,7))*dl2pij(i2p(5,
 
690
     +  7),i2p(6,7))*dl3qi(7)/del2pi(1))
 
691
        cs(14)=-4*cc0i(1)*DBLE(dl2pij(i2p(5,7),i2p(6,7))*dl3qi(i2p(5,
 
692
     +  7)))
 
693
        cs(15)=-2*cc0i(2)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*dl2pij(i2p(5,
 
694
     +  7),i2p(7,8))*dl3qi(i2p(7,8))/del2pi(2))
 
695
        cs(16)=+4*cc0i(2)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*dl3qi(i2p(5,
 
696
     +  7)))
 
697
        cs(17)=-2*cc0i(3)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*dl2pij(i2p(5,
 
698
     +  7),i2p(5,8))*dl3qi(i2p(5,8))/del2pi(3))
 
699
        cs(18)=+4*cc0i(3)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*dl3qi(i2p(5,
 
700
     +  7)))
 
701
        cs(19)=-2*cc0i(4)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*dl2pij(i2p(5,
 
702
     +  6),i2p(5,7))*dl3qi(i2p(5,6))/del2pi(4))
 
703
        cs(20)=+4*cc0i(4)*DBLE(dl2pij(i2p(5,6),i2p(5,7))*dl3qi(i2p(5,
 
704
     +  7)))
 
705
        cs(21)=+4*cd0*DBLE(dl3qi(i2p(5,7))*dl3qi(i2p(5,7)))
 
706
*
 
707
        cd2pp(2,2) = 0
 
708
        xmax = 0
 
709
        do 250 i=1,21
 
710
            cd2pp(2,2) = cd2pp(2,2) + cs(i)
 
711
            a = abs(cs(i))
 
712
            xmax = max(xmax,a)
 
713
  250   continue
 
714
        if ( lwarn .and. abs(cd2pp(2,2)) .lt. xloss*xmax ) then
 
715
            a = abs(cd2pp(2,2))
 
716
            ier0 = ier
 
717
            call ffwarn(190,ier0,a,xmax)
 
718
            ier1 = max(ier1,ier0)
 
719
            if ( lwrite ) print *,'cs,cd2pp(2,2) = ',(cs(i),i=1,21),
 
720
     +          cd2pp(2,2)
 
721
        endif
 
722
        cd2pp(2,2) = cd2pp(2,2)*DBLE(1/(4*Del3p**2))
 
723
*
 
724
*-      #] D2pp(2,2):
 
725
*-      #[ D2pp(2,3):
 
726
*
 
727
        if ( lwrite ) print *,'D2pp(2,3)'
 
728
*
 
729
        cs(1)=+cb0ij(1,2)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*piDpj(5,
 
730
     +  5)*del3p/del2pi(3))
 
731
        cs(2)=+cb0ij(1,2)*DBLE(piDpj(5,5)*del3p)
 
732
        cs(3)=+cb0ij(1,3)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*piDpj(7,
 
733
     +  9)*del3p/del2pi(2))
 
734
        cs(4)=+cb0ij(1,3)*DBLE(piDpj(5,9)*del3p)
 
735
        cs(5)=+cb0ij(1,4)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*piDpj(5,
 
736
     +  8)*del3p/del2pi(3))
 
737
        cs(6)=-cb0ij(1,4)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*piDpj(7,
 
738
     +  8)*del3p/del2pi(2))
 
739
        cs(7)=+cb0ij(2,3)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*piDpj(6,
 
740
     +  7)*del3p/del2pi(1))
 
741
        cs(8)=+cb0ij(2,3)*DBLE(piDpj(5,6)*del3p)
 
742
        cs(9)=+cb0ij(2,4)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*piDpj(5,
 
743
     +  10)*del3p/del2pi(3))
 
744
        cs(10)=-cb0ij(2,4)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*piDpj(7,
 
745
     +  10)*del3p/del2pi(1))
 
746
        cs(11)=+cb0ij(3,4)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*piDpj(7,
 
747
     +  7)*del3p/del2pi(1))
 
748
        cs(12)=-cb0ij(3,4)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*piDpj(7,
 
749
     +  7)*del3p/del2pi(2))
 
750
        cs(13)=-2*cc0i(1)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*dl2pij(i2p(5,
 
751
     +  7),i2p(6,7))*dl3qi(7)/del2pi(1))
 
752
        cs(14)=+2*cc0i(1)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*dl3qi(i2p(5,
 
753
     +  7)))
 
754
        cs(15)=+2*cc0i(1)*DBLE(dl2pij(i2p(5,7),i2p(6,7))*dl3qi(i2p(5,
 
755
     +  6)))
 
756
        cs(16)=+2*cc0i(2)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*dl2pij(i2p(5,
 
757
     +  7),i2p(7,8))*dl3qi(i2p(7,8))/del2pi(2))
 
758
        cs(17)=-2*cc0i(2)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*dl3qi(i2p(5,
 
759
     +  7)))
 
760
        cs(18)=-2*cc0i(2)*DBLE(dl2pij(i2p(5,7),i2p(7,8))*dl3qi(i2p(5,
 
761
     +  6)))
 
762
        cs(19)=+2*cc0i(3)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*dl2pij(i2p(5,
 
763
     +  7),i2p(5,8))*dl3qi(i2p(5,8))/del2pi(3))
 
764
        cs(20)=-2*cc0i(3)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*dl3qi(i2p(5,
 
765
     +  7)))
 
766
        cs(21)=-2*cc0i(3)*DBLE(dl2pij(i2p(5,7),i2p(5,8))*dl3qi(i2p(5,
 
767
     +  6)))
 
768
        cs(22)=-2*cc0i(4)*DBLE(dl3qi(i2p(5,7))*del2pi(4))
 
769
        cs(23)=-4*cd0*DBLE(dl3qi(i2p(5,6))*dl3qi(i2p(5,7)))
 
770
*
 
771
        cd2pp(2,3) = 0
 
772
        xmax = 0
 
773
        do 260 i=1,23
 
774
            cd2pp(2,3) = cd2pp(2,3) + cs(i)
 
775
            a = abs(cs(i))
 
776
            xmax = max(xmax,a)
 
777
  260   continue
 
778
        if ( lwarn .and. abs(cd2pp(2,3)) .lt. xloss*xmax ) then
 
779
            a = abs(cd2pp(2,3))
 
780
            ier = ier0
 
781
            call ffwarn(190,ier0,a,xmax)
 
782
            ier1 = max(ier1,ier0)
 
783
            if ( lwrite ) print *,'cs,cd2pp(2,3) = ',(cs(i),i=1,23),
 
784
     +          cd2pp(2,3)
 
785
        endif
 
786
        cd2pp(2,3) = cd2pp(2,3)*DBLE(1/(4*Del3p**2))
 
787
        cd2pp(3,2) = cd2pp(2,3)
 
788
*
 
789
*-      #] D2pp(2,3):
 
790
*-      #[ D2pp(3,3):
 
791
*
 
792
        if ( lwrite ) print *,'D2pp(3,3)'
 
793
        cs(1)=+cb0ij(1,2)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*piDpj(5,
 
794
     +  5)*del3p/del2pi(3))
 
795
        cs(2)=+cb0ij(1,3)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*piDpj(9,
 
796
     +  9)*del3p/del2pi(2))
 
797
        cs(3)=+cb0ij(1,4)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*piDpj(5,
 
798
     +  8)*del3p/del2pi(3))
 
799
        cs(4)=-cb0ij(1,4)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*piDpj(8,
 
800
     +  9)*del3p/del2pi(2))
 
801
        cs(5)=-cb0ij(2,3)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*piDpj(6,
 
802
     +  6)*del3p/del2pi(1))
 
803
        cs(6)=+cb0ij(2,4)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*piDpj(5,
 
804
     +  10)*del3p/del2pi(3))
 
805
        cs(7)=+cb0ij(2,4)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*piDpj(6,
 
806
     +  10)*del3p/del2pi(1))
 
807
        cs(8)=-cb0ij(3,4)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*piDpj(6,
 
808
     +  7)*del3p/del2pi(1))
 
809
        cs(9)=-cb0ij(3,4)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*piDpj(7,
 
810
     +  9)*del3p/del2pi(2))
 
811
        cs(10)=+2*cc0i(1)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*dl2pij(i2p(5,
 
812
     +  6),i2p(6,7))*dl3qi(7)/del2pi(1))
 
813
        cs(11)=-4*cc0i(1)*DBLE(dl2pij(i2p(5,6),i2p(6,7))*dl3qi(i2p(5,
 
814
     +  6)))
 
815
        cs(12)=-2*cc0i(2)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*dl2pij(i2p(5,
 
816
     +  6),i2p(7,8))*dl3qi(i2p(7,8))/del2pi(2))
 
817
        cs(13)=+4*cc0i(2)*DBLE(dl2pij(i2p(5,6),i2p(7,8))*dl3qi(i2p(5,
 
818
     +  6)))
 
819
        cs(14)=-2*cc0i(3)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*dl2pij(i2p(5,
 
820
     +  6),i2p(5,8))*dl3qi(i2p(5,8))/del2pi(3))
 
821
        cs(15)=+4*cc0i(3)*DBLE(dl2pij(i2p(5,6),i2p(5,8))*dl3qi(i2p(5,
 
822
     +  6)))
 
823
        cs(16)=+2*cc0i(4)*DBLE(dl3qi(i2p(5,6))*del2pi(4))
 
824
        cs(17)=+4*cd0*DBLE(dl3qi(i2p(5,6))*dl3qi(i2p(5,6)))
 
825
*
 
826
        cd2pp(3,3) = 0
 
827
        xmax = 0
 
828
        do 270 i=1,17
 
829
            cd2pp(3,3) = cd2pp(3,3) + cs(i)
 
830
            a = abs(cs(i))
 
831
            xmax = max(xmax,a)
 
832
  270   continue
 
833
        if ( lwarn .and. abs(cd2pp(3,3)) .lt. xloss*xmax ) then
 
834
            a = abs(cd2pp(3,3))
 
835
            ier0 = ier
 
836
            call ffwarn(190,ier0,a,xmax)
 
837
            ier1 = max(ier1,ier0)
 
838
            if ( lwrite ) print *,'cs,cd2pp(3,3) = ',(cs(i),i=1,17),
 
839
     +          cd2pp(3,3)
 
840
        endif
 
841
        cd2pp(3,3) = cd2pp(3,3)*DBLE(1/(4*Del3p**2))
 
842
*
 
843
*-      #] D2pp(3,3):
 
844
*-      #[ print output:
 
845
        if ( lwrite ) then
 
846
            print '(a,2e20.13)','cd2del = ',cd2del
 
847
            print '(a)','cd2pp  = '
 
848
            print '(6e20.13)',cd2pp
 
849
            print *,'ier    = ',ier1
 
850
        endif
 
851
        if ( ltest ) then
 
852
            xlosn = xloss*DBLE(10)**(-2-mod(ier1,50))
 
853
            cs(1) =   DBLE(piDpj(5,5))*cd2pp(1,1)
 
854
            cs(2) = 2*DBLE(piDpj(5,6))*cd2pp(1,2)
 
855
            cs(3) = 2*DBLE(piDpj(5,7))*cd2pp(1,3)
 
856
            cs(4) =   DBLE(piDpj(6,6))*cd2pp(2,2)
 
857
            cs(5) = 2*DBLE(piDpj(6,7))*cd2pp(2,3)
 
858
            cs(6) =   DBLE(piDpj(7,7))*cd2pp(3,3)
 
859
            cs(7) =   DBLE(del3p)*cd2del
 
860
            cs(8) = - cc0i(1)
 
861
            cs(9) = - DBLE(piDpj(1,1))*cd0
 
862
            cnul = 0
 
863
            xmax = 0
 
864
            do 910 i=1,9
 
865
                cnul = cnul + cs(i)
 
866
                a = abs(cs(i))
 
867
                xmax = max(xmax,a)
 
868
  910       continue
 
869
            if ( lwrite ) print *,'ffxdi: checking D2.gmumu= ',cnul,xmax
 
870
            if ( xlosn*abs(cnul) .gt. precc*xmax ) print *,'ffxdi: ',
 
871
     +          'error: D2(mu,mu) not correct ',cnul,xmax,ier1
 
872
            cs(1) = 4*DBLE(piDpj(5,5)*piDpj(7,5))*cd2pp(1,1)
 
873
            cs(2) = 4*DBLE(piDpj(5,5)*piDpj(7,6))*cd2pp(1,2)
 
874
            cs(3) = 4*DBLE(piDpj(5,6)*piDpj(7,5))*cd2pp(1,2)
 
875
            cs(4) = 4*DBLE(piDpj(5,5)*piDpj(7,7))*cd2pp(1,3)
 
876
            cs(5) = 4*DBLE(piDpj(5,7)*piDpj(7,5))*cd2pp(1,3)
 
877
            cs(6) = 4*DBLE(piDpj(5,6)*piDpj(7,6))*cd2pp(2,2)
 
878
            cs(7) = 4*DBLE(piDpj(5,6)*piDpj(7,7))*cd2pp(2,3)
 
879
            cs(8) = 4*DBLE(piDpj(5,7)*piDpj(7,6))*cd2pp(2,3)
 
880
            cs(9) = 4*DBLE(piDpj(5,7)*piDpj(7,7))*cd2pp(3,3)
 
881
            cs(10)= - cb0ij(1,3)
 
882
            cs(11)= + cb0ij(1,4)
 
883
            cs(12)= + cb0ij(2,3)
 
884
            cs(13)= - cb0ij(2,4)
 
885
            cs(14)= - 2*DBLE(piDpj(1,7))*cc0i(2)
 
886
            cs(15)= + 2*DBLE(piDpj(1,7))*cc0i(1)
 
887
            cs(16)= - 2*DBLE(piDpj(1,5))*cc0i(4)
 
888
            cs(17)= + 2*DBLE(piDpj(1,5))*cc0i(3)
 
889
            cs(18)= - 4*DBLE(piDpj(1,5)*piDpj(1,7))*cd0
 
890
            cnul = 0
 
891
            xmax = 0
 
892
            do 920 i=1,18
 
893
                cnul = cnul + cs(i)
 
894
                a = abs(cs(i))
 
895
                xmax = max(xmax,a)
 
896
  920       continue
 
897
            if ( lwrite ) print *,'ffxdi: checking D2.p1p3 = ',cnul,xmax
 
898
            if ( xlosn*abs(cnul) .gt. precc*xmax ) print *,'ffxdi :',
 
899
     +          'error: D2(p1,p3) not correct ',cnul,xmax,ier1
 
900
            cs(1) = 4*DBLE(piDpj(6,5)*piDpj(8,5))*cd2pp(1,1)
 
901
            cs(2) = 4*DBLE(piDpj(6,5)*piDpj(8,6))*cd2pp(1,2)
 
902
            cs(3) = 4*DBLE(piDpj(6,6)*piDpj(8,5))*cd2pp(1,2)
 
903
            cs(4) = 4*DBLE(piDpj(6,5)*piDpj(8,7))*cd2pp(1,3)
 
904
            cs(5) = 4*DBLE(piDpj(6,7)*piDpj(8,5))*cd2pp(1,3)
 
905
            cs(6) = 4*DBLE(piDpj(6,6)*piDpj(8,6))*cd2pp(2,2)
 
906
            cs(7) = 4*DBLE(piDpj(6,6)*piDpj(8,7))*cd2pp(2,3)
 
907
            cs(8) = 4*DBLE(piDpj(6,7)*piDpj(8,6))*cd2pp(2,3)
 
908
            cs(9) = 4*DBLE(piDpj(6,7)*piDpj(8,7))*cd2pp(3,3)
 
909
            cs(10)= - cb0ij(2,4)
 
910
            cs(11)= + cb0ij(1,2)
 
911
            cs(12)= + cb0ij(3,4)
 
912
            cs(13)= - cb0ij(1,3)
 
913
            cs(14)= - 2*DBLE(piDpj(1,8))*cc0i(3)
 
914
            cs(15)= + 2*DBLE(piDpj(1,8))*cc0i(2)
 
915
            cs(16)= - 2*DBLE(piDpj(1,6))*cc0i(1)
 
916
            cs(17)= + 2*DBLE(piDpj(1,6))*cc0i(4)
 
917
            cs(18)= - 4*DBLE(piDpj(1,6)*piDpj(1,8))*cd0
 
918
            cnul = 0
 
919
            xmax = 0
 
920
            do 930 i=1,18
 
921
                cnul = cnul + cs(i)
 
922
                a = abs(cs(i))
 
923
                xmax = max(xmax,a)
 
924
  930       continue
 
925
            if ( lwrite ) print *,'ffxdi: checking D2.p2p4 = ',cnul,xmax
 
926
            if ( xlosn*abs(cnul) .gt. precc*xmax ) print *,'ffxdi :',
 
927
     +          'error: D2(p2,p4) not correct ',cnul,xmax,ier1
 
928
        endif
 
929
*-      #] print output:
 
930
        if ( degree .eq. 2 ) then
 
931
            ier = ier1
 
932
            return
 
933
        endif
 
934
*  #] D2:
 
935
        print *,'ffxdi: error: D3 not ready'
 
936
        stop
 
937
*###] ffxdi:
 
938
        end