← Back to branch summary

~alifson/chiralityflow/trunk

« back to all changes in this revision

Viewing changes to vendor/IREGI/src/oneloop/src/avh_olo_dilog.f90

  • 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:
vendor/IREGI/src/oneloop/src/avh_olo_dilog.f90
 
1
!!
 
2
!! Copyright (C) 2014 Andreas van Hameren. 
 
3
!!
 
4
!! This file is part of OneLOop-3.4.
 
5
!!
 
6
!! OneLOop-3.4 is free software: you can redistribute it and/or modify
 
7
!! it under the terms of the GNU General Public License as published by
 
8
!! the Free Software Foundation, either version 3 of the License, or
 
9
!! (at your option) any later version.
 
10
!!
 
11
!! OneLOop-3.4 is distributed in the hope that it will be useful,
 
12
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
!! GNU General Public License for more details.
 
15
!!
 
16
!! You should have received a copy of the GNU General Public License
 
17
!! along with OneLOop-3.4.  If not, see <http://www.gnu.org/licenses/>.
 
18
!!
 
19
 
 
20
 
 
21
module avh_olo_forIREGI_dilog
 
22
!***********************************************************************
 
23
!                     /1    ln(1-zz*t)
 
24
!   dilog(xx,iph) = - |  dt ---------- 
 
25
!                     /0        t
 
26
! with  zz = 1 - xx*exp(imag*pi*iph)  [pi, NOT 2*pi]
 
27
!
 
28
!   dilog(x1,i1,x2,i2) = ( dilog(x1,i1)-dilog(x2,i2) )/( x1-x2 )
 
29
!
 
30
! Arguments xx,x1,x2, may be all real or all complex,
 
31
! arguments iph,i1,i2 must be all integer.
 
32
!***********************************************************************
 
33
  use avh_olo_forIREGI_units
 
34
  use avh_olo_forIREGI_prec
 
35
  use avh_olo_forIREGI_print
 
36
  use avh_olo_forIREGI_auxfun
 
37
  use avh_olo_forIREGI_arrays
 
38
  implicit none
 
39
  private
 
40
  public :: update_dilog,dilog
 
41
 
 
42
  include 'avh_olo_real.h90'
 
43
         ,allocatable,save :: coeff(:)
 
44
  include 'avh_olo_real.h90'
 
45
         ,allocatable,save :: thrs(:,:)
 
46
  integer,allocatable,save :: ntrm(:,:)
 
47
  integer,parameter :: nStp=6
 
48
 
 
49
  include 'avh_olo_real.h90'
 
50
         ,allocatable :: bern(:),fact(:)
 
51
 
 
52
  interface dilog
 
53
    module procedure dilog_c,dilog_r,dilog2_c,dilog2_r
 
54
  end interface
 
55
 
 
56
contains
 
57
 
 
58
  subroutine update_dilog
 
59
!***********************************************************************
 
60
!***********************************************************************
 
61
  include 'avh_olo_real.h90'
 
62
    :: tt
 
63
  integer :: nn,ii,jj
 
64
  logical :: highestSoFar
 
65
!  real(kind(1d0)) :: xx(6) !DEBUG
 
66
!
 
67
  if (allocated(thrs)) then
 
68
    call shift2( thrs ,prcpar )
 
69
    call shift2( ntrm ,prcpar )
 
70
  else
 
71
    allocate(thrs(1:nStp,1:1))
 
72
    allocate(ntrm(1:nStp,1:1))
 
73
    if (prcpar.ne.1) then
 
74
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop update_dilog'
 
75
      stop
 
76
    endif
 
77
  endif
 
78
!
 
79
  highestSoFar = prcpar.eq.ubound(ntrm,2)
 
80
  if (highestSoFar) then
 
81
    if (allocated(coeff)) deallocate(coeff)
 
82
    allocate(coeff(0:-1)) ! allocate at size=0
 
83
  endif
 
84
!
 
85
  if (prcpar.gt.1) then ;nn=ntrm(nStp,prcpar-1)-1
 
86
                   else ;nn=2
 
87
  endif
 
88
!
 
89
  do
 
90
    nn = nn+1
 
91
    if (nn.gt.ubound(coeff,1)) call update_coeff( 2*nn )
 
92
    tt = 1
 
93
    tt = (EPSN/abs(coeff(nn)))**(tt/(2*nn))
 
94
! expansion parameter is smaller than 1.05
 
95
    if (100*tt.gt.105*RONE) exit
 
96
  enddo
 
97
!
 
98
  if (highestSoFar) call resize( coeff ,0,nn )
 
99
!
 
100
  ntrm(nStp,prcpar) = nn
 
101
  thrs(nStp,prcpar) = tt
 
102
  nn = max(1,nint(nn*1d0/nStp))
 
103
  do ii=nStp-1,1,-1
 
104
    ntrm(ii,prcpar) = ntrm(ii+1,prcpar)-nn
 
105
    if (ntrm(ii,prcpar).le.2) then
 
106
      do jj=1,ii
 
107
        ntrm(jj,prcpar) = max(2,ntrm(ii,prcpar))
 
108
        thrs(jj,prcpar) = 0 
 
109
      enddo
 
110
      exit
 
111
    endif
 
112
    jj = ntrm(ii,prcpar)
 
113
    tt = 1
 
114
    tt = (EPSN/abs(coeff(jj)))**(tt/(2*jj))
 
115
    thrs(ii,prcpar) = tt
 
116
  enddo
 
117
!
 
118
  if (allocated(bern)) deallocate(bern)
 
119
  if (allocated(fact)) deallocate(fact)
 
120
!
 
121
!  do ii=lbound(thrs,2),ubound(thrs,2) !DEBUG
 
122
!    do jj=1,nStp                      !DEBUG
 
123
!      xx(jj) = thrs(jj,ii)            !DEBUG
 
124
!    enddo                             !DEBUG
 
125
!    write(*,'(99e10.3)') xx(:)        !DEBUG
 
126
!    write(*,'(99i10)'  ) ntrm(:,ii)   !DEBUG
 
127
!  enddo                               !DEBUG
 
128
  end subroutine
 
129
 
 
130
 
 
131
  subroutine update_coeff( ncf )
 
132
!*******************************************************************
 
133
!   coeff(0)=-1/4
 
134
!   coeff(n)=bern(2*n)/(2*n+1)
 
135
!    bern(n)=bernoulli(n)/n!
 
136
!    fact(n)=n!
 
137
! DO NOT SKIP THE ODD bern IN THE RECURSIVE LOOP
 
138
! DO NOT PUT THE ODD bern TO ZERO
 
139
!*******************************************************************
 
140
  integer ,intent(in) :: ncf
 
141
  integer :: ii,jj,nbern,nold
 
142
!
 
143
  if (allocated(bern)) then ;nold=ubound(bern,1)
 
144
                       else ;nold=0
 
145
  endif
 
146
!
 
147
  nbern = 2*ncf
 
148
!
 
149
  call enlarge( bern  ,1,nbern   )
 
150
  call enlarge( fact  ,0,nbern+1 )
 
151
  call enlarge( coeff ,0,ncf     )
 
152
!
 
153
  fact(0) = 1
 
154
  do ii=nold+1,nbern+1
 
155
    fact(ii) = fact(ii-1)*ii
 
156
  enddo
 
157
!
 
158
  do ii=nold+1,nbern
 
159
    bern(ii) = -1/fact(ii+1)
 
160
    do jj=1,ii-1
 
161
      bern(ii) = bern(ii) - bern(jj)/fact(ii+1-jj)
 
162
    enddo
 
163
  enddo
 
164
!
 
165
  coeff(0) = 1
 
166
  coeff(0) =-coeff(0)/4
 
167
  do ii=nold+2,nbern,2
 
168
    coeff(ii/2) = bern(ii)/(ii+1)
 
169
  enddo
 
170
!
 
171
  end subroutine
 
172
 
 
173
 
 
174
  function dilog_c(xx,iph) result(rslt)
 
175
!*******************************************************************
 
176
!*******************************************************************
 
177
  include 'avh_olo_complex.h90'
 
178
          ,intent(in) :: xx
 
179
  integer ,intent(in) :: iph
 
180
  include 'avh_olo_complex.h90'
 
181
    :: rslt ,yy,lyy,loy,zz,z2
 
182
  include 'avh_olo_real.h90'
 
183
    :: rex,imx,az
 
184
  integer :: ii,jj,ntwo,odd,nn
 
185
  logical :: r_gt_1 , y_lt_h
 
186
!
 
187
  rex = areal(xx)
 
188
  imx = aimag(xx)
 
189
!
 
190
  if (abs(imx).le.EPSN*abs(rex)) then
 
191
    if (rex.ge.RZRO) then
 
192
      rslt = dilog_r( rex, iph )
 
193
    else
 
194
      rslt = dilog_r(-rex, iph+sgnRe(imx) )
 
195
    endif
 
196
    return
 
197
  endif
 
198
!
 
199
  if (rex.gt.RZRO) then ;yy= xx ;jj=iph
 
200
                   else ;yy=-xx ;jj=iph+sgnRe(imx)
 
201
  endif
 
202
!
 
203
  odd = mod(jj,2)
 
204
  ntwo = jj-odd
 
205
! 
 
206
  r_gt_1 = (rex*rex+imx*imx.gt.RONE)
 
207
  lyy = log(yy)
 
208
  if (odd.ne.0) yy = -yy
 
209
!
 
210
  if (r_gt_1) then
 
211
    yy   = 1/yy
 
212
    lyy  =-lyy
 
213
    ntwo =-ntwo
 
214
    odd  =-odd
 
215
  endif
 
216
  loy = log(1-yy)
 
217
!
 
218
  y_lt_h = (2*areal(yy).lt.RONE)
 
219
  if (y_lt_h) then ;zz=-loy
 
220
              else ;zz=-lyy
 
221
  endif
 
222
!
 
223
  az = abs(zz)
 
224
! if (az.gt.thrs(6,prcpar)) ERROR az to big 
 
225
  if     (az.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
226
  elseif (az.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
227
  elseif (az.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
228
  elseif (az.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
229
  elseif (az.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
230
                                else ;nn=ntrm(1,prcpar)
 
231
  endif
 
232
  z2 = zz*zz
 
233
  rslt = coeff(nn)
 
234
  do ii=nn,2,-1
 
235
    rslt = coeff(ii-1) + z2*rslt
 
236
  enddo
 
237
  rslt = zz*( 1 + zz*( coeff(0) + zz*rslt ) )
 
238
!
 
239
  if (y_lt_h) then
 
240
    rslt = 4*PISQo24 - rslt - loy*(lyy+IPI*(ntwo+odd))
 
241
  else
 
242
    rslt = rslt - loy*IPI*ntwo
 
243
  endif
 
244
!
 
245
  if (r_gt_1) rslt = -rslt - (lyy+IPI*(ntwo+odd))**2/2
 
246
  end function
 
247
 
 
248
 
 
249
 
 
250
  function dilog_r(xx,iph) result(rslt)
 
251
!*******************************************************************
 
252
!*******************************************************************
 
253
  include 'avh_olo_real.h90'
 
254
          ,intent(in) :: xx
 
255
  integer ,intent(in) :: iph
 
256
  include 'avh_olo_complex.h90'
 
257
    :: rslt
 
258
  include 'avh_olo_real.h90'
 
259
    :: yy,lyy,loy,zz,z2,liox,az
 
260
  integer :: jj,ii,ntwo,odd,nn
 
261
  logical :: r_gt_1 , y_lt_h
 
262
!
 
263
  if (xx.eq.RZRO) then
 
264
    rslt = 4*PISQo24
 
265
    return
 
266
  elseif (xx.gt.RZRO) then ;yy= xx ;jj=iph
 
267
                      else ;yy=-xx ;jj=iph+1 ! log(-1)=i*pi
 
268
  endif
 
269
!
 
270
  odd = mod(jj,2)
 
271
  ntwo = jj-odd
 
272
! 
 
273
  if (yy.eq.RONE.and.odd.eq.0) then
 
274
    if (ntwo.ne.0) then
 
275
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog_r: ' &
 
276
        ,'|x|,iph = ',trim(myprint(yy)),',',jj,', returning 0'
 
277
    endif
 
278
    rslt = 0
 
279
    return
 
280
  endif
 
281
!
 
282
  r_gt_1 = (yy.gt.RONE)
 
283
  lyy = log(yy)
 
284
  if (odd.ne.0) yy = -yy
 
285
!
 
286
  if (r_gt_1) then
 
287
    yy   = 1/yy
 
288
    lyy  =-lyy
 
289
    ntwo =-ntwo
 
290
    odd  =-odd
 
291
  endif
 
292
  loy = log(1-yy) ! log(1-yy) is always real
 
293
!
 
294
  y_lt_h = (2*yy.lt.RONE)
 
295
  if (y_lt_h) then
 
296
    zz = -loy ! log(1-yy) is real
 
297
  else
 
298
    zz = -lyy ! yy>0.5 => log(yy) is real
 
299
  endif
 
300
!
 
301
  az = abs(zz)
 
302
! if (az.gt.thrs(6,prcpar)) ERROR az to big 
 
303
  if     (az.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
304
  elseif (az.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
305
  elseif (az.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
306
  elseif (az.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
307
  elseif (az.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
308
                                else ;nn=ntrm(1,prcpar)
 
309
  endif
 
310
  z2 = zz*zz
 
311
  liox = coeff(nn)
 
312
  do ii=nn,2,-1
 
313
    liox = coeff(ii-1) + z2*liox
 
314
  enddo
 
315
  liox = zz*( 1 + zz*( coeff(0) + zz*liox ) )
 
316
!
 
317
  rslt = acmplx(liox)
 
318
!
 
319
  if (y_lt_h) then
 
320
    rslt = 4*PISQo24 - rslt - acmplx(loy*lyy,loy*ONEPI*(ntwo+odd))
 
321
  else
 
322
    rslt = rslt + acmplx( 0 ,-loy*ONEPI*ntwo )
 
323
  endif
 
324
!
 
325
  if (r_gt_1) rslt = -rslt - acmplx(lyy,ONEPI*(ntwo+odd))**2/2
 
326
  end function
 
327
 
 
328
 
 
329
  function dilog2_c( x1,i1 ,x2,i2 ) result(rslt)
 
330
!*******************************************************************
 
331
!*******************************************************************
 
332
  use avh_olo_forIREGI_olog
 
333
  include 'avh_olo_complex.h90'
 
334
          ,intent(in) :: x1,x2
 
335
  integer ,intent(in) :: i1,i2
 
336
  include 'avh_olo_complex.h90'
 
337
    :: rslt ,y1,y2 ,ff,gg,logr1,logr2,logo1,logo2,r1,r2,rr
 
338
  include 'avh_olo_real.h90'
 
339
    :: eps ,re1,im1,re2,im2,a1,a2,aa,ao1,ao2
 
340
  integer :: j1,j2,ii,nn,oo
 
341
  integer,parameter :: pp(-1:1,-1:1)=&
 
342
                      reshape((/-2,-2,2 ,-2,0,2 ,-2,2,2/),(/3,3/))
 
343
!
 
344
  re1=areal(x1) ;re2=areal(x2)
 
345
  im1=aimag(x1) ;im2=aimag(x2)
 
346
!
 
347
  if (abs(im1).le.EPSN*abs(re1).and.abs(im2).le.EPSN*abs(re2)) then
 
348
    if (re1.ge.RZRO) then
 
349
      if (re2.ge.RZRO) then
 
350
        rslt = dilog2_r( re1,i1 , re2,i2 )
 
351
      else
 
352
        rslt = dilog2_r( re1,i1 ,-re2,i2+sgnRe(im2) )
 
353
      endif
 
354
    elseif (re2.ge.RZRO) then
 
355
      rslt = dilog2_r(-re1,i1+sgnRe(im1) , re2,i2 )
 
356
    else
 
357
      rslt = dilog2_r(-re1,i1+sgnRe(im1) ,-re2,i2+sgnRe(im2) )
 
358
    endif
 
359
    return
 
360
  endif
 
361
!
 
362
  if (re1.ge.RZRO) then ;r1= x1 ;j1=i1
 
363
                   else ;r1=-x1 ;j1=i1+sgnRe(im1,1)
 
364
  endif
 
365
  if (re2.ge.RZRO) then ;r2= x2 ;j2=i2
 
366
                   else ;r2=-x2 ;j2=i2+sgnRe(im2,1)
 
367
  endif
 
368
!
 
369
  a1=abs(r1) ;a2=abs(r2)
 
370
  if (a1.gt.a2) then
 
371
    aa=a1;a1=a2;a2=aa
 
372
    rr=r1;r1=r2;r2=rr
 
373
    ii=j1;j1=j2;j2=ii
 
374
  endif
 
375
!
 
376
  oo=mod(j1,2) ;nn=j1-oo ;y1=r1 ;if (oo.ne.0) y1=-y1
 
377
  oo=mod(j2,2) ;nn=j2-oo ;y2=r2 ;if (oo.ne.0) y2=-y2
 
378
!
 
379
  eps = 10*EPSN
 
380
!
 
381
  if (j1.ne.j2) then
 
382
    if (r1.eq.r2) then
 
383
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
384
        ,'j1,j2,r1-r2',j1,j2,',',trim(myprint(r1-r2)),', returning 0'
 
385
      rslt = 0
 
386
!      write(*,*) 'dilog2_c j1=/=j2,r1=r2' !DEBUG
 
387
      return
 
388
    else
 
389
      rslt = ( dilog_c(r1,j1)-dilog_c(r2,j2) )/(y1-y2)
 
390
!      write(*,*) 'dilog2_c j1=/=j2' !DEBUG
 
391
      return
 
392
    endif
 
393
  endif
 
394
!
 
395
  if (a1.lt.eps) then
 
396
    if (a2.lt.eps) then
 
397
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
398
        ,'r1,r2 =',trim(myprint(r1)),',',trim(myprint(r2)),', returning 0'
 
399
      rslt = 0
 
400
!      write(*,*) 'dilog2_c r1<eps,r2<eps' !DEBUG
 
401
      return
 
402
    else
 
403
      rslt = (dilog_c(r2,j2)-4*PISQo24)/y2
 
404
!      write(*,*) 'dilog2_c r1<eps' !DEBUG
 
405
      return
 
406
    endif
 
407
  endif
 
408
!
 
409
  logr1=log(r1) ;logr2=log(r2)
 
410
!
 
411
  ao1=abs(1-y1) ;ao2=abs(1-y2)
 
412
  if (10*ao1.lt.RONE.or.10*ao2.lt.RONE) then
 
413
    aa = abs(r1/r2-1)
 
414
    if (10*aa.gt.RONE) then
 
415
      rslt = (dilog_c(r1,j1)-dilog_c(r2,j2))/(y1-y2)
 
416
!      write(*,*) 'dilog2_c ||1-y1|/|1-y2|-1|>0.1' !DEBUG
 
417
      return
 
418
    elseif (oo.eq.0.and.ao1.lt.eps) then
 
419
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
420
        ,'r1,oo,nn =',trim(myprint(r1)),',',oo,nn,', putting nn=0'
 
421
      if (ao2.lt.eps) then
 
422
        rslt = -1
 
423
!        write(*,*) 'dilog2_c |1-y1|' !DEBUG
 
424
        return
 
425
      else
 
426
        y1=1-eps ;nn=0 ;logr1=0 ;r1=1-eps
 
427
      endif
 
428
    elseif (oo.eq.0.and.ao2.lt.eps) then
 
429
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
430
        ,'r2,oo,nn =',trim(myprint(r2)),',',oo,nn,', putting nn=0'
 
431
      y2=1-eps ;nn=0 ;logr2=0 ;r2=1-eps
 
432
    endif
 
433
  else
 
434
    aa = abs((logr1+oo*IPI)/(logr2+oo*IPI)-1)
 
435
    if (10*aa.gt.RONE) then
 
436
      rslt = (dilog_c(r1,j1)-dilog_c(r2,j2))/(y1-y2)
 
437
!      write(*,*) 'dilog2_c |logr1/logr2-1|>0.1',logr1,logr2 !DEBUG
 
438
      return
 
439
    elseif (aa.lt.eps) then
 
440
      ii = 0
 
441
      if (a1.gt.RONE) ii = ii + (nn+pp(oo,sgnIm(y2)))
 
442
      if (a2.gt.RONE) ii = ii - (nn+pp(oo,sgnIm(y2)))
 
443
      ii = nn*ii
 
444
      if (ii.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
445
        ,'r1,r2,nn =',trim(myprint(r1)),',',trim(myprint(r2)),',',nn &
 
446
        ,', putting nn=0'
 
447
      rslt = -olog2(y2,0)
 
448
!      write(*,*) 'dilog2_c |logr1/lorg2|<eps' !DEBUG
 
449
      return
 
450
    endif
 
451
  endif
 
452
!
 
453
  if (a1.gt.RONE) then
 
454
    y1=1/y1 ;logr1=-logr1
 
455
    y2=1/y2 ;logr2=-logr2
 
456
    nn=-nn ;oo=-oo
 
457
  endif
 
458
!
 
459
  ff=y1/y2         ;ff=-olog2(ff,0)/y2
 
460
  gg=(1-y1)/(1-y2) ;gg=-olog2(gg,0)/(1-y2)
 
461
!
 
462
  if (2*areal(y1).ge.RONE) then
 
463
!    write(*,*) 'dilog2_c re>1/2' !DEBUG
 
464
    rslt = ff*sumterms_c(-logr1,-logr2) - nn*IPI*gg
 
465
  else
 
466
!    write(*,*) 'dilog2_c re<1/2' !DEBUG
 
467
    logo1 = log(1-y1)
 
468
    logo2 = log(1-y2)
 
469
    rslt = gg*( sumterms_c(-logo1,-logo2) - (nn+oo)*IPI - logr2 ) + ff*logo1
 
470
  endif
 
471
!
 
472
  if (a1.gt.RONE) then !implies also r2>1
 
473
!    write(*,*) 'dilog2_c r1>1,r2>1' !DEBUG
 
474
    rslt = y1*y2*( rslt - ff*((logr1+logr2)/2 + (nn+oo)*IPI) )
 
475
  elseif (a2.gt.RONE.and.nn.ne.0) then
 
476
!    write(*,*) 'dilog2_c r1<1,r2>1',oo,sgnIm(y2)!DEBUG
 
477
    rslt = rslt - 12*nn*( nn + pp(oo,sgnIm(y2)) )*PISQo24/(y1-y2)
 
478
  endif
 
479
!
 
480
  end function
 
481
 
 
482
 
 
483
  function dilog2_r( x1,i1 ,x2,i2 ) result(rslt)
 
484
!*******************************************************************
 
485
!*******************************************************************
 
486
  use avh_olo_forIREGI_olog
 
487
  include 'avh_olo_real.h90'
 
488
          ,intent(in) :: x1,x2
 
489
  integer ,intent(in) :: i1,i2
 
490
  include 'avh_olo_complex.h90'
 
491
    :: rslt
 
492
  include 'avh_olo_real.h90'
 
493
    :: y1,y2 ,ff,gg,logr1,logr2,logo1,logo2
 
494
  include 'avh_olo_real.h90'
 
495
    :: eps,r1,r2,rr,ro1,ro2
 
496
  integer :: j1,j2,ii,nn,oo
 
497
!
 
498
  if (x1.ge.RZRO) then ;r1= x1 ;j1=i1
 
499
                  else ;r1=-x1 ;j1=i1+1 ! log(-1)=i*pi
 
500
  endif
 
501
  if (x2.ge.RZRO) then ;r2= x2 ;j2=i2
 
502
                  else ;r2=-x2 ;j2=i2+1 ! log(-1)=i*pi
 
503
  endif
 
504
!
 
505
  if (r1.gt.r2) then
 
506
    rr=r1;r1=r2;r2=rr
 
507
    ii=j1;j1=j2;j2=ii
 
508
  endif
 
509
!
 
510
  oo=mod(j1,2) ;nn=j1-oo ;y1=r1 ;if (oo.ne.0) y1=-y1
 
511
  oo=mod(j2,2) ;nn=j2-oo ;y2=r2 ;if (oo.ne.0) y2=-y2
 
512
!
 
513
  eps = 10*EPSN
 
514
!
 
515
  if (j1.ne.j2) then
 
516
    if (r1.eq.r2) then
 
517
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
518
        ,'j1,j2,r1-r2',j1,j2,',',trim(myprint(r1-r2)),', returning 0'
 
519
      rslt = 0
 
520
!      write(*,*) 'dilog2_r j1=/=j2,r1=r2' !DEBUG
 
521
      return
 
522
    else
 
523
      rslt = ( dilog_r(r1,j1)-dilog_r(r2,j2) )/(y1-y2)
 
524
!      write(*,*) 'dilog2_r j1=/=j2' !DEBUG
 
525
      return
 
526
    endif
 
527
  endif
 
528
!
 
529
  if (r1.lt.eps) then
 
530
    if (r2.lt.eps) then
 
531
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
532
        ,'r1,r2 =',trim(myprint(r1)),',',trim(myprint(r2)),', returning 0'
 
533
      rslt = 0
 
534
!      write(*,*) 'dilog2_r r1<eps,r2<eps' !DEBUG
 
535
      return
 
536
    else
 
537
      rslt = (dilog_r(r2,j2)-4*PISQo24)/y2
 
538
!      write(*,*) 'dilog2_r r1<eps' !DEBUG
 
539
      return
 
540
    endif
 
541
  endif
 
542
!
 
543
  logr1=log(r1) ;logr2=log(r2)
 
544
!
 
545
  ro1=abs(1-y1) ;ro2=abs(1-y2)
 
546
  if (10*ro1.lt.RONE.or.10*ro2.lt.RONE) then
 
547
    rr = abs(r1/r2-1)
 
548
    if (10*rr.gt.RONE) then
 
549
      rslt = (dilog_r(r1,j1)-dilog_r(r2,j2))/(y1-y2)
 
550
!      write(*,*) 'dilog2_r ||1-y1|/|1-y2|-1|>0.1' !DEBUG
 
551
      return
 
552
    elseif (oo.eq.0.and.ro1.lt.eps) then
 
553
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
554
        ,'r1,oo,nn =',trim(myprint(r1)),',',oo,nn,', putting nn=0'
 
555
      if (ro2.lt.eps) then
 
556
        rslt = -1
 
557
!        write(*,*) 'dilog2_r |1-y1|' !DEBUG
 
558
        return
 
559
      else
 
560
        y1=1-eps ;nn=0 ;logr1=0 ;r1=1-eps
 
561
      endif
 
562
    elseif (oo.eq.0.and.ro2.lt.eps) then
 
563
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
564
        ,'r2,oo,nn =',trim(myprint(r2)),',',oo,nn,', putting nn=0'
 
565
      y2=1-eps ;nn=0 ;logr2=0 ;r2=1-eps
 
566
    endif
 
567
  else
 
568
    rr = abs((logr1+oo*IPI)/(logr2+oo*IPI)-1)
 
569
    if (10*rr.gt.RONE) then
 
570
      rslt = (dilog_r(r1,j1)-dilog_r(r2,j2))/(y1-y2)
 
571
!      write(*,*) 'dilog2_r |logr1/logr2-1|>0.1',logr1,logr2 !DEBUG
 
572
      return
 
573
    elseif (rr.lt.eps) then
 
574
      ii = 0
 
575
      if (r1.gt.RONE) ii = ii + (nn+2*oo)
 
576
      if (r2.gt.RONE) ii = ii - (nn+2*oo)
 
577
      ii = nn*ii
 
578
      if (ii.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
579
        ,'r1,r2,nn =',trim(myprint(r1)),',',trim(myprint(r2)),',',nn &
 
580
        ,', putting nn=0'
 
581
      rslt = -olog2(y2,0)
 
582
!      write(*,*) 'dilog2_r |logr1/lorg2|<eps' !DEBUG
 
583
      return
 
584
    endif
 
585
  endif
 
586
!
 
587
  if (r1.gt.RONE) then
 
588
    y1=1/y1 ;logr1=-logr1
 
589
    y2=1/y2 ;logr2=-logr2
 
590
    nn=-nn ;oo=-oo
 
591
  endif
 
592
!
 
593
  ff=y1/y2         ;ff=-olog2(ff,0)/y2
 
594
  gg=(1-y1)/(1-y2) ;gg=-olog2(gg,0)/(1-y2)
 
595
!
 
596
  if (2*y1.ge.RONE) then
 
597
!    write(*,*) 'dilog2_r re>1/2' !DEBUG
 
598
    rslt = ff*sumterms_r(-logr1,-logr2) - nn*IPI*gg
 
599
  else
 
600
!    write(*,*) 'dilog2_r re<1/2' !DEBUG
 
601
    logo1 = log(1-y1)
 
602
    logo2 = log(1-y2)
 
603
    rslt = gg*( sumterms_r(-logo1,-logo2) - (nn+oo)*IPI - logr2 ) + ff*logo1
 
604
  endif
 
605
!
 
606
  if (r1.gt.RONE) then !implies also r2>1
 
607
!    write(*,*) 'dilog2_r r1>1,r2>1' !DEBUG
 
608
    rslt = y1*y2*( rslt - ff*((logr1+logr2)/2 + (nn+oo)*IPI) )
 
609
  elseif (r2.gt.RONE.and.nn.ne.0) then
 
610
!    write(*,*) 'dilog2_r r1<1,r2>1' !DEBUG
 
611
    rslt = rslt - 12*nn*PISQo24*(nn+2*oo)/(y1-y2)
 
612
  endif
 
613
!
 
614
  end function
 
615
 
 
616
 
 
617
  function sumterms_c( z1,z2 ) result(rslt)
 
618
!***********************************************************************
 
619
! ( f(z1)-f(z2) )/( z1-z2 ), where
 
620
! f(z)= z + c0*z^2 + c1*z^3 + c2*z^5 + c3*z^7 + ...
 
621
!***********************************************************************
 
622
  include 'avh_olo_complex.h90'
 
623
    ,intent(in) :: z1,z2
 
624
  include 'avh_olo_complex.h90'
 
625
    :: rslt,yy,zz
 
626
  include 'avh_olo_real.h90'
 
627
    :: az
 
628
  integer :: ii,nn
 
629
  az = max(abs(z1),abs(z2))
 
630
  if     (az.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
631
  elseif (az.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
632
  elseif (az.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
633
  elseif (az.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
634
  elseif (az.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
635
                                else ;nn=ntrm(1,prcpar)
 
636
  endif
 
637
! calculates all z(i)=(z1^i-z2^i)/(z1-z2) numerically stable
 
638
!  zz(1) = 1
 
639
!  yy    = 1
 
640
!  do ii=2,2*nn+1
 
641
!    yy = z2*yy
 
642
!    zz(ii) = z1*zz(ii-1) + yy
 
643
!  enddo
 
644
  zz = 1
 
645
  yy = 1
 
646
  rslt = zz
 
647
  yy = z2*yy
 
648
  zz = z1*zz+yy
 
649
  rslt = rslt + coeff(0)*zz
 
650
  do ii=1,nn
 
651
    yy = z2*yy
 
652
    zz = z1*zz+yy
 
653
    rslt = rslt + coeff(ii)*zz
 
654
    yy = z2*yy
 
655
    zz = z1*zz+yy
 
656
  enddo
 
657
  end function  
 
658
 
 
659
 
 
660
  function sumterms_r( z1,z2 ) result(rslt)
 
661
!***********************************************************************
 
662
! ( f(z1)-f(z2) )/( z1-z2 ), where
 
663
! f(z)= z + c0*z^2 + c1*z^3 + c2*z^5 + c3*z^7 + ...
 
664
!***********************************************************************
 
665
  include 'avh_olo_real.h90'
 
666
    ,intent(in) :: z1,z2
 
667
  include 'avh_olo_real.h90'
 
668
    :: rslt,yy,zz
 
669
  include 'avh_olo_real.h90'
 
670
    :: az
 
671
  integer :: ii,nn
 
672
  az = max(abs(z1),abs(z2))
 
673
  if     (az.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
674
  elseif (az.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
675
  elseif (az.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
676
  elseif (az.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
677
  elseif (az.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
678
                                else ;nn=ntrm(1,prcpar)
 
679
  endif
 
680
  zz = 1
 
681
  yy = 1
 
682
  rslt = zz
 
683
  yy = z2*yy
 
684
  zz = z1*zz+yy
 
685
  rslt = rslt + coeff(0)*zz
 
686
  do ii=1,nn
 
687
    yy = z2*yy
 
688
    zz = z1*zz+yy
 
689
    rslt = rslt + coeff(ii)*zz
 
690
    yy = z2*yy
 
691
    zz = z1*zz+yy
 
692
  enddo
 
693
  end function  
 
694
 
 
695
end module