← Back to branch summary

~alifson/chiralityflow/trunk

« back to all changes in this revision

Viewing changes to tests/input_files/IOTestsComparison/short_ML_SMQCD_default/gg_ttx/%..%..%Source%MODEL%model_functions.f

  • Committer: andrew.lifson at lu
  • Date: 2021-09-01 15:34:39 UTC
  • Revision ID: andrew.lifson@thep.lu.se-20210901153439-7fasjhav4cp4m88r
testing a new repository of a madgraph folder

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/input_files/IOTestsComparison/short_ML_SMQCD_default/gg_ttx/%..%..%Source%MODEL%model_functions.f
 
1
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
2
c      written by the UFO converter
 
3
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
4
 
 
5
      DOUBLE COMPLEX FUNCTION COND(CONDITION,TRUECASE,FALSECASE)
 
6
      IMPLICIT NONE
 
7
      DOUBLE COMPLEX CONDITION,TRUECASE,FALSECASE
 
8
      IF(CONDITION.EQ.(0.0D0,0.0D0)) THEN
 
9
        COND=TRUECASE
 
10
      ELSE
 
11
        COND=FALSECASE
 
12
      ENDIF
 
13
      END
 
14
 
 
15
      DOUBLE COMPLEX FUNCTION CONDIF(CONDITION,TRUECASE,FALSECASE)
 
16
      IMPLICIT NONE
 
17
      LOGICAL CONDITION
 
18
      DOUBLE COMPLEX TRUECASE,FALSECASE
 
19
      IF(CONDITION) THEN
 
20
        CONDIF=TRUECASE
 
21
      ELSE
 
22
        CONDIF=FALSECASE
 
23
      ENDIF
 
24
      END
 
25
 
 
26
      DOUBLE COMPLEX FUNCTION RECMS(CONDITION,EXPR)
 
27
      IMPLICIT NONE
 
28
      LOGICAL CONDITION
 
29
      DOUBLE COMPLEX EXPR
 
30
      IF(CONDITION)THEN
 
31
        RECMS=EXPR
 
32
      ELSE
 
33
        RECMS=DCMPLX(DBLE(EXPR))
 
34
      ENDIF
 
35
      END
 
36
 
 
37
      DOUBLE COMPLEX FUNCTION REGLOG(ARG_IN)
 
38
      IMPLICIT NONE
 
39
      DOUBLE COMPLEX TWOPII
 
40
      PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
 
41
      DOUBLE COMPLEX ARG_IN
 
42
      DOUBLE COMPLEX ARG
 
43
      ARG=ARG_IN
 
44
      IF(DABS(DIMAG(ARG)).EQ.0.0D0)THEN
 
45
        ARG=DCMPLX(DBLE(ARG),0.0D0)
 
46
      ENDIF
 
47
      IF(DABS(DBLE(ARG)).EQ.0.0D0)THEN
 
48
        ARG=DCMPLX(0.0D0,DIMAG(ARG))
 
49
      ENDIF
 
50
      IF(ARG.EQ.(0.0D0,0.0D0)) THEN
 
51
        REGLOG=(0.0D0,0.0D0)
 
52
      ELSE
 
53
        REGLOG=LOG(ARG)
 
54
      ENDIF
 
55
      END
 
56
 
 
57
      DOUBLE COMPLEX FUNCTION REGLOGP(ARG_IN)
 
58
      IMPLICIT NONE
 
59
      DOUBLE COMPLEX TWOPII
 
60
      PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
 
61
      DOUBLE COMPLEX ARG_IN
 
62
      DOUBLE COMPLEX ARG
 
63
      ARG=ARG_IN
 
64
      IF(DABS(DIMAG(ARG)).EQ.0.0D0)THEN
 
65
        ARG=DCMPLX(DBLE(ARG),0.0D0)
 
66
      ENDIF
 
67
      IF(DABS(DBLE(ARG)).EQ.0.0D0)THEN
 
68
        ARG=DCMPLX(0.0D0,DIMAG(ARG))
 
69
      ENDIF
 
70
      IF(ARG.EQ.(0.0D0,0.0D0))THEN
 
71
        REGLOGP=(0.0D0,0.0D0)
 
72
      ELSE
 
73
        IF(DBLE(ARG).LT.0.0D0.AND.DIMAG(ARG).LT.0.0D0)THEN
 
74
          REGLOGP=LOG(ARG) + TWOPII
 
75
        ELSE
 
76
          REGLOGP=LOG(ARG)
 
77
        ENDIF
 
78
      ENDIF
 
79
      END
 
80
 
 
81
      DOUBLE COMPLEX FUNCTION REGLOGM(ARG_IN)
 
82
      IMPLICIT NONE
 
83
      DOUBLE COMPLEX TWOPII
 
84
      PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
 
85
      DOUBLE COMPLEX ARG_IN
 
86
      DOUBLE COMPLEX ARG
 
87
      ARG=ARG_IN
 
88
      IF(DABS(DIMAG(ARG)).EQ.0.0D0)THEN
 
89
        ARG=DCMPLX(DBLE(ARG),0.0D0)
 
90
      ENDIF
 
91
      IF(DABS(DBLE(ARG)).EQ.0.0D0)THEN
 
92
        ARG=DCMPLX(0.0D0,DIMAG(ARG))
 
93
      ENDIF
 
94
      IF(ARG.EQ.(0.0D0,0.0D0))THEN
 
95
        REGLOGM=(0.0D0,0.0D0)
 
96
      ELSE
 
97
        IF(DBLE(ARG).LT.0.0D0.AND.DIMAG(ARG).GT.0.0D0)THEN
 
98
          REGLOGM=LOG(ARG) - TWOPII
 
99
        ELSE
 
100
          REGLOGM=LOG(ARG)
 
101
        ENDIF
 
102
      ENDIF
 
103
      END
 
104
 
 
105
      DOUBLE COMPLEX FUNCTION REGSQRT(ARG_IN)
 
106
      IMPLICIT NONE
 
107
      DOUBLE COMPLEX ARG_IN
 
108
      DOUBLE COMPLEX ARG
 
109
      ARG=ARG_IN
 
110
      IF(DABS(DIMAG(ARG)).EQ.0.0D0)THEN
 
111
        ARG=DCMPLX(DBLE(ARG),0.0D0)
 
112
      ENDIF
 
113
      IF(DABS(DBLE(ARG)).EQ.0.0D0)THEN
 
114
        ARG=DCMPLX(0.0D0,DIMAG(ARG))
 
115
      ENDIF
 
116
      REGSQRT=SQRT(ARG)
 
117
      END
 
118
 
 
119
      DOUBLE COMPLEX FUNCTION GRREGLOG(LOGSW,EXPR1_IN,EXPR2_IN)
 
120
      IMPLICIT NONE
 
121
      DOUBLE COMPLEX TWOPII
 
122
      PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
 
123
      DOUBLE COMPLEX EXPR1_IN,EXPR2_IN
 
124
      DOUBLE COMPLEX EXPR1,EXPR2
 
125
      DOUBLE PRECISION LOGSW
 
126
      DOUBLE PRECISION IMAGEXPR
 
127
      LOGICAL FIRSTSHEET
 
128
      EXPR1=EXPR1_IN
 
129
      EXPR2=EXPR2_IN
 
130
      IF(DABS(DIMAG(EXPR1)).EQ.0.0D0)THEN
 
131
        EXPR1=DCMPLX(DBLE(EXPR1),0.0D0)
 
132
      ENDIF
 
133
      IF(DABS(DBLE(EXPR1)).EQ.0.0D0)THEN
 
134
        EXPR1=DCMPLX(0.0D0,DIMAG(EXPR1))
 
135
      ENDIF
 
136
      IF(DABS(DIMAG(EXPR2)).EQ.0.0D0)THEN
 
137
        EXPR2=DCMPLX(DBLE(EXPR2),0.0D0)
 
138
      ENDIF
 
139
      IF(DABS(DBLE(EXPR2)).EQ.0.0D0)THEN
 
140
        EXPR2=DCMPLX(0.0D0,DIMAG(EXPR2))
 
141
      ENDIF
 
142
      IF(EXPR1.EQ.(0.0D0,0.0D0))THEN
 
143
        GRREGLOG=(0.0D0,0.0D0)
 
144
      ELSE
 
145
        IMAGEXPR=DIMAG(EXPR1)*DIMAG(EXPR2)
 
146
        FIRSTSHEET=IMAGEXPR.GE.0.0D0
 
147
        FIRSTSHEET=FIRSTSHEET.OR.DBLE(EXPR1).GE.0.0D0
 
148
        FIRSTSHEET=FIRSTSHEET.OR.DBLE(EXPR2).GE.0.0D0
 
149
        IF(FIRSTSHEET)THEN
 
150
          GRREGLOG=LOG(EXPR1)
 
151
        ELSE
 
152
          IF(DIMAG(EXPR1).GT.0.0D0)THEN
 
153
            GRREGLOG=LOG(EXPR1) - LOGSW*TWOPII
 
154
          ELSE
 
155
            GRREGLOG=LOG(EXPR1) + LOGSW*TWOPII
 
156
          ENDIF
 
157
        ENDIF
 
158
      ENDIF
 
159
      END
 
160
 
 
161
      MODULE B0F_CACHING
 
162
 
 
163
      TYPE B0F_NODE
 
164
        DOUBLE COMPLEX P2,M12,M22
 
165
        DOUBLE COMPLEX VALUE
 
166
        TYPE(B0F_NODE),POINTER::PARENT
 
167
        TYPE(B0F_NODE),POINTER::LEFT
 
168
        TYPE(B0F_NODE),POINTER::RIGHT
 
169
        END TYPE B0F_NODE
 
170
 
 
171
        CONTAINS
 
172
 
 
173
        SUBROUTINE B0F_SEARCH(ITEM, HEAD, FIND)
 
174
        IMPLICIT NONE
 
175
        TYPE(B0F_NODE),POINTER,INTENT(INOUT)::HEAD,ITEM
 
176
        LOGICAL,INTENT(OUT)::FIND
 
177
        TYPE(B0F_NODE),POINTER::ITEM1
 
178
        INTEGER::ICOMP
 
179
        FIND=.FALSE.
 
180
        NULLIFY(ITEM%%PARENT)
 
181
        NULLIFY(ITEM%%LEFT)
 
182
        NULLIFY(ITEM%%RIGHT)
 
183
        IF(.NOT.ASSOCIATED(HEAD))THEN
 
184
          HEAD => ITEM
 
185
          RETURN
 
186
        ENDIF
 
187
        ITEM1 => HEAD
 
188
        DO
 
189
        ICOMP=B0F_NODE_COMPARE(ITEM,ITEM1)
 
190
        IF(ICOMP.LT.0)THEN
 
191
          IF(.NOT.ASSOCIATED(ITEM1%%LEFT))THEN
 
192
            ITEM1%%LEFT => ITEM
 
193
            ITEM%%PARENT => ITEM1
 
194
            EXIT
 
195
          ELSE
 
196
            ITEM1 => ITEM1%%LEFT
 
197
          ENDIF
 
198
        ELSEIF(ICOMP.GT.0)THEN
 
199
          IF(.NOT.ASSOCIATED(ITEM1%%RIGHT))THEN
 
200
            ITEM1%%RIGHT => ITEM
 
201
            ITEM%%PARENT => ITEM1
 
202
            EXIT
 
203
          ELSE
 
204
            ITEM1 => ITEM1%%RIGHT
 
205
          ENDIF
 
206
        ELSE
 
207
          FIND=.TRUE.
 
208
          ITEM%%VALUE=ITEM1%%VALUE
 
209
          EXIT
 
210
        ENDIF
 
211
        ENDDO
 
212
        RETURN
 
213
        END
 
214
 
 
215
        INTEGER FUNCTION B0F_NODE_COMPARE(ITEM1,ITEM2) RESULT(RES)
 
216
        IMPLICIT NONE
 
217
        TYPE(B0F_NODE),POINTER,INTENT(IN)::ITEM1,ITEM2
 
218
        RES=COMPLEX_COMPARE(ITEM1%%P2,ITEM2%%P2)
 
219
        IF(RES.NE.0)RETURN
 
220
        RES=COMPLEX_COMPARE(ITEM1%%M22,ITEM2%%M22)
 
221
        IF(RES.NE.0)RETURN
 
222
        RES=COMPLEX_COMPARE(ITEM1%%M12,ITEM2%%M12)
 
223
        RETURN
 
224
        END
 
225
 
 
226
        INTEGER FUNCTION REAL_COMPARE(R1,R2) RESULT(RES)
 
227
        IMPLICIT NONE
 
228
        DOUBLE PRECISION R1,R2
 
229
        DOUBLE PRECISION MAXR,DIFF
 
230
        DOUBLE PRECISION TINY
 
231
        PARAMETER (TINY=-1D-14)
 
232
        MAXR=MAX(ABS(R1),ABS(R2))
 
233
        DIFF=R1-R2
 
234
        IF(MAXR.LE.1D-99.OR.ABS(DIFF)/MAX(MAXR,1D-99).LE.ABS(TINY))THEN
 
235
          RES=0
 
236
          RETURN
 
237
        ENDIF
 
238
        IF(DIFF.GT.0D0)THEN
 
239
          RES=1
 
240
          RETURN
 
241
        ELSE
 
242
          RES=-1
 
243
          RETURN
 
244
        ENDIF
 
245
        END
 
246
 
 
247
        INTEGER FUNCTION COMPLEX_COMPARE(C1,C2) RESULT(RES)
 
248
        IMPLICIT NONE
 
249
        DOUBLE COMPLEX C1,C2
 
250
        DOUBLE PRECISION R1,R2
 
251
        R1=DBLE(C1)
 
252
        R2=DBLE(C2)
 
253
        RES=REAL_COMPARE(R1,R2)
 
254
        IF(RES.NE.0)RETURN
 
255
        R1=DIMAG(C1)
 
256
        R2=DIMAG(C2)
 
257
        RES=REAL_COMPARE(R1,R2)
 
258
        RETURN
 
259
        END
 
260
 
 
261
        END MODULE B0F_CACHING
 
262
 
 
263
        DOUBLE COMPLEX FUNCTION B0F(P2,M12,M22)
 
264
        USE B0F_CACHING
 
265
        IMPLICIT NONE
 
266
        DOUBLE COMPLEX P2,M12,M22
 
267
        DOUBLE COMPLEX ZERO,TWOPII
 
268
        PARAMETER (ZERO=(0.0D0,0.0D0))
 
269
        PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
 
270
        DOUBLE PRECISION M,M2,GA,GA2
 
271
        DOUBLE PRECISION TINY
 
272
        PARAMETER (TINY=-1D-14)
 
273
        DOUBLE COMPLEX LOGTERMS
 
274
        DOUBLE COMPLEX LOG_TRAJECTORY
 
275
        LOGICAL USE_CACHING
 
276
        PARAMETER (USE_CACHING=.TRUE.)
 
277
        TYPE(B0F_NODE),POINTER::ITEM
 
278
        TYPE(B0F_NODE),POINTER,SAVE::B0F_BT
 
279
        INTEGER INIT
 
280
        SAVE INIT
 
281
        DATA INIT /0/
 
282
        LOGICAL FIND
 
283
        IF(M12.EQ.ZERO)THEN
 
284
C         it is a special case
 
285
C         refer to Eq.(5.48) in arXiv:1804.10017
 
286
          M=DBLE(P2)  ! M^2
 
287
          M2=DBLE(M22)  ! M2^2
 
288
          IF(M.LT.TINY.OR.M2.LT.TINY)THEN
 
289
            WRITE(*,*)'ERROR:B0F is not well defined when M^2,M2^2<0'
 
290
            STOP
 
291
          ENDIF
 
292
          M=DSQRT(DABS(M))
 
293
          M2=DSQRT(DABS(M2))
 
294
          IF(M.EQ.0D0)THEN
 
295
            GA=0D0
 
296
          ELSE
 
297
            GA=-DIMAG(P2)/M
 
298
          ENDIF
 
299
          IF(M2.EQ.0D0)THEN
 
300
            GA2=0D0
 
301
          ELSE
 
302
            GA2=-DIMAG(M22)/M2
 
303
          ENDIF
 
304
          IF(P2.NE.M22.AND.P2.NE.ZERO.AND.M22.NE.ZERO)THEN
 
305
            B0F=(M22-P2)/P2*LOG((M22-P2)/M22)
 
306
            IF(M.GT.M2.AND.GA*M2.GT.GA2*M)THEN
 
307
              B0F=B0F-TWOPII
 
308
            ENDIF
 
309
            RETURN
 
310
          ELSE
 
311
            WRITE(*,*)'ERROR:B0F is not supported for a simple form'
 
312
            STOP
 
313
          ENDIF
 
314
        ENDIF
 
315
C       the general case
 
316
C       trajectory method as advocated in arXiv:1804.10017 (Eq.(E.47))
 
317
        IF(USE_CACHING)THEN
 
318
          IF(INIT.EQ.0)THEN
 
319
            NULLIFY(B0F_BT)
 
320
            INIT=1
 
321
          ENDIF
 
322
          ALLOCATE(ITEM)
 
323
          ITEM%%P2=P2
 
324
          ITEM%%M12=M12
 
325
          ITEM%%M22=M22
 
326
          FIND=.FALSE.
 
327
          CALL B0F_SEARCH(ITEM,B0F_BT,FIND)
 
328
          IF(FIND)THEN
 
329
            B0F=ITEM%%VALUE
 
330
            DEALLOCATE(ITEM)
 
331
            RETURN
 
332
          ELSE
 
333
            LOGTERMS=LOG_TRAJECTORY(100,P2,M12,M22)
 
334
            B0F=-LOG(P2/M22)+LOGTERMS
 
335
            ITEM%%VALUE=B0F
 
336
            RETURN
 
337
          ENDIF
 
338
        ELSE
 
339
          LOGTERMS=LOG_TRAJECTORY(100,P2,M12,M22)
 
340
          B0F=-LOG(P2/M22)+LOGTERMS
 
341
        ENDIF
 
342
        RETURN
 
343
        END
 
344
 
 
345
        DOUBLE COMPLEX FUNCTION SQRT_TRAJECTORY(N_SEG,P2,M12,M22)
 
346
C       only needed when p2*m12*m22=\=0
 
347
        IMPLICIT NONE
 
348
        INTEGER N_SEG  ! number of segments
 
349
        DOUBLE COMPLEX P2,M12,M22
 
350
        DOUBLE COMPLEX ZERO,ONE
 
351
        PARAMETER (ZERO=(0.0D0,0.0D0),ONE=(1.0D0,0.0D0))
 
352
        DOUBLE COMPLEX GAMMA0,GAMMA1
 
353
        DOUBLE PRECISION M,GA,DGA,GA_START
 
354
        DOUBLE PRECISION GAI,INTERSECTION
 
355
        DOUBLE COMPLEX ARGIM1,ARGI,P2I
 
356
        DOUBLE COMPLEX GAMMA0I,GAMMA1I
 
357
        DOUBLE PRECISION TINY
 
358
        PARAMETER (TINY=-1D-24)
 
359
        INTEGER I
 
360
        DOUBLE PRECISION PREFACTOR
 
361
        IF(ABS(P2*M12*M22).EQ.0D0)THEN
 
362
          WRITE(*,*)'ERROR:sqrt_trajectory works when p2*m12*m22/=0'
 
363
          STOP
 
364
        ENDIF
 
365
        M=DBLE(P2)  ! M^2
 
366
        M=DSQRT(DABS(M))
 
367
        IF(M.EQ.0D0)THEN
 
368
          GA=0D0
 
369
        ELSE
 
370
          GA=-DIMAG(P2)/M
 
371
        ENDIF
 
372
C       Eq.(5.37) in arXiv:1804.10017
 
373
        GAMMA0=ONE+M12/P2-M22/P2
 
374
        GAMMA1=M12/P2-DCMPLX(0D0,1D0)*ABS(TINY)/P2
 
375
        IF(ABS(GA).EQ.0D0)THEN
 
376
          SQRT_TRAJECTORY=SQRT(GAMMA0**2-4D0*GAMMA1)
 
377
          RETURN
 
378
        ENDIF
 
379
C       segments from -DABS(tiny*Ga) to Ga
 
380
        GA_START=-DABS(TINY*GA)
 
381
        DGA=(GA-GA_START)/N_SEG
 
382
        PREFACTOR=1D0
 
383
        GAI=GA_START
 
384
        P2I=DCMPLX(M**2,-GAI*M)
 
385
        GAMMA0I=ONE+M12/P2I-M22/P2I
 
386
        GAMMA1I=M12/P2I-DCMPLX(0D0,1D0)*ABS(TINY)/P2I
 
387
        ARGIM1=GAMMA0I**2-4D0*GAMMA1I
 
388
        DO I=1,N_SEG
 
389
          GAI=DGA*I+GA_START
 
390
          P2I=DCMPLX(M**2,-GAI*M)
 
391
          GAMMA0I=ONE+M12/P2I-M22/P2I
 
392
          GAMMA1I=M12/P2I-DCMPLX(0D0,1D0)*ABS(TINY)/P2I
 
393
          ARGI=GAMMA0I**2-4D0*GAMMA1I
 
394
          IF(DIMAG(ARGI)*DIMAG(ARGIM1).LT.0D0)THEN
 
395
            INTERSECTION=DIMAG(ARGIM1)*(DBLE(ARGI)-DBLE(ARGIM1))
 
396
            INTERSECTION=INTERSECTION/(DIMAG(ARGI)-DIMAG(ARGIM1))
 
397
            INTERSECTION=INTERSECTION-DBLE(ARGIM1)
 
398
            IF(INTERSECTION.GT.0D0)THEN
 
399
              PREFACTOR=-PREFACTOR
 
400
            ENDIF
 
401
          ENDIF
 
402
          ARGIM1=ARGI
 
403
        ENDDO
 
404
        SQRT_TRAJECTORY=SQRT(GAMMA0**2-4D0*GAMMA1)*PREFACTOR
 
405
        RETURN
 
406
        END
 
407
 
 
408
        DOUBLE COMPLEX FUNCTION LOG_TRAJECTORY(N_SEG,P2,M12,M22)
 
409
C       sum of log terms appearing in Eq.(5.35) of arXiv:1804.10017
 
410
C       only needed when p2*m12*m22=\=0
 
411
        IMPLICIT NONE
 
412
C       4 possible logarithms appearing in Eq.(5.35) of
 
413
C        arXiv:1804.10017
 
414
C       log(arg(i)) with arg(i) for i=1 to 4
 
415
C       i=1: (ga_{+}-1)
 
416
C       i=2: (ga_{-}-1)
 
417
C       i=3: (ga_{+}-1)/ga_{+}
 
418
C       i=4: (ga_{-}-1)/ga_{-}
 
419
        INTEGER N_SEG  ! number of segments
 
420
        DOUBLE COMPLEX P2,M12,M22
 
421
        DOUBLE COMPLEX ZERO,ONE,HALF,TWOPII
 
422
        PARAMETER (ZERO=(0.0D0,0.0D0),ONE=(1.0D0,0.0D0))
 
423
        PARAMETER (HALF=(0.5D0,0.0D0))
 
424
        PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
 
425
        DOUBLE COMPLEX GAMMA0,GAMMAP,GAMMAM,SQRTTERM
 
426
        DOUBLE PRECISION M,GA,DGA,GA_START
 
427
        DOUBLE PRECISION GAI,INTERSECTION
 
428
        DOUBLE COMPLEX ARGIM1(4),ARGI(4),P2I,SQRTTERMI
 
429
        DOUBLE COMPLEX GAMMA0I,GAMMAPI,GAMMAMI
 
430
        DOUBLE PRECISION TINY
 
431
        PARAMETER (TINY=-1D-14)
 
432
        INTEGER I,J
 
433
        DOUBLE COMPLEX ADDFACTOR(4)
 
434
        DOUBLE COMPLEX SQRT_TRAJECTORY
 
435
        IF(ABS(P2*M12*M22).EQ.0D0)THEN
 
436
          WRITE(*,*)'ERROR:log_trajectory works when p2*m12*m22/=0'
 
437
          STOP
 
438
        ENDIF
 
439
        M=DBLE(P2)  ! M^2
 
440
        M=DSQRT(DABS(M))
 
441
        IF(M.EQ.0D0)THEN
 
442
          GA=0D0
 
443
        ELSE
 
444
          GA=-DIMAG(P2)/M
 
445
        ENDIF
 
446
C       Eq.(5.36-5.38) in arXiv:1804.10017
 
447
        SQRTTERM=SQRT_TRAJECTORY(N_SEG,P2,M12,M22)
 
448
        GAMMA0=ONE+M12/P2-M22/P2
 
449
        GAMMAP=HALF*(GAMMA0+SQRTTERM)
 
450
        GAMMAM=HALF*(GAMMA0-SQRTTERM)
 
451
        IF(ABS(GA).EQ.0D0)THEN
 
452
          LOG_TRAJECTORY=-LOG(GAMMAP-ONE)-LOG(GAMMAM-ONE)+GAMMAP
 
453
     $     *LOG((GAMMAP-ONE)/GAMMAP)+GAMMAM*LOG((GAMMAM-ONE)/GAMMAM)
 
454
          RETURN
 
455
        ENDIF
 
456
C       segments from -DABS(tiny*Ga) to Ga
 
457
        GA_START=-DABS(TINY*GA)
 
458
        DGA=(GA-GA_START)/N_SEG
 
459
        ADDFACTOR(1:4)=ZERO
 
460
        GAI=GA_START
 
461
        P2I=DCMPLX(M**2,-GAI*M)
 
462
        SQRTTERMI=SQRT_TRAJECTORY(N_SEG,P2I,M12,M22)
 
463
        GAMMA0I=ONE+M12/P2I-M22/P2I
 
464
        GAMMAPI=HALF*(GAMMA0I+SQRTTERMI)
 
465
        GAMMAMI=HALF*(GAMMA0I-SQRTTERMI)
 
466
        ARGIM1(1)=GAMMAPI-ONE
 
467
        ARGIM1(2)=GAMMAMI-ONE
 
468
        ARGIM1(3)=(GAMMAPI-ONE)/GAMMAPI
 
469
        ARGIM1(4)=(GAMMAMI-ONE)/GAMMAMI
 
470
        DO I=1,N_SEG
 
471
          GAI=DGA*I+GA_START
 
472
          P2I=DCMPLX(M**2,-GAI*M)
 
473
          SQRTTERMI=SQRT_TRAJECTORY(N_SEG,P2I,M12,M22)
 
474
          GAMMA0I=ONE+M12/P2I-M22/P2I
 
475
          GAMMAPI=HALF*(GAMMA0I+SQRTTERMI)
 
476
          GAMMAMI=HALF*(GAMMA0I-SQRTTERMI)
 
477
          ARGI(1)=GAMMAPI-ONE
 
478
          ARGI(2)=GAMMAMI-ONE
 
479
          ARGI(3)=(GAMMAPI-ONE)/GAMMAPI
 
480
          ARGI(4)=(GAMMAMI-ONE)/GAMMAMI
 
481
          DO J=1,4
 
482
            IF(DIMAG(ARGI(J))*DIMAG(ARGIM1(J)).LT.0D0)THEN
 
483
              INTERSECTION=DIMAG(ARGIM1(J))*(DBLE(ARGI(J))
 
484
     $         -DBLE(ARGIM1(J)))
 
485
              INTERSECTION=INTERSECTION/(DIMAG(ARGI(J))-DIMAG(ARGIM1(J)
 
486
     $         ))
 
487
              INTERSECTION=INTERSECTION-DBLE(ARGIM1(J))
 
488
              IF(INTERSECTION.GT.0D0)THEN
 
489
                IF(DIMAG(ARGIM1(J)).LT.0)THEN
 
490
                  ADDFACTOR(J)=ADDFACTOR(J)-TWOPII
 
491
                ELSE
 
492
                  ADDFACTOR(J)=ADDFACTOR(J)+TWOPII
 
493
                ENDIF
 
494
              ENDIF
 
495
            ENDIF
 
496
            ARGIM1(J)=ARGI(J)
 
497
          ENDDO
 
498
        ENDDO
 
499
        LOG_TRAJECTORY=-(LOG(GAMMAP-ONE)+ADDFACTOR(1))-(LOG(GAMMAM-ONE)
 
500
     $   +ADDFACTOR(2))
 
501
        LOG_TRAJECTORY=LOG_TRAJECTORY+GAMMAP*(LOG((GAMMAP-ONE)/GAMMAP)
 
502
     $   +ADDFACTOR(3))
 
503
        LOG_TRAJECTORY=LOG_TRAJECTORY+GAMMAM*(LOG((GAMMAM-ONE)/GAMMAM)
 
504
     $   +ADDFACTOR(4))
 
505
        RETURN
 
506
        END
 
507
 
 
508
        DOUBLE COMPLEX FUNCTION ARG(COMNUM)
 
509
        IMPLICIT NONE
 
510
        DOUBLE COMPLEX COMNUM
 
511
        DOUBLE COMPLEX IIM
 
512
        IIM = (0.0D0,1.0D0)
 
513
        IF(COMNUM.EQ.(0.0D0,0.0D0)) THEN
 
514
          ARG=(0.0D0,0.0D0)
 
515
        ELSE
 
516
          ARG=LOG(COMNUM/ABS(COMNUM))/IIM
 
517
        ENDIF
 
518
        END
 
519
 
 
520
 
 
521
        COMPLEX*32 FUNCTION MP_COND(CONDITION,TRUECASE,FALSECASE)
 
522
        IMPLICIT NONE
 
523
        COMPLEX*32 CONDITION,TRUECASE,FALSECASE
 
524
        IF(CONDITION.EQ.(0.0E0_16,0.0E0_16)) THEN
 
525
          MP_COND=TRUECASE
 
526
        ELSE
 
527
          MP_COND=FALSECASE
 
528
        ENDIF
 
529
        END
 
530
 
 
531
        COMPLEX*32 FUNCTION MP_CONDIF(CONDITION,TRUECASE,FALSECASE)
 
532
        IMPLICIT NONE
 
533
        LOGICAL CONDITION
 
534
        COMPLEX*32 TRUECASE,FALSECASE
 
535
        IF(CONDITION) THEN
 
536
          MP_CONDIF=TRUECASE
 
537
        ELSE
 
538
          MP_CONDIF=FALSECASE
 
539
        ENDIF
 
540
        END
 
541
 
 
542
        COMPLEX*32 FUNCTION MP_RECMS(CONDITION,EXPR)
 
543
        IMPLICIT NONE
 
544
        LOGICAL CONDITION
 
545
        COMPLEX*32 EXPR
 
546
        IF(CONDITION)THEN
 
547
          MP_RECMS=EXPR
 
548
        ELSE
 
549
          MP_RECMS=CMPLX(REAL(EXPR),KIND=16)
 
550
        ENDIF
 
551
        END
 
552
 
 
553
 
 
554
        COMPLEX*32 FUNCTION MP_REGLOG(ARG_IN)
 
555
        IMPLICIT NONE
 
556
        COMPLEX*32 TWOPII
 
557
        PARAMETER (TWOPII=2.0E0_16
 
558
     $   *3.14169258478796109557151794433593750E0_16*(0.0E0_16
 
559
     $   ,1.0E0_16))
 
560
        COMPLEX*32 ARG_IN
 
561
        COMPLEX*32 ARG
 
562
        ARG=ARG_IN
 
563
        IF(ABS(IMAGPART(ARG)).EQ.0.0E0_16)THEN
 
564
          ARG=CMPLX(REAL(ARG,KIND=16),0.0E0_16)
 
565
        ENDIF
 
566
        IF(ABS(REAL(ARG,KIND=16)).EQ.0.0E0_16)THEN
 
567
          ARG=CMPLX(0.0E0_16,IMAGPART(ARG))
 
568
        ENDIF
 
569
        IF(ARG.EQ.(0.0E0_16,0.0E0_16)) THEN
 
570
          MP_REGLOG=(0.0E0_16,0.0E0_16)
 
571
        ELSE
 
572
          MP_REGLOG=LOG(ARG)
 
573
        ENDIF
 
574
        END
 
575
 
 
576
        COMPLEX*32 FUNCTION MP_REGLOGP(ARG_IN)
 
577
        IMPLICIT NONE
 
578
        COMPLEX*32 TWOPII
 
579
        PARAMETER (TWOPII=2.0E0_16
 
580
     $   *3.14169258478796109557151794433593750E0_16*(0.0E0_16
 
581
     $   ,1.0E0_16))
 
582
        COMPLEX*32 ARG_IN
 
583
        COMPLEX*32 ARG
 
584
        ARG=ARG_IN
 
585
        IF(ABS(IMAGPART(ARG)).EQ.0.0E0_16)THEN
 
586
          ARG=CMPLX(REAL(ARG,KIND=16),0.0E0_16)
 
587
        ENDIF
 
588
        IF(ABS(REAL(ARG,KIND=16)).EQ.0.0E0_16)THEN
 
589
          ARG=CMPLX(0.0E0_16,IMAGPART(ARG))
 
590
        ENDIF
 
591
        IF(ARG.EQ.(0.0E0_16,0.0E0_16))THEN
 
592
          MP_REGLOGP=(0.0E0_16,0.0E0_16)
 
593
        ELSE
 
594
          IF(REAL(ARG,KIND=16).LT.0.0E0_16.AND.IMAGPART(ARG)
 
595
     $     .LT.0.0E0_16)THEN
 
596
            MP_REGLOGP=LOG(ARG) + TWOPII
 
597
          ELSE
 
598
            MP_REGLOGP=LOG(ARG)
 
599
          ENDIF
 
600
        ENDIF
 
601
        END
 
602
 
 
603
        COMPLEX*32 FUNCTION MP_REGLOGM(ARG_IN)
 
604
        IMPLICIT NONE
 
605
        COMPLEX*32 TWOPII
 
606
        PARAMETER (TWOPII=2.0E0_16
 
607
     $   *3.14169258478796109557151794433593750E0_16*(0.0E0_16
 
608
     $   ,1.0E0_16))
 
609
        COMPLEX*32 ARG_IN
 
610
        COMPLEX*32 ARG
 
611
        ARG=ARG_IN
 
612
        IF(ABS(IMAGPART(ARG)).EQ.0.0E0_16)THEN
 
613
          ARG=CMPLX(REAL(ARG,KIND=16),0.0E0_16)
 
614
        ENDIF
 
615
        IF(ABS(REAL(ARG,KIND=16)).EQ.0.0E0_16)THEN
 
616
          ARG=CMPLX(0.0E0_16,IMAGPART(ARG))
 
617
        ENDIF
 
618
        IF(ARG.EQ.(0.0E0_16,0.0E0_16))THEN
 
619
          MP_REGLOGM=(0.0E0_16,0.0E0_16)
 
620
        ELSE
 
621
          IF(REAL(ARG,KIND=16).LT.0.0E0_16.AND.IMAGPART(ARG)
 
622
     $     .GT.0.0E0_16)THEN
 
623
            MP_REGLOGM=LOG(ARG) - TWOPII
 
624
          ELSE
 
625
            MP_REGLOGM=LOG(ARG)
 
626
          ENDIF
 
627
        ENDIF
 
628
        END
 
629
 
 
630
        COMPLEX*32 FUNCTION MP_REGSQRT(ARG_IN)
 
631
        IMPLICIT NONE
 
632
        COMPLEX*32 ARG_IN
 
633
        COMPLEX*32 ARG
 
634
        ARG=ARG_IN
 
635
        IF(ABS(IMAGPART(ARG)).EQ.0.0E0_16)THEN
 
636
          ARG=CMPLX(REAL(ARG,KIND=16),0.0E0_16)
 
637
        ENDIF
 
638
        IF(ABS(REAL(ARG,KIND=16)).EQ.0.0E0_16)THEN
 
639
          ARG=CMPLX(0.0E0_16,IMAGPART(ARG))
 
640
        ENDIF
 
641
        MP_REGSQRT=SQRT(ARG)
 
642
        END
 
643
 
 
644
        COMPLEX*32 FUNCTION MP_GRREGLOG(LOGSW,EXPR1_IN,EXPR2_IN)
 
645
        IMPLICIT NONE
 
646
        COMPLEX*32 TWOPII
 
647
        PARAMETER (TWOPII=2.0E0_16
 
648
     $   *3.14169258478796109557151794433593750E0_16*(0.0E0_16
 
649
     $   ,1.0E0_16))
 
650
        COMPLEX*32 EXPR1_IN,EXPR2_IN
 
651
        COMPLEX*32 EXPR1,EXPR2
 
652
        REAL*16 LOGSW
 
653
        REAL*16 IMAGEXPR
 
654
        LOGICAL FIRSTSHEET
 
655
        EXPR1=EXPR1_IN
 
656
        EXPR2=EXPR2_IN
 
657
        IF(ABS(IMAGPART(EXPR1)).EQ.0.0E0_16)THEN
 
658
          EXPR1=CMPLX(REAL(EXPR1,KIND=16),0.0E0_16)
 
659
        ENDIF
 
660
        IF(ABS(REAL(EXPR1,KIND=16)).EQ.0.0E0_16)THEN
 
661
          EXPR1=CMPLX(0.0E0_16,IMAGPART(EXPR1))
 
662
        ENDIF
 
663
        IF(ABS(IMAGPART(EXPR2)).EQ.0.0E0_16)THEN
 
664
          EXPR2=CMPLX(REAL(EXPR2,KIND=16),0.0E0_16)
 
665
        ENDIF
 
666
        IF(ABS(REAL(EXPR2,KIND=16)).EQ.0.0E0_16)THEN
 
667
          EXPR2=CMPLX(0.0E0_16,IMAGPART(EXPR2))
 
668
        ENDIF
 
669
        IF(EXPR1.EQ.(0.0E0_16,0.0E0_16))THEN
 
670
          MP_GRREGLOG=(0.0E0_16,0.0E0_16)
 
671
        ELSE
 
672
          IMAGEXPR=IMAGPART(EXPR1)*IMAGPART(EXPR2)
 
673
          FIRSTSHEET=IMAGEXPR.GE.0.0E0_16
 
674
          FIRSTSHEET=FIRSTSHEET.OR.REAL(EXPR1,KIND=16).GE.0.0E0_16
 
675
          FIRSTSHEET=FIRSTSHEET.OR.REAL(EXPR2,KIND=16).GE.0.0E0_16
 
676
          IF(FIRSTSHEET)THEN
 
677
            MP_GRREGLOG=LOG(EXPR1)
 
678
          ELSE
 
679
            IF(IMAGPART(EXPR1).GT.0.0E0_16)THEN
 
680
              MP_GRREGLOG=LOG(EXPR1) - LOGSW*TWOPII
 
681
            ELSE
 
682
              MP_GRREGLOG=LOG(EXPR1) + LOGSW*TWOPII
 
683
            ENDIF
 
684
          ENDIF
 
685
        ENDIF
 
686
        END
 
687
 
 
688
        MODULE MP_B0F_CACHING
 
689
 
 
690
        TYPE MP_B0F_NODE
 
691
          COMPLEX*32 P2,M12,M22
 
692
          COMPLEX*32 VALUE
 
693
          TYPE(MP_B0F_NODE),POINTER::PARENT
 
694
          TYPE(MP_B0F_NODE),POINTER::LEFT
 
695
          TYPE(MP_B0F_NODE),POINTER::RIGHT
 
696
          END TYPE MP_B0F_NODE
 
697
 
 
698
          CONTAINS
 
699
 
 
700
          SUBROUTINE MP_B0F_SEARCH(ITEM, HEAD, FIND)
 
701
          IMPLICIT NONE
 
702
          TYPE(MP_B0F_NODE),POINTER,INTENT(INOUT)::HEAD,ITEM
 
703
          LOGICAL,INTENT(OUT)::FIND
 
704
          TYPE(MP_B0F_NODE),POINTER::ITEM1
 
705
          INTEGER::ICOMP
 
706
          FIND=.FALSE.
 
707
          NULLIFY(ITEM%%PARENT)
 
708
          NULLIFY(ITEM%%LEFT)
 
709
          NULLIFY(ITEM%%RIGHT)
 
710
          IF(.NOT.ASSOCIATED(HEAD))THEN
 
711
            HEAD => ITEM
 
712
            RETURN
 
713
          ENDIF
 
714
          ITEM1 => HEAD
 
715
          DO
 
716
          ICOMP=MP_B0F_NODE_COMPARE(ITEM,ITEM1)
 
717
          IF(ICOMP.LT.0)THEN
 
718
            IF(.NOT.ASSOCIATED(ITEM1%%LEFT))THEN
 
719
              ITEM1%%LEFT => ITEM
 
720
              ITEM%%PARENT => ITEM1
 
721
              EXIT
 
722
            ELSE
 
723
              ITEM1 => ITEM1%%LEFT
 
724
            ENDIF
 
725
          ELSEIF(ICOMP.GT.0)THEN
 
726
            IF(.NOT.ASSOCIATED(ITEM1%%RIGHT))THEN
 
727
              ITEM1%%RIGHT => ITEM
 
728
              ITEM%%PARENT => ITEM1
 
729
              EXIT
 
730
            ELSE
 
731
              ITEM1 => ITEM1%%RIGHT
 
732
            ENDIF
 
733
          ELSE
 
734
            FIND=.TRUE.
 
735
            ITEM%%VALUE=ITEM1%%VALUE
 
736
            EXIT
 
737
          ENDIF
 
738
          ENDDO
 
739
          RETURN
 
740
          END
 
741
 
 
742
          INTEGER FUNCTION MP_B0F_NODE_COMPARE(ITEM1,ITEM2) RESULT(RES)
 
743
          IMPLICIT NONE
 
744
          TYPE(MP_B0F_NODE),POINTER,INTENT(IN)::ITEM1,ITEM2
 
745
          RES=MP_COMPLEX_COMPARE(ITEM1%%P2,ITEM2%%P2)
 
746
          IF(RES.NE.0)RETURN
 
747
          RES=MP_COMPLEX_COMPARE(ITEM1%%M22,ITEM2%%M22)
 
748
          IF(RES.NE.0)RETURN
 
749
          RES=MP_COMPLEX_COMPARE(ITEM1%%M12,ITEM2%%M12)
 
750
          RETURN
 
751
          END
 
752
 
 
753
          INTEGER FUNCTION MP_REAL_COMPARE(R1,R2) RESULT(RES)
 
754
          IMPLICIT NONE
 
755
          REAL*16 R1,R2
 
756
          REAL*16 MAXR,DIFF
 
757
          REAL*16 TINY
 
758
          PARAMETER (TINY=-1.0E-14_16)
 
759
          MAXR=MAX(ABS(R1),ABS(R2))
 
760
          DIFF=R1-R2
 
761
          IF(MAXR.LE.1.0E-99_16.OR.ABS(DIFF)/MAX(MAXR,1.0E-99_16)
 
762
     $     .LE.ABS(TINY))THEN
 
763
            RES=0
 
764
            RETURN
 
765
          ENDIF
 
766
          IF(DIFF.GT.0.0E0_16)THEN
 
767
            RES=1
 
768
            RETURN
 
769
          ELSE
 
770
            RES=-1
 
771
            RETURN
 
772
          ENDIF
 
773
          END
 
774
 
 
775
          INTEGER FUNCTION MP_COMPLEX_COMPARE(C1,C2) RESULT(RES)
 
776
          IMPLICIT NONE
 
777
          COMPLEX*32 C1,C2
 
778
          REAL*16 R1,R2
 
779
          R1=REAL(C1,KIND=16)
 
780
          R2=REAL(C2,KIND=16)
 
781
          RES=MP_REAL_COMPARE(R1,R2)
 
782
          IF(RES.NE.0)RETURN
 
783
          R1=IMAGPART(C1)
 
784
          R2=IMAGPART(C2)
 
785
          RES=MP_REAL_COMPARE(R1,R2)
 
786
          RETURN
 
787
          END
 
788
 
 
789
          END MODULE MP_B0F_CACHING
 
790
 
 
791
          COMPLEX*32 FUNCTION MP_B0F(P2,M12,M22)
 
792
          USE MP_B0F_CACHING
 
793
          IMPLICIT NONE
 
794
          COMPLEX*32 P2,M12,M22
 
795
          COMPLEX*32 ZERO,TWOPII
 
796
          PARAMETER (ZERO=(0.0E0_16,0.0E0_16))
 
797
          PARAMETER (TWOPII=2.0E0_16
 
798
     $     *3.14169258478796109557151794433593750E0_16*(0.0E0_16
 
799
     $     ,1.0E0_16))
 
800
          REAL*16 M,M2,GA,GA2
 
801
          REAL*16 TINY
 
802
          PARAMETER (TINY=-1.0E-14_16)
 
803
          COMPLEX*32 LOGTERMS
 
804
          COMPLEX*32 MP_LOG_TRAJECTORY
 
805
          LOGICAL USE_CACHING
 
806
          PARAMETER (USE_CACHING=.TRUE.)
 
807
          TYPE(MP_B0F_NODE),POINTER::ITEM
 
808
          TYPE(MP_B0F_NODE),POINTER,SAVE::B0F_BT
 
809
          INTEGER INIT
 
810
          SAVE INIT
 
811
          DATA INIT /0/
 
812
          LOGICAL FIND
 
813
          IF(M12.EQ.ZERO)THEN
 
814
            M=REAL(P2,KIND=16)
 
815
            M2=REAL(M22,KIND=16)
 
816
            IF(M.LT.TINY.OR.M2.LT.TINY)THEN
 
817
              WRITE(*,*)'ERROR:MP_B0F is not well defined when M^2'
 
818
     $         //',M2^2<0'
 
819
              STOP
 
820
            ENDIF
 
821
            M=SQRT(ABS(M))
 
822
            M2=SQRT(ABS(M2))
 
823
            IF(M.EQ.0.0E0_16)THEN
 
824
              GA=0.0E0_16
 
825
            ELSE
 
826
              GA=-IMAGPART(P2)/M
 
827
            ENDIF
 
828
            IF(M2.EQ.0.0E0_16)THEN
 
829
              GA2=0.0E0_16
 
830
            ELSE
 
831
              GA2=-IMAGPART(M22)/M2
 
832
            ENDIF
 
833
            IF(P2.NE.M22.AND.P2.NE.ZERO.AND.M22.NE.ZERO)THEN
 
834
              MP_B0F=(M22-P2)/P2*LOG((M22-P2)/M22)
 
835
              IF(M.GT.M2.AND.GA*M2.GT.GA2*M)THEN
 
836
                MP_B0F=MP_B0F-TWOPII
 
837
              ENDIF
 
838
              RETURN
 
839
            ELSE
 
840
              WRITE(*,*)'ERROR:MP_B0F is not supported for a simple'
 
841
     $         //' form'
 
842
              STOP
 
843
            ENDIF
 
844
          ENDIF
 
845
          IF(USE_CACHING)THEN
 
846
            IF(INIT.EQ.0)THEN
 
847
              NULLIFY(B0F_BT)
 
848
              INIT=1
 
849
            ENDIF
 
850
            ALLOCATE(ITEM)
 
851
            ITEM%%P2=P2
 
852
            ITEM%%M12=M12
 
853
            ITEM%%M22=M22
 
854
            FIND=.FALSE.
 
855
            CALL MP_B0F_SEARCH(ITEM, B0F_BT, FIND)
 
856
            IF(FIND)THEN
 
857
              MP_B0F=ITEM%%VALUE
 
858
              DEALLOCATE(ITEM)
 
859
              RETURN
 
860
            ELSE
 
861
              LOGTERMS=MP_LOG_TRAJECTORY(100,P2,M12,M22)
 
862
              MP_B0F=-LOG(P2/M22)+LOGTERMS
 
863
              ITEM%%VALUE=MP_B0F
 
864
              RETURN
 
865
            ENDIF
 
866
          ELSE
 
867
            LOGTERMS=MP_LOG_TRAJECTORY(100,P2,M12,M22)
 
868
            MP_B0F=-LOG(P2/M22)+LOGTERMS
 
869
          ENDIF
 
870
          RETURN
 
871
          END
 
872
 
 
873
          COMPLEX*32 FUNCTION MP_SQRT_TRAJECTORY(N_SEG,P2,M12,M22)
 
874
          IMPLICIT NONE
 
875
          INTEGER N_SEG
 
876
          COMPLEX*32 P2,M12,M22
 
877
          COMPLEX*32 ZERO,ONE
 
878
          PARAMETER (ZERO=(0.0E0_16,0.0E0_16),ONE=(1.0E0_16,0.0E0_16))
 
879
          COMPLEX*32 GAMMA0,GAMMA1
 
880
          REAL*16 M,GA,DGA,GA_START
 
881
          REAL*16 GAI,INTERSECTION
 
882
          COMPLEX*32 ARGIM1,ARGI,P2I
 
883
          COMPLEX*32 GAMMA0I,GAMMA1I
 
884
          REAL*16 TINY
 
885
          PARAMETER (TINY=-1.0E-24_16)
 
886
          INTEGER I
 
887
          REAL*16 PREFACTOR
 
888
          IF(ABS(P2*M12*M22).EQ.0.0E0_16)THEN
 
889
            WRITE(*,*)'ERROR:mp_sqrt_trajectory works when p2*m12*m22'
 
890
     $       //'/=0'
 
891
            STOP
 
892
          ENDIF
 
893
          M=REAL(P2,KIND=16)
 
894
          M=SQRT(ABS(M))
 
895
          IF(M.EQ.0.0E0_16)THEN
 
896
            GA=0.0E0_16
 
897
          ELSE
 
898
            GA=-IMAGPART(P2)/M
 
899
          ENDIF
 
900
          GAMMA0=ONE+M12/P2-M22/P2
 
901
          GAMMA1=M12/P2-CMPLX(0.0E0_16,1.0E0_16)*ABS(TINY)/P2
 
902
          IF(ABS(GA).EQ.0.0E0_16)THEN
 
903
            MP_SQRT_TRAJECTORY=SQRT(GAMMA0**2-4.0E0_16*GAMMA1)
 
904
            RETURN
 
905
          ENDIF
 
906
          GA_START=-ABS(TINY*GA)
 
907
          DGA=(GA-GA_START)/N_SEG
 
908
          PREFACTOR=1.0E0_16
 
909
          GAI=GA_START
 
910
          P2I=CMPLX(M**2,-GAI*M)
 
911
          GAMMA0I=ONE+M12/P2I-M22/P2I
 
912
          GAMMA1I=M12/P2I-CMPLX(0.0E0_16,1.0E0_16)*ABS(TINY)/P2I
 
913
          ARGIM1=GAMMA0I**2-4.0E0_16*GAMMA1I
 
914
          DO I=1,N_SEG
 
915
            GAI=DGA*I+GA_START
 
916
            P2I=CMPLX(M**2,-GAI*M)
 
917
            GAMMA0I=ONE+M12/P2I-M22/P2I
 
918
            GAMMA1I=M12/P2I-CMPLX(0.0E0_16,1.0E0_16)*ABS(TINY)/P2I
 
919
            ARGI=GAMMA0I**2-4.0E0_16*GAMMA1I
 
920
            IF(IMAGPART(ARGI)*IMAGPART(ARGIM1).LT.0.0E0_16)THEN
 
921
              INTERSECTION=IMAGPART(ARGIM1)*(REAL(ARGI,KIND=16)
 
922
     $         -REAL(ARGIM1,KIND=16))
 
923
              INTERSECTION=INTERSECTION/(IMAGPART(ARGI)
 
924
     $         -IMAGPART(ARGIM1))
 
925
              INTERSECTION=INTERSECTION-REAL(ARGIM1,KIND=16)
 
926
              IF(INTERSECTION.GT.0.0E0_16)THEN
 
927
                PREFACTOR=-PREFACTOR
 
928
              ENDIF
 
929
            ENDIF
 
930
            ARGIM1=ARGI
 
931
          ENDDO
 
932
          MP_SQRT_TRAJECTORY=SQRT(GAMMA0**2-4.0E0_16*GAMMA1)*PREFACTOR
 
933
          RETURN
 
934
          END
 
935
 
 
936
          COMPLEX*32 FUNCTION MP_LOG_TRAJECTORY(N_SEG,P2,M12,M22)
 
937
          IMPLICIT NONE
 
938
          INTEGER N_SEG
 
939
          COMPLEX*32 P2,M12,M22
 
940
          COMPLEX*32 ZERO,ONE,HALF,TWOPII
 
941
          PARAMETER (ZERO=(0.0E0_16,0.0E0_16),ONE=(1.0E0_16,0.0E0_16))
 
942
          PARAMETER (HALF=(0.5E0_16,0.0E0_16))
 
943
          PARAMETER (TWOPII=2.0E0_16
 
944
     $     *3.14169258478796109557151794433593750E0_16*(0.0E0_16
 
945
     $     ,1.0E0_16))
 
946
          COMPLEX*32 GAMMA0,GAMMAP,GAMMAM,SQRTTERM
 
947
          REAL*16 M,GA,DGA,GA_START
 
948
          REAL*16 GAI,INTERSECTION
 
949
          COMPLEX*32 ARGIM1(4),ARGI(4),P2I,SQRTTERMI
 
950
          COMPLEX*32 GAMMA0I,GAMMAPI,GAMMAMI
 
951
          REAL*16 TINY
 
952
          PARAMETER (TINY=-1.0E-14_16)
 
953
          INTEGER I,J
 
954
          COMPLEX*32 ADDFACTOR(4)
 
955
          COMPLEX*32 MP_SQRT_TRAJECTORY
 
956
          IF(ABS(P2*M12*M22).EQ.0.0E0_16)THEN
 
957
            WRITE(*,*)'ERROR:mp_log_trajectory works when p2*m12*m22'
 
958
     $       //'/=0'
 
959
            STOP
 
960
          ENDIF
 
961
          M=REAL(P2,KIND=16)
 
962
          M=SQRT(ABS(M))
 
963
          IF(M.EQ.0.0E0_16)THEN
 
964
            GA=0.0E0_16
 
965
          ELSE
 
966
            GA=-IMAGPART(P2)/M
 
967
          ENDIF
 
968
          SQRTTERM=MP_SQRT_TRAJECTORY(N_SEG,P2,M12,M22)
 
969
          GAMMA0=ONE+M12/P2-M22/P2
 
970
          GAMMAP=HALF*(GAMMA0+SQRTTERM)
 
971
          GAMMAM=HALF*(GAMMA0-SQRTTERM)
 
972
          IF(ABS(GA).EQ.0.0E0_16)THEN
 
973
            MP_LOG_TRAJECTORY=-LOG(GAMMAP-ONE)-LOG(GAMMAM-ONE)+GAMMAP
 
974
     $       *LOG((GAMMAP-ONE)/GAMMAP)+GAMMAM*LOG((GAMMAM-ONE)/GAMMAM)
 
975
            RETURN
 
976
          ENDIF
 
977
          GA_START=-ABS(TINY*GA)
 
978
          DGA=(GA-GA_START)/N_SEG
 
979
          ADDFACTOR(1:4)=ZERO
 
980
          GAI=GA_START
 
981
          P2I=CMPLX(M**2,-GAI*M)
 
982
          SQRTTERMI=MP_SQRT_TRAJECTORY(N_SEG,P2I,M12,M22)
 
983
          GAMMA0I=ONE+M12/P2I-M22/P2I
 
984
          GAMMAPI=HALF*(GAMMA0I+SQRTTERMI)
 
985
          GAMMAMI=HALF*(GAMMA0I-SQRTTERMI)
 
986
          ARGIM1(1)=GAMMAPI-ONE
 
987
          ARGIM1(2)=GAMMAMI-ONE
 
988
          ARGIM1(3)=(GAMMAPI-ONE)/GAMMAPI
 
989
          ARGIM1(4)=(GAMMAMI-ONE)/GAMMAMI
 
990
          DO I=1,N_SEG
 
991
            GAI=DGA*I+GA_START
 
992
            P2I=CMPLX(M**2,-GAI*M)
 
993
            SQRTTERMI=MP_SQRT_TRAJECTORY(N_SEG,P2I,M12,M22)
 
994
            GAMMA0I=ONE+M12/P2I-M22/P2I
 
995
            GAMMAPI=HALF*(GAMMA0I+SQRTTERMI)
 
996
            GAMMAMI=HALF*(GAMMA0I-SQRTTERMI)
 
997
            ARGI(1)=GAMMAPI-ONE
 
998
            ARGI(2)=GAMMAMI-ONE
 
999
            ARGI(3)=(GAMMAPI-ONE)/GAMMAPI
 
1000
            ARGI(4)=(GAMMAMI-ONE)/GAMMAMI
 
1001
            DO J=1,4
 
1002
              IF(IMAGPART(ARGI(J))*IMAGPART(ARGIM1(J)).LT.0.0E0_16)THEN
 
1003
                INTERSECTION=IMAGPART(ARGIM1(J))*(REAL(ARGI(J),KIND=16)
 
1004
     $           -REAL(ARGIM1(J),KIND=16))
 
1005
                INTERSECTION=INTERSECTION/(IMAGPART(ARGI(J))
 
1006
     $           -IMAGPART(ARGIM1(J)))
 
1007
                INTERSECTION=INTERSECTION-REAL(ARGIM1(J),KIND=16)
 
1008
                IF(INTERSECTION.GT.0.0E0_16)THEN
 
1009
                  IF(IMAGPART(ARGIM1(J)).LT.0.0E0_16)THEN
 
1010
                    ADDFACTOR(J)=ADDFACTOR(J)-TWOPII
 
1011
                  ELSE
 
1012
                    ADDFACTOR(J)=ADDFACTOR(J)+TWOPII
 
1013
                  ENDIF
 
1014
                ENDIF
 
1015
              ENDIF
 
1016
              ARGIM1(J)=ARGI(J)
 
1017
            ENDDO
 
1018
          ENDDO
 
1019
          MP_LOG_TRAJECTORY=-(LOG(GAMMAP-ONE)+ADDFACTOR(1))
 
1020
     $     -(LOG(GAMMAM-ONE)+ADDFACTOR(2))
 
1021
          MP_LOG_TRAJECTORY=MP_LOG_TRAJECTORY+GAMMAP*(LOG((GAMMAP-ONE)
 
1022
     $     /GAMMAP)+ADDFACTOR(3))
 
1023
          MP_LOG_TRAJECTORY=MP_LOG_TRAJECTORY+GAMMAM*(LOG((GAMMAM-ONE)
 
1024
     $     /GAMMAM)+ADDFACTOR(4))
 
1025
          RETURN
 
1026
          END
 
1027
 
 
1028
          COMPLEX*32 FUNCTION MP_ARG(COMNUM)
 
1029
          IMPLICIT NONE
 
1030
          COMPLEX*32 COMNUM
 
1031
          COMPLEX*32 IMM
 
1032
          IMM = (0.0E0_16,1.0E0_16)
 
1033
          IF(COMNUM.EQ.(0.0E0_16,0.0E0_16)) THEN
 
1034
            MP_ARG=(0.0E0_16,0.0E0_16)
 
1035
          ELSE
 
1036
            MP_ARG=LOG(COMNUM/ABS(COMNUM))/IMM
 
1037
          ENDIF
 
1038
          END