← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

Viewing changes to vendor/CutTools/src/avh/avh_olo.f90

merged with 2.3 rev286

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/CutTools/src/avh/avh_olo.f90
1
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
 
2
! Copyright (C) 2015 Andreas van Hameren. 
 
3
!
 
4
! This file is part of OneLOop-3.6.
 
5
!
 
6
! OneLOop-3.6 is free software: you can redistribute it and/or modify
7
7
! it under the terms of the GNU General Public License as published by
8
8
! the Free Software Foundation, either version 3 of the License, or
9
9
! (at your option) any later version.
10
10
!
11
 
! OneLOop-3.4 is distributed in the hope that it will be useful,
 
11
! OneLOop-3.6 is distributed in the hope that it will be useful,
12
12
! but WITHOUT ANY WARRANTY; without even the implied warranty of
13
13
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14
14
! GNU General Public License for more details.
15
15
!
16
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/>.
 
17
! along with OneLOop-3.6.  If not, see <http://www.gnu.org/licenses/>.
18
18
!
19
19
 
20
20
 
28
28
  if (done) return ;done=.true.
29
29
  write(*,'(a72)') '########################################################################'
30
30
  write(*,'(a72)') '#                                                                      #'
31
 
  write(*,'(a72)') '#                      You are using OneLOop-3.4                       #'
 
31
  write(*,'(a72)') '#                      You are using OneLOop-3.6                       #'
32
32
  write(*,'(a72)') '#                                                                      #'
33
33
  write(*,'(a72)') '# for the evaluation of 1-loop scalar 1-, 2-, 3- and 4-point functions #'
34
34
  write(*,'(a72)') '#                                                                      #'
35
35
  write(*,'(a72)') '# author: Andreas van Hameren <hamerenREMOVETHIS@ifj.edu.pl>           #'
36
 
  write(*,'(a72)') '#   date: 02-01-2014                                                   #'
 
36
  write(*,'(a72)') '#   date: 18-02-2015                                                   #'
37
37
  write(*,'(a72)') '#                                                                      #'
38
38
  write(*,'(a72)') '# Please cite                                                          #'
39
39
  write(*,'(a72)') '#    A. van Hameren,                                                   #'
50
50
module avh_olo_units
51
51
  implicit none
52
52
  integer :: eunit=6
53
 
  integer :: dunit=-1
54
53
  integer :: wunit=6
55
54
  integer :: munit=6
56
55
  integer :: punit=0 ! print all
781
780
      gg=xd1*pq1 ;hh=yd1*uv1
782
781
      rx2 = gg+hh
783
782
      if (abs(rx2).lt.neglig(prcpar)*max(abs(gg),abs(hh))) rx2 = 0
784
 
    else
 
783
    elseif (abs(pq2).gt.abs(pq1)) then
785
784
      rx2 = pq2
786
785
      gg=xd2*pq2 ;hh=yd2*uv2
787
786
      rx1 = gg+hh
788
787
      if (abs(rx1).lt.neglig(prcpar)*max(abs(gg),abs(hh))) rx1 = 0
 
788
    else
 
789
      rx1 = pq1
 
790
      rx2 = pq2
789
791
    endif
790
792
    if (abs(uv1).gt.abs(uv2)) then
791
793
      ix1 = uv1
792
794
      gg=yd1*pq1 ;hh=xd1*uv1
793
795
      ix2 = gg-hh
794
796
      if (abs(ix2).lt.neglig(prcpar)*max(abs(gg),abs(hh))) ix2 = 0
795
 
    else
 
797
    elseif (abs(uv2).gt.abs(uv1)) then
796
798
      ix2 = uv2
797
799
      gg=yd2*pq2 ;hh=xd2*uv2
798
800
      ix1 = gg-hh
799
801
      if (abs(ix1).lt.neglig(prcpar)*max(abs(gg),abs(hh))) ix1 = 0
 
802
    else
 
803
      ix1 = uv1
 
804
      ix2 = uv2
800
805
    endif
801
806
    x1 = acmplx(rx1,ix1)
802
807
    x2 = acmplx(rx2,ix2)
1138
1143
!***********************************************************************
1139
1144
! Provides the functions
1140
1145
!   olog(x,n) = log(x) + n*pi*imag  
1141
 
!   olog2(x,n) = olog(x,n)/(x-1)
 
1146
!   olog1(x,n) = olog(x,n)/(x-1)
 
1147
!   olog2(x,n) = ( olog1(x,n) - 1 )/(x-1)
 
1148
!   olog3(x,n) = ( olog2(x,n) + 1/2 )/(x-1)
1142
1149
! In the vicinity of x=1,n=0, the logarithm of complex argument is
1143
1150
! evaluated with a series expansion.
1144
1151
!***********************************************************************
1148
1155
  use avh_olo_dp_auxfun
1149
1156
  implicit none
1150
1157
  private
1151
 
  public :: update_olog,olog,olog2
 
1158
  public :: update_olog,olog,olog1,olog2,olog3
1152
1159
 
1153
1160
  real(kindr2) &  
1154
1161
         ,allocatable,save :: thrs(:,:)
1158
1165
  interface olog
1159
1166
    module procedure log_c,log_r
1160
1167
  end interface
 
1168
  interface olog1
 
1169
    module procedure log1_c,log1_r
 
1170
  end interface
1161
1171
  interface olog2
1162
1172
    module procedure log2_c,log2_r
1163
1173
  end interface
 
1174
  interface olog3
 
1175
    module procedure log3_c,log3_r
 
1176
  end interface
1164
1177
 
1165
1178
contains
1166
1179
 
1271
1284
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
1272
1285
                                else ;nn=ntrm(1,prcpar)
1273
1286
  endif
 
1287
! convergence acceleration using  z=(y-1)/(y+1)
 
1288
! rslt = 2 * ( z + z^3/3 + z^5/5 + z^7/7 + ... )  
1274
1289
  zz = zz/(yy+1)
1275
1290
  z2 = zz*zz
1276
1291
  aa = 2
1308
1323
  end function
1309
1324
 
1310
1325
 
1311
 
  function log2_c(xx,iph) result(rslt)
 
1326
  function log1_c(xx,iph) result(rslt)
1312
1327
!***********************************************************************
1313
1328
!***********************************************************************
1314
1329
  complex(kindr2) &   
1325
1340
!
1326
1341
  if (abs(imx).le.EPSN*abs(rex)) then
1327
1342
    if (rex.ge.RZRO) then
1328
 
      rslt = log2_r( rex, iph )
 
1343
      rslt = log1_r( rex, iph )
1329
1344
    else
1330
 
      rslt = log2_r(-rex, iph+sgnRe(imx) )
 
1345
      rslt = log1_r(-rex, iph+sgnRe(imx) )
1331
1346
    endif
1332
1347
    return
1333
1348
  endif
1353
1368
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
1354
1369
                                else ;nn=ntrm(1,prcpar)
1355
1370
  endif
 
1371
! convergence acceleration using  z=(y-1)/(y+1)
 
1372
! rslt = 2/(y+1) * ( 1 + z^2/3 + z^4/5 + z^6/7 + ... )  
1356
1373
  zz = zz/(yy+1)
1357
1374
  z2 = zz*zz
1358
1375
  aa = 2
1365
1382
  end function
1366
1383
 
1367
1384
 
1368
 
  function log2_r(xx,iph) result(rslt)
 
1385
  function log1_r(xx,iph) result(rslt)
1369
1386
!***********************************************************************
1370
1387
!***********************************************************************
1371
1388
  real(kindr2) &  
1381
1398
!  integer :: nn,ii
1382
1399
!
1383
1400
  if (xx.eq.RZRO) then
1384
 
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_r: ' &
 
1401
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log1_r: ' &
1385
1402
       ,'xx =',trim(myprint(xx)),', returning 0'
1386
1403
    rslt = 0
1387
1404
    return
1393
1410
!
1394
1411
  if (abs(yy-1).le.10*EPSN) then
1395
1412
    if (jj.ne.0) then
1396
 
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_r: ' &
 
1413
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log1_r: ' &
1397
1414
        ,'rr,jj =',trim(myprint(rr)),jj,', putting jj to 0'
1398
1415
    endif
1399
1416
    rslt = 1 - (yy-1)/2
1403
1420
  rslt = ( log(rr) + IPI*jj )/(yy-1)
1404
1421
  end function
1405
1422
 
 
1423
 
 
1424
  function log2_r(xx,iph) result(rslt)
 
1425
!***********************************************************************
 
1426
!***********************************************************************
 
1427
  real(kindr2) &  
 
1428
          ,intent(in) :: xx
 
1429
  integer ,intent(in) :: iph
 
1430
  complex(kindr2) &   
 
1431
    :: rslt
 
1432
  rslt = log2_c(xx*CONE,iph)
 
1433
  end function
 
1434
 
 
1435
 
 
1436
  function log2_c(xx,iph) result(rslt)
 
1437
!***********************************************************************
 
1438
!***********************************************************************
 
1439
  complex(kindr2) &   
 
1440
          ,intent(in) :: xx
 
1441
  integer ,intent(in) :: iph
 
1442
  complex(kindr2) &   
 
1443
    :: rslt ,yy,zz,z2
 
1444
  real(kindr2) &  
 
1445
    :: aa,rex,imx
 
1446
  integer :: nn,ii,jj
 
1447
!
 
1448
  rex = areal(xx)
 
1449
  imx = aimag(xx)
 
1450
!
 
1451
  if (rex.eq.RZRO.and.imx.eq.RZRO) then
 
1452
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_c: ' &
 
1453
       ,'xx = 0, returning 0'
 
1454
    rslt = 0
 
1455
    return
 
1456
  endif
 
1457
!
 
1458
  if (mod(iph,2).eq.0) then ;yy= xx ;jj=iph
 
1459
                       else ;yy=-xx ;jj=iph+sgnRe(imx)
 
1460
  endif
 
1461
!
 
1462
  if (jj.ne.0) then
 
1463
    rslt = ( olog1(yy,jj) - 1 )/(yy-1)
 
1464
    return
 
1465
  endif
 
1466
!
 
1467
  zz = yy-1
 
1468
  aa = abs(zz)
 
1469
  if     (aa.ge.thrs(6,prcpar)) then
 
1470
    rslt = (log(yy)/zz-1)/zz
 
1471
    return
 
1472
  elseif (aa.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
1473
  elseif (aa.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
1474
  elseif (aa.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
1475
  elseif (aa.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
1476
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
1477
                                else ;nn=ntrm(1,prcpar)
 
1478
  endif
 
1479
! convergence acceleration using  z=(y-1)/(y+1)
 
1480
! rslt = -1/(y+1) + 2/(y+1)^2 * ( z/3 + z^3/5 + z^5/7 + ... )  
 
1481
  zz = zz/(yy+1)
 
1482
  z2 = zz*zz
 
1483
  aa = 2
 
1484
  nn = 2*nn-1
 
1485
  rslt = aa/nn
 
1486
  do ii=nn-2,3,-2
 
1487
    rslt = aa/ii + z2*rslt
 
1488
  enddo
 
1489
  rslt = ( -1 + zz*rslt/(yy+1) )/(yy+1)
 
1490
  end function
 
1491
 
 
1492
 
 
1493
  function log3_r(xx,iph) result(rslt)
 
1494
!***********************************************************************
 
1495
!***********************************************************************
 
1496
  real(kindr2) &  
 
1497
          ,intent(in) :: xx
 
1498
  integer ,intent(in) :: iph
 
1499
  complex(kindr2) &   
 
1500
    :: rslt
 
1501
  rslt = log3_c(xx*CONE,iph)
 
1502
  end function
 
1503
 
 
1504
 
 
1505
  function log3_c(xx,iph) result(rslt)
 
1506
!***********************************************************************
 
1507
!***********************************************************************
 
1508
  complex(kindr2) &   
 
1509
          ,intent(in) :: xx
 
1510
  integer ,intent(in) :: iph
 
1511
  complex(kindr2) &   
 
1512
    :: rslt ,yy,zz,z2,HLF
 
1513
  real(kindr2) &  
 
1514
    :: aa,rex,imx
 
1515
  integer :: nn,ii,jj
 
1516
!
 
1517
  rex = areal(xx)
 
1518
  imx = aimag(xx)
 
1519
!
 
1520
  if (rex.eq.RZRO.and.imx.eq.RZRO) then
 
1521
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_c: ' &
 
1522
       ,'xx = 0, returning 0'
 
1523
    rslt = 0
 
1524
    return
 
1525
  endif
 
1526
!
 
1527
  HLF = CONE/2
 
1528
!
 
1529
  if (mod(iph,2).eq.0) then ;yy= xx ;jj=iph
 
1530
                       else ;yy=-xx ;jj=iph+sgnRe(imx)
 
1531
  endif
 
1532
!
 
1533
  if (jj.ne.0) then
 
1534
    rslt = ( olog2(xx,jj) + HLF )/(yy-1)
 
1535
    return
 
1536
  endif
 
1537
!
 
1538
  zz = yy-1
 
1539
  aa = abs(zz)
 
1540
  if     (aa.ge.thrs(6,prcpar)) then
 
1541
    rslt = ((log(yy)/zz-1)/zz+HLF)/zz
 
1542
    return
 
1543
  elseif (aa.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
1544
  elseif (aa.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
1545
  elseif (aa.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
1546
  elseif (aa.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
1547
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
1548
                                else ;nn=ntrm(1,prcpar)
 
1549
  endif
 
1550
! convergence acceleration using  z=(y-1)/(y+1)
 
1551
! rslt = 1/(2*(y+1)) + 2/(y+1)^3 * ( 1/3 + z^2/5 + z^4/7 + ... )
 
1552
  zz = zz/(yy+1)
 
1553
  z2 = zz*zz
 
1554
  aa = 2
 
1555
  nn = 2*nn-1
 
1556
  rslt = aa/nn
 
1557
  do ii=nn-2,3,-2
 
1558
    rslt = aa/ii + z2*rslt
 
1559
  enddo
 
1560
  rslt = ( HLF + rslt/(yy+1)**2 )/(yy+1)
 
1561
  end function
 
1562
 
1406
1563
end module
1407
1564
 
1408
1565
 
 
1566
 
 
1567
 
1409
1568
module avh_olo_dp_dilog
1410
1569
!***********************************************************************
1411
1570
!                     /1    ln(1-zz*t)
1660
1819
! 
1661
1820
  if (yy.eq.RONE.and.odd.eq.0) then
1662
1821
    if (ntwo.ne.0) then
1663
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog_r: ' &
 
1822
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog_r: ' &
1664
1823
        ,'|x|,iph = ',trim(myprint(yy)),',',jj,', returning 0'
1665
1824
    endif
1666
1825
    rslt = 0
1768
1927
!
1769
1928
  if (j1.ne.j2) then
1770
1929
    if (r1.eq.r2) then
1771
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
1930
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
1772
1931
        ,'j1,j2,r1-r2',j1,j2,',',trim(myprint(r1-r2)),', returning 0'
1773
1932
      rslt = 0
1774
1933
!      write(*,*) 'dilog2_c j1=/=j2,r1=r2' !DEBUG
1782
1941
!
1783
1942
  if (a1.lt.eps) then
1784
1943
    if (a2.lt.eps) then
1785
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
1944
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
1786
1945
        ,'r1,r2 =',trim(myprint(r1)),',',trim(myprint(r2)),', returning 0'
1787
1946
      rslt = 0
1788
1947
!      write(*,*) 'dilog2_c r1<eps,r2<eps' !DEBUG
1804
1963
!      write(*,*) 'dilog2_c ||1-y1|/|1-y2|-1|>0.1' !DEBUG
1805
1964
      return
1806
1965
    elseif (oo.eq.0.and.ao1.lt.eps) then
1807
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
1966
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
1808
1967
        ,'r1,oo,nn =',trim(myprint(r1)),',',oo,nn,', putting nn=0'
1809
1968
      if (ao2.lt.eps) then
1810
1969
        rslt = -1
1814
1973
        y1=1-eps ;nn=0 ;logr1=0 ;r1=1-eps
1815
1974
      endif
1816
1975
    elseif (oo.eq.0.and.ao2.lt.eps) then
1817
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
1976
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
1818
1977
        ,'r2,oo,nn =',trim(myprint(r2)),',',oo,nn,', putting nn=0'
1819
1978
      y2=1-eps ;nn=0 ;logr2=0 ;r2=1-eps
1820
1979
    endif
1829
1988
      if (a1.gt.RONE) ii = ii + (nn+pp(oo,sgnIm(y2)))
1830
1989
      if (a2.gt.RONE) ii = ii - (nn+pp(oo,sgnIm(y2)))
1831
1990
      ii = nn*ii
1832
 
      if (ii.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
1991
      if (ii.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
1833
1992
        ,'r1,r2,nn =',trim(myprint(r1)),',',trim(myprint(r2)),',',nn &
1834
1993
        ,', putting nn=0'
1835
 
      rslt = -olog2(y2,0)
 
1994
      rslt = -olog1(y2,0)
1836
1995
!      write(*,*) 'dilog2_c |logr1/lorg2|<eps' !DEBUG
1837
1996
      return
1838
1997
    endif
1844
2003
    nn=-nn ;oo=-oo
1845
2004
  endif
1846
2005
!
1847
 
  ff=y1/y2         ;ff=-olog2(ff,0)/y2
1848
 
  gg=(1-y1)/(1-y2) ;gg=-olog2(gg,0)/(1-y2)
 
2006
  ff=y1/y2         ;ff=-olog1(ff,0)/y2
 
2007
  gg=(1-y1)/(1-y2) ;gg=-olog1(gg,0)/(1-y2)
1849
2008
!
1850
2009
  if (2*areal(y1).ge.RONE) then
1851
2010
!    write(*,*) 'dilog2_c re>1/2' !DEBUG
1902
2061
!
1903
2062
  if (j1.ne.j2) then
1904
2063
    if (r1.eq.r2) then
1905
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
2064
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
1906
2065
        ,'j1,j2,r1-r2',j1,j2,',',trim(myprint(r1-r2)),', returning 0'
1907
2066
      rslt = 0
1908
2067
!      write(*,*) 'dilog2_r j1=/=j2,r1=r2' !DEBUG
1916
2075
!
1917
2076
  if (r1.lt.eps) then
1918
2077
    if (r2.lt.eps) then
1919
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
2078
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
1920
2079
        ,'r1,r2 =',trim(myprint(r1)),',',trim(myprint(r2)),', returning 0'
1921
2080
      rslt = 0
1922
2081
!      write(*,*) 'dilog2_r r1<eps,r2<eps' !DEBUG
1938
2097
!      write(*,*) 'dilog2_r ||1-y1|/|1-y2|-1|>0.1' !DEBUG
1939
2098
      return
1940
2099
    elseif (oo.eq.0.and.ro1.lt.eps) then
1941
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
2100
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
1942
2101
        ,'r1,oo,nn =',trim(myprint(r1)),',',oo,nn,', putting nn=0'
1943
2102
      if (ro2.lt.eps) then
1944
2103
        rslt = -1
1948
2107
        y1=1-eps ;nn=0 ;logr1=0 ;r1=1-eps
1949
2108
      endif
1950
2109
    elseif (oo.eq.0.and.ro2.lt.eps) then
1951
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
2110
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
1952
2111
        ,'r2,oo,nn =',trim(myprint(r2)),',',oo,nn,', putting nn=0'
1953
2112
      y2=1-eps ;nn=0 ;logr2=0 ;r2=1-eps
1954
2113
    endif
1963
2122
      if (r1.gt.RONE) ii = ii + (nn+2*oo)
1964
2123
      if (r2.gt.RONE) ii = ii - (nn+2*oo)
1965
2124
      ii = nn*ii
1966
 
      if (ii.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
2125
      if (ii.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
1967
2126
        ,'r1,r2,nn =',trim(myprint(r1)),',',trim(myprint(r2)),',',nn &
1968
2127
        ,', putting nn=0'
1969
 
      rslt = -olog2(y2,0)
 
2128
      rslt = -olog1(y2,0)
1970
2129
!      write(*,*) 'dilog2_r |logr1/lorg2|<eps' !DEBUG
1971
2130
      return
1972
2131
    endif
1978
2137
    nn=-nn ;oo=-oo
1979
2138
  endif
1980
2139
!
1981
 
  ff=y1/y2         ;ff=-olog2(ff,0)/y2
1982
 
  gg=(1-y1)/(1-y2) ;gg=-olog2(gg,0)/(1-y2)
 
2140
  ff=y1/y2         ;ff=-olog1(ff,0)/y2
 
2141
  gg=(1-y1)/(1-y2) ;gg=-olog1(gg,0)/(1-y2)
1983
2142
!
1984
2143
  if (2*y1.ge.RONE) then
1985
2144
!    write(*,*) 'dilog2_r re>1/2' !DEBUG
2409
2568
 
2410
2569
  implicit none
2411
2570
  private
2412
 
  public :: qmplx_type,qonv,directly,sheet,logc,logc2,li2c,li2c2
 
2571
  public :: qmplx_type,qonv,directly,sheet,logc,logc2,logc3,li2c,li2c2
2413
2572
  public :: operator (*) ,operator (/)
2414
2573
 
2415
2574
  type :: qmplx_type
2656
2815
  type(qmplx_type) ,intent(in) :: xx
2657
2816
  complex(kindr2) &   
2658
2817
    :: rslt
2659
 
!  rslt = -olog2(acmplx(xx%c),xx%p)
2660
 
  rslt = -olog2(xx%c,xx%p)
 
2818
!  rslt = -olog1(acmplx(xx%c),xx%p)
 
2819
  rslt = -olog1(xx%c,xx%p)
 
2820
  end function
 
2821
 
 
2822
  function logc3(xx) result(rslt)
 
2823
!*******************************************************************
 
2824
!  ( log(xx)/(1-xx) + 1 )/(1-xx)
 
2825
!*******************************************************************
 
2826
  type(qmplx_type) ,intent(in) :: xx
 
2827
  complex(kindr2) &   
 
2828
    :: rslt
 
2829
!  rslt = olog2(acmplx(xx%c),xx%p)
 
2830
  rslt = olog2(xx%c,xx%p)
2661
2831
  end function
2662
2832
 
2663
2833
  function li2c(xx) result(rslt)
2697
2867
  use avh_olo_dp_auxfun
2698
2868
  use avh_olo_dp_bnlog
2699
2869
  use avh_olo_dp_qmplx
 
2870
  use avh_olo_dp_olog
2700
2871
  implicit none
2701
2872
  private
2702
 
  public :: tadp ,tadpn ,bub0 ,bub1 ,bub11 ,bub111 ,bub1111
 
2873
  public :: tadp ,tadpn ,bub0 ,dbub0 ,bub1 ,bub11 ,bub111 ,bub1111
2703
2874
 
2704
2875
contains
2705
2876
 
3290
3461
  b0011(0) = ( a0(0,0) - x1*b111(0) + x2*b11(0) + 4*b0011(1) )/10
3291
3462
  end subroutine
3292
3463
 
 
3464
 
 
3465
!*******************************************************************
 
3466
! Derivative of B0
 
3467
! expects  m0<m1
 
3468
! only finite case, so input must not be  m0=0 & m1=pp
 
3469
!*******************************************************************
 
3470
 
 
3471
  subroutine dbub0( rslt &
 
3472
                   ,pp,m0,m1 ,app,am0,am1 )
 
3473
  complex(kindr2) &   
 
3474
    ,intent(out) :: rslt
 
3475
  complex(kindr2) &   
 
3476
    ,intent(in)  :: pp,m0,m1
 
3477
  real(kindr2) &  
 
3478
    ,intent(in)  :: app,am0,am1
 
3479
  complex(kindr2) &   
 
3480
    :: ch,x1,x2,lambda
 
3481
  real(kindr2) &  
 
3482
    :: ax1,ax2,ax1x2,maxa
 
3483
  type(qmplx_type) :: q1,q2,q1o,q2o
 
3484
  integer :: sgn
 
3485
!
 
3486
  if (am1.eq.RZRO) then
 
3487
    if (app.eq.RZRO) then
 
3488
      rslt = 0
 
3489
      return
 
3490
    endif
 
3491
  endif
 
3492
!
 
3493
  if (app.eq.RZRO) then
 
3494
    if (abs(m0-m1).le.am1*EPSN*10) then
 
3495
      rslt = 1/(6*m1)
 
3496
    else
 
3497
      ch = m0/m1
 
3498
      rslt = ( CONE/2 - ch*olog3(ch,0) )/m1 
 
3499
    endif
 
3500
  elseif (am1.eq.RZRO) then
 
3501
    rslt =-1/pp
 
3502
  else
 
3503
    call solabc( x1,x2 ,lambda ,pp ,(m0-m1)-pp ,m1 ,0 )
 
3504
    sgn =-sgnRe(pp)*sgnRe(x2-x1)
 
3505
    q1  = qonv(x1  , sgn)
 
3506
    q1o = qonv(x1-1, sgn)
 
3507
    q2  = qonv(x2  ,-sgn)
 
3508
    q2o = qonv(x2-1,-sgn)
 
3509
    ax1 = abs(x1)
 
3510
    ax2 = abs(x2)
 
3511
    ax1x2 = abs(x1-x2)
 
3512
    maxa = max(ax1,ax2)
 
3513
    if (ax1x2.lt.maxa*EPSN*10) then
 
3514
      rslt = ( (x1+x2-1)*logc(q2/q2o) - 2 )/pp
 
3515
    elseif (ax1x2*2.lt.maxa) then
 
3516
      if     (x1.eq.CZRO.or.x1.eq.CONE) then
 
3517
        rslt = ( (x1+x2-1)*logc(q2/q2o) - 1 )/pp
 
3518
      elseif (x2.eq.CZRO.or.x2.eq.CONE) then
 
3519
        rslt = ( (x1+x2-1)*logc(q1/q1o) - 1 )/pp
 
3520
      else
 
3521
        rslt = x1*(x1-1)*( logc2(q1o/q2o)/(x2-1) - logc2(q1/q2)/x2 ) &
 
3522
             + (x1+x2-1)*logc(q2/q2o) - 1
 
3523
        rslt = rslt/pp
 
3524
      endif
 
3525
    else
 
3526
      rslt = 0
 
3527
      if (ax1.ne.RZRO) then
 
3528
        if (ax1.lt.2*RONE) then
 
3529
          rslt = rslt - x1
 
3530
          if (x1.ne.CONE) rslt = rslt - x1*logc2(q1/q1o)
 
3531
        else
 
3532
          rslt = rslt + x1/(x1-1)*logc3(q1/q1o)
 
3533
        endif
 
3534
      endif
 
3535
      if (ax2.ne.RZRO) then
 
3536
        if (ax2.lt.2*RONE) then
 
3537
          rslt = rslt + x2
 
3538
          if (x2.ne.CONE) rslt = rslt + x2*logc2(q2/q2o)
 
3539
        else
 
3540
          rslt = rslt - x2/(x2-1)*logc3(q2/q2o)
 
3541
        endif
 
3542
      endif
 
3543
      rslt = rslt/lambda
 
3544
    endif
 
3545
  endif
 
3546
!
 
3547
  end subroutine
 
3548
 
 
3549
 
3293
3550
end module
3294
3551
 
3295
3552
 
3523
3780
    log2 = olog( abs(rp2/rmu2) ,i2 )
3524
3781
    log3 = olog( abs(rp3/rmu2) ,i3 )
3525
3782
    rslt(2) = 0
3526
 
    rslt(1) = -olog2( abs(rp3/rp2) ,i3-i2 )/rp2
 
3783
    rslt(1) = -olog1( abs(rp3/rp2) ,i3-i2 )/rp2
3527
3784
    rslt(0) = -rslt(1)*(log3+log2)/2
3528
3785
  elseif (icase.eq.3) then
3529
3786
! 3 masses non-zero
6096
6353
     ,intent(in) :: y1,y2
6097
6354
  complex(kindr2) &   
6098
6355
     :: rslt ,oy1,oy2
 
6356
!
6099
6357
   oy1 = 1-y1
6100
6358
   oy2 = 1-y2
6101
6359
   rslt = logc2( qonv(-y2)/qonv(-y1) )/y1 &
6164
6422
  private
6165
6423
  public :: olo_unit ,olo_scale ,olo_onshell ,olo_setting
6166
6424
  public :: olo_precision
6167
 
  public :: olo_a0 ,olo_b0 ,olo_b11 ,olo_c0 ,olo_d0
 
6425
  public :: olo_a0 ,olo_b0 ,olo_db0 ,olo_b11 ,olo_c0 ,olo_d0
6168
6426
  public :: olo_an ,olo_bn
6169
6427
  public :: olo
6170
6428
  public :: olo_get_scale ,olo_get_onshell ,olo_get_precision
6171
 
  public :: a0_r,a0rr,a0_c,a0cr
6172
 
  public :: an_r,anrr,an_c,ancr
6173
 
  public :: b0rr,b0rrr,b0rc,b0rcr,b0cc,b0ccr
6174
 
  public :: b11rr,b11rrr,b11rc,b11rcr,b11cc,b11ccr
6175
 
  public :: bnrr,bnrrr,bnrc,bnrcr,bncc,bnccr
6176
 
  public :: c0rr,c0rrr,c0rc,c0rcr,c0cc,c0ccr
6177
 
  public :: d0rr,d0rrr,d0rc,d0rcr,d0cc,d0ccr
6178
6429
!
6179
6430
  integer ,public ,parameter :: olo_kind=kindr2    
6180
6431
!
6200
6451
  interface olo_b0
6201
6452
    module procedure b0rr,b0rrr,b0rc,b0rcr,b0cc,b0ccr
6202
6453
  end interface 
 
6454
  interface olo_db0
 
6455
    module procedure db0rr,db0rrr,db0rc,db0rcr,db0cc,db0ccr
 
6456
  end interface 
6203
6457
  interface olo_b11
6204
6458
    module procedure b11rr,b11rrr,b11rc,b11rcr,b11cc,b11ccr
6205
6459
  end interface 
7186
7440
  endif
7187
7441
  end subroutine
7188
7442
 
 
7443
!*******************************************************************
 
7444
! Derivative of B0
 
7445
!*******************************************************************
 
7446
  subroutine db0cc( rslt ,pp,m1,m2 )
 
7447
!
 
7448
  use avh_olo_dp_bub ,only: dbub0
 
7449
  use avh_olo_dp_olog ,only: olog
 
7450
!
 
7451
  complex(kindr2) &   
 
7452
    ,intent(out) :: rslt(0:2)
 
7453
  complex(kindr2) &   
 
7454
    ,intent(in)  :: pp
 
7455
  complex(kindr2) &   
 
7456
    ,intent(in)  :: m1,m2
 
7457
!
 
7458
  complex(kindr2) &   
 
7459
    :: ss,r1,r2,ch
 
7460
  real(kindr2) &  
 
7461
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
7462
  character(26+99) ,parameter :: warning=&
 
7463
                     'WARNING from OneLOop db0: '//warnonshell
 
7464
  if (initz) call init
 
7465
  ss = pp
 
7466
  r1 = m1
 
7467
  r2 = m2
 
7468
!
 
7469
  app = areal(ss)
 
7470
  if (aimag(ss).ne.RZRO) then
 
7471
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
7472
      ,'ss has non-zero imaginary part, putting it to zero.'
 
7473
    ss = acmplx( app )
 
7474
  endif
 
7475
  app = abs(app)
 
7476
!
 
7477
  am1 = areal(r1)
 
7478
  hh  = aimag(r1)
 
7479
  if (hh.gt.RZRO) then
 
7480
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
7481
      ,'r1 has positive imaginary part, switching its sign.'
 
7482
    r1 = acmplx( am1 ,-hh )
 
7483
  endif
 
7484
  am1 = abs(am1) + abs(hh)
 
7485
!
 
7486
  am2 = areal(r2)
 
7487
  hh  = aimag(r2)
 
7488
  if (hh.gt.RZRO) then
 
7489
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
7490
      ,'r2 has positive imaginary part, switching its sign.'
 
7491
    r2 = acmplx( am2 ,-hh )
 
7492
  endif
 
7493
  am2 = abs(am2) + abs(hh)
 
7494
!
 
7495
  mulocal = muscale 
 
7496
!
 
7497
  mulocal2 = mulocal*mulocal
 
7498
!
 
7499
  if (am1.gt.am2) then
 
7500
    ch=r1 ; r1=r2 ; r2=ch
 
7501
    hh=am1;am1=am2;am2=hh
 
7502
  endif
 
7503
  ssr2 = abs(ss-r2)
 
7504
!
 
7505
  if (nonzerothrs) then
 
7506
    hh = onshellthrs
 
7507
    if (app.lt.hh) app = 0
 
7508
    if (am1.lt.hh) am1 = 0
 
7509
    if (am2.lt.hh) am2 = 0
 
7510
    if (ssr2.lt.hh) ssr2 = 0
 
7511
  elseif (wunit.gt.0) then
 
7512
    hh = onshellthrs*max(app,max(am1,am2))
 
7513
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
7514
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
7515
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
7516
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
7517
  endif
 
7518
!
 
7519
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
7520
!
 
7521
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
7522
    rslt(1) =-1/(2*ss)
 
7523
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
7524
  else
 
7525
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
7526
  endif
 
7527
!
 
7528
  if (punit.gt.0) then
 
7529
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
7530
    write(punit,*) ' pp:',trim(myprint(pp))
 
7531
    write(punit,*) ' m1:',trim(myprint(m1))
 
7532
    write(punit,*) ' m2:',trim(myprint(m2))
 
7533
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
7534
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
7535
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
7536
  endif
 
7537
  end subroutine
 
7538
 
 
7539
  subroutine db0ccr( rslt ,pp,m1,m2 ,rmu )
 
7540
!
 
7541
  use avh_olo_dp_bub ,only: dbub0
 
7542
  use avh_olo_dp_olog ,only: olog
 
7543
!
 
7544
  complex(kindr2) &   
 
7545
    ,intent(out) :: rslt(0:2)
 
7546
  complex(kindr2) &   
 
7547
    ,intent(in)  :: pp
 
7548
  complex(kindr2) &   
 
7549
    ,intent(in)  :: m1,m2
 
7550
  real(kindr2) &  
 
7551
   ,intent(in)  :: rmu       
 
7552
!
 
7553
  complex(kindr2) &   
 
7554
    :: ss,r1,r2,ch
 
7555
  real(kindr2) &  
 
7556
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
7557
  character(26+99) ,parameter :: warning=&
 
7558
                     'WARNING from OneLOop db0: '//warnonshell
 
7559
  if (initz) call init
 
7560
  ss = pp
 
7561
  r1 = m1
 
7562
  r2 = m2
 
7563
!
 
7564
  app = areal(ss)
 
7565
  if (aimag(ss).ne.RZRO) then
 
7566
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
7567
      ,'ss has non-zero imaginary part, putting it to zero.'
 
7568
    ss = acmplx( app )
 
7569
  endif
 
7570
  app = abs(app)
 
7571
!
 
7572
  am1 = areal(r1)
 
7573
  hh  = aimag(r1)
 
7574
  if (hh.gt.RZRO) then
 
7575
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
7576
      ,'r1 has positive imaginary part, switching its sign.'
 
7577
    r1 = acmplx( am1 ,-hh )
 
7578
  endif
 
7579
  am1 = abs(am1) + abs(hh)
 
7580
!
 
7581
  am2 = areal(r2)
 
7582
  hh  = aimag(r2)
 
7583
  if (hh.gt.RZRO) then
 
7584
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
7585
      ,'r2 has positive imaginary part, switching its sign.'
 
7586
    r2 = acmplx( am2 ,-hh )
 
7587
  endif
 
7588
  am2 = abs(am2) + abs(hh)
 
7589
!
 
7590
  mulocal = rmu     
 
7591
!
 
7592
  mulocal2 = mulocal*mulocal
 
7593
!
 
7594
  if (am1.gt.am2) then
 
7595
    ch=r1 ; r1=r2 ; r2=ch
 
7596
    hh=am1;am1=am2;am2=hh
 
7597
  endif
 
7598
  ssr2 = abs(ss-r2)
 
7599
!
 
7600
  if (nonzerothrs) then
 
7601
    hh = onshellthrs
 
7602
    if (app.lt.hh) app = 0
 
7603
    if (am1.lt.hh) am1 = 0
 
7604
    if (am2.lt.hh) am2 = 0
 
7605
    if (ssr2.lt.hh) ssr2 = 0
 
7606
  elseif (wunit.gt.0) then
 
7607
    hh = onshellthrs*max(app,max(am1,am2))
 
7608
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
7609
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
7610
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
7611
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
7612
  endif
 
7613
!
 
7614
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
7615
!
 
7616
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
7617
    rslt(1) =-1/(2*ss)
 
7618
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
7619
  else
 
7620
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
7621
  endif
 
7622
!
 
7623
  if (punit.gt.0) then
 
7624
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
7625
    write(punit,*) ' pp:',trim(myprint(pp))
 
7626
    write(punit,*) ' m1:',trim(myprint(m1))
 
7627
    write(punit,*) ' m2:',trim(myprint(m2))
 
7628
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
7629
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
7630
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
7631
  endif
 
7632
  end subroutine
 
7633
 
 
7634
  subroutine db0rc( rslt ,pp ,m1,m2 )
 
7635
!
 
7636
  use avh_olo_dp_bub ,only: dbub0
 
7637
  use avh_olo_dp_olog ,only: olog
 
7638
!
 
7639
  complex(kindr2) &   
 
7640
    ,intent(out) :: rslt(0:2)
 
7641
  real(kindr2) &  
 
7642
    ,intent(in)  :: pp
 
7643
  complex(kindr2) &   
 
7644
    ,intent(in)  :: m1,m2
 
7645
!
 
7646
  complex(kindr2) &   
 
7647
    :: ss,r1,r2,ch
 
7648
  real(kindr2) &  
 
7649
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
7650
  character(26+99) ,parameter :: warning=&
 
7651
                     'WARNING from OneLOop db0: '//warnonshell
 
7652
  if (initz) call init
 
7653
  ss = pp
 
7654
  r1 = m1
 
7655
  r2 = m2
 
7656
!
 
7657
  app = abs(pp)
 
7658
!
 
7659
  am1 = areal(r1)
 
7660
  hh  = aimag(r1)
 
7661
  if (hh.gt.RZRO) then
 
7662
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
7663
      ,'r1 has positive imaginary part, switching its sign.'
 
7664
    r1 = acmplx( am1 ,-hh )
 
7665
  endif
 
7666
  am1 = abs(am1) + abs(hh)
 
7667
!
 
7668
  am2 = areal(r2)
 
7669
  hh  = aimag(r2)
 
7670
  if (hh.gt.RZRO) then
 
7671
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
7672
      ,'r2 has positive imaginary part, switching its sign.'
 
7673
    r2 = acmplx( am2 ,-hh )
 
7674
  endif
 
7675
  am2 = abs(am2) + abs(hh)
 
7676
!
 
7677
  mulocal = muscale 
 
7678
!
 
7679
  mulocal2 = mulocal*mulocal
 
7680
!
 
7681
  if (am1.gt.am2) then
 
7682
    ch=r1 ; r1=r2 ; r2=ch
 
7683
    hh=am1;am1=am2;am2=hh
 
7684
  endif
 
7685
  ssr2 = abs(ss-r2)
 
7686
!
 
7687
  if (nonzerothrs) then
 
7688
    hh = onshellthrs
 
7689
    if (app.lt.hh) app = 0
 
7690
    if (am1.lt.hh) am1 = 0
 
7691
    if (am2.lt.hh) am2 = 0
 
7692
    if (ssr2.lt.hh) ssr2 = 0
 
7693
  elseif (wunit.gt.0) then
 
7694
    hh = onshellthrs*max(app,max(am1,am2))
 
7695
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
7696
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
7697
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
7698
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
7699
  endif
 
7700
!
 
7701
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
7702
!
 
7703
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
7704
    rslt(1) =-1/(2*ss)
 
7705
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
7706
  else
 
7707
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
7708
  endif
 
7709
!
 
7710
  if (punit.gt.0) then
 
7711
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
7712
    write(punit,*) ' pp:',trim(myprint(pp))
 
7713
    write(punit,*) ' m1:',trim(myprint(m1))
 
7714
    write(punit,*) ' m2:',trim(myprint(m2))
 
7715
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
7716
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
7717
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
7718
  endif
 
7719
  end subroutine
 
7720
 
 
7721
  subroutine db0rcr( rslt ,pp,m1,m2 ,rmu )
 
7722
!
 
7723
  use avh_olo_dp_bub ,only: dbub0
 
7724
  use avh_olo_dp_olog ,only: olog
 
7725
!
 
7726
  complex(kindr2) &   
 
7727
    ,intent(out) :: rslt(0:2)
 
7728
  real(kindr2) &  
 
7729
    ,intent(in)  :: pp
 
7730
  complex(kindr2) &   
 
7731
    ,intent(in)  :: m1,m2
 
7732
  real(kindr2) &  
 
7733
   ,intent(in)  :: rmu       
 
7734
!
 
7735
  complex(kindr2) &   
 
7736
    :: ss,r1,r2,ch
 
7737
  real(kindr2) &  
 
7738
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
7739
  character(26+99) ,parameter :: warning=&
 
7740
                     'WARNING from OneLOop db0: '//warnonshell
 
7741
  if (initz) call init
 
7742
  ss = pp
 
7743
  r1 = m1
 
7744
  r2 = m2
 
7745
!
 
7746
  app = abs(pp)
 
7747
!
 
7748
  am1 = areal(r1)
 
7749
  hh  = aimag(r1)
 
7750
  if (hh.gt.RZRO) then
 
7751
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
7752
      ,'r1 has positive imaginary part, switching its sign.'
 
7753
    r1 = acmplx( am1 ,-hh )
 
7754
  endif
 
7755
  am1 = abs(am1) + abs(hh)
 
7756
!
 
7757
  am2 = areal(r2)
 
7758
  hh  = aimag(r2)
 
7759
  if (hh.gt.RZRO) then
 
7760
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
7761
      ,'r2 has positive imaginary part, switching its sign.'
 
7762
    r2 = acmplx( am2 ,-hh )
 
7763
  endif
 
7764
  am2 = abs(am2) + abs(hh)
 
7765
!
 
7766
  mulocal = rmu     
 
7767
!
 
7768
  mulocal2 = mulocal*mulocal
 
7769
!
 
7770
  if (am1.gt.am2) then
 
7771
    ch=r1 ; r1=r2 ; r2=ch
 
7772
    hh=am1;am1=am2;am2=hh
 
7773
  endif
 
7774
  ssr2 = abs(ss-r2)
 
7775
!
 
7776
  if (nonzerothrs) then
 
7777
    hh = onshellthrs
 
7778
    if (app.lt.hh) app = 0
 
7779
    if (am1.lt.hh) am1 = 0
 
7780
    if (am2.lt.hh) am2 = 0
 
7781
    if (ssr2.lt.hh) ssr2 = 0
 
7782
  elseif (wunit.gt.0) then
 
7783
    hh = onshellthrs*max(app,max(am1,am2))
 
7784
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
7785
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
7786
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
7787
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
7788
  endif
 
7789
!
 
7790
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
7791
!
 
7792
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
7793
    rslt(1) =-1/(2*ss)
 
7794
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
7795
  else
 
7796
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
7797
  endif
 
7798
!
 
7799
  if (punit.gt.0) then
 
7800
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
7801
    write(punit,*) ' pp:',trim(myprint(pp))
 
7802
    write(punit,*) ' m1:',trim(myprint(m1))
 
7803
    write(punit,*) ' m2:',trim(myprint(m2))
 
7804
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
7805
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
7806
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
7807
  endif
 
7808
  end subroutine
 
7809
 
 
7810
  subroutine db0rr( rslt ,pp ,m1,m2 )
 
7811
!
 
7812
  use avh_olo_dp_bub ,only: dbub0
 
7813
  use avh_olo_dp_olog ,only: olog
 
7814
!
 
7815
  complex(kindr2) &   
 
7816
    ,intent(out) :: rslt(0:2)
 
7817
  real(kindr2) &  
 
7818
    ,intent(in)  :: pp
 
7819
  real(kindr2) &  
 
7820
    ,intent(in)  :: m1,m2
 
7821
!
 
7822
  complex(kindr2) &   
 
7823
    :: ss,r1,r2,ch
 
7824
  real(kindr2) &  
 
7825
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
7826
  character(26+99) ,parameter :: warning=&
 
7827
                     'WARNING from OneLOop db0: '//warnonshell
 
7828
  if (initz) call init
 
7829
  ss = pp
 
7830
  r1 = m1
 
7831
  r2 = m2
 
7832
!
 
7833
  app = abs(pp)
 
7834
!
 
7835
  am1 = abs(m1)
 
7836
  am2 = abs(m2)
 
7837
!
 
7838
  mulocal = muscale 
 
7839
!
 
7840
  mulocal2 = mulocal*mulocal
 
7841
!
 
7842
  if (am1.gt.am2) then
 
7843
    ch=r1 ; r1=r2 ; r2=ch
 
7844
    hh=am1;am1=am2;am2=hh
 
7845
  endif
 
7846
  ssr2 = abs(ss-r2)
 
7847
!
 
7848
  if (nonzerothrs) then
 
7849
    hh = onshellthrs
 
7850
    if (app.lt.hh) app = 0
 
7851
    if (am1.lt.hh) am1 = 0
 
7852
    if (am2.lt.hh) am2 = 0
 
7853
    if (ssr2.lt.hh) ssr2 = 0
 
7854
  elseif (wunit.gt.0) then
 
7855
    hh = onshellthrs*max(app,max(am1,am2))
 
7856
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
7857
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
7858
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
7859
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
7860
  endif
 
7861
!
 
7862
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
7863
!
 
7864
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
7865
    rslt(1) =-1/(2*ss)
 
7866
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
7867
  else
 
7868
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
7869
  endif
 
7870
!
 
7871
  if (punit.gt.0) then
 
7872
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
7873
    write(punit,*) ' pp:',trim(myprint(pp))
 
7874
    write(punit,*) ' m1:',trim(myprint(m1))
 
7875
    write(punit,*) ' m2:',trim(myprint(m2))
 
7876
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
7877
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
7878
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
7879
  endif
 
7880
  end subroutine
 
7881
 
 
7882
  subroutine db0rrr( rslt ,pp ,m1,m2 ,rmu )
 
7883
!
 
7884
  use avh_olo_dp_bub ,only: dbub0
 
7885
  use avh_olo_dp_olog ,only: olog
 
7886
!
 
7887
  complex(kindr2) &   
 
7888
    ,intent(out) :: rslt(0:2)
 
7889
  real(kindr2) &  
 
7890
    ,intent(in)  :: pp
 
7891
  real(kindr2) &  
 
7892
    ,intent(in)  :: m1,m2
 
7893
  real(kindr2) &  
 
7894
   ,intent(in)  :: rmu       
 
7895
!
 
7896
  complex(kindr2) &   
 
7897
    :: ss,r1,r2,ch
 
7898
  real(kindr2) &  
 
7899
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
7900
  character(26+99) ,parameter :: warning=&
 
7901
                     'WARNING from OneLOop db0: '//warnonshell
 
7902
  if (initz) call init
 
7903
  ss = pp
 
7904
  r1 = m1
 
7905
  r2 = m2
 
7906
!
 
7907
  app = abs(pp)
 
7908
!
 
7909
  am1 = abs(m1)
 
7910
  am2 = abs(m2)
 
7911
!
 
7912
  mulocal = rmu     
 
7913
!
 
7914
  mulocal2 = mulocal*mulocal
 
7915
!
 
7916
  if (am1.gt.am2) then
 
7917
    ch=r1 ; r1=r2 ; r2=ch
 
7918
    hh=am1;am1=am2;am2=hh
 
7919
  endif
 
7920
  ssr2 = abs(ss-r2)
 
7921
!
 
7922
  if (nonzerothrs) then
 
7923
    hh = onshellthrs
 
7924
    if (app.lt.hh) app = 0
 
7925
    if (am1.lt.hh) am1 = 0
 
7926
    if (am2.lt.hh) am2 = 0
 
7927
    if (ssr2.lt.hh) ssr2 = 0
 
7928
  elseif (wunit.gt.0) then
 
7929
    hh = onshellthrs*max(app,max(am1,am2))
 
7930
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
7931
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
7932
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
7933
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
7934
  endif
 
7935
!
 
7936
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
7937
!
 
7938
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
7939
    rslt(1) =-1/(2*ss)
 
7940
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
7941
  else
 
7942
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
7943
  endif
 
7944
!
 
7945
  if (punit.gt.0) then
 
7946
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
7947
    write(punit,*) ' pp:',trim(myprint(pp))
 
7948
    write(punit,*) ' m1:',trim(myprint(m1))
 
7949
    write(punit,*) ' m2:',trim(myprint(m2))
 
7950
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
7951
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
7952
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
7953
  endif
 
7954
  end subroutine
 
7955
 
 
7956
 
7189
7957
 
7190
7958
!*******************************************************************
7191
7959
! Return the Papparino-Veltman functions b11,b00,b1,b0 , for
9557
10325
               .or.(     areal(ss(1)).ge.-small  &
9558
10326
                    .and.areal(ss(2)).ge.-small  &
9559
10327
                    .and.areal(ss(3)).ge.-small  &
9560
 
                    .and.areal(ss(4)).ge.-small) )
 
10328
                    .and.areal(ss(4)).ge.-small) &
 
10329
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
9561
10330
    if (useboxc) then
9562
10331
      call boxc( rslt ,ss,rr ,as ,smax )
9563
10332
    else
9573
10342
                 .or.(     areal(ss(1)).ge.-small  &
9574
10343
                      .and.areal(ss(2)).ge.-small  &
9575
10344
                      .and.areal(ss(3)).ge.-small  &
9576
 
!OLD                      .and.areal(ss(4)).ge.-small) )
9577
 
                      .and.areal(ss(4)).ge.-small) & !NEW
9578
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
10345
                      .and.areal(ss(4)).ge.-small) &
 
10346
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
9579
10347
      if (useboxc) then
9580
10348
        call boxc( rslt ,ss,rr ,as ,smax )
9581
10349
      else
9832
10600
               .or.(     areal(ss(1)).ge.-small  &
9833
10601
                    .and.areal(ss(2)).ge.-small  &
9834
10602
                    .and.areal(ss(3)).ge.-small  &
9835
 
                    .and.areal(ss(4)).ge.-small) )
 
10603
                    .and.areal(ss(4)).ge.-small) &
 
10604
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
9836
10605
    if (useboxc) then
9837
10606
      call boxc( rslt ,ss,rr ,as ,smax )
9838
10607
    else
9848
10617
                 .or.(     areal(ss(1)).ge.-small  &
9849
10618
                      .and.areal(ss(2)).ge.-small  &
9850
10619
                      .and.areal(ss(3)).ge.-small  &
9851
 
!OLD                      .and.areal(ss(4)).ge.-small) )
9852
 
                      .and.areal(ss(4)).ge.-small) & !NEW
9853
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
10620
                      .and.areal(ss(4)).ge.-small) &
 
10621
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
9854
10622
      if (useboxc) then
9855
10623
        call boxc( rslt ,ss,rr ,as ,smax )
9856
10624
      else
10099
10867
               .or.(     areal(ss(1)).ge.-small  &
10100
10868
                    .and.areal(ss(2)).ge.-small  &
10101
10869
                    .and.areal(ss(3)).ge.-small  &
10102
 
                    .and.areal(ss(4)).ge.-small) )
 
10870
                    .and.areal(ss(4)).ge.-small) &
 
10871
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
10103
10872
    if (useboxc) then
10104
10873
      call boxc( rslt ,ss,rr ,as ,smax )
10105
10874
    else
10115
10884
                 .or.(     areal(ss(1)).ge.-small  &
10116
10885
                      .and.areal(ss(2)).ge.-small  &
10117
10886
                      .and.areal(ss(3)).ge.-small  &
10118
 
!OLD                      .and.areal(ss(4)).ge.-small) )
10119
 
                      .and.areal(ss(4)).ge.-small) & !NEW
10120
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
10887
                      .and.areal(ss(4)).ge.-small) &
 
10888
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
10121
10889
      if (useboxc) then
10122
10890
        call boxc( rslt ,ss,rr ,as ,smax )
10123
10891
      else
10368
11136
               .or.(     areal(ss(1)).ge.-small  &
10369
11137
                    .and.areal(ss(2)).ge.-small  &
10370
11138
                    .and.areal(ss(3)).ge.-small  &
10371
 
                    .and.areal(ss(4)).ge.-small) )
 
11139
                    .and.areal(ss(4)).ge.-small) &
 
11140
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
10372
11141
    if (useboxc) then
10373
11142
      call boxc( rslt ,ss,rr ,as ,smax )
10374
11143
    else
10384
11153
                 .or.(     areal(ss(1)).ge.-small  &
10385
11154
                      .and.areal(ss(2)).ge.-small  &
10386
11155
                      .and.areal(ss(3)).ge.-small  &
10387
 
!OLD                      .and.areal(ss(4)).ge.-small) )
10388
 
                      .and.areal(ss(4)).ge.-small) & !NEW
10389
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
11156
                      .and.areal(ss(4)).ge.-small) &
 
11157
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
10390
11158
      if (useboxc) then
10391
11159
        call boxc( rslt ,ss,rr ,as ,smax )
10392
11160
      else
10628
11396
               .or.(     areal(ss(1)).ge.-small  &
10629
11397
                    .and.areal(ss(2)).ge.-small  &
10630
11398
                    .and.areal(ss(3)).ge.-small  &
10631
 
                    .and.areal(ss(4)).ge.-small) )
 
11399
                    .and.areal(ss(4)).ge.-small) &
 
11400
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
10632
11401
    if (useboxc) then
10633
11402
      call boxc( rslt ,ss,rr ,as ,smax )
10634
11403
    else
10644
11413
                 .or.(     areal(ss(1)).ge.-small  &
10645
11414
                      .and.areal(ss(2)).ge.-small  &
10646
11415
                      .and.areal(ss(3)).ge.-small  &
10647
 
!OLD                      .and.areal(ss(4)).ge.-small) )
10648
 
                      .and.areal(ss(4)).ge.-small) & !NEW
10649
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
11416
                      .and.areal(ss(4)).ge.-small) &
 
11417
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
10650
11418
      if (useboxc) then
10651
11419
        call boxc( rslt ,ss,rr ,as ,smax )
10652
11420
      else
10890
11658
               .or.(     areal(ss(1)).ge.-small  &
10891
11659
                    .and.areal(ss(2)).ge.-small  &
10892
11660
                    .and.areal(ss(3)).ge.-small  &
10893
 
                    .and.areal(ss(4)).ge.-small) )
 
11661
                    .and.areal(ss(4)).ge.-small) &
 
11662
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
10894
11663
    if (useboxc) then
10895
11664
      call boxc( rslt ,ss,rr ,as ,smax )
10896
11665
    else
10906
11675
                 .or.(     areal(ss(1)).ge.-small  &
10907
11676
                      .and.areal(ss(2)).ge.-small  &
10908
11677
                      .and.areal(ss(3)).ge.-small  &
10909
 
!OLD                      .and.areal(ss(4)).ge.-small) )
10910
 
                      .and.areal(ss(4)).ge.-small) & !NEW
10911
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
11678
                      .and.areal(ss(4)).ge.-small) &
 
11679
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
10912
11680
      if (useboxc) then
10913
11681
        call boxc( rslt ,ss,rr ,as ,smax )
10914
11682
      else
11482
12250
  write(aa,'(i10)') min(len(cc),ndec+novh+1) ;aa=adjustl(aa)
11483
12251
  write(bb,'(i10)') min(len(cc),ndec       ) ;bb=adjustl(bb)
11484
12252
  aa = '(e'//trim(aa)//'.'//trim(bb)//')'
11485
 
  write(cc,aa) dble(xx)  ;cc=adjustl(cc)
 
12253
  write(cc,aa) xx  ;cc=adjustl(cc)
11486
12254
  if (cc(1:2).eq.'-0') then ;rslt = '-'//cc(3:len(cc))
11487
12255
  else                      ;rslt = ' '//cc(2:len(cc))
11488
12256
  endif
11709
12477
      gg=xd1*pq1 ;hh=yd1*uv1
11710
12478
      rx2 = gg+hh
11711
12479
      if (abs(rx2).lt.neglig(prcpar)*max(abs(gg),abs(hh))) rx2 = 0
11712
 
    else
 
12480
    elseif (abs(pq2).gt.abs(pq1)) then
11713
12481
      rx2 = pq2
11714
12482
      gg=xd2*pq2 ;hh=yd2*uv2
11715
12483
      rx1 = gg+hh
11716
12484
      if (abs(rx1).lt.neglig(prcpar)*max(abs(gg),abs(hh))) rx1 = 0
 
12485
    else
 
12486
      rx1 = pq1
 
12487
      rx2 = pq2
11717
12488
    endif
11718
12489
    if (abs(uv1).gt.abs(uv2)) then
11719
12490
      ix1 = uv1
11720
12491
      gg=yd1*pq1 ;hh=xd1*uv1
11721
12492
      ix2 = gg-hh
11722
12493
      if (abs(ix2).lt.neglig(prcpar)*max(abs(gg),abs(hh))) ix2 = 0
11723
 
    else
 
12494
    elseif (abs(uv2).gt.abs(uv1)) then
11724
12495
      ix2 = uv2
11725
12496
      gg=yd2*pq2 ;hh=xd2*uv2
11726
12497
      ix1 = gg-hh
11727
12498
      if (abs(ix1).lt.neglig(prcpar)*max(abs(gg),abs(hh))) ix1 = 0
 
12499
    else
 
12500
      ix1 = uv1
 
12501
      ix2 = uv2
11728
12502
    endif
11729
12503
    x1 = acmplx(rx1,ix1)
11730
12504
    x2 = acmplx(rx2,ix2)
12066
12840
!***********************************************************************
12067
12841
! Provides the functions
12068
12842
!   olog(x,n) = log(x) + n*pi*imag  
12069
 
!   olog2(x,n) = olog(x,n)/(x-1)
 
12843
!   olog1(x,n) = olog(x,n)/(x-1)
 
12844
!   olog2(x,n) = ( olog1(x,n) - 1 )/(x-1)
 
12845
!   olog3(x,n) = ( olog2(x,n) + 1/2 )/(x-1)
12070
12846
! In the vicinity of x=1,n=0, the logarithm of complex argument is
12071
12847
! evaluated with a series expansion.
12072
12848
!***********************************************************************
12076
12852
  use avh_olo_qp_auxfun
12077
12853
  implicit none
12078
12854
  private
12079
 
  public :: update_olog,olog,olog2
 
12855
  public :: update_olog,olog,olog1,olog2,olog3
12080
12856
 
12081
12857
  real(kindr2) &  
12082
12858
         ,allocatable,save :: thrs(:,:)
12086
12862
  interface olog
12087
12863
    module procedure log_c,log_r
12088
12864
  end interface
 
12865
  interface olog1
 
12866
    module procedure log1_c,log1_r
 
12867
  end interface
12089
12868
  interface olog2
12090
12869
    module procedure log2_c,log2_r
12091
12870
  end interface
 
12871
  interface olog3
 
12872
    module procedure log3_c,log3_r
 
12873
  end interface
12092
12874
 
12093
12875
contains
12094
12876
 
12199
12981
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
12200
12982
                                else ;nn=ntrm(1,prcpar)
12201
12983
  endif
 
12984
! convergence acceleration using  z=(y-1)/(y+1)
 
12985
! rslt = 2 * ( z + z^3/3 + z^5/5 + z^7/7 + ... )  
12202
12986
  zz = zz/(yy+1)
12203
12987
  z2 = zz*zz
12204
12988
  aa = 2
12236
13020
  end function
12237
13021
 
12238
13022
 
12239
 
  function log2_c(xx,iph) result(rslt)
 
13023
  function log1_c(xx,iph) result(rslt)
12240
13024
!***********************************************************************
12241
13025
!***********************************************************************
12242
13026
  complex(kindr2) &   
12253
13037
!
12254
13038
  if (abs(imx).le.EPSN*abs(rex)) then
12255
13039
    if (rex.ge.RZRO) then
12256
 
      rslt = log2_r( rex, iph )
 
13040
      rslt = log1_r( rex, iph )
12257
13041
    else
12258
 
      rslt = log2_r(-rex, iph+sgnRe(imx) )
 
13042
      rslt = log1_r(-rex, iph+sgnRe(imx) )
12259
13043
    endif
12260
13044
    return
12261
13045
  endif
12281
13065
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
12282
13066
                                else ;nn=ntrm(1,prcpar)
12283
13067
  endif
 
13068
! convergence acceleration using  z=(y-1)/(y+1)
 
13069
! rslt = 2/(y+1) * ( 1 + z^2/3 + z^4/5 + z^6/7 + ... )  
12284
13070
  zz = zz/(yy+1)
12285
13071
  z2 = zz*zz
12286
13072
  aa = 2
12293
13079
  end function
12294
13080
 
12295
13081
 
12296
 
  function log2_r(xx,iph) result(rslt)
 
13082
  function log1_r(xx,iph) result(rslt)
12297
13083
!***********************************************************************
12298
13084
!***********************************************************************
12299
13085
  real(kindr2) &  
12309
13095
!  integer :: nn,ii
12310
13096
!
12311
13097
  if (xx.eq.RZRO) then
12312
 
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_r: ' &
 
13098
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log1_r: ' &
12313
13099
       ,'xx =',trim(myprint(xx)),', returning 0'
12314
13100
    rslt = 0
12315
13101
    return
12321
13107
!
12322
13108
  if (abs(yy-1).le.10*EPSN) then
12323
13109
    if (jj.ne.0) then
12324
 
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_r: ' &
 
13110
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log1_r: ' &
12325
13111
        ,'rr,jj =',trim(myprint(rr)),jj,', putting jj to 0'
12326
13112
    endif
12327
13113
    rslt = 1 - (yy-1)/2
12331
13117
  rslt = ( log(rr) + IPI*jj )/(yy-1)
12332
13118
  end function
12333
13119
 
 
13120
 
 
13121
  function log2_r(xx,iph) result(rslt)
 
13122
!***********************************************************************
 
13123
!***********************************************************************
 
13124
  real(kindr2) &  
 
13125
          ,intent(in) :: xx
 
13126
  integer ,intent(in) :: iph
 
13127
  complex(kindr2) &   
 
13128
    :: rslt
 
13129
  rslt = log2_c(xx*CONE,iph)
 
13130
  end function
 
13131
 
 
13132
 
 
13133
  function log2_c(xx,iph) result(rslt)
 
13134
!***********************************************************************
 
13135
!***********************************************************************
 
13136
  complex(kindr2) &   
 
13137
          ,intent(in) :: xx
 
13138
  integer ,intent(in) :: iph
 
13139
  complex(kindr2) &   
 
13140
    :: rslt ,yy,zz,z2
 
13141
  real(kindr2) &  
 
13142
    :: aa,rex,imx
 
13143
  integer :: nn,ii,jj
 
13144
!
 
13145
  rex = areal(xx)
 
13146
  imx = aimag(xx)
 
13147
!
 
13148
  if (rex.eq.RZRO.and.imx.eq.RZRO) then
 
13149
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_c: ' &
 
13150
       ,'xx = 0, returning 0'
 
13151
    rslt = 0
 
13152
    return
 
13153
  endif
 
13154
!
 
13155
  if (mod(iph,2).eq.0) then ;yy= xx ;jj=iph
 
13156
                       else ;yy=-xx ;jj=iph+sgnRe(imx)
 
13157
  endif
 
13158
!
 
13159
  if (jj.ne.0) then
 
13160
    rslt = ( olog1(yy,jj) - 1 )/(yy-1)
 
13161
    return
 
13162
  endif
 
13163
!
 
13164
  zz = yy-1
 
13165
  aa = abs(zz)
 
13166
  if     (aa.ge.thrs(6,prcpar)) then
 
13167
    rslt = (log(yy)/zz-1)/zz
 
13168
    return
 
13169
  elseif (aa.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
13170
  elseif (aa.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
13171
  elseif (aa.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
13172
  elseif (aa.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
13173
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
13174
                                else ;nn=ntrm(1,prcpar)
 
13175
  endif
 
13176
! convergence acceleration using  z=(y-1)/(y+1)
 
13177
! rslt = -1/(y+1) + 2/(y+1)^2 * ( z/3 + z^3/5 + z^5/7 + ... )  
 
13178
  zz = zz/(yy+1)
 
13179
  z2 = zz*zz
 
13180
  aa = 2
 
13181
  nn = 2*nn-1
 
13182
  rslt = aa/nn
 
13183
  do ii=nn-2,3,-2
 
13184
    rslt = aa/ii + z2*rslt
 
13185
  enddo
 
13186
  rslt = ( -1 + zz*rslt/(yy+1) )/(yy+1)
 
13187
  end function
 
13188
 
 
13189
 
 
13190
  function log3_r(xx,iph) result(rslt)
 
13191
!***********************************************************************
 
13192
!***********************************************************************
 
13193
  real(kindr2) &  
 
13194
          ,intent(in) :: xx
 
13195
  integer ,intent(in) :: iph
 
13196
  complex(kindr2) &   
 
13197
    :: rslt
 
13198
  rslt = log3_c(xx*CONE,iph)
 
13199
  end function
 
13200
 
 
13201
 
 
13202
  function log3_c(xx,iph) result(rslt)
 
13203
!***********************************************************************
 
13204
!***********************************************************************
 
13205
  complex(kindr2) &   
 
13206
          ,intent(in) :: xx
 
13207
  integer ,intent(in) :: iph
 
13208
  complex(kindr2) &   
 
13209
    :: rslt ,yy,zz,z2,HLF
 
13210
  real(kindr2) &  
 
13211
    :: aa,rex,imx
 
13212
  integer :: nn,ii,jj
 
13213
!
 
13214
  rex = areal(xx)
 
13215
  imx = aimag(xx)
 
13216
!
 
13217
  if (rex.eq.RZRO.and.imx.eq.RZRO) then
 
13218
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_c: ' &
 
13219
       ,'xx = 0, returning 0'
 
13220
    rslt = 0
 
13221
    return
 
13222
  endif
 
13223
!
 
13224
  HLF = CONE/2
 
13225
!
 
13226
  if (mod(iph,2).eq.0) then ;yy= xx ;jj=iph
 
13227
                       else ;yy=-xx ;jj=iph+sgnRe(imx)
 
13228
  endif
 
13229
!
 
13230
  if (jj.ne.0) then
 
13231
    rslt = ( olog2(xx,jj) + HLF )/(yy-1)
 
13232
    return
 
13233
  endif
 
13234
!
 
13235
  zz = yy-1
 
13236
  aa = abs(zz)
 
13237
  if     (aa.ge.thrs(6,prcpar)) then
 
13238
    rslt = ((log(yy)/zz-1)/zz+HLF)/zz
 
13239
    return
 
13240
  elseif (aa.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
13241
  elseif (aa.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
13242
  elseif (aa.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
13243
  elseif (aa.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
13244
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
13245
                                else ;nn=ntrm(1,prcpar)
 
13246
  endif
 
13247
! convergence acceleration using  z=(y-1)/(y+1)
 
13248
! rslt = 1/(2*(y+1)) + 2/(y+1)^3 * ( 1/3 + z^2/5 + z^4/7 + ... )
 
13249
  zz = zz/(yy+1)
 
13250
  z2 = zz*zz
 
13251
  aa = 2
 
13252
  nn = 2*nn-1
 
13253
  rslt = aa/nn
 
13254
  do ii=nn-2,3,-2
 
13255
    rslt = aa/ii + z2*rslt
 
13256
  enddo
 
13257
  rslt = ( HLF + rslt/(yy+1)**2 )/(yy+1)
 
13258
  end function
 
13259
 
12334
13260
end module
12335
13261
 
12336
13262
 
 
13263
 
 
13264
 
12337
13265
module avh_olo_qp_dilog
12338
13266
!***********************************************************************
12339
13267
!                     /1    ln(1-zz*t)
12588
13516
! 
12589
13517
  if (yy.eq.RONE.and.odd.eq.0) then
12590
13518
    if (ntwo.ne.0) then
12591
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog_r: ' &
 
13519
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog_r: ' &
12592
13520
        ,'|x|,iph = ',trim(myprint(yy)),',',jj,', returning 0'
12593
13521
    endif
12594
13522
    rslt = 0
12696
13624
!
12697
13625
  if (j1.ne.j2) then
12698
13626
    if (r1.eq.r2) then
12699
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
13627
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
12700
13628
        ,'j1,j2,r1-r2',j1,j2,',',trim(myprint(r1-r2)),', returning 0'
12701
13629
      rslt = 0
12702
13630
!      write(*,*) 'dilog2_c j1=/=j2,r1=r2' !DEBUG
12710
13638
!
12711
13639
  if (a1.lt.eps) then
12712
13640
    if (a2.lt.eps) then
12713
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
13641
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
12714
13642
        ,'r1,r2 =',trim(myprint(r1)),',',trim(myprint(r2)),', returning 0'
12715
13643
      rslt = 0
12716
13644
!      write(*,*) 'dilog2_c r1<eps,r2<eps' !DEBUG
12732
13660
!      write(*,*) 'dilog2_c ||1-y1|/|1-y2|-1|>0.1' !DEBUG
12733
13661
      return
12734
13662
    elseif (oo.eq.0.and.ao1.lt.eps) then
12735
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
13663
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
12736
13664
        ,'r1,oo,nn =',trim(myprint(r1)),',',oo,nn,', putting nn=0'
12737
13665
      if (ao2.lt.eps) then
12738
13666
        rslt = -1
12742
13670
        y1=1-eps ;nn=0 ;logr1=0 ;r1=1-eps
12743
13671
      endif
12744
13672
    elseif (oo.eq.0.and.ao2.lt.eps) then
12745
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
13673
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
12746
13674
        ,'r2,oo,nn =',trim(myprint(r2)),',',oo,nn,', putting nn=0'
12747
13675
      y2=1-eps ;nn=0 ;logr2=0 ;r2=1-eps
12748
13676
    endif
12757
13685
      if (a1.gt.RONE) ii = ii + (nn+pp(oo,sgnIm(y2)))
12758
13686
      if (a2.gt.RONE) ii = ii - (nn+pp(oo,sgnIm(y2)))
12759
13687
      ii = nn*ii
12760
 
      if (ii.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
13688
      if (ii.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
12761
13689
        ,'r1,r2,nn =',trim(myprint(r1)),',',trim(myprint(r2)),',',nn &
12762
13690
        ,', putting nn=0'
12763
 
      rslt = -olog2(y2,0)
 
13691
      rslt = -olog1(y2,0)
12764
13692
!      write(*,*) 'dilog2_c |logr1/lorg2|<eps' !DEBUG
12765
13693
      return
12766
13694
    endif
12772
13700
    nn=-nn ;oo=-oo
12773
13701
  endif
12774
13702
!
12775
 
  ff=y1/y2         ;ff=-olog2(ff,0)/y2
12776
 
  gg=(1-y1)/(1-y2) ;gg=-olog2(gg,0)/(1-y2)
 
13703
  ff=y1/y2         ;ff=-olog1(ff,0)/y2
 
13704
  gg=(1-y1)/(1-y2) ;gg=-olog1(gg,0)/(1-y2)
12777
13705
!
12778
13706
  if (2*areal(y1).ge.RONE) then
12779
13707
!    write(*,*) 'dilog2_c re>1/2' !DEBUG
12830
13758
!
12831
13759
  if (j1.ne.j2) then
12832
13760
    if (r1.eq.r2) then
12833
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
13761
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
12834
13762
        ,'j1,j2,r1-r2',j1,j2,',',trim(myprint(r1-r2)),', returning 0'
12835
13763
      rslt = 0
12836
13764
!      write(*,*) 'dilog2_r j1=/=j2,r1=r2' !DEBUG
12844
13772
!
12845
13773
  if (r1.lt.eps) then
12846
13774
    if (r2.lt.eps) then
12847
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
13775
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
12848
13776
        ,'r1,r2 =',trim(myprint(r1)),',',trim(myprint(r2)),', returning 0'
12849
13777
      rslt = 0
12850
13778
!      write(*,*) 'dilog2_r r1<eps,r2<eps' !DEBUG
12866
13794
!      write(*,*) 'dilog2_r ||1-y1|/|1-y2|-1|>0.1' !DEBUG
12867
13795
      return
12868
13796
    elseif (oo.eq.0.and.ro1.lt.eps) then
12869
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
13797
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
12870
13798
        ,'r1,oo,nn =',trim(myprint(r1)),',',oo,nn,', putting nn=0'
12871
13799
      if (ro2.lt.eps) then
12872
13800
        rslt = -1
12876
13804
        y1=1-eps ;nn=0 ;logr1=0 ;r1=1-eps
12877
13805
      endif
12878
13806
    elseif (oo.eq.0.and.ro2.lt.eps) then
12879
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
13807
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
12880
13808
        ,'r2,oo,nn =',trim(myprint(r2)),',',oo,nn,', putting nn=0'
12881
13809
      y2=1-eps ;nn=0 ;logr2=0 ;r2=1-eps
12882
13810
    endif
12891
13819
      if (r1.gt.RONE) ii = ii + (nn+2*oo)
12892
13820
      if (r2.gt.RONE) ii = ii - (nn+2*oo)
12893
13821
      ii = nn*ii
12894
 
      if (ii.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
13822
      if (ii.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
12895
13823
        ,'r1,r2,nn =',trim(myprint(r1)),',',trim(myprint(r2)),',',nn &
12896
13824
        ,', putting nn=0'
12897
 
      rslt = -olog2(y2,0)
 
13825
      rslt = -olog1(y2,0)
12898
13826
!      write(*,*) 'dilog2_r |logr1/lorg2|<eps' !DEBUG
12899
13827
      return
12900
13828
    endif
12906
13834
    nn=-nn ;oo=-oo
12907
13835
  endif
12908
13836
!
12909
 
  ff=y1/y2         ;ff=-olog2(ff,0)/y2
12910
 
  gg=(1-y1)/(1-y2) ;gg=-olog2(gg,0)/(1-y2)
 
13837
  ff=y1/y2         ;ff=-olog1(ff,0)/y2
 
13838
  gg=(1-y1)/(1-y2) ;gg=-olog1(gg,0)/(1-y2)
12911
13839
!
12912
13840
  if (2*y1.ge.RONE) then
12913
13841
!    write(*,*) 'dilog2_r re>1/2' !DEBUG
13337
14265
 
13338
14266
  implicit none
13339
14267
  private
13340
 
  public :: qmplx_type,qonv,directly,sheet,logc,logc2,li2c,li2c2
 
14268
  public :: qmplx_type,qonv,directly,sheet,logc,logc2,logc3,li2c,li2c2
13341
14269
  public :: operator (*) ,operator (/)
13342
14270
 
13343
14271
  type :: qmplx_type
13584
14512
  type(qmplx_type) ,intent(in) :: xx
13585
14513
  complex(kindr2) &   
13586
14514
    :: rslt
13587
 
!  rslt = -olog2(acmplx(xx%c),xx%p)
13588
 
  rslt = -olog2(xx%c,xx%p)
 
14515
!  rslt = -olog1(acmplx(xx%c),xx%p)
 
14516
  rslt = -olog1(xx%c,xx%p)
 
14517
  end function
 
14518
 
 
14519
  function logc3(xx) result(rslt)
 
14520
!*******************************************************************
 
14521
!  ( log(xx)/(1-xx) + 1 )/(1-xx)
 
14522
!*******************************************************************
 
14523
  type(qmplx_type) ,intent(in) :: xx
 
14524
  complex(kindr2) &   
 
14525
    :: rslt
 
14526
!  rslt = olog2(acmplx(xx%c),xx%p)
 
14527
  rslt = olog2(xx%c,xx%p)
13589
14528
  end function
13590
14529
 
13591
14530
  function li2c(xx) result(rslt)
13625
14564
  use avh_olo_qp_auxfun
13626
14565
  use avh_olo_qp_bnlog
13627
14566
  use avh_olo_qp_qmplx
 
14567
  use avh_olo_qp_olog
13628
14568
  implicit none
13629
14569
  private
13630
 
  public :: tadp ,tadpn ,bub0 ,bub1 ,bub11 ,bub111 ,bub1111
 
14570
  public :: tadp ,tadpn ,bub0 ,dbub0 ,bub1 ,bub11 ,bub111 ,bub1111
13631
14571
 
13632
14572
contains
13633
14573
 
14218
15158
  b0011(0) = ( a0(0,0) - x1*b111(0) + x2*b11(0) + 4*b0011(1) )/10
14219
15159
  end subroutine
14220
15160
 
 
15161
 
 
15162
!*******************************************************************
 
15163
! Derivative of B0
 
15164
! expects  m0<m1
 
15165
! only finite case, so input must not be  m0=0 & m1=pp
 
15166
!*******************************************************************
 
15167
 
 
15168
  subroutine dbub0( rslt &
 
15169
                   ,pp,m0,m1 ,app,am0,am1 )
 
15170
  complex(kindr2) &   
 
15171
    ,intent(out) :: rslt
 
15172
  complex(kindr2) &   
 
15173
    ,intent(in)  :: pp,m0,m1
 
15174
  real(kindr2) &  
 
15175
    ,intent(in)  :: app,am0,am1
 
15176
  complex(kindr2) &   
 
15177
    :: ch,x1,x2,lambda
 
15178
  real(kindr2) &  
 
15179
    :: ax1,ax2,ax1x2,maxa
 
15180
  type(qmplx_type) :: q1,q2,q1o,q2o
 
15181
  integer :: sgn
 
15182
!
 
15183
  if (am1.eq.RZRO) then
 
15184
    if (app.eq.RZRO) then
 
15185
      rslt = 0
 
15186
      return
 
15187
    endif
 
15188
  endif
 
15189
!
 
15190
  if (app.eq.RZRO) then
 
15191
    if (abs(m0-m1).le.am1*EPSN*10) then
 
15192
      rslt = 1/(6*m1)
 
15193
    else
 
15194
      ch = m0/m1
 
15195
      rslt = ( CONE/2 - ch*olog3(ch,0) )/m1 
 
15196
    endif
 
15197
  elseif (am1.eq.RZRO) then
 
15198
    rslt =-1/pp
 
15199
  else
 
15200
    call solabc( x1,x2 ,lambda ,pp ,(m0-m1)-pp ,m1 ,0 )
 
15201
    sgn =-sgnRe(pp)*sgnRe(x2-x1)
 
15202
    q1  = qonv(x1  , sgn)
 
15203
    q1o = qonv(x1-1, sgn)
 
15204
    q2  = qonv(x2  ,-sgn)
 
15205
    q2o = qonv(x2-1,-sgn)
 
15206
    ax1 = abs(x1)
 
15207
    ax2 = abs(x2)
 
15208
    ax1x2 = abs(x1-x2)
 
15209
    maxa = max(ax1,ax2)
 
15210
    if (ax1x2.lt.maxa*EPSN*10) then
 
15211
      rslt = ( (x1+x2-1)*logc(q2/q2o) - 2 )/pp
 
15212
    elseif (ax1x2*2.lt.maxa) then
 
15213
      if     (x1.eq.CZRO.or.x1.eq.CONE) then
 
15214
        rslt = ( (x1+x2-1)*logc(q2/q2o) - 1 )/pp
 
15215
      elseif (x2.eq.CZRO.or.x2.eq.CONE) then
 
15216
        rslt = ( (x1+x2-1)*logc(q1/q1o) - 1 )/pp
 
15217
      else
 
15218
        rslt = x1*(x1-1)*( logc2(q1o/q2o)/(x2-1) - logc2(q1/q2)/x2 ) &
 
15219
             + (x1+x2-1)*logc(q2/q2o) - 1
 
15220
        rslt = rslt/pp
 
15221
      endif
 
15222
    else
 
15223
      rslt = 0
 
15224
      if (ax1.ne.RZRO) then
 
15225
        if (ax1.lt.2*RONE) then
 
15226
          rslt = rslt - x1
 
15227
          if (x1.ne.CONE) rslt = rslt - x1*logc2(q1/q1o)
 
15228
        else
 
15229
          rslt = rslt + x1/(x1-1)*logc3(q1/q1o)
 
15230
        endif
 
15231
      endif
 
15232
      if (ax2.ne.RZRO) then
 
15233
        if (ax2.lt.2*RONE) then
 
15234
          rslt = rslt + x2
 
15235
          if (x2.ne.CONE) rslt = rslt + x2*logc2(q2/q2o)
 
15236
        else
 
15237
          rslt = rslt - x2/(x2-1)*logc3(q2/q2o)
 
15238
        endif
 
15239
      endif
 
15240
      rslt = rslt/lambda
 
15241
    endif
 
15242
  endif
 
15243
!
 
15244
  end subroutine
 
15245
 
 
15246
 
14221
15247
end module
14222
15248
 
14223
15249
 
14451
15477
    log2 = olog( abs(rp2/rmu2) ,i2 )
14452
15478
    log3 = olog( abs(rp3/rmu2) ,i3 )
14453
15479
    rslt(2) = 0
14454
 
    rslt(1) = -olog2( abs(rp3/rp2) ,i3-i2 )/rp2
 
15480
    rslt(1) = -olog1( abs(rp3/rp2) ,i3-i2 )/rp2
14455
15481
    rslt(0) = -rslt(1)*(log3+log2)/2
14456
15482
  elseif (icase.eq.3) then
14457
15483
! 3 masses non-zero
17024
18050
     ,intent(in) :: y1,y2
17025
18051
  complex(kindr2) &   
17026
18052
     :: rslt ,oy1,oy2
 
18053
!
17027
18054
   oy1 = 1-y1
17028
18055
   oy2 = 1-y2
17029
18056
   rslt = logc2( qonv(-y2)/qonv(-y1) )/y1 &
17092
18119
  private
17093
18120
  public :: olo_unit ,olo_scale ,olo_onshell ,olo_setting
17094
18121
  public :: olo_precision
17095
 
  public :: olo_a0 ,olo_b0 ,olo_b11 ,olo_c0 ,olo_d0
 
18122
  public :: olo_a0 ,olo_b0 ,olo_db0 ,olo_b11 ,olo_c0 ,olo_d0
17096
18123
  public :: olo_an ,olo_bn
17097
18124
  public :: olo
17098
18125
  public :: olo_get_scale ,olo_get_onshell ,olo_get_precision
17099
 
  public :: a0_r,a0rr,a0_c,a0cr
17100
 
  public :: an_r,anrr,an_c,ancr
17101
 
  public :: b0rr,b0rrr,b0rc,b0rcr,b0cc,b0ccr
17102
 
  public :: b11rr,b11rrr,b11rc,b11rcr,b11cc,b11ccr
17103
 
  public :: bnrr,bnrrr,bnrc,bnrcr,bncc,bnccr
17104
 
  public :: c0rr,c0rrr,c0rc,c0rcr,c0cc,c0ccr
17105
 
  public :: d0rr,d0rrr,d0rc,d0rcr,d0cc,d0ccr
17106
18126
!
17107
18127
  integer ,public ,parameter :: olo_kind=kindr2    
17108
18128
!
17128
18148
  interface olo_b0
17129
18149
    module procedure b0rr,b0rrr,b0rc,b0rcr,b0cc,b0ccr
17130
18150
  end interface 
 
18151
  interface olo_db0
 
18152
    module procedure db0rr,db0rrr,db0rc,db0rcr,db0cc,db0ccr
 
18153
  end interface 
17131
18154
  interface olo_b11
17132
18155
    module procedure b11rr,b11rrr,b11rc,b11rcr,b11cc,b11ccr
17133
18156
  end interface 
18114
19137
  endif
18115
19138
  end subroutine
18116
19139
 
 
19140
!*******************************************************************
 
19141
! Derivative of B0
 
19142
!*******************************************************************
 
19143
  subroutine db0cc( rslt ,pp,m1,m2 )
 
19144
!
 
19145
  use avh_olo_qp_bub ,only: dbub0
 
19146
  use avh_olo_qp_olog ,only: olog
 
19147
!
 
19148
  complex(kindr2) &   
 
19149
    ,intent(out) :: rslt(0:2)
 
19150
  complex(kindr2) &   
 
19151
    ,intent(in)  :: pp
 
19152
  complex(kindr2) &   
 
19153
    ,intent(in)  :: m1,m2
 
19154
!
 
19155
  complex(kindr2) &   
 
19156
    :: ss,r1,r2,ch
 
19157
  real(kindr2) &  
 
19158
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
19159
  character(26+99) ,parameter :: warning=&
 
19160
                     'WARNING from OneLOop db0: '//warnonshell
 
19161
  if (initz) call init
 
19162
  ss = pp
 
19163
  r1 = m1
 
19164
  r2 = m2
 
19165
!
 
19166
  app = areal(ss)
 
19167
  if (aimag(ss).ne.RZRO) then
 
19168
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
19169
      ,'ss has non-zero imaginary part, putting it to zero.'
 
19170
    ss = acmplx( app )
 
19171
  endif
 
19172
  app = abs(app)
 
19173
!
 
19174
  am1 = areal(r1)
 
19175
  hh  = aimag(r1)
 
19176
  if (hh.gt.RZRO) then
 
19177
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
19178
      ,'r1 has positive imaginary part, switching its sign.'
 
19179
    r1 = acmplx( am1 ,-hh )
 
19180
  endif
 
19181
  am1 = abs(am1) + abs(hh)
 
19182
!
 
19183
  am2 = areal(r2)
 
19184
  hh  = aimag(r2)
 
19185
  if (hh.gt.RZRO) then
 
19186
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
19187
      ,'r2 has positive imaginary part, switching its sign.'
 
19188
    r2 = acmplx( am2 ,-hh )
 
19189
  endif
 
19190
  am2 = abs(am2) + abs(hh)
 
19191
!
 
19192
  mulocal = muscale 
 
19193
!
 
19194
  mulocal2 = mulocal*mulocal
 
19195
!
 
19196
  if (am1.gt.am2) then
 
19197
    ch=r1 ; r1=r2 ; r2=ch
 
19198
    hh=am1;am1=am2;am2=hh
 
19199
  endif
 
19200
  ssr2 = abs(ss-r2)
 
19201
!
 
19202
  if (nonzerothrs) then
 
19203
    hh = onshellthrs
 
19204
    if (app.lt.hh) app = 0
 
19205
    if (am1.lt.hh) am1 = 0
 
19206
    if (am2.lt.hh) am2 = 0
 
19207
    if (ssr2.lt.hh) ssr2 = 0
 
19208
  elseif (wunit.gt.0) then
 
19209
    hh = onshellthrs*max(app,max(am1,am2))
 
19210
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
19211
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
19212
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
19213
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
19214
  endif
 
19215
!
 
19216
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
19217
!
 
19218
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
19219
    rslt(1) =-1/(2*ss)
 
19220
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
19221
  else
 
19222
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
19223
  endif
 
19224
!
 
19225
  if (punit.gt.0) then
 
19226
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
19227
    write(punit,*) ' pp:',trim(myprint(pp))
 
19228
    write(punit,*) ' m1:',trim(myprint(m1))
 
19229
    write(punit,*) ' m2:',trim(myprint(m2))
 
19230
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
19231
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
19232
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
19233
  endif
 
19234
  end subroutine
 
19235
 
 
19236
  subroutine db0ccr( rslt ,pp,m1,m2 ,rmu )
 
19237
!
 
19238
  use avh_olo_qp_bub ,only: dbub0
 
19239
  use avh_olo_qp_olog ,only: olog
 
19240
!
 
19241
  complex(kindr2) &   
 
19242
    ,intent(out) :: rslt(0:2)
 
19243
  complex(kindr2) &   
 
19244
    ,intent(in)  :: pp
 
19245
  complex(kindr2) &   
 
19246
    ,intent(in)  :: m1,m2
 
19247
  real(kindr2) &  
 
19248
   ,intent(in)  :: rmu       
 
19249
!
 
19250
  complex(kindr2) &   
 
19251
    :: ss,r1,r2,ch
 
19252
  real(kindr2) &  
 
19253
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
19254
  character(26+99) ,parameter :: warning=&
 
19255
                     'WARNING from OneLOop db0: '//warnonshell
 
19256
  if (initz) call init
 
19257
  ss = pp
 
19258
  r1 = m1
 
19259
  r2 = m2
 
19260
!
 
19261
  app = areal(ss)
 
19262
  if (aimag(ss).ne.RZRO) then
 
19263
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
19264
      ,'ss has non-zero imaginary part, putting it to zero.'
 
19265
    ss = acmplx( app )
 
19266
  endif
 
19267
  app = abs(app)
 
19268
!
 
19269
  am1 = areal(r1)
 
19270
  hh  = aimag(r1)
 
19271
  if (hh.gt.RZRO) then
 
19272
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
19273
      ,'r1 has positive imaginary part, switching its sign.'
 
19274
    r1 = acmplx( am1 ,-hh )
 
19275
  endif
 
19276
  am1 = abs(am1) + abs(hh)
 
19277
!
 
19278
  am2 = areal(r2)
 
19279
  hh  = aimag(r2)
 
19280
  if (hh.gt.RZRO) then
 
19281
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
19282
      ,'r2 has positive imaginary part, switching its sign.'
 
19283
    r2 = acmplx( am2 ,-hh )
 
19284
  endif
 
19285
  am2 = abs(am2) + abs(hh)
 
19286
!
 
19287
  mulocal = rmu     
 
19288
!
 
19289
  mulocal2 = mulocal*mulocal
 
19290
!
 
19291
  if (am1.gt.am2) then
 
19292
    ch=r1 ; r1=r2 ; r2=ch
 
19293
    hh=am1;am1=am2;am2=hh
 
19294
  endif
 
19295
  ssr2 = abs(ss-r2)
 
19296
!
 
19297
  if (nonzerothrs) then
 
19298
    hh = onshellthrs
 
19299
    if (app.lt.hh) app = 0
 
19300
    if (am1.lt.hh) am1 = 0
 
19301
    if (am2.lt.hh) am2 = 0
 
19302
    if (ssr2.lt.hh) ssr2 = 0
 
19303
  elseif (wunit.gt.0) then
 
19304
    hh = onshellthrs*max(app,max(am1,am2))
 
19305
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
19306
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
19307
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
19308
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
19309
  endif
 
19310
!
 
19311
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
19312
!
 
19313
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
19314
    rslt(1) =-1/(2*ss)
 
19315
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
19316
  else
 
19317
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
19318
  endif
 
19319
!
 
19320
  if (punit.gt.0) then
 
19321
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
19322
    write(punit,*) ' pp:',trim(myprint(pp))
 
19323
    write(punit,*) ' m1:',trim(myprint(m1))
 
19324
    write(punit,*) ' m2:',trim(myprint(m2))
 
19325
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
19326
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
19327
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
19328
  endif
 
19329
  end subroutine
 
19330
 
 
19331
  subroutine db0rc( rslt ,pp ,m1,m2 )
 
19332
!
 
19333
  use avh_olo_qp_bub ,only: dbub0
 
19334
  use avh_olo_qp_olog ,only: olog
 
19335
!
 
19336
  complex(kindr2) &   
 
19337
    ,intent(out) :: rslt(0:2)
 
19338
  real(kindr2) &  
 
19339
    ,intent(in)  :: pp
 
19340
  complex(kindr2) &   
 
19341
    ,intent(in)  :: m1,m2
 
19342
!
 
19343
  complex(kindr2) &   
 
19344
    :: ss,r1,r2,ch
 
19345
  real(kindr2) &  
 
19346
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
19347
  character(26+99) ,parameter :: warning=&
 
19348
                     'WARNING from OneLOop db0: '//warnonshell
 
19349
  if (initz) call init
 
19350
  ss = pp
 
19351
  r1 = m1
 
19352
  r2 = m2
 
19353
!
 
19354
  app = abs(pp)
 
19355
!
 
19356
  am1 = areal(r1)
 
19357
  hh  = aimag(r1)
 
19358
  if (hh.gt.RZRO) then
 
19359
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
19360
      ,'r1 has positive imaginary part, switching its sign.'
 
19361
    r1 = acmplx( am1 ,-hh )
 
19362
  endif
 
19363
  am1 = abs(am1) + abs(hh)
 
19364
!
 
19365
  am2 = areal(r2)
 
19366
  hh  = aimag(r2)
 
19367
  if (hh.gt.RZRO) then
 
19368
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
19369
      ,'r2 has positive imaginary part, switching its sign.'
 
19370
    r2 = acmplx( am2 ,-hh )
 
19371
  endif
 
19372
  am2 = abs(am2) + abs(hh)
 
19373
!
 
19374
  mulocal = muscale 
 
19375
!
 
19376
  mulocal2 = mulocal*mulocal
 
19377
!
 
19378
  if (am1.gt.am2) then
 
19379
    ch=r1 ; r1=r2 ; r2=ch
 
19380
    hh=am1;am1=am2;am2=hh
 
19381
  endif
 
19382
  ssr2 = abs(ss-r2)
 
19383
!
 
19384
  if (nonzerothrs) then
 
19385
    hh = onshellthrs
 
19386
    if (app.lt.hh) app = 0
 
19387
    if (am1.lt.hh) am1 = 0
 
19388
    if (am2.lt.hh) am2 = 0
 
19389
    if (ssr2.lt.hh) ssr2 = 0
 
19390
  elseif (wunit.gt.0) then
 
19391
    hh = onshellthrs*max(app,max(am1,am2))
 
19392
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
19393
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
19394
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
19395
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
19396
  endif
 
19397
!
 
19398
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
19399
!
 
19400
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
19401
    rslt(1) =-1/(2*ss)
 
19402
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
19403
  else
 
19404
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
19405
  endif
 
19406
!
 
19407
  if (punit.gt.0) then
 
19408
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
19409
    write(punit,*) ' pp:',trim(myprint(pp))
 
19410
    write(punit,*) ' m1:',trim(myprint(m1))
 
19411
    write(punit,*) ' m2:',trim(myprint(m2))
 
19412
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
19413
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
19414
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
19415
  endif
 
19416
  end subroutine
 
19417
 
 
19418
  subroutine db0rcr( rslt ,pp,m1,m2 ,rmu )
 
19419
!
 
19420
  use avh_olo_qp_bub ,only: dbub0
 
19421
  use avh_olo_qp_olog ,only: olog
 
19422
!
 
19423
  complex(kindr2) &   
 
19424
    ,intent(out) :: rslt(0:2)
 
19425
  real(kindr2) &  
 
19426
    ,intent(in)  :: pp
 
19427
  complex(kindr2) &   
 
19428
    ,intent(in)  :: m1,m2
 
19429
  real(kindr2) &  
 
19430
   ,intent(in)  :: rmu       
 
19431
!
 
19432
  complex(kindr2) &   
 
19433
    :: ss,r1,r2,ch
 
19434
  real(kindr2) &  
 
19435
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
19436
  character(26+99) ,parameter :: warning=&
 
19437
                     'WARNING from OneLOop db0: '//warnonshell
 
19438
  if (initz) call init
 
19439
  ss = pp
 
19440
  r1 = m1
 
19441
  r2 = m2
 
19442
!
 
19443
  app = abs(pp)
 
19444
!
 
19445
  am1 = areal(r1)
 
19446
  hh  = aimag(r1)
 
19447
  if (hh.gt.RZRO) then
 
19448
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
19449
      ,'r1 has positive imaginary part, switching its sign.'
 
19450
    r1 = acmplx( am1 ,-hh )
 
19451
  endif
 
19452
  am1 = abs(am1) + abs(hh)
 
19453
!
 
19454
  am2 = areal(r2)
 
19455
  hh  = aimag(r2)
 
19456
  if (hh.gt.RZRO) then
 
19457
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
19458
      ,'r2 has positive imaginary part, switching its sign.'
 
19459
    r2 = acmplx( am2 ,-hh )
 
19460
  endif
 
19461
  am2 = abs(am2) + abs(hh)
 
19462
!
 
19463
  mulocal = rmu     
 
19464
!
 
19465
  mulocal2 = mulocal*mulocal
 
19466
!
 
19467
  if (am1.gt.am2) then
 
19468
    ch=r1 ; r1=r2 ; r2=ch
 
19469
    hh=am1;am1=am2;am2=hh
 
19470
  endif
 
19471
  ssr2 = abs(ss-r2)
 
19472
!
 
19473
  if (nonzerothrs) then
 
19474
    hh = onshellthrs
 
19475
    if (app.lt.hh) app = 0
 
19476
    if (am1.lt.hh) am1 = 0
 
19477
    if (am2.lt.hh) am2 = 0
 
19478
    if (ssr2.lt.hh) ssr2 = 0
 
19479
  elseif (wunit.gt.0) then
 
19480
    hh = onshellthrs*max(app,max(am1,am2))
 
19481
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
19482
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
19483
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
19484
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
19485
  endif
 
19486
!
 
19487
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
19488
!
 
19489
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
19490
    rslt(1) =-1/(2*ss)
 
19491
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
19492
  else
 
19493
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
19494
  endif
 
19495
!
 
19496
  if (punit.gt.0) then
 
19497
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
19498
    write(punit,*) ' pp:',trim(myprint(pp))
 
19499
    write(punit,*) ' m1:',trim(myprint(m1))
 
19500
    write(punit,*) ' m2:',trim(myprint(m2))
 
19501
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
19502
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
19503
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
19504
  endif
 
19505
  end subroutine
 
19506
 
 
19507
  subroutine db0rr( rslt ,pp ,m1,m2 )
 
19508
!
 
19509
  use avh_olo_qp_bub ,only: dbub0
 
19510
  use avh_olo_qp_olog ,only: olog
 
19511
!
 
19512
  complex(kindr2) &   
 
19513
    ,intent(out) :: rslt(0:2)
 
19514
  real(kindr2) &  
 
19515
    ,intent(in)  :: pp
 
19516
  real(kindr2) &  
 
19517
    ,intent(in)  :: m1,m2
 
19518
!
 
19519
  complex(kindr2) &   
 
19520
    :: ss,r1,r2,ch
 
19521
  real(kindr2) &  
 
19522
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
19523
  character(26+99) ,parameter :: warning=&
 
19524
                     'WARNING from OneLOop db0: '//warnonshell
 
19525
  if (initz) call init
 
19526
  ss = pp
 
19527
  r1 = m1
 
19528
  r2 = m2
 
19529
!
 
19530
  app = abs(pp)
 
19531
!
 
19532
  am1 = abs(m1)
 
19533
  am2 = abs(m2)
 
19534
!
 
19535
  mulocal = muscale 
 
19536
!
 
19537
  mulocal2 = mulocal*mulocal
 
19538
!
 
19539
  if (am1.gt.am2) then
 
19540
    ch=r1 ; r1=r2 ; r2=ch
 
19541
    hh=am1;am1=am2;am2=hh
 
19542
  endif
 
19543
  ssr2 = abs(ss-r2)
 
19544
!
 
19545
  if (nonzerothrs) then
 
19546
    hh = onshellthrs
 
19547
    if (app.lt.hh) app = 0
 
19548
    if (am1.lt.hh) am1 = 0
 
19549
    if (am2.lt.hh) am2 = 0
 
19550
    if (ssr2.lt.hh) ssr2 = 0
 
19551
  elseif (wunit.gt.0) then
 
19552
    hh = onshellthrs*max(app,max(am1,am2))
 
19553
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
19554
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
19555
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
19556
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
19557
  endif
 
19558
!
 
19559
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
19560
!
 
19561
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
19562
    rslt(1) =-1/(2*ss)
 
19563
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
19564
  else
 
19565
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
19566
  endif
 
19567
!
 
19568
  if (punit.gt.0) then
 
19569
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
19570
    write(punit,*) ' pp:',trim(myprint(pp))
 
19571
    write(punit,*) ' m1:',trim(myprint(m1))
 
19572
    write(punit,*) ' m2:',trim(myprint(m2))
 
19573
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
19574
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
19575
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
19576
  endif
 
19577
  end subroutine
 
19578
 
 
19579
  subroutine db0rrr( rslt ,pp ,m1,m2 ,rmu )
 
19580
!
 
19581
  use avh_olo_qp_bub ,only: dbub0
 
19582
  use avh_olo_qp_olog ,only: olog
 
19583
!
 
19584
  complex(kindr2) &   
 
19585
    ,intent(out) :: rslt(0:2)
 
19586
  real(kindr2) &  
 
19587
    ,intent(in)  :: pp
 
19588
  real(kindr2) &  
 
19589
    ,intent(in)  :: m1,m2
 
19590
  real(kindr2) &  
 
19591
   ,intent(in)  :: rmu       
 
19592
!
 
19593
  complex(kindr2) &   
 
19594
    :: ss,r1,r2,ch
 
19595
  real(kindr2) &  
 
19596
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
19597
  character(26+99) ,parameter :: warning=&
 
19598
                     'WARNING from OneLOop db0: '//warnonshell
 
19599
  if (initz) call init
 
19600
  ss = pp
 
19601
  r1 = m1
 
19602
  r2 = m2
 
19603
!
 
19604
  app = abs(pp)
 
19605
!
 
19606
  am1 = abs(m1)
 
19607
  am2 = abs(m2)
 
19608
!
 
19609
  mulocal = rmu     
 
19610
!
 
19611
  mulocal2 = mulocal*mulocal
 
19612
!
 
19613
  if (am1.gt.am2) then
 
19614
    ch=r1 ; r1=r2 ; r2=ch
 
19615
    hh=am1;am1=am2;am2=hh
 
19616
  endif
 
19617
  ssr2 = abs(ss-r2)
 
19618
!
 
19619
  if (nonzerothrs) then
 
19620
    hh = onshellthrs
 
19621
    if (app.lt.hh) app = 0
 
19622
    if (am1.lt.hh) am1 = 0
 
19623
    if (am2.lt.hh) am2 = 0
 
19624
    if (ssr2.lt.hh) ssr2 = 0
 
19625
  elseif (wunit.gt.0) then
 
19626
    hh = onshellthrs*max(app,max(am1,am2))
 
19627
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
19628
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
19629
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
19630
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
19631
  endif
 
19632
!
 
19633
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
19634
!
 
19635
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
19636
    rslt(1) =-1/(2*ss)
 
19637
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
19638
  else
 
19639
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
19640
  endif
 
19641
!
 
19642
  if (punit.gt.0) then
 
19643
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
19644
    write(punit,*) ' pp:',trim(myprint(pp))
 
19645
    write(punit,*) ' m1:',trim(myprint(m1))
 
19646
    write(punit,*) ' m2:',trim(myprint(m2))
 
19647
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
19648
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
19649
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
19650
  endif
 
19651
  end subroutine
 
19652
 
 
19653
 
18117
19654
 
18118
19655
!*******************************************************************
18119
19656
! Return the Papparino-Veltman functions b11,b00,b1,b0 , for
20485
22022
               .or.(     areal(ss(1)).ge.-small  &
20486
22023
                    .and.areal(ss(2)).ge.-small  &
20487
22024
                    .and.areal(ss(3)).ge.-small  &
20488
 
                    .and.areal(ss(4)).ge.-small) )
 
22025
                    .and.areal(ss(4)).ge.-small) &
 
22026
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
20489
22027
    if (useboxc) then
20490
22028
      call boxc( rslt ,ss,rr ,as ,smax )
20491
22029
    else
20501
22039
                 .or.(     areal(ss(1)).ge.-small  &
20502
22040
                      .and.areal(ss(2)).ge.-small  &
20503
22041
                      .and.areal(ss(3)).ge.-small  &
20504
 
!OLD                      .and.areal(ss(4)).ge.-small) )
20505
 
                      .and.areal(ss(4)).ge.-small) & !NEW
20506
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
22042
                      .and.areal(ss(4)).ge.-small) &
 
22043
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
20507
22044
      if (useboxc) then
20508
22045
        call boxc( rslt ,ss,rr ,as ,smax )
20509
22046
      else
20760
22297
               .or.(     areal(ss(1)).ge.-small  &
20761
22298
                    .and.areal(ss(2)).ge.-small  &
20762
22299
                    .and.areal(ss(3)).ge.-small  &
20763
 
                    .and.areal(ss(4)).ge.-small) )
 
22300
                    .and.areal(ss(4)).ge.-small) &
 
22301
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
20764
22302
    if (useboxc) then
20765
22303
      call boxc( rslt ,ss,rr ,as ,smax )
20766
22304
    else
20776
22314
                 .or.(     areal(ss(1)).ge.-small  &
20777
22315
                      .and.areal(ss(2)).ge.-small  &
20778
22316
                      .and.areal(ss(3)).ge.-small  &
20779
 
!OLD                      .and.areal(ss(4)).ge.-small) )
20780
 
                      .and.areal(ss(4)).ge.-small) & !NEW
20781
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
22317
                      .and.areal(ss(4)).ge.-small) &
 
22318
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
20782
22319
      if (useboxc) then
20783
22320
        call boxc( rslt ,ss,rr ,as ,smax )
20784
22321
      else
21027
22564
               .or.(     areal(ss(1)).ge.-small  &
21028
22565
                    .and.areal(ss(2)).ge.-small  &
21029
22566
                    .and.areal(ss(3)).ge.-small  &
21030
 
                    .and.areal(ss(4)).ge.-small) )
 
22567
                    .and.areal(ss(4)).ge.-small) &
 
22568
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
21031
22569
    if (useboxc) then
21032
22570
      call boxc( rslt ,ss,rr ,as ,smax )
21033
22571
    else
21043
22581
                 .or.(     areal(ss(1)).ge.-small  &
21044
22582
                      .and.areal(ss(2)).ge.-small  &
21045
22583
                      .and.areal(ss(3)).ge.-small  &
21046
 
!OLD                      .and.areal(ss(4)).ge.-small) )
21047
 
                      .and.areal(ss(4)).ge.-small) & !NEW
21048
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
22584
                      .and.areal(ss(4)).ge.-small) &
 
22585
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
21049
22586
      if (useboxc) then
21050
22587
        call boxc( rslt ,ss,rr ,as ,smax )
21051
22588
      else
21296
22833
               .or.(     areal(ss(1)).ge.-small  &
21297
22834
                    .and.areal(ss(2)).ge.-small  &
21298
22835
                    .and.areal(ss(3)).ge.-small  &
21299
 
                    .and.areal(ss(4)).ge.-small) )
 
22836
                    .and.areal(ss(4)).ge.-small) &
 
22837
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
21300
22838
    if (useboxc) then
21301
22839
      call boxc( rslt ,ss,rr ,as ,smax )
21302
22840
    else
21312
22850
                 .or.(     areal(ss(1)).ge.-small  &
21313
22851
                      .and.areal(ss(2)).ge.-small  &
21314
22852
                      .and.areal(ss(3)).ge.-small  &
21315
 
!OLD                      .and.areal(ss(4)).ge.-small) )
21316
 
                      .and.areal(ss(4)).ge.-small) & !NEW
21317
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
22853
                      .and.areal(ss(4)).ge.-small) &
 
22854
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
21318
22855
      if (useboxc) then
21319
22856
        call boxc( rslt ,ss,rr ,as ,smax )
21320
22857
      else
21556
23093
               .or.(     areal(ss(1)).ge.-small  &
21557
23094
                    .and.areal(ss(2)).ge.-small  &
21558
23095
                    .and.areal(ss(3)).ge.-small  &
21559
 
                    .and.areal(ss(4)).ge.-small) )
 
23096
                    .and.areal(ss(4)).ge.-small) &
 
23097
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
21560
23098
    if (useboxc) then
21561
23099
      call boxc( rslt ,ss,rr ,as ,smax )
21562
23100
    else
21572
23110
                 .or.(     areal(ss(1)).ge.-small  &
21573
23111
                      .and.areal(ss(2)).ge.-small  &
21574
23112
                      .and.areal(ss(3)).ge.-small  &
21575
 
!OLD                      .and.areal(ss(4)).ge.-small) )
21576
 
                      .and.areal(ss(4)).ge.-small) & !NEW
21577
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
23113
                      .and.areal(ss(4)).ge.-small) &
 
23114
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
21578
23115
      if (useboxc) then
21579
23116
        call boxc( rslt ,ss,rr ,as ,smax )
21580
23117
      else
21818
23355
               .or.(     areal(ss(1)).ge.-small  &
21819
23356
                    .and.areal(ss(2)).ge.-small  &
21820
23357
                    .and.areal(ss(3)).ge.-small  &
21821
 
                    .and.areal(ss(4)).ge.-small) )
 
23358
                    .and.areal(ss(4)).ge.-small) &
 
23359
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
21822
23360
    if (useboxc) then
21823
23361
      call boxc( rslt ,ss,rr ,as ,smax )
21824
23362
    else
21834
23372
                 .or.(     areal(ss(1)).ge.-small  &
21835
23373
                      .and.areal(ss(2)).ge.-small  &
21836
23374
                      .and.areal(ss(3)).ge.-small  &
21837
 
!OLD                      .and.areal(ss(4)).ge.-small) )
21838
 
                      .and.areal(ss(4)).ge.-small) & !NEW
21839
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
23375
                      .and.areal(ss(4)).ge.-small) &
 
23376
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
21840
23377
      if (useboxc) then
21841
23378
        call boxc( rslt ,ss,rr ,as ,smax )
21842
23379
      else
22660
24197
      gg=xd1*pq1 ;hh=yd1*uv1
22661
24198
      rx2 = gg+hh
22662
24199
      if (abs(rx2).lt.neglig(prcpar)*max(abs(gg),abs(hh))) rx2 = 0
22663
 
    else
 
24200
    elseif (abs(pq2).gt.abs(pq1)) then
22664
24201
      rx2 = pq2
22665
24202
      gg=xd2*pq2 ;hh=yd2*uv2
22666
24203
      rx1 = gg+hh
22667
24204
      if (abs(rx1).lt.neglig(prcpar)*max(abs(gg),abs(hh))) rx1 = 0
 
24205
    else
 
24206
      rx1 = pq1
 
24207
      rx2 = pq2
22668
24208
    endif
22669
24209
    if (abs(uv1).gt.abs(uv2)) then
22670
24210
      ix1 = uv1
22671
24211
      gg=yd1*pq1 ;hh=xd1*uv1
22672
24212
      ix2 = gg-hh
22673
24213
      if (abs(ix2).lt.neglig(prcpar)*max(abs(gg),abs(hh))) ix2 = 0
22674
 
    else
 
24214
    elseif (abs(uv2).gt.abs(uv1)) then
22675
24215
      ix2 = uv2
22676
24216
      gg=yd2*pq2 ;hh=xd2*uv2
22677
24217
      ix1 = gg-hh
22678
24218
      if (abs(ix1).lt.neglig(prcpar)*max(abs(gg),abs(hh))) ix1 = 0
 
24219
    else
 
24220
      ix1 = uv1
 
24221
      ix2 = uv2
22679
24222
    endif
22680
24223
    x1 = acmplx(rx1,ix1)
22681
24224
    x2 = acmplx(rx2,ix2)
23017
24560
!***********************************************************************
23018
24561
! Provides the functions
23019
24562
!   olog(x,n) = log(x) + n*pi*imag  
23020
 
!   olog2(x,n) = olog(x,n)/(x-1)
 
24563
!   olog1(x,n) = olog(x,n)/(x-1)
 
24564
!   olog2(x,n) = ( olog1(x,n) - 1 )/(x-1)
 
24565
!   olog3(x,n) = ( olog2(x,n) + 1/2 )/(x-1)
23021
24566
! In the vicinity of x=1,n=0, the logarithm of complex argument is
23022
24567
! evaluated with a series expansion.
23023
24568
!***********************************************************************
23027
24572
  use avh_olo_mp_auxfun
23028
24573
  implicit none
23029
24574
  private
23030
 
  public :: update_olog,olog,olog2
 
24575
  public :: update_olog,olog,olog1,olog2,olog3
23031
24576
 
23032
24577
  type(mp_real) & 
23033
24578
         ,allocatable,save :: thrs(:,:)
23037
24582
  interface olog
23038
24583
    module procedure log_c,log_r
23039
24584
  end interface
 
24585
  interface olog1
 
24586
    module procedure log1_c,log1_r
 
24587
  end interface
23040
24588
  interface olog2
23041
24589
    module procedure log2_c,log2_r
23042
24590
  end interface
 
24591
  interface olog3
 
24592
    module procedure log3_c,log3_r
 
24593
  end interface
23043
24594
 
23044
24595
contains
23045
24596
 
23150
24701
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
23151
24702
                                else ;nn=ntrm(1,prcpar)
23152
24703
  endif
 
24704
! convergence acceleration using  z=(y-1)/(y+1)
 
24705
! rslt = 2 * ( z + z^3/3 + z^5/5 + z^7/7 + ... )  
23153
24706
  zz = zz/(yy+1)
23154
24707
  z2 = zz*zz
23155
24708
  aa = 2
23187
24740
  end function
23188
24741
 
23189
24742
 
23190
 
  function log2_c(xx,iph) result(rslt)
 
24743
  function log1_c(xx,iph) result(rslt)
23191
24744
!***********************************************************************
23192
24745
!***********************************************************************
23193
24746
  type(mp_complex) &  
23204
24757
!
23205
24758
  if (abs(imx).le.EPSN*abs(rex)) then
23206
24759
    if (rex.ge.RZRO) then
23207
 
      rslt = log2_r( rex, iph )
 
24760
      rslt = log1_r( rex, iph )
23208
24761
    else
23209
 
      rslt = log2_r(-rex, iph+sgnRe(imx) )
 
24762
      rslt = log1_r(-rex, iph+sgnRe(imx) )
23210
24763
    endif
23211
24764
    return
23212
24765
  endif
23232
24785
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
23233
24786
                                else ;nn=ntrm(1,prcpar)
23234
24787
  endif
 
24788
! convergence acceleration using  z=(y-1)/(y+1)
 
24789
! rslt = 2/(y+1) * ( 1 + z^2/3 + z^4/5 + z^6/7 + ... )  
23235
24790
  zz = zz/(yy+1)
23236
24791
  z2 = zz*zz
23237
24792
  aa = 2
23244
24799
  end function
23245
24800
 
23246
24801
 
23247
 
  function log2_r(xx,iph) result(rslt)
 
24802
  function log1_r(xx,iph) result(rslt)
23248
24803
!***********************************************************************
23249
24804
!***********************************************************************
23250
24805
  type(mp_real) & 
23260
24815
!  integer :: nn,ii
23261
24816
!
23262
24817
  if (xx.eq.RZRO) then
23263
 
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_r: ' &
 
24818
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log1_r: ' &
23264
24819
       ,'xx =',trim(myprint(xx)),', returning 0'
23265
24820
    rslt = 0
23266
24821
    return
23272
24827
!
23273
24828
  if (abs(yy-1).le.10*EPSN) then
23274
24829
    if (jj.ne.0) then
23275
 
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_r: ' &
 
24830
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log1_r: ' &
23276
24831
        ,'rr,jj =',trim(myprint(rr)),jj,', putting jj to 0'
23277
24832
    endif
23278
24833
    rslt = 1 - (yy-1)/2
23282
24837
  rslt = ( log(rr) + IPI*jj )/(yy-1)
23283
24838
  end function
23284
24839
 
 
24840
 
 
24841
  function log2_r(xx,iph) result(rslt)
 
24842
!***********************************************************************
 
24843
!***********************************************************************
 
24844
  type(mp_real) & 
 
24845
          ,intent(in) :: xx
 
24846
  integer ,intent(in) :: iph
 
24847
  type(mp_complex) &  
 
24848
    :: rslt
 
24849
  rslt = log2_c(xx*CONE,iph)
 
24850
  end function
 
24851
 
 
24852
 
 
24853
  function log2_c(xx,iph) result(rslt)
 
24854
!***********************************************************************
 
24855
!***********************************************************************
 
24856
  type(mp_complex) &  
 
24857
          ,intent(in) :: xx
 
24858
  integer ,intent(in) :: iph
 
24859
  type(mp_complex) &  
 
24860
    :: rslt ,yy,zz,z2
 
24861
  type(mp_real) & 
 
24862
    :: aa,rex,imx
 
24863
  integer :: nn,ii,jj
 
24864
!
 
24865
  rex = areal(xx)
 
24866
  imx = aimag(xx)
 
24867
!
 
24868
  if (rex.eq.RZRO.and.imx.eq.RZRO) then
 
24869
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_c: ' &
 
24870
       ,'xx = 0, returning 0'
 
24871
    rslt = 0
 
24872
    return
 
24873
  endif
 
24874
!
 
24875
  if (mod(iph,2).eq.0) then ;yy= xx ;jj=iph
 
24876
                       else ;yy=-xx ;jj=iph+sgnRe(imx)
 
24877
  endif
 
24878
!
 
24879
  if (jj.ne.0) then
 
24880
    rslt = ( olog1(yy,jj) - 1 )/(yy-1)
 
24881
    return
 
24882
  endif
 
24883
!
 
24884
  zz = yy-1
 
24885
  aa = abs(zz)
 
24886
  if     (aa.ge.thrs(6,prcpar)) then
 
24887
    rslt = (log(yy)/zz-1)/zz
 
24888
    return
 
24889
  elseif (aa.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
24890
  elseif (aa.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
24891
  elseif (aa.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
24892
  elseif (aa.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
24893
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
24894
                                else ;nn=ntrm(1,prcpar)
 
24895
  endif
 
24896
! convergence acceleration using  z=(y-1)/(y+1)
 
24897
! rslt = -1/(y+1) + 2/(y+1)^2 * ( z/3 + z^3/5 + z^5/7 + ... )  
 
24898
  zz = zz/(yy+1)
 
24899
  z2 = zz*zz
 
24900
  aa = 2
 
24901
  nn = 2*nn-1
 
24902
  rslt = aa/nn
 
24903
  do ii=nn-2,3,-2
 
24904
    rslt = aa/ii + z2*rslt
 
24905
  enddo
 
24906
  rslt = ( -1 + zz*rslt/(yy+1) )/(yy+1)
 
24907
  end function
 
24908
 
 
24909
 
 
24910
  function log3_r(xx,iph) result(rslt)
 
24911
!***********************************************************************
 
24912
!***********************************************************************
 
24913
  type(mp_real) & 
 
24914
          ,intent(in) :: xx
 
24915
  integer ,intent(in) :: iph
 
24916
  type(mp_complex) &  
 
24917
    :: rslt
 
24918
  rslt = log3_c(xx*CONE,iph)
 
24919
  end function
 
24920
 
 
24921
 
 
24922
  function log3_c(xx,iph) result(rslt)
 
24923
!***********************************************************************
 
24924
!***********************************************************************
 
24925
  type(mp_complex) &  
 
24926
          ,intent(in) :: xx
 
24927
  integer ,intent(in) :: iph
 
24928
  type(mp_complex) &  
 
24929
    :: rslt ,yy,zz,z2,HLF
 
24930
  type(mp_real) & 
 
24931
    :: aa,rex,imx
 
24932
  integer :: nn,ii,jj
 
24933
!
 
24934
  rex = areal(xx)
 
24935
  imx = aimag(xx)
 
24936
!
 
24937
  if (rex.eq.RZRO.and.imx.eq.RZRO) then
 
24938
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop log2_c: ' &
 
24939
       ,'xx = 0, returning 0'
 
24940
    rslt = 0
 
24941
    return
 
24942
  endif
 
24943
!
 
24944
  HLF = CONE/2
 
24945
!
 
24946
  if (mod(iph,2).eq.0) then ;yy= xx ;jj=iph
 
24947
                       else ;yy=-xx ;jj=iph+sgnRe(imx)
 
24948
  endif
 
24949
!
 
24950
  if (jj.ne.0) then
 
24951
    rslt = ( olog2(xx,jj) + HLF )/(yy-1)
 
24952
    return
 
24953
  endif
 
24954
!
 
24955
  zz = yy-1
 
24956
  aa = abs(zz)
 
24957
  if     (aa.ge.thrs(6,prcpar)) then
 
24958
    rslt = ((log(yy)/zz-1)/zz+HLF)/zz
 
24959
    return
 
24960
  elseif (aa.ge.thrs(5,prcpar)) then ;nn=ntrm(6,prcpar)
 
24961
  elseif (aa.ge.thrs(4,prcpar)) then ;nn=ntrm(5,prcpar)
 
24962
  elseif (aa.ge.thrs(3,prcpar)) then ;nn=ntrm(4,prcpar)
 
24963
  elseif (aa.ge.thrs(2,prcpar)) then ;nn=ntrm(3,prcpar)
 
24964
  elseif (aa.ge.thrs(1,prcpar)) then ;nn=ntrm(2,prcpar)
 
24965
                                else ;nn=ntrm(1,prcpar)
 
24966
  endif
 
24967
! convergence acceleration using  z=(y-1)/(y+1)
 
24968
! rslt = 1/(2*(y+1)) + 2/(y+1)^3 * ( 1/3 + z^2/5 + z^4/7 + ... )
 
24969
  zz = zz/(yy+1)
 
24970
  z2 = zz*zz
 
24971
  aa = 2
 
24972
  nn = 2*nn-1
 
24973
  rslt = aa/nn
 
24974
  do ii=nn-2,3,-2
 
24975
    rslt = aa/ii + z2*rslt
 
24976
  enddo
 
24977
  rslt = ( HLF + rslt/(yy+1)**2 )/(yy+1)
 
24978
  end function
 
24979
 
23285
24980
end module
23286
24981
 
23287
24982
 
 
24983
 
 
24984
 
23288
24985
module avh_olo_mp_dilog
23289
24986
!***********************************************************************
23290
24987
!                     /1    ln(1-zz*t)
23539
25236
! 
23540
25237
  if (yy.eq.RONE.and.odd.eq.0) then
23541
25238
    if (ntwo.ne.0) then
23542
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog_r: ' &
 
25239
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog_r: ' &
23543
25240
        ,'|x|,iph = ',trim(myprint(yy)),',',jj,', returning 0'
23544
25241
    endif
23545
25242
    rslt = 0
23647
25344
!
23648
25345
  if (j1.ne.j2) then
23649
25346
    if (r1.eq.r2) then
23650
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
25347
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
23651
25348
        ,'j1,j2,r1-r2',j1,j2,',',trim(myprint(r1-r2)),', returning 0'
23652
25349
      rslt = 0
23653
25350
!      write(*,*) 'dilog2_c j1=/=j2,r1=r2' !DEBUG
23661
25358
!
23662
25359
  if (a1.lt.eps) then
23663
25360
    if (a2.lt.eps) then
23664
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
25361
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
23665
25362
        ,'r1,r2 =',trim(myprint(r1)),',',trim(myprint(r2)),', returning 0'
23666
25363
      rslt = 0
23667
25364
!      write(*,*) 'dilog2_c r1<eps,r2<eps' !DEBUG
23683
25380
!      write(*,*) 'dilog2_c ||1-y1|/|1-y2|-1|>0.1' !DEBUG
23684
25381
      return
23685
25382
    elseif (oo.eq.0.and.ao1.lt.eps) then
23686
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
25383
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
23687
25384
        ,'r1,oo,nn =',trim(myprint(r1)),',',oo,nn,', putting nn=0'
23688
25385
      if (ao2.lt.eps) then
23689
25386
        rslt = -1
23693
25390
        y1=1-eps ;nn=0 ;logr1=0 ;r1=1-eps
23694
25391
      endif
23695
25392
    elseif (oo.eq.0.and.ao2.lt.eps) then
23696
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
25393
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
23697
25394
        ,'r2,oo,nn =',trim(myprint(r2)),',',oo,nn,', putting nn=0'
23698
25395
      y2=1-eps ;nn=0 ;logr2=0 ;r2=1-eps
23699
25396
    endif
23708
25405
      if (a1.gt.RONE) ii = ii + (nn+pp(oo,sgnIm(y2)))
23709
25406
      if (a2.gt.RONE) ii = ii - (nn+pp(oo,sgnIm(y2)))
23710
25407
      ii = nn*ii
23711
 
      if (ii.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
 
25408
      if (ii.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_c: ' &
23712
25409
        ,'r1,r2,nn =',trim(myprint(r1)),',',trim(myprint(r2)),',',nn &
23713
25410
        ,', putting nn=0'
23714
 
      rslt = -olog2(y2,0)
 
25411
      rslt = -olog1(y2,0)
23715
25412
!      write(*,*) 'dilog2_c |logr1/lorg2|<eps' !DEBUG
23716
25413
      return
23717
25414
    endif
23723
25420
    nn=-nn ;oo=-oo
23724
25421
  endif
23725
25422
!
23726
 
  ff=y1/y2         ;ff=-olog2(ff,0)/y2
23727
 
  gg=(1-y1)/(1-y2) ;gg=-olog2(gg,0)/(1-y2)
 
25423
  ff=y1/y2         ;ff=-olog1(ff,0)/y2
 
25424
  gg=(1-y1)/(1-y2) ;gg=-olog1(gg,0)/(1-y2)
23728
25425
!
23729
25426
  if (2*areal(y1).ge.RONE) then
23730
25427
!    write(*,*) 'dilog2_c re>1/2' !DEBUG
23781
25478
!
23782
25479
  if (j1.ne.j2) then
23783
25480
    if (r1.eq.r2) then
23784
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
25481
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
23785
25482
        ,'j1,j2,r1-r2',j1,j2,',',trim(myprint(r1-r2)),', returning 0'
23786
25483
      rslt = 0
23787
25484
!      write(*,*) 'dilog2_r j1=/=j2,r1=r2' !DEBUG
23795
25492
!
23796
25493
  if (r1.lt.eps) then
23797
25494
    if (r2.lt.eps) then
23798
 
      if (dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
25495
      if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
23799
25496
        ,'r1,r2 =',trim(myprint(r1)),',',trim(myprint(r2)),', returning 0'
23800
25497
      rslt = 0
23801
25498
!      write(*,*) 'dilog2_r r1<eps,r2<eps' !DEBUG
23817
25514
!      write(*,*) 'dilog2_r ||1-y1|/|1-y2|-1|>0.1' !DEBUG
23818
25515
      return
23819
25516
    elseif (oo.eq.0.and.ro1.lt.eps) then
23820
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
25517
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
23821
25518
        ,'r1,oo,nn =',trim(myprint(r1)),',',oo,nn,', putting nn=0'
23822
25519
      if (ro2.lt.eps) then
23823
25520
        rslt = -1
23827
25524
        y1=1-eps ;nn=0 ;logr1=0 ;r1=1-eps
23828
25525
      endif
23829
25526
    elseif (oo.eq.0.and.ro2.lt.eps) then
23830
 
      if (nn.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
25527
      if (nn.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
23831
25528
        ,'r2,oo,nn =',trim(myprint(r2)),',',oo,nn,', putting nn=0'
23832
25529
      y2=1-eps ;nn=0 ;logr2=0 ;r2=1-eps
23833
25530
    endif
23842
25539
      if (r1.gt.RONE) ii = ii + (nn+2*oo)
23843
25540
      if (r2.gt.RONE) ii = ii - (nn+2*oo)
23844
25541
      ii = nn*ii
23845
 
      if (ii.ne.0.and.dunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
 
25542
      if (ii.ne.0.and.eunit.gt.0) write(eunit,*) 'ERROR in OneLOop dilog2_r: ' &
23846
25543
        ,'r1,r2,nn =',trim(myprint(r1)),',',trim(myprint(r2)),',',nn &
23847
25544
        ,', putting nn=0'
23848
 
      rslt = -olog2(y2,0)
 
25545
      rslt = -olog1(y2,0)
23849
25546
!      write(*,*) 'dilog2_r |logr1/lorg2|<eps' !DEBUG
23850
25547
      return
23851
25548
    endif
23857
25554
    nn=-nn ;oo=-oo
23858
25555
  endif
23859
25556
!
23860
 
  ff=y1/y2         ;ff=-olog2(ff,0)/y2
23861
 
  gg=(1-y1)/(1-y2) ;gg=-olog2(gg,0)/(1-y2)
 
25557
  ff=y1/y2         ;ff=-olog1(ff,0)/y2
 
25558
  gg=(1-y1)/(1-y2) ;gg=-olog1(gg,0)/(1-y2)
23862
25559
!
23863
25560
  if (2*y1.ge.RONE) then
23864
25561
!    write(*,*) 'dilog2_r re>1/2' !DEBUG
24288
25985
 
24289
25986
  implicit none
24290
25987
  private
24291
 
  public :: qmplx_type,qonv,directly,sheet,logc,logc2,li2c,li2c2
 
25988
  public :: qmplx_type,qonv,directly,sheet,logc,logc2,logc3,li2c,li2c2
24292
25989
  public :: operator (*) ,operator (/)
24293
25990
 
24294
25991
  type :: qmplx_type
24535
26232
  type(qmplx_type) ,intent(in) :: xx
24536
26233
  type(mp_complex) &  
24537
26234
    :: rslt
24538
 
!  rslt = -olog2(acmplx(xx%c),xx%p)
24539
 
  rslt = -olog2(xx%c,xx%p)
 
26235
!  rslt = -olog1(acmplx(xx%c),xx%p)
 
26236
  rslt = -olog1(xx%c,xx%p)
 
26237
  end function
 
26238
 
 
26239
  function logc3(xx) result(rslt)
 
26240
!*******************************************************************
 
26241
!  ( log(xx)/(1-xx) + 1 )/(1-xx)
 
26242
!*******************************************************************
 
26243
  type(qmplx_type) ,intent(in) :: xx
 
26244
  type(mp_complex) &  
 
26245
    :: rslt
 
26246
!  rslt = olog2(acmplx(xx%c),xx%p)
 
26247
  rslt = olog2(xx%c,xx%p)
24540
26248
  end function
24541
26249
 
24542
26250
  function li2c(xx) result(rslt)
24576
26284
  use avh_olo_mp_auxfun
24577
26285
  use avh_olo_mp_bnlog
24578
26286
  use avh_olo_mp_qmplx
 
26287
  use avh_olo_mp_olog
24579
26288
  implicit none
24580
26289
  private
24581
 
  public :: tadp ,tadpn ,bub0 ,bub1 ,bub11 ,bub111 ,bub1111
 
26290
  public :: tadp ,tadpn ,bub0 ,dbub0 ,bub1 ,bub11 ,bub111 ,bub1111
24582
26291
 
24583
26292
contains
24584
26293
 
25169
26878
  b0011(0) = ( a0(0,0) - x1*b111(0) + x2*b11(0) + 4*b0011(1) )/10
25170
26879
  end subroutine
25171
26880
 
 
26881
 
 
26882
!*******************************************************************
 
26883
! Derivative of B0
 
26884
! expects  m0<m1
 
26885
! only finite case, so input must not be  m0=0 & m1=pp
 
26886
!*******************************************************************
 
26887
 
 
26888
  subroutine dbub0( rslt &
 
26889
                   ,pp,m0,m1 ,app,am0,am1 )
 
26890
  type(mp_complex) &  
 
26891
    ,intent(out) :: rslt
 
26892
  type(mp_complex) &  
 
26893
    ,intent(in)  :: pp,m0,m1
 
26894
  type(mp_real) & 
 
26895
    ,intent(in)  :: app,am0,am1
 
26896
  type(mp_complex) &  
 
26897
    :: ch,x1,x2,lambda
 
26898
  type(mp_real) & 
 
26899
    :: ax1,ax2,ax1x2,maxa
 
26900
  type(qmplx_type) :: q1,q2,q1o,q2o
 
26901
  integer :: sgn
 
26902
!
 
26903
  if (am1.eq.RZRO) then
 
26904
    if (app.eq.RZRO) then
 
26905
      rslt = 0
 
26906
      return
 
26907
    endif
 
26908
  endif
 
26909
!
 
26910
  if (app.eq.RZRO) then
 
26911
    if (abs(m0-m1).le.am1*EPSN*10) then
 
26912
      rslt = 1/(6*m1)
 
26913
    else
 
26914
      ch = m0/m1
 
26915
      rslt = ( CONE/2 - ch*olog3(ch,0) )/m1 
 
26916
    endif
 
26917
  elseif (am1.eq.RZRO) then
 
26918
    rslt =-1/pp
 
26919
  else
 
26920
    call solabc( x1,x2 ,lambda ,pp ,(m0-m1)-pp ,m1 ,0 )
 
26921
    sgn =-sgnRe(pp)*sgnRe(x2-x1)
 
26922
    q1  = qonv(x1  , sgn)
 
26923
    q1o = qonv(x1-1, sgn)
 
26924
    q2  = qonv(x2  ,-sgn)
 
26925
    q2o = qonv(x2-1,-sgn)
 
26926
    ax1 = abs(x1)
 
26927
    ax2 = abs(x2)
 
26928
    ax1x2 = abs(x1-x2)
 
26929
    maxa = max(ax1,ax2)
 
26930
    if (ax1x2.lt.maxa*EPSN*10) then
 
26931
      rslt = ( (x1+x2-1)*logc(q2/q2o) - 2 )/pp
 
26932
    elseif (ax1x2*2.lt.maxa) then
 
26933
      if     (x1.eq.CZRO.or.x1.eq.CONE) then
 
26934
        rslt = ( (x1+x2-1)*logc(q2/q2o) - 1 )/pp
 
26935
      elseif (x2.eq.CZRO.or.x2.eq.CONE) then
 
26936
        rslt = ( (x1+x2-1)*logc(q1/q1o) - 1 )/pp
 
26937
      else
 
26938
        rslt = x1*(x1-1)*( logc2(q1o/q2o)/(x2-1) - logc2(q1/q2)/x2 ) &
 
26939
             + (x1+x2-1)*logc(q2/q2o) - 1
 
26940
        rslt = rslt/pp
 
26941
      endif
 
26942
    else
 
26943
      rslt = 0
 
26944
      if (ax1.ne.RZRO) then
 
26945
        if (ax1.lt.2*RONE) then
 
26946
          rslt = rslt - x1
 
26947
          if (x1.ne.CONE) rslt = rslt - x1*logc2(q1/q1o)
 
26948
        else
 
26949
          rslt = rslt + x1/(x1-1)*logc3(q1/q1o)
 
26950
        endif
 
26951
      endif
 
26952
      if (ax2.ne.RZRO) then
 
26953
        if (ax2.lt.2*RONE) then
 
26954
          rslt = rslt + x2
 
26955
          if (x2.ne.CONE) rslt = rslt + x2*logc2(q2/q2o)
 
26956
        else
 
26957
          rslt = rslt - x2/(x2-1)*logc3(q2/q2o)
 
26958
        endif
 
26959
      endif
 
26960
      rslt = rslt/lambda
 
26961
    endif
 
26962
  endif
 
26963
!
 
26964
  end subroutine
 
26965
 
 
26966
 
25172
26967
end module
25173
26968
 
25174
26969
 
25402
27197
    log2 = olog( abs(rp2/rmu2) ,i2 )
25403
27198
    log3 = olog( abs(rp3/rmu2) ,i3 )
25404
27199
    rslt(2) = 0
25405
 
    rslt(1) = -olog2( abs(rp3/rp2) ,i3-i2 )/rp2
 
27200
    rslt(1) = -olog1( abs(rp3/rp2) ,i3-i2 )/rp2
25406
27201
    rslt(0) = -rslt(1)*(log3+log2)/2
25407
27202
  elseif (icase.eq.3) then
25408
27203
! 3 masses non-zero
27975
29770
     ,intent(in) :: y1,y2
27976
29771
  type(mp_complex) &  
27977
29772
     :: rslt ,oy1,oy2
 
29773
!
27978
29774
   oy1 = 1-y1
27979
29775
   oy2 = 1-y2
27980
29776
   rslt = logc2( qonv(-y2)/qonv(-y1) )/y1 &
28043
29839
  private
28044
29840
  public :: olo_unit ,olo_scale ,olo_onshell ,olo_setting
28045
29841
  public :: olo_precision
28046
 
  public :: olo_a0 ,olo_b0 ,olo_b11 ,olo_c0 ,olo_d0
 
29842
  public :: olo_a0 ,olo_b0 ,olo_db0 ,olo_b11 ,olo_c0 ,olo_d0
28047
29843
  public :: olo_an ,olo_bn
28048
29844
  public :: olo
28049
29845
  public :: olo_get_scale ,olo_get_onshell ,olo_get_precision
28050
 
  public :: a0_r,a0rr,a0_c,a0cr
28051
 
  public :: an_r,anrr,an_c,ancr
28052
 
  public :: b0rr,b0rrr,b0rc,b0rcr,b0cc,b0ccr
28053
 
  public :: b11rr,b11rrr,b11rc,b11rcr,b11cc,b11ccr
28054
 
  public :: bnrr,bnrrr,bnrc,bnrcr,bncc,bnccr
28055
 
  public :: c0rr,c0rrr,c0rc,c0rcr,c0cc,c0ccr
28056
 
  public :: d0rr,d0rrr,d0rc,d0rcr,d0cc,d0ccr
28057
29846
!
28058
29847
 integer ,public ,parameter :: olo_kind=kind(1d0) 
28059
29848
!
28079
29868
  interface olo_b0
28080
29869
    module procedure b0rr,b0rrr,b0rc,b0rcr,b0cc,b0ccr
28081
29870
  end interface 
 
29871
  interface olo_db0
 
29872
    module procedure db0rr,db0rrr,db0rc,db0rcr,db0cc,db0ccr
 
29873
  end interface 
28082
29874
  interface olo_b11
28083
29875
    module procedure b11rr,b11rrr,b11rc,b11rcr,b11cc,b11ccr
28084
29876
  end interface 
29065
30857
  endif
29066
30858
  end subroutine
29067
30859
 
 
30860
!*******************************************************************
 
30861
! Derivative of B0
 
30862
!*******************************************************************
 
30863
  subroutine db0cc( rslt ,pp,m1,m2 )
 
30864
!
 
30865
  use avh_olo_mp_bub ,only: dbub0
 
30866
  use avh_olo_mp_olog ,only: olog
 
30867
!
 
30868
  type(mp_complex) &  
 
30869
    ,intent(out) :: rslt(0:2)
 
30870
  type(mp_complex) &  
 
30871
    ,intent(in)  :: pp
 
30872
  type(mp_complex) &  
 
30873
    ,intent(in)  :: m1,m2
 
30874
!
 
30875
  type(mp_complex) &  
 
30876
    :: ss,r1,r2,ch
 
30877
  type(mp_real) & 
 
30878
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
30879
  character(26+99) ,parameter :: warning=&
 
30880
                     'WARNING from OneLOop db0: '//warnonshell
 
30881
  if (initz) call init
 
30882
  ss = pp
 
30883
  r1 = m1
 
30884
  r2 = m2
 
30885
!
 
30886
  app = areal(ss)
 
30887
  if (aimag(ss).ne.RZRO) then
 
30888
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
30889
      ,'ss has non-zero imaginary part, putting it to zero.'
 
30890
    ss = acmplx( app )
 
30891
  endif
 
30892
  app = abs(app)
 
30893
!
 
30894
  am1 = areal(r1)
 
30895
  hh  = aimag(r1)
 
30896
  if (hh.gt.RZRO) then
 
30897
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
30898
      ,'r1 has positive imaginary part, switching its sign.'
 
30899
    r1 = acmplx( am1 ,-hh )
 
30900
  endif
 
30901
  am1 = abs(am1) + abs(hh)
 
30902
!
 
30903
  am2 = areal(r2)
 
30904
  hh  = aimag(r2)
 
30905
  if (hh.gt.RZRO) then
 
30906
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
30907
      ,'r2 has positive imaginary part, switching its sign.'
 
30908
    r2 = acmplx( am2 ,-hh )
 
30909
  endif
 
30910
  am2 = abs(am2) + abs(hh)
 
30911
!
 
30912
  mulocal = muscale 
 
30913
!
 
30914
  mulocal2 = mulocal*mulocal
 
30915
!
 
30916
  if (am1.gt.am2) then
 
30917
    ch=r1 ; r1=r2 ; r2=ch
 
30918
    hh=am1;am1=am2;am2=hh
 
30919
  endif
 
30920
  ssr2 = abs(ss-r2)
 
30921
!
 
30922
  if (nonzerothrs) then
 
30923
    hh = onshellthrs
 
30924
    if (app.lt.hh) app = 0
 
30925
    if (am1.lt.hh) am1 = 0
 
30926
    if (am2.lt.hh) am2 = 0
 
30927
    if (ssr2.lt.hh) ssr2 = 0
 
30928
  elseif (wunit.gt.0) then
 
30929
    hh = onshellthrs*max(app,max(am1,am2))
 
30930
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
30931
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
30932
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
30933
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
30934
  endif
 
30935
!
 
30936
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
30937
!
 
30938
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
30939
    rslt(1) =-1/(2*ss)
 
30940
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
30941
  else
 
30942
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
30943
  endif
 
30944
!
 
30945
  if (punit.gt.0) then
 
30946
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
30947
    write(punit,*) ' pp:',trim(myprint(pp))
 
30948
    write(punit,*) ' m1:',trim(myprint(m1))
 
30949
    write(punit,*) ' m2:',trim(myprint(m2))
 
30950
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
30951
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
30952
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
30953
  endif
 
30954
  end subroutine
 
30955
 
 
30956
  subroutine db0ccr( rslt ,pp,m1,m2 ,rmu )
 
30957
!
 
30958
  use avh_olo_mp_bub ,only: dbub0
 
30959
  use avh_olo_mp_olog ,only: olog
 
30960
!
 
30961
  type(mp_complex) &  
 
30962
    ,intent(out) :: rslt(0:2)
 
30963
  type(mp_complex) &  
 
30964
    ,intent(in)  :: pp
 
30965
  type(mp_complex) &  
 
30966
    ,intent(in)  :: m1,m2
 
30967
  type(mp_real) & 
 
30968
   ,intent(in)  :: rmu       
 
30969
!
 
30970
  type(mp_complex) &  
 
30971
    :: ss,r1,r2,ch
 
30972
  type(mp_real) & 
 
30973
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
30974
  character(26+99) ,parameter :: warning=&
 
30975
                     'WARNING from OneLOop db0: '//warnonshell
 
30976
  if (initz) call init
 
30977
  ss = pp
 
30978
  r1 = m1
 
30979
  r2 = m2
 
30980
!
 
30981
  app = areal(ss)
 
30982
  if (aimag(ss).ne.RZRO) then
 
30983
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
30984
      ,'ss has non-zero imaginary part, putting it to zero.'
 
30985
    ss = acmplx( app )
 
30986
  endif
 
30987
  app = abs(app)
 
30988
!
 
30989
  am1 = areal(r1)
 
30990
  hh  = aimag(r1)
 
30991
  if (hh.gt.RZRO) then
 
30992
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
30993
      ,'r1 has positive imaginary part, switching its sign.'
 
30994
    r1 = acmplx( am1 ,-hh )
 
30995
  endif
 
30996
  am1 = abs(am1) + abs(hh)
 
30997
!
 
30998
  am2 = areal(r2)
 
30999
  hh  = aimag(r2)
 
31000
  if (hh.gt.RZRO) then
 
31001
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
31002
      ,'r2 has positive imaginary part, switching its sign.'
 
31003
    r2 = acmplx( am2 ,-hh )
 
31004
  endif
 
31005
  am2 = abs(am2) + abs(hh)
 
31006
!
 
31007
  mulocal = rmu     
 
31008
!
 
31009
  mulocal2 = mulocal*mulocal
 
31010
!
 
31011
  if (am1.gt.am2) then
 
31012
    ch=r1 ; r1=r2 ; r2=ch
 
31013
    hh=am1;am1=am2;am2=hh
 
31014
  endif
 
31015
  ssr2 = abs(ss-r2)
 
31016
!
 
31017
  if (nonzerothrs) then
 
31018
    hh = onshellthrs
 
31019
    if (app.lt.hh) app = 0
 
31020
    if (am1.lt.hh) am1 = 0
 
31021
    if (am2.lt.hh) am2 = 0
 
31022
    if (ssr2.lt.hh) ssr2 = 0
 
31023
  elseif (wunit.gt.0) then
 
31024
    hh = onshellthrs*max(app,max(am1,am2))
 
31025
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
31026
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
31027
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
31028
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
31029
  endif
 
31030
!
 
31031
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
31032
!
 
31033
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
31034
    rslt(1) =-1/(2*ss)
 
31035
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
31036
  else
 
31037
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
31038
  endif
 
31039
!
 
31040
  if (punit.gt.0) then
 
31041
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
31042
    write(punit,*) ' pp:',trim(myprint(pp))
 
31043
    write(punit,*) ' m1:',trim(myprint(m1))
 
31044
    write(punit,*) ' m2:',trim(myprint(m2))
 
31045
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
31046
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
31047
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
31048
  endif
 
31049
  end subroutine
 
31050
 
 
31051
  subroutine db0rc( rslt ,pp ,m1,m2 )
 
31052
!
 
31053
  use avh_olo_mp_bub ,only: dbub0
 
31054
  use avh_olo_mp_olog ,only: olog
 
31055
!
 
31056
  type(mp_complex) &  
 
31057
    ,intent(out) :: rslt(0:2)
 
31058
  type(mp_real) & 
 
31059
    ,intent(in)  :: pp
 
31060
  type(mp_complex) &  
 
31061
    ,intent(in)  :: m1,m2
 
31062
!
 
31063
  type(mp_complex) &  
 
31064
    :: ss,r1,r2,ch
 
31065
  type(mp_real) & 
 
31066
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
31067
  character(26+99) ,parameter :: warning=&
 
31068
                     'WARNING from OneLOop db0: '//warnonshell
 
31069
  if (initz) call init
 
31070
  ss = pp
 
31071
  r1 = m1
 
31072
  r2 = m2
 
31073
!
 
31074
  app = abs(pp)
 
31075
!
 
31076
  am1 = areal(r1)
 
31077
  hh  = aimag(r1)
 
31078
  if (hh.gt.RZRO) then
 
31079
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
31080
      ,'r1 has positive imaginary part, switching its sign.'
 
31081
    r1 = acmplx( am1 ,-hh )
 
31082
  endif
 
31083
  am1 = abs(am1) + abs(hh)
 
31084
!
 
31085
  am2 = areal(r2)
 
31086
  hh  = aimag(r2)
 
31087
  if (hh.gt.RZRO) then
 
31088
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
31089
      ,'r2 has positive imaginary part, switching its sign.'
 
31090
    r2 = acmplx( am2 ,-hh )
 
31091
  endif
 
31092
  am2 = abs(am2) + abs(hh)
 
31093
!
 
31094
  mulocal = muscale 
 
31095
!
 
31096
  mulocal2 = mulocal*mulocal
 
31097
!
 
31098
  if (am1.gt.am2) then
 
31099
    ch=r1 ; r1=r2 ; r2=ch
 
31100
    hh=am1;am1=am2;am2=hh
 
31101
  endif
 
31102
  ssr2 = abs(ss-r2)
 
31103
!
 
31104
  if (nonzerothrs) then
 
31105
    hh = onshellthrs
 
31106
    if (app.lt.hh) app = 0
 
31107
    if (am1.lt.hh) am1 = 0
 
31108
    if (am2.lt.hh) am2 = 0
 
31109
    if (ssr2.lt.hh) ssr2 = 0
 
31110
  elseif (wunit.gt.0) then
 
31111
    hh = onshellthrs*max(app,max(am1,am2))
 
31112
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
31113
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
31114
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
31115
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
31116
  endif
 
31117
!
 
31118
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
31119
!
 
31120
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
31121
    rslt(1) =-1/(2*ss)
 
31122
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
31123
  else
 
31124
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
31125
  endif
 
31126
!
 
31127
  if (punit.gt.0) then
 
31128
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
31129
    write(punit,*) ' pp:',trim(myprint(pp))
 
31130
    write(punit,*) ' m1:',trim(myprint(m1))
 
31131
    write(punit,*) ' m2:',trim(myprint(m2))
 
31132
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
31133
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
31134
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
31135
  endif
 
31136
  end subroutine
 
31137
 
 
31138
  subroutine db0rcr( rslt ,pp,m1,m2 ,rmu )
 
31139
!
 
31140
  use avh_olo_mp_bub ,only: dbub0
 
31141
  use avh_olo_mp_olog ,only: olog
 
31142
!
 
31143
  type(mp_complex) &  
 
31144
    ,intent(out) :: rslt(0:2)
 
31145
  type(mp_real) & 
 
31146
    ,intent(in)  :: pp
 
31147
  type(mp_complex) &  
 
31148
    ,intent(in)  :: m1,m2
 
31149
  type(mp_real) & 
 
31150
   ,intent(in)  :: rmu       
 
31151
!
 
31152
  type(mp_complex) &  
 
31153
    :: ss,r1,r2,ch
 
31154
  type(mp_real) & 
 
31155
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
31156
  character(26+99) ,parameter :: warning=&
 
31157
                     'WARNING from OneLOop db0: '//warnonshell
 
31158
  if (initz) call init
 
31159
  ss = pp
 
31160
  r1 = m1
 
31161
  r2 = m2
 
31162
!
 
31163
  app = abs(pp)
 
31164
!
 
31165
  am1 = areal(r1)
 
31166
  hh  = aimag(r1)
 
31167
  if (hh.gt.RZRO) then
 
31168
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop db0: ' &
 
31169
      ,'r1 has positive imaginary part, switching its sign.'
 
31170
    r1 = acmplx( am1 ,-hh )
 
31171
  endif
 
31172
  am1 = abs(am1) + abs(hh)
 
31173
!
 
31174
  am2 = areal(r2)
 
31175
  hh  = aimag(r2)
 
31176
  if (hh.gt.RZRO) then
 
31177
    if (eunit.gt.0) write(eunit,*) 'ERROR in OneLOop b0: ' &
 
31178
      ,'r2 has positive imaginary part, switching its sign.'
 
31179
    r2 = acmplx( am2 ,-hh )
 
31180
  endif
 
31181
  am2 = abs(am2) + abs(hh)
 
31182
!
 
31183
  mulocal = rmu     
 
31184
!
 
31185
  mulocal2 = mulocal*mulocal
 
31186
!
 
31187
  if (am1.gt.am2) then
 
31188
    ch=r1 ; r1=r2 ; r2=ch
 
31189
    hh=am1;am1=am2;am2=hh
 
31190
  endif
 
31191
  ssr2 = abs(ss-r2)
 
31192
!
 
31193
  if (nonzerothrs) then
 
31194
    hh = onshellthrs
 
31195
    if (app.lt.hh) app = 0
 
31196
    if (am1.lt.hh) am1 = 0
 
31197
    if (am2.lt.hh) am2 = 0
 
31198
    if (ssr2.lt.hh) ssr2 = 0
 
31199
  elseif (wunit.gt.0) then
 
31200
    hh = onshellthrs*max(app,max(am1,am2))
 
31201
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
31202
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
31203
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
31204
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
31205
  endif
 
31206
!
 
31207
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
31208
!
 
31209
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
31210
    rslt(1) =-1/(2*ss)
 
31211
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
31212
  else
 
31213
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
31214
  endif
 
31215
!
 
31216
  if (punit.gt.0) then
 
31217
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
31218
    write(punit,*) ' pp:',trim(myprint(pp))
 
31219
    write(punit,*) ' m1:',trim(myprint(m1))
 
31220
    write(punit,*) ' m2:',trim(myprint(m2))
 
31221
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
31222
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
31223
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
31224
  endif
 
31225
  end subroutine
 
31226
 
 
31227
  subroutine db0rr( rslt ,pp ,m1,m2 )
 
31228
!
 
31229
  use avh_olo_mp_bub ,only: dbub0
 
31230
  use avh_olo_mp_olog ,only: olog
 
31231
!
 
31232
  type(mp_complex) &  
 
31233
    ,intent(out) :: rslt(0:2)
 
31234
  type(mp_real) & 
 
31235
    ,intent(in)  :: pp
 
31236
  type(mp_real) & 
 
31237
    ,intent(in)  :: m1,m2
 
31238
!
 
31239
  type(mp_complex) &  
 
31240
    :: ss,r1,r2,ch
 
31241
  type(mp_real) & 
 
31242
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
31243
  character(26+99) ,parameter :: warning=&
 
31244
                     'WARNING from OneLOop db0: '//warnonshell
 
31245
  if (initz) call init
 
31246
  ss = pp
 
31247
  r1 = m1
 
31248
  r2 = m2
 
31249
!
 
31250
  app = abs(pp)
 
31251
!
 
31252
  am1 = abs(m1)
 
31253
  am2 = abs(m2)
 
31254
!
 
31255
  mulocal = muscale 
 
31256
!
 
31257
  mulocal2 = mulocal*mulocal
 
31258
!
 
31259
  if (am1.gt.am2) then
 
31260
    ch=r1 ; r1=r2 ; r2=ch
 
31261
    hh=am1;am1=am2;am2=hh
 
31262
  endif
 
31263
  ssr2 = abs(ss-r2)
 
31264
!
 
31265
  if (nonzerothrs) then
 
31266
    hh = onshellthrs
 
31267
    if (app.lt.hh) app = 0
 
31268
    if (am1.lt.hh) am1 = 0
 
31269
    if (am2.lt.hh) am2 = 0
 
31270
    if (ssr2.lt.hh) ssr2 = 0
 
31271
  elseif (wunit.gt.0) then
 
31272
    hh = onshellthrs*max(app,max(am1,am2))
 
31273
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
31274
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
31275
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
31276
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
31277
  endif
 
31278
!
 
31279
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
31280
!
 
31281
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
31282
    rslt(1) =-1/(2*ss)
 
31283
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
31284
  else
 
31285
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
31286
  endif
 
31287
!
 
31288
  if (punit.gt.0) then
 
31289
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
31290
    write(punit,*) ' pp:',trim(myprint(pp))
 
31291
    write(punit,*) ' m1:',trim(myprint(m1))
 
31292
    write(punit,*) ' m2:',trim(myprint(m2))
 
31293
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
31294
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
31295
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
31296
  endif
 
31297
  end subroutine
 
31298
 
 
31299
  subroutine db0rrr( rslt ,pp ,m1,m2 ,rmu )
 
31300
!
 
31301
  use avh_olo_mp_bub ,only: dbub0
 
31302
  use avh_olo_mp_olog ,only: olog
 
31303
!
 
31304
  type(mp_complex) &  
 
31305
    ,intent(out) :: rslt(0:2)
 
31306
  type(mp_real) & 
 
31307
    ,intent(in)  :: pp
 
31308
  type(mp_real) & 
 
31309
    ,intent(in)  :: m1,m2
 
31310
  type(mp_real) & 
 
31311
   ,intent(in)  :: rmu       
 
31312
!
 
31313
  type(mp_complex) &  
 
31314
    :: ss,r1,r2,ch
 
31315
  type(mp_real) & 
 
31316
    :: app,am1,am2,hh,mulocal,mulocal2,ssr2
 
31317
  character(26+99) ,parameter :: warning=&
 
31318
                     'WARNING from OneLOop db0: '//warnonshell
 
31319
  if (initz) call init
 
31320
  ss = pp
 
31321
  r1 = m1
 
31322
  r2 = m2
 
31323
!
 
31324
  app = abs(pp)
 
31325
!
 
31326
  am1 = abs(m1)
 
31327
  am2 = abs(m2)
 
31328
!
 
31329
  mulocal = rmu     
 
31330
!
 
31331
  mulocal2 = mulocal*mulocal
 
31332
!
 
31333
  if (am1.gt.am2) then
 
31334
    ch=r1 ; r1=r2 ; r2=ch
 
31335
    hh=am1;am1=am2;am2=hh
 
31336
  endif
 
31337
  ssr2 = abs(ss-r2)
 
31338
!
 
31339
  if (nonzerothrs) then
 
31340
    hh = onshellthrs
 
31341
    if (app.lt.hh) app = 0
 
31342
    if (am1.lt.hh) am1 = 0
 
31343
    if (am2.lt.hh) am2 = 0
 
31344
    if (ssr2.lt.hh) ssr2 = 0
 
31345
  elseif (wunit.gt.0) then
 
31346
    hh = onshellthrs*max(app,max(am1,am2))
 
31347
    if (RZRO.lt.app.and.app.lt.hh) write(wunit,*) warning
 
31348
    if (RZRO.lt.am1.and.am1.lt.hh) write(wunit,*) warning
 
31349
    if (RZRO.lt.am2.and.am2.lt.hh) write(wunit,*) warning
 
31350
    if (RZRO.lt.ssr2.and.ssr2.lt.hh) write(wunit,*) warning
 
31351
  endif
 
31352
!
 
31353
  rslt(0)=0;rslt(1)=0;rslt(2)=0
 
31354
!
 
31355
  if (am1.eq.RZRO.and.ssr2.eq.RZRO) then
 
31356
    rslt(1) =-1/(2*ss)
 
31357
    rslt(0) =-( 1 + olog(mulocal2/ss,0)/2 )/ss
 
31358
  else
 
31359
    call dbub0( rslt(0) ,ss,r1,r2 ,app,am1,am2 )
 
31360
  endif
 
31361
!
 
31362
  if (punit.gt.0) then
 
31363
    if (nonzerothrs) write(punit,*) 'onshell:',trim(myprint(onshellthrs))
 
31364
    write(punit,*) ' pp:',trim(myprint(pp))
 
31365
    write(punit,*) ' m1:',trim(myprint(m1))
 
31366
    write(punit,*) ' m2:',trim(myprint(m2))
 
31367
    write(punit,*) 'db0(2):',trim(myprint(rslt(2)))
 
31368
    write(punit,*) 'db0(1):',trim(myprint(rslt(1)))
 
31369
    write(punit,*) 'db0(0):',trim(myprint(rslt(0)))
 
31370
  endif
 
31371
  end subroutine
 
31372
 
 
31373
 
29068
31374
 
29069
31375
!*******************************************************************
29070
31376
! Return the Papparino-Veltman functions b11,b00,b1,b0 , for
31436
33742
               .or.(     areal(ss(1)).ge.-small  &
31437
33743
                    .and.areal(ss(2)).ge.-small  &
31438
33744
                    .and.areal(ss(3)).ge.-small  &
31439
 
                    .and.areal(ss(4)).ge.-small) )
 
33745
                    .and.areal(ss(4)).ge.-small) &
 
33746
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
31440
33747
    if (useboxc) then
31441
33748
      call boxc( rslt ,ss,rr ,as ,smax )
31442
33749
    else
31452
33759
                 .or.(     areal(ss(1)).ge.-small  &
31453
33760
                      .and.areal(ss(2)).ge.-small  &
31454
33761
                      .and.areal(ss(3)).ge.-small  &
31455
 
!OLD                      .and.areal(ss(4)).ge.-small) )
31456
 
                      .and.areal(ss(4)).ge.-small) & !NEW
31457
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
33762
                      .and.areal(ss(4)).ge.-small) &
 
33763
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
31458
33764
      if (useboxc) then
31459
33765
        call boxc( rslt ,ss,rr ,as ,smax )
31460
33766
      else
31711
34017
               .or.(     areal(ss(1)).ge.-small  &
31712
34018
                    .and.areal(ss(2)).ge.-small  &
31713
34019
                    .and.areal(ss(3)).ge.-small  &
31714
 
                    .and.areal(ss(4)).ge.-small) )
 
34020
                    .and.areal(ss(4)).ge.-small) &
 
34021
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
31715
34022
    if (useboxc) then
31716
34023
      call boxc( rslt ,ss,rr ,as ,smax )
31717
34024
    else
31727
34034
                 .or.(     areal(ss(1)).ge.-small  &
31728
34035
                      .and.areal(ss(2)).ge.-small  &
31729
34036
                      .and.areal(ss(3)).ge.-small  &
31730
 
!OLD                      .and.areal(ss(4)).ge.-small) )
31731
 
                      .and.areal(ss(4)).ge.-small) & !NEW
31732
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
34037
                      .and.areal(ss(4)).ge.-small) &
 
34038
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
31733
34039
      if (useboxc) then
31734
34040
        call boxc( rslt ,ss,rr ,as ,smax )
31735
34041
      else
31978
34284
               .or.(     areal(ss(1)).ge.-small  &
31979
34285
                    .and.areal(ss(2)).ge.-small  &
31980
34286
                    .and.areal(ss(3)).ge.-small  &
31981
 
                    .and.areal(ss(4)).ge.-small) )
 
34287
                    .and.areal(ss(4)).ge.-small) &
 
34288
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
31982
34289
    if (useboxc) then
31983
34290
      call boxc( rslt ,ss,rr ,as ,smax )
31984
34291
    else
31994
34301
                 .or.(     areal(ss(1)).ge.-small  &
31995
34302
                      .and.areal(ss(2)).ge.-small  &
31996
34303
                      .and.areal(ss(3)).ge.-small  &
31997
 
!OLD                      .and.areal(ss(4)).ge.-small) )
31998
 
                      .and.areal(ss(4)).ge.-small) & !NEW
31999
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
34304
                      .and.areal(ss(4)).ge.-small) &
 
34305
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
32000
34306
      if (useboxc) then
32001
34307
        call boxc( rslt ,ss,rr ,as ,smax )
32002
34308
      else
32247
34553
               .or.(     areal(ss(1)).ge.-small  &
32248
34554
                    .and.areal(ss(2)).ge.-small  &
32249
34555
                    .and.areal(ss(3)).ge.-small  &
32250
 
                    .and.areal(ss(4)).ge.-small) )
 
34556
                    .and.areal(ss(4)).ge.-small) &
 
34557
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
32251
34558
    if (useboxc) then
32252
34559
      call boxc( rslt ,ss,rr ,as ,smax )
32253
34560
    else
32263
34570
                 .or.(     areal(ss(1)).ge.-small  &
32264
34571
                      .and.areal(ss(2)).ge.-small  &
32265
34572
                      .and.areal(ss(3)).ge.-small  &
32266
 
!OLD                      .and.areal(ss(4)).ge.-small) )
32267
 
                      .and.areal(ss(4)).ge.-small) & !NEW
32268
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
34573
                      .and.areal(ss(4)).ge.-small) &
 
34574
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
32269
34575
      if (useboxc) then
32270
34576
        call boxc( rslt ,ss,rr ,as ,smax )
32271
34577
      else
32507
34813
               .or.(     areal(ss(1)).ge.-small  &
32508
34814
                    .and.areal(ss(2)).ge.-small  &
32509
34815
                    .and.areal(ss(3)).ge.-small  &
32510
 
                    .and.areal(ss(4)).ge.-small) )
 
34816
                    .and.areal(ss(4)).ge.-small) &
 
34817
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
32511
34818
    if (useboxc) then
32512
34819
      call boxc( rslt ,ss,rr ,as ,smax )
32513
34820
    else
32523
34830
                 .or.(     areal(ss(1)).ge.-small  &
32524
34831
                      .and.areal(ss(2)).ge.-small  &
32525
34832
                      .and.areal(ss(3)).ge.-small  &
32526
 
!OLD                      .and.areal(ss(4)).ge.-small) )
32527
 
                      .and.areal(ss(4)).ge.-small) & !NEW
32528
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
34833
                      .and.areal(ss(4)).ge.-small) &
 
34834
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
32529
34835
      if (useboxc) then
32530
34836
        call boxc( rslt ,ss,rr ,as ,smax )
32531
34837
      else
32769
35075
               .or.(     areal(ss(1)).ge.-small  &
32770
35076
                    .and.areal(ss(2)).ge.-small  &
32771
35077
                    .and.areal(ss(3)).ge.-small  &
32772
 
                    .and.areal(ss(4)).ge.-small) )
 
35078
                    .and.areal(ss(4)).ge.-small) &
 
35079
               .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
32773
35080
    if (useboxc) then
32774
35081
      call boxc( rslt ,ss,rr ,as ,smax )
32775
35082
    else
32785
35092
                 .or.(     areal(ss(1)).ge.-small  &
32786
35093
                      .and.areal(ss(2)).ge.-small  &
32787
35094
                      .and.areal(ss(3)).ge.-small  &
32788
 
!OLD                      .and.areal(ss(4)).ge.-small) )
32789
 
                      .and.areal(ss(4)).ge.-small) & !NEW
32790
 
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small)) !NEW
 
35095
                      .and.areal(ss(4)).ge.-small) &
 
35096
                 .or.(areal(ss(5)).ge.-small.and.areal(ss(6)).ge.-small))
32791
35097
      if (useboxc) then
32792
35098
        call boxc( rslt ,ss,rr ,as ,smax )
32793
35099
      else
32891
35197
    ,olo_dp_scale=>olo_get_scale &
32892
35198
    ,olo_dp_onshell=>olo_get_onshell &
32893
35199
    ,olo_dp_precision=>olo_get_precision &
32894
 
    ,olo,olo_a0,olo_an,olo_b0,olo_b11,olo_bn,olo_c0,olo_d0
 
35200
    ,olo,olo_a0,olo_an,olo_b0,olo_db0,olo_b11,olo_bn,olo_c0,olo_d0
32895
35201
  use avh_olo_qp ,only: &
32896
35202
     olo_qp_kind=>olo_kind &
32897
35203
    ,olo_qp_scale=>olo_get_scale &
32898
35204
    ,olo_qp_onshell=>olo_get_onshell &
32899
35205
    ,olo_qp_precision=>olo_get_precision &
32900
 
    ,olo,olo_a0,olo_an,olo_b0,olo_b11,olo_bn,olo_c0,olo_d0
 
35206
    ,olo,olo_a0,olo_an,olo_b0,olo_db0,olo_b11,olo_bn,olo_c0,olo_d0
32901
35207
  use avh_olo_mp ,only: &
32902
35208
     olo_mp_kind=>olo_kind &
32903
35209
    ,olo_mp_scale=>olo_get_scale &
32904
35210
    ,olo_mp_onshell=>olo_get_onshell &
32905
35211
    ,olo_mp_precision=>olo_get_precision &
32906
 
    ,olo,olo_a0,olo_an,olo_b0,olo_b11,olo_bn,olo_c0,olo_d0
 
35212
    ,olo,olo_a0,olo_an,olo_b0,olo_db0,olo_b11,olo_bn,olo_c0,olo_d0
32907
35213
 
32908
35214
  implicit none
32909
35215