JHUGen MELA  v2.4.1
Matrix element calculations as used in JHUGen. MELA is an important tool that was used for the Higgs boson discovery and for precise measurements of its structure and interactions. Please see the website https://spin.pha.jhu.edu/ and papers cited there for more details, and kindly cite those papers when using this code.
Functions/Subroutines
hto_sp_fun Module Reference

Functions/Subroutines

recursive real *8 function, dimension(6, 2) hto_s_niels_up4 (x)
 
real *8 function, dimension(2) hto_li2_srsz (x, y, unit)
 
real *8 function, dimension(2) hto_li3_srsz (x, y, unit)
 
real *8 function, dimension(3, 2) hto_poly_unit (theta)
 
real *8 function, dimension(2) hto_li2 (x)
 
real *8 function, dimension(size(x)) hto_li3 (x)
 
subroutine hto_init_niels
 
real *8 function, dimension(2) hto_cqlnx (arg)
 
real *8 function, dimension(2) hto_cqlnomx (arg, omarg)
 

Function/Subroutine Documentation

◆ hto_cqlnomx()

real*8 function, dimension(2) hto_sp_fun::hto_cqlnomx ( real*8, dimension(2), intent(in)  arg,
real*8, dimension(2), intent(in)  omarg 
)

Definition at line 1970 of file CALLING_cpHTO.f.

1970 *
1971  IMPLICIT NONE
1972 *
1973  INTEGER n,k
1974  real*8 zr,zi,zm2,zm
1975  real*8, dimension(2) :: arg,omarg,res,ares
1976  real*8, dimension(10) :: ct,sn
1977  INTENT(IN) arg,omarg
1978 *
1979  zr= arg(1)
1980  zi= arg(2)
1981  zm2= zr*zr+zi*zi
1982  zm= sqrt(zm2)
1983  IF(zm < 1.d-7) THEN
1984  ct(1)= zr/zm
1985  sn(1)= zi/zm
1986  DO n=2,10
1987  ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1)
1988  sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1)
1989  ENDDO
1990  ares(1)= ct(10)/10.d0
1991  ares(2)= sn(10)/10.d0
1992  DO k=9,1,-1
1993  ares(1)= ares(1)*zm+ct(k)/k
1994  ares(2)= ares(2)*zm+sn(k)/k
1995  ENDDO
1996  ares(1)= -ares(1)*zm
1997  ares(2)= -ares(2)*zm
1998  ELSE
1999  ares= hto_cqlnx(omarg)
2000  ENDIF
2001  res= ares
2002 *
2003  RETURN
2004 *

◆ hto_cqlnx()

real*8 function, dimension(2) hto_sp_fun::hto_cqlnx ( real*8, dimension(2), intent(in)  arg)

Definition at line 1912 of file CALLING_cpHTO.f.

1912  USE hto_riemann
1913  USE hto_units
1914  IMPLICIT NONE
1915 *
1916  real*8 teta,zm,zm2,tnteta,sr,si
1917  real*8, dimension(2) :: arg,aarg,res
1918  INTENT(IN) arg
1919 *
1920  IF((abs(arg(2))-eps).eq.0.d0) THEN
1921  res(1)= log(abs(arg(1)))
1922  IF(arg(1) > 0.d0) THEN
1923  res(2)= 0.d0
1924  ELSE
1925  res(2)= pi*sign(one,arg(2))
1926  ENDIF
1927  RETURN
1928  ENDIF
1929 *
1930  aarg= abs(arg)
1931  zm2= (arg(1))**2+(arg(2))**2
1932  zm= sqrt(zm2)
1933  res(1)= log(zm)
1934  IF(arg(1).eq.0.d0) THEN
1935  IF(arg(2) > 0.d0) THEN
1936  teta= pi/2.d0
1937  ELSE
1938  teta= -pi/2.d0
1939  ENDIF
1940  res(2)= teta
1941  RETURN
1942  ELSE IF(arg(2).eq.0.d0) THEN
1943  IF(arg(1) > 0.d0) THEN
1944  teta= 0.d0
1945  ELSE
1946  teta= pi
1947  ENDIF
1948  res(2)= teta
1949  RETURN
1950  ELSE
1951  tnteta= aarg(2)/aarg(1)
1952  teta= atan(tnteta)
1953  sr= arg(1)/aarg(1)
1954  si= arg(2)/aarg(2)
1955  IF(sr > 0.d0) THEN
1956  res(2)= si*teta
1957  ELSE
1958  res(2)= si*(pi-teta)
1959  ENDIF
1960  RETURN
1961  ENDIF
1962 *

◆ hto_init_niels()

subroutine hto_sp_fun::hto_init_niels

Definition at line 1822 of file CALLING_cpHTO.f.

1822  USE hto_common_niels
1823 *
1824  plr= 0.d0
1825 *
1826  plr(1,1)= 1.d0
1827  plr(1,2)= -2.5d-1
1828  plr(1,3)= 2.77777777777777777d-2
1829  plr(1,5)= -2.77777777777777777d-4
1830  plr(1,7)= 4.72411186696900982d-6
1831  plr(1,9)= -9.18577307466196355d-8
1832  plr(1,11)= 1.89788699889709990d-9
1833  plr(1,13)= -4.06476164514422552d-11
1834  plr(1,15)= 8.92169102045645255d-13
1835 *
1836  plr(2,1)= 1.d0
1837  plr(2,2)= -3.75d-1
1838  plr(2,3)= 7.87037037037037037d-2
1839  plr(2,4)= -8.68055555555555555d-3
1840  plr(2,5)= 1.29629629629629629d-4
1841  plr(2,6)= 8.10185185185185185d-5
1842  plr(2,7)= -3.41935716085375949d-6
1843  plr(2,8)= -1.32865646258503401d-6
1844  plr(2,9)= 8.66087175610985134d-8
1845  plr(2,10)= 2.52608759553203997d-8
1846  plr(2,11)= -2.14469446836406476d-9
1847  plr(2,12)= -5.14011062201297891d-10
1848  plr(2,13)= 5.24958211460082943d-11
1849  plr(2,14)= 1.08877544066363183d-11
1850  plr(2,15)= -1.27793960944936953d-12
1851 *
1852  plr(3,2)= 2.5d-1
1853  plr(3,3)= -8.33333333333333333d-2
1854  plr(3,4)= 1.04166666666666666d-2
1855  plr(3,6)= -1.15740740740740740d-4
1856  plr(3,8)= 2.06679894179894179d-6
1857  plr(3,10)= -4.13359788359788359d-8
1858  plr(3,12)= 8.69864874494504124d-10
1859  plr(3,14)= -1.88721076381696185d-11
1860 *
1861  plr_4= 0.d0
1862 *
1863  plr_4(1,1)= 1.d0
1864  plr_4(1,2)= -4.375d-1
1865  plr_4(1,3)= 1.16512345679012345d-1
1866  plr_4(1,4)= -1.98206018518518518d-2
1867  plr_4(1,5)= 1.92793209876543209d-3
1868  plr_4(1,6)= -3.10570987654320987d-5
1869  plr_4(1,7)= -1.56240091148578352d-5
1870  plr_4(1,8)= 8.48512354677320663d-7
1871  plr_4(1,9)= 2.29096166031897114d-7
1872  plr_4(1,10)= -2.18326142185269169d-8
1873  plr_4(1,11)= -3.88282487917201557d-9
1874  plr_4(1,12)= 5.44629210322033211d-10
1875  plr_4(1,13)= 6.96080521068272540d-11
1876  plr_4(1,14)= -1.33757376864452151d-11
1877  plr_4(1,15)= -1.27848526852665716d-12
1878 *
1879  plr_4(2,2)= 2.5d-1
1880  plr_4(2,3)= -1.25d-1
1881  plr_4(2,4)= 2.95138888888888888d-2
1882  plr_4(2,5)= -3.47222222222222222d-3
1883  plr_4(2,6)= 5.40123456790123456d-5
1884  plr_4(2,7)= 3.47222222222222222d-5
1885  plr_4(2,8)= -1.49596875787351977d-6
1886  plr_4(2,9)= -5.90513983371126228d-7
1887  plr_4(2,10)= 3.89739229024943310d-8
1888  plr_4(2,11)= 1.14822163433274544d-8
1889  plr_4(2,12)= -9.82984964666863015d-10
1890  plr_4(2,13)= -2.37235874862137488d-10
1891  plr_4(2,14)= 2.43730598177895652d-11
1892  plr_4(2,15)= 5.08095205643028190d-12
1893 *
1894  plr_4(3,2)= 5.55555555555555555d-2
1895  plr_4(3,3)= -2.08333333333333333d-2
1896  plr_4(3,4)= 2.77777777777777777d-3
1897  plr_4(3,6)= -3.30687830687830687d-5
1898  plr_4(3,8)= 6.12384871644130903d-7
1899  plr_4(3,10)= -1.25260541927208593d-8
1900  plr_4(3,12)= 2.67650730613693576d-10
1901  plr_4(3,14)= -5.87132237631943687d-12
1902 *
1903  RETURN
1904 *

◆ hto_li2()

real*8 function, dimension(2) hto_sp_fun::hto_li2 ( real*8, dimension(2), intent(in)  x)

Definition at line 1302 of file CALLING_cpHTO.f.

1302  USE hto_acmplx_pro
1303  USE hto_acmplx_rat
1304  USE hto_full_ln
1305  USE hto_bernoulli
1306  USE hto_riemann
1307  USE hto_units
1308 *
1309  IMPLICIT NONE
1310 *
1311  real*8, dimension(2) ::x,value
1312  INTENT(IN) x
1313 *
1314  INTEGER n
1315  real*8 ym,ym2,sign1,sign2,sign3,fact, parr,pari,parm,parm2,zr
1316  real*8, dimension(2) :: clnx,clnomx,clnoy,clnz,clnomz,
1317  # add1,add2,add3,par,res
1318  real*8, dimension(0:14) :: bf
1319  real*8, dimension(15) :: ct,sn
1320  real*8, dimension(2) :: omx,y,oy,omy,z,omz,t,omt
1321 *
1322  omx= co-x
1323  IF(x(1) < 0.d0) THEN
1324  y= omx
1325  sign1= -1.d0
1326  clnx= x(1).fln.x(2)
1327  clnomx= omx(1).fln.omx(2)
1328  add1= pis/6.d0*co-(clnx.cp.clnomx)
1329  ELSE
1330  y= x
1331  sign1= 1.d0
1332  add1= 0.d0
1333  ENDIF
1334  omy= co-y
1335  ym2= y(1)*y(1)+y(2)*y(2)
1336  ym= sqrt(ym2)
1337  IF(ym > 1.d0) THEN
1338  z(1)= y(1)/ym2
1339  z(2)= -y(2)/ym2
1340  sign2= -1.d0
1341  oy= -y
1342  clnoy= oy(1).fln.oy(2)
1343  add2= -pis/6.d0*co-0.5d0*(clnoy.cp.clnoy)
1344  ELSE
1345  z= y
1346  sign2= 1.d0
1347  add2= 0.d0
1348  ENDIF
1349  omz= co-z
1350  zr= z(1)
1351  IF(zr > 0.5d0) THEN
1352  t= co-z
1353  omt= co-t
1354  sign3= -1.d0
1355  clnz= z(1).fln.z(2)
1356  clnomz= omz(1).fln.omz(2)
1357  add3= pis/6.d0*co-(clnz.cp.clnomz)
1358  ELSE
1359  t= z
1360  omt= co-t
1361  sign3= 1.d0
1362  add3= 0.d0
1363  ENDIF
1364  par= omt(1).fln.omt(2)
1365  fact= 1.d0
1366  DO n=0,14
1367  bf(n)= b_num(n)/fact
1368  fact= fact*(n+2.d0)
1369  ENDDO
1370  parr= par(1)
1371  pari= par(2)
1372  parm2= parr*parr+pari*pari
1373  parm= sqrt(parm2)
1374  ct(1)= parr/parm
1375  sn(1)= pari/parm
1376  DO n=2,15
1377  ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1)
1378  sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1)
1379  ENDDO
1380 *
1381  res(1)= -((((((((bf(14)*ct(15)*parm2+bf(12)*ct(13))*parm2+
1382  # bf(10)*ct(11))*parm2+bf(8)*ct(9))*parm2+
1383  # bf(6)*ct(7))*parm2+bf(4)*ct(5))*parm2+
1384  # bf(2)*ct(3))*(-parm)+bf(1)*ct(2))*(-parm)+
1385  # bf(0)*ct(1))*parm
1386  res(2)= -((((((((bf(14)*sn(15)*parm2+bf(12)*sn(13))*parm2+
1387  # bf(10)*sn(11))*parm2+bf(8)*sn(9))*parm2+
1388  # bf(6)*sn(7))*parm2+bf(4)*sn(5))*parm2+
1389  # bf(2)*sn(3))*(-parm)+bf(1)*sn(2))*(-parm)+
1390  # bf(0)*sn(1))*parm
1391 *
1392  value= sign1*(sign2*(sign3*res+add3)+add2)+add1
1393 *
1394  RETURN
1395 *

◆ hto_li2_srsz()

real*8 function, dimension(2) hto_sp_fun::hto_li2_srsz ( real*8, dimension (:), intent(in)  x,
real*8, dimension (:), intent(in)  y,
integer  unit 
)

Definition at line 1065 of file CALLING_cpHTO.f.

1065  USE hto_acmplx_pro
1066  USE hto_full_ln
1067  USE hto_riemann
1068  USE hto_units
1069  USE hto_kountac
1070 *
1071  IMPLICIT NONE
1072 *
1073  INTEGER unit
1074  real*8 xms,atheta,theta,sr,si
1075  real*8, dimension (:) :: x,y
1076  real*8, dimension(3,2) :: aux
1077  real*8, dimension(2) :: lnc,add,value
1078  INTENT(IN) x,y
1079 *
1080  IF(abs(y(2)).ne.1.d0) THEN
1081  IF(unit.eq.1) THEN
1082  xms= x(1)*x(1)+x(2)*x(2)
1083  IF(abs(1.d0-sqrt(xms)).gt.1.d-12) THEN
1084  print*,' apparent inconsistency '
1085  ENDIF
1086  atheta= atan(abs(x(2)/x(1)))
1087  sr= x(1)/abs(x(1))
1088  si= x(2)/abs(x(2))
1089  IF(sr > 0.d0.and.si > 0.d0) THEN
1090  theta= atheta
1091  ELSEIF(sr > 0.d0.and.si < 0.d0) THEN
1092  theta= 2.d0*pi-atheta
1093  ELSEIF(sr < 0.d0.and.si > 0.d0) THEN
1094  theta= pi-atheta
1095  ELSEIF(sr < 0.d0.and.si < 0.d0) THEN
1096  theta= pi+atheta
1097  ENDIF
1098  aux= hto_poly_unit(theta)
1099  value(1:2)= aux(1,1:2)
1100  ELSE
1101  value= hto_li2(x)
1102  ENDIF
1103  RETURN
1104  ENDIF
1105 *
1106  IF(x(1) > 1.d0) THEN
1107  IF(y(1) < 1.d0) THEN
1108  value= hto_li2(x)
1109  print*,'+++++++++++++++++++++++++++++++++++'
1110  print*,' anomaly Li2 '
1111  print*,x
1112  print*,y
1113  print*,'+++++++++++++++++++++++++++++++++++'
1114  ENDIF
1115  lnc= x(1).fln.x(2)
1116  add(1)= -lnc(2)
1117  add(2)= lnc(1)
1118  IF(y(2) < 0.d0.and.x(2) > 0.d0) THEN
1119  value= hto_li2(x)-2.d0*pi*add
1120  km= km+1
1121  ELSEIF(y(2) > 0.d0.and.x(2) < 0.d0) THEN
1122  value= hto_li2(x)+2.d0*pi*add
1123  kp= kp+1
1124  ELSE
1125  value= hto_li2(x)
1126  ENDIF
1127  ELSE
1128  value= hto_li2(x)
1129  ENDIF
1130 *
1131  RETURN
1132 *

◆ hto_li3()

real*8 function, dimension(size(x)) hto_sp_fun::hto_li3 ( real*8, dimension (:), intent(in)  x)

Definition at line 1401 of file CALLING_cpHTO.f.

1401  USE hto_riemann
1402  USE hto_units
1403 *
1404  IMPLICIT NONE
1405 *
1406  real*8, dimension (:) :: x
1407  real*8 hto_li3(size(x))
1408  INTENT(IN) x
1409  INTEGER n
1410  real*8 xm,xm2,xtst,pr,pj,p2,pm,pr1,pj1,pm1,p12,pr2,pj2,p22,pm2,
1411  # rln2_x,iln2_x,tm2,y2r,ym2,ytst
1412  real*8, dimension(0:14) :: bf
1413  real*8, dimension(15) :: ct,sn,ct1,sn1,ct2,sn2
1414  real*8, dimension(2) :: y,addx,ox,clnx,par,res,u1,u2,ln_y,omy,
1415  # ln_omy,addy,par1,par2,res1,res2,t,resa,resb,ln_t,res3,res4,
1416  # ln_omt,addt,addt2,omt,omu1,omu2
1417  real*8 :: b(0:14)=(/1.d0,-0.75d0,
1418  # 0.236111111111111111111111111111111d0,
1419  # -3.472222222222222222222222222222222d-2,
1420  # 6.481481481481481481481481481481482d-4,
1421  # 4.861111111111111111111111111111111d-4,
1422  # -2.393550012597631645250692869740488d-5,
1423  # -1.062925170068027210884353741496599d-5,
1424  # 7.794784580498866213151927437641723d-7,
1425  # 2.526087595532039976484420928865373d-7,
1426  # -2.359163915200471237027273583310139d-8,
1427  # -6.168132746415574698402981231264060d-9,
1428  # 6.824456748981078267312315451125495d-10,
1429  # 1.524285616929084572552216019859487d-10,
1430  # -1.916909414174054295837274763110831d-11/)
1431 *
1432  FORALL (n=0:14) bf(n)= b(n)/(n+1.d0)
1433 *
1434  xm2= x(1)*x(1)+x(2)*x(2)
1435  xm= sqrt(xm2)
1436 *
1437 *-----the modulus of x is checked
1438 *
1439  xtst= xm-1.d0
1440  IF(xtst <= 1.d-20) THEN
1441  y= x
1442  addx= 0.d0
1443  ELSE IF(xm > 1.d-20) THEN
1444  y(1)= x(1)/xm2
1445  y(2)= -x(2)/xm2
1446  ox= -x
1447  clnx= hto_cqlnx(ox)
1448  rln2_x= clnx(1)*clnx(1)
1449  iln2_x= clnx(2)*clnx(2)
1450  addx(1)= -clnx(1)*(rz2+1.d0/6.d0*(rln2_x-3.d0*iln2_x))
1451  addx(2)= -clnx(2)*(rz2+1.d0/6.d0*(3.d0*rln2_x-iln2_x))
1452  ENDIF
1453 *
1454 *-----once x --> y, |y|<1 the sign of re(y) is checked
1455 * if re(y)>0 a transformation is required for re(y)>1/2
1456 *
1457  y2r= y(1)*y(1)-y(2)*y(2)
1458  IF(y(1) >= 0.d0.or.y2r < 0.d0) THEN
1459  ytst= y(1)-0.5d0
1460  IF(ytst <= 1.d-20) THEN
1461 *
1462 *-----li_3(y) is computed
1463 *
1464  omy= co-y
1465  par= hto_cqlnomx(y,omy)
1466  pr= -par(1)
1467  pj= -par(2)
1468  p2= pr*pr+pj*pj
1469  pm= sqrt(p2)
1470  ct(1)= pr/pm
1471  sn(1)= pj/pm
1472  DO n=2,15
1473  ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1)
1474  sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1)
1475  ENDDO
1476  res(1)= pm*(bf(0)*ct(1)+pm*(bf(1)*ct(2)+pm*
1477  # (bf(2)*ct(3)+pm*(bf(3)*ct(4)+pm*
1478  # (bf(4)*ct(5)+pm*(bf(5)*ct(6)+pm*
1479  # (bf(6)*ct(7)+pm*(bf(7)*ct(8)+pm*
1480  # (bf(8)*ct(9)+pm*(bf(9)*ct(10)+pm*
1481  # (bf(10)*ct(11)+pm*(bf(11)*ct(12)+pm*
1482  # (bf(12)*ct(13)+pm*(bf(13)*ct(14)+pm*
1483  # (bf(14)*ct(15))))))))))))))))
1484  res(2)= pm*(bf(0)*sn(1)+pm*(bf(1)*sn(2)+pm*
1485  # (bf(2)*sn(3)+pm*(bf(3)*sn(4)+pm*
1486  # (bf(4)*sn(5)+pm*(bf(5)*sn(6)+pm*
1487  # (bf(6)*sn(7)+pm*(bf(7)*sn(8)+pm*
1488  # (bf(8)*sn(9)+pm*(bf(9)*sn(10)+pm*
1489  # (bf(10)*sn(11)+pm*(bf(11)*sn(12)+pm*
1490  # (bf(12)*sn(13)+pm*(bf(13)*sn(14)+pm*
1491  # (bf(14)*sn(15))))))))))))))))
1492  hto_li3= res+addx
1493  RETURN
1494  ELSE IF(ytst > 1.d-20) THEN
1495  ym2= y(1)*y(1)+y(2)*y(2)
1496  u1(1)= 1.d0-y(1)/ym2
1497  u1(2)= y(2)/ym2
1498  u2= co-y
1499  ln_y= hto_cqlnx(y)
1500  omy= co-y
1501  ln_omy= hto_cqlnomx(y,omy)
1502  addy(1)= rz3+rz2*ln_y(1)+1.d0/6.d0*ln_y(1)*
1503  # (ln_y(1)*ln_y(1)-3.d0*ln_y(2)*ln_y(2))-
1504  # 0.5d0*ln_omy(1)*(ln_y(1)*ln_y(1)-ln_y(2)*
1505  # ln_y(2))+ln_y(1)*ln_y(2)*ln_omy(2)
1506  addy(2)= rz2*ln_y(2)+1.d0/6.d0*ln_y(2)*(3.d0*
1507  # ln_y(1)*ln_y(1)-ln_y(2)*ln_y(2))-0.5d0*
1508  # ln_omy(2)*(ln_y(1)*ln_y(1)-ln_y(2)*ln_y(2))-
1509  # ln_y(1)*ln_omy(1)*ln_y(2)
1510 *
1511 *-----li_3(1-1/y) is computed
1512 *
1513  omu1= co-u1
1514  par1= hto_cqlnomx(u1,omu1)
1515  pr1= -par1(1)
1516  pj1= -par1(2)
1517  p12= pr1*pr1+pj1*pj1
1518  pm1= sqrt(p12)
1519  ct1(1)= pr1/pm1
1520  sn1(1)= pj1/pm1
1521  DO n=2,15
1522  ct1(n)= ct1(1)*ct1(n-1)-sn1(1)*sn1(n-1)
1523  sn1(n)= sn1(1)*ct1(n-1)+ct1(1)*sn1(n-1)
1524  ENDDO
1525  res1(1)= pm1*(bf(0)*ct1(1)+pm1*(bf(1)*ct1(2)+pm1*
1526  # (bf(2)*ct1(3)+pm1*(bf(3)*ct1(4)+pm1*
1527  # (bf(4)*ct1(5)+pm1*(bf(5)*ct1(6)+pm1*
1528  # (bf(6)*ct1(7)+pm1*(bf(7)*ct1(8)+pm1*
1529  # (bf(8)*ct1(9)+pm1*(bf(9)*ct1(10)+pm1*
1530  # (bf(10)*ct1(11)+pm1*(bf(11)*ct1(12)+pm1*
1531  # (bf(12)*ct1(13)+pm1*(bf(13)*ct1(14)+pm1*
1532  # (bf(14)*ct1(15))))))))))))))))
1533  res1(2)= pm1*(bf(0)*sn1(1)+pm1*(bf(1)*sn1(2)+pm1*
1534  # (bf(2)*sn1(3)+pm1*(bf(3)*sn1(4)+pm1*
1535  # (bf(4)*sn1(5)+pm1*(bf(5)*sn1(6)+pm1*
1536  # (bf(6)*sn1(7)+pm1*(bf(7)*sn1(8)+pm1*
1537  # (bf(8)*sn1(9)+pm1*(bf(9)*sn1(10)+pm1*
1538  # (bf(10)*sn1(11)+pm1*(bf(11)*sn1(12)+pm1*
1539  # (bf(12)*sn1(13)+pm1*(bf(13)*sn1(14)+pm1*
1540  # (bf(14)*sn1(15))))))))))))))))
1541 *
1542 *-----li_3(1-y) is computed
1543 *
1544  omu2= co-u2
1545  par2= hto_cqlnomx(u2,omu2)
1546  pr2= -par2(1)
1547  pj2= -par2(2)
1548  p22= pr2*pr2+pj2*pj2
1549  pm2= sqrt(p22)
1550  ct2(1)= pr2/pm2
1551  sn2(1)= pj2/pm2
1552  DO n=2,15
1553  ct2(n)= ct2(1)*ct2(n-1)-sn2(1)*sn2(n-1)
1554  sn2(n)= sn2(1)*ct2(n-1)+ct2(1)*sn2(n-1)
1555  ENDDO
1556  res2(1)= pm2*(bf(0)*ct2(1)+pm2*(bf(1)*ct2(2)+pm2*
1557  # (bf(2)*ct2(3)+pm2*(bf(3)*ct2(4)+pm2*
1558  # (bf(4)*ct2(5)+pm2*(bf(5)*ct2(6)+pm2*
1559  # (bf(6)*ct2(7)+pm2*(bf(7)*ct2(8)+pm2*
1560  # (bf(8)*ct2(9)+pm2*(bf(9)*ct2(10)+pm2*
1561  # (bf(10)*ct2(11)+pm2*(bf(11)*ct2(12)+pm2*
1562  # (bf(12)*ct2(13)+pm2*(bf(13)*ct2(14)+pm2*
1563  # (bf(14)*ct2(15))))))))))))))))
1564  res2(2)= pm2*(bf(0)*sn2(1)+pm2*(bf(1)*sn2(2)+pm2*
1565  # (bf(2)*sn2(3)+pm2*(bf(3)*sn2(4)+pm2*
1566  # (bf(4)*sn2(5)+pm2*(bf(5)*sn2(6)+pm2*
1567  # (bf(6)*sn2(7)+pm2*(bf(7)*sn2(8)+pm2*
1568  # (bf(8)*sn2(9)+pm2*(bf(9)*sn2(10)+pm2*
1569  # (bf(10)*sn2(11)+pm2*(bf(11)*sn2(12)+pm2*
1570  # (bf(12)*sn2(13)+pm2*(bf(13)*sn2(14)+pm2*
1571  # (bf(14)*sn2(15))))))))))))))))
1572  hto_li3= -res1-res2+addx+addy
1573  RETURN
1574  ENDIF
1575 *
1576 *-----if re(y)<0 a transformation is required in terms of t = -y
1577 * and of t^2
1578 *
1579  ELSE IF(y(1) < 0.d0) THEN
1580 *
1581 *-----first t
1582 *
1583  t= -y
1584  IF(t(1) <= 0.5d0) THEN
1585 *
1586 *-----li_3(t) is computed
1587 *
1588  omt= co-t
1589  par= hto_cqlnomx(t,omt)
1590  pr= -par(1)
1591  pj= -par(2)
1592  p2= pr*pr+pj*pj
1593  pm= sqrt(p2)
1594  ct(1)= pr/pm
1595  sn(1)= pj/pm
1596  DO n=2,15
1597  ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1)
1598  sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1)
1599  ENDDO
1600  resa(1)= pm*(bf(0)*ct(1)+pm*(bf(1)*ct(2)+pm*
1601  # (bf(2)*ct(3)+pm*(bf(3)*ct(4)+pm*
1602  # (bf(4)*ct(5)+pm*(bf(5)*ct(6)+pm*
1603  # (bf(6)*ct(7)+pm*(bf(7)*ct(8)+pm*
1604  # (bf(8)*ct(9)+pm*(bf(9)*ct(10)+pm*
1605  # (bf(10)*ct(11)+pm*(bf(11)*ct(12)+pm*
1606  # (bf(12)*ct(13)+pm*(bf(13)*ct(14)+pm*
1607  # (bf(14)*ct(15))))))))))))))))
1608  resa(2)= pm*(bf(0)*sn(1)+pm*(bf(1)*sn(2)+pm*
1609  # (bf(2)*sn(3)+pm*(bf(3)*sn(4)+pm*
1610  # (bf(4)*sn(5)+pm*(bf(5)*sn(6)+pm*
1611  # (bf(6)*sn(7)+pm*(bf(7)*sn(8)+pm*
1612  # (bf(8)*sn(9)+pm*(bf(9)*sn(10)+pm*
1613  # (bf(10)*sn(11)+pm*(bf(11)*sn(12)+pm*
1614  # (bf(12)*sn(13)+pm*(bf(13)*sn(14)+pm*
1615  # (bf(14)*sn(15))))))))))))))))
1616  ELSE IF(t(1) > 0.5d0) THEN
1617  tm2= t(1)*t(1)+t(2)*t(2)
1618  u1(1)= 1.d0-t(1)/tm2
1619  u1(2)= t(2)/tm2
1620  u2= co-t
1621  ln_t= hto_cqlnx(t)
1622  omt= co-t
1623  ln_omt= hto_cqlnomx(t,omt)
1624  addt(1)= rz3+rz2*ln_t(1)+1.d0/6.d0*ln_t(1)*
1625  # (ln_t(1)*ln_t(1)-3.d0*ln_t(2)*ln_t(2))-
1626  # 0.5d0*ln_omt(1)*(ln_t(1)*ln_t(1)-ln_t(2)*
1627  # ln_t(2))+ln_t(1)*ln_t(2)*ln_omt(2)
1628  addt(2)= rz2*ln_t(2)+1.d0/6.d0*ln_t(2)*(3.d0*
1629  # ln_t(1)*ln_t(1)-ln_t(2)*ln_t(2))-0.5d0*
1630  # ln_omt(2)*(ln_t(1)*ln_t(1)-ln_t(2)*ln_t(2))-
1631  # ln_t(1)*ln_omt(1)*ln_t(2)
1632 *
1633 *-----li3(1-1/t) is computed
1634 *
1635  omu1= co-u1
1636  par1= hto_cqlnomx(u1,omu1)
1637  pr1= -par1(1)
1638  pj1= -par1(2)
1639  p12= pr1*pr1+pj1*pj1
1640  pm1= sqrt(p12)
1641  ct1(1)= pr1/pm1
1642  sn1(1)= pj1/pm1
1643  DO n=2,15
1644  ct1(n)= ct1(1)*ct1(n-1)-sn1(1)*sn1(n-1)
1645  sn1(n)= sn1(1)*ct1(n-1)+ct1(1)*sn1(n-1)
1646  ENDDO
1647  res1(1)= pm1*(bf(0)*ct1(1)+pm1*(bf(1)*ct1(2)+pm1*
1648  # (bf(2)*ct1(3)+pm1*(bf(3)*ct1(4)+pm1*
1649  # (bf(4)*ct1(5)+pm1*(bf(5)*ct1(6)+pm1*
1650  # (bf(6)*ct1(7)+pm1*(bf(7)*ct1(8)+pm1*
1651  # (bf(8)*ct1(9)+pm1*(bf(9)*ct1(10)+pm1*
1652  # (bf(10)*ct1(11)+pm1*(bf(11)*ct1(12)+pm1*
1653  # (bf(12)*ct1(13)+pm1*(bf(13)*ct1(14)+pm1*
1654  # (bf(14)*ct1(15))))))))))))))))
1655  res1(2)= pm1*(bf(0)*sn1(1)+pm1*(bf(1)*sn1(2)+pm1*
1656  # (bf(2)*sn1(3)+pm1*(bf(3)*sn1(4)+pm1*
1657  # (bf(4)*sn1(5)+pm1*(bf(5)*sn1(6)+pm1*
1658  # (bf(6)*sn1(7)+pm1*(bf(7)*sn1(8)+pm1*
1659  # (bf(8)*sn1(9)+pm1*(bf(9)*sn1(10)+pm1*
1660  # (bf(10)*sn1(11)+pm1*(bf(11)*sn1(12)+pm1*
1661  # (bf(12)*sn1(13)+pm1*(bf(13)*sn1(14)+pm1*
1662  # (bf(14)*sn1(15))))))))))))))))
1663 *
1664 *-----li3(1-t) is computed
1665 *
1666  omu2= co-u2
1667  par2= hto_cqlnomx(u2,omu2)
1668  pr2= -par2(1)
1669  pj2= -par2(2)
1670  p22= pr2*pr2+pj2*pj2
1671  pm2= sqrt(p22)
1672  ct2(1)= pr2/pm2
1673  sn2(1)= pj2/pm2
1674  DO n=2,15
1675  ct2(n)= ct2(1)*ct2(n-1)-sn2(1)*sn2(n-1)
1676  sn2(n)= sn2(1)*ct2(n-1)+ct2(1)*sn2(n-1)
1677  ENDDO
1678  res2(1)= pm2*(bf(0)*ct2(1)+pm2*(bf(1)*ct2(2)+pm2*
1679  # (bf(2)*ct2(3)+pm2*(bf(3)*ct2(4)+pm2*
1680  # (bf(4)*ct2(5)+pm2*(bf(5)*ct2(6)+pm2*
1681  # (bf(6)*ct2(7)+pm2*(bf(7)*ct2(8)+pm2*
1682  # (bf(8)*ct2(9)+pm2*(bf(9)*ct2(10)+pm2*
1683  # (bf(10)*ct2(11)+pm2*(bf(11)*ct2(12)+pm2*
1684  # (bf(12)*ct2(13)+pm2*(bf(13)*ct2(14)+pm2*
1685  # (bf(14)*ct2(15))))))))))))))))
1686  res2(2)= pm2*(bf(0)*sn2(1)+pm2*(bf(1)*sn2(2)+pm2*
1687  # (bf(2)*sn2(3)+pm2*(bf(3)*sn2(4)+pm2*
1688  # (bf(4)*sn2(5)+pm2*(bf(5)*sn2(6)+pm2*
1689  # (bf(6)*sn2(7)+pm2*(bf(7)*sn2(8)+pm2*
1690  # (bf(8)*sn2(9)+pm2*(bf(9)*sn2(10)+pm2*
1691  # (bf(10)*sn2(11)+pm2*(bf(11)*sn2(12)+pm2*
1692  # (bf(12)*sn2(13)+pm2*(bf(13)*sn2(14)+pm2*
1693  # (bf(14)*sn2(15))))))))))))))))
1694  resa= -res1-res2+addt
1695  ENDIF
1696 *
1697 *-----THEN t^2
1698 *
1699  t(1)= y(1)*y(1)-y(2)*y(2)
1700  t(2)= 2.d0*y(1)*y(2)
1701  IF(t(1) <= 0.5d0) THEN
1702 *
1703 *-----li_3(t^2) is computed
1704 *
1705  omt= co-t
1706  par= hto_cqlnomx(t,omt)
1707  pr= -par(1)
1708  pj= -par(2)
1709  p2= pr*pr+pj*pj
1710  pm= sqrt(p2)
1711  ct(1)= pr/pm
1712  sn(1)= pj/pm
1713  DO n=2,15
1714  ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1)
1715  sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1)
1716  ENDDO
1717  resb(1)= pm*(bf(0)*ct(1)+pm*(bf(1)*ct(2)+pm*
1718  # (bf(2)*ct(3)+pm*(bf(3)*ct(4)+pm*
1719  # (bf(4)*ct(5)+pm*(bf(5)*ct(6)+pm*
1720  # (bf(6)*ct(7)+pm*(bf(7)*ct(8)+pm*
1721  # (bf(8)*ct(9)+pm*(bf(9)*ct(10)+pm*
1722  # (bf(10)*ct(11)+pm*(bf(11)*ct(12)+pm*
1723  # (bf(12)*ct(13)+pm*(bf(13)*ct(14)+pm*
1724  # (bf(14)*ct(15))))))))))))))))
1725  resb(2)= pm*(bf(0)*sn(1)+pm*(bf(1)*sn(2)+pm*
1726  # (bf(2)*sn(3)+pm*(bf(3)*sn(4)+pm*
1727  # (bf(4)*sn(5)+pm*(bf(5)*sn(6)+pm*
1728  # (bf(6)*sn(7)+pm*(bf(7)*sn(8)+pm*
1729  # (bf(8)*sn(9)+pm*(bf(9)*sn(10)+pm*
1730  # (bf(10)*sn(11)+pm*(bf(11)*sn(12)+pm*
1731  # (bf(12)*sn(13)+pm*(bf(13)*sn(14)+pm*
1732  # (bf(14)*sn(15))))))))))))))))
1733  ELSE IF(t(1) > 0.5d0) THEN
1734  tm2= t(1)*t(1)+t(2)*t(2)
1735  u1(1)= 1.d0-t(1)/tm2
1736  u1(2)= t(2)/tm2
1737  u2= co-t
1738  ln_t= hto_cqlnx(t)
1739  omt= co-t
1740  ln_omt= hto_cqlnomx(t,omt)
1741  addt2(1)= rz3+rz2*ln_t(1)+1.d0/6.d0*ln_t(1)*
1742  # (ln_t(1)*ln_t(1)-3.d0*ln_t(2)*ln_t(2))-
1743  # 0.5d0*ln_omt(1)*(ln_t(1)*ln_t(1)-ln_t(2)*
1744  # ln_t(2))+ln_t(1)*ln_t(2)*ln_omt(2)
1745  addt2(2)= rz2*ln_t(2)+1.d0/6.d0*ln_t(2)*(3.d0*
1746  # ln_t(1)*ln_t(1)-ln_t(2)*ln_t(2))-0.5d0*
1747  # ln_omt(2)*(ln_t(1)*ln_t(1)-ln_t(2)*ln_t(2))-
1748  # ln_t(1)*ln_omt(1)*ln_t(2)
1749 *
1750 *-----li_3(1-1/t^2) is computed
1751 *
1752  omu1= co-u1
1753  par1= hto_cqlnomx(u1,omu1)
1754  pr1= -par1(1)
1755  pj1= -par1(2)
1756  p12= pr1*pr1+pj1*pj1
1757  pm1= sqrt(p12)
1758  ct1(1)= pr1/pm1
1759  sn1(1)= pj1/pm1
1760  DO n=2,15
1761  ct1(n)= ct1(1)*ct1(n-1)-sn1(1)*sn1(n-1)
1762  sn1(n)= sn1(1)*ct1(n-1)+ct1(1)*sn1(n-1)
1763  ENDDO
1764  res3(1)= pm1*(bf(0)*ct1(1)+pm1*(bf(1)*ct1(2)+pm1*
1765  # (bf(2)*ct1(3)+pm1*(bf(3)*ct1(4)+pm1*
1766  # (bf(4)*ct1(5)+pm1*(bf(5)*ct1(6)+pm1*
1767  # (bf(6)*ct1(7)+pm1*(bf(7)*ct1(8)+pm1*
1768  # (bf(8)*ct1(9)+pm1*(bf(9)*ct1(10)+pm1*
1769  # (bf(10)*ct1(11)+pm1*(bf(11)*ct1(12)+pm1*
1770  # (bf(12)*ct1(13)+pm1*(bf(13)*ct1(14)+pm1*
1771  # (bf(14)*ct1(15))))))))))))))))
1772  res3(2)= pm1*(bf(0)*sn1(1)+pm1*(bf(1)*sn1(2)+pm1*
1773  # (bf(2)*sn1(3)+pm1*(bf(3)*sn1(4)+pm1*
1774  # (bf(4)*sn1(5)+pm1*(bf(5)*sn1(6)+pm1*
1775  # (bf(6)*sn1(7)+pm1*(bf(7)*sn1(8)+pm1*
1776  # (bf(8)*sn1(9)+pm1*(bf(9)*sn1(10)+pm1*
1777  # (bf(10)*sn1(11)+pm1*(bf(11)*sn1(12)+pm1*
1778  # (bf(12)*sn1(13)+pm1*(bf(13)*sn1(14)+pm1*
1779  # (bf(14)*sn1(15))))))))))))))))
1780 *
1781 *-----li_3(1-t^2) is computed
1782 *
1783  omu2= co-u2
1784  par2= hto_cqlnomx(u2,omu2)
1785  pr2= -par2(1)
1786  pj2= -par2(2)
1787  p22= pr2*pr2+pj2*pj2
1788  pm2= sqrt(p22)
1789  ct2(1)= pr2/pm2
1790  sn2(1)= pj2/pm2
1791  DO n=2,15
1792  ct2(n)= ct2(1)*ct2(n-1)-sn2(1)*sn2(n-1)
1793  sn2(n)= sn2(1)*ct2(n-1)+ct2(1)*sn2(n-1)
1794  ENDDO
1795  res4(1)= pm2*(bf(0)*ct2(1)+pm2*(bf(1)*ct2(2)+pm2*
1796  # (bf(2)*ct2(3)+pm2*(bf(3)*ct2(4)+pm2*
1797  # (bf(4)*ct2(5)+pm2*(bf(5)*ct2(6)+pm2*
1798  # (bf(6)*ct2(7)+pm2*(bf(7)*ct2(8)+pm2*
1799  # (bf(8)*ct2(9)+pm2*(bf(9)*ct2(10)+pm2*
1800  # (bf(10)*ct2(11)+pm2*(bf(11)*ct2(12)+pm2*
1801  # (bf(12)*ct2(13)+pm2*(bf(13)*ct2(14)+pm2*
1802  # (bf(14)*ct2(15))))))))))))))))
1803  res4(2)= pm2*(bf(0)*sn2(1)+pm2*(bf(1)*sn2(2)+pm2*
1804  # (bf(2)*sn2(3)+pm2*(bf(3)*sn2(4)+pm2*
1805  # (bf(4)*sn2(5)+pm2*(bf(5)*sn2(6)+pm2*
1806  # (bf(6)*sn2(7)+pm2*(bf(7)*sn2(8)+pm2*
1807  # (bf(8)*sn2(9)+pm2*(bf(9)*sn2(10)+pm2*
1808  # (bf(10)*sn2(11)+pm2*(bf(11)*sn2(12)+pm2*
1809  # (bf(12)*sn2(13)+pm2*(bf(13)*sn2(14)+pm2*
1810  # (bf(14)*sn2(15))))))))))))))))
1811  resb= -res3-res4+addt2
1812  ENDIF
1813  hto_li3= -resa+0.25d0*resb+addx
1814  RETURN
1815  ENDIF
1816 *

◆ hto_li3_srsz()

real*8 function, dimension(2) hto_sp_fun::hto_li3_srsz ( real*8, dimension (:), intent(in)  x,
real*8, dimension (:), intent(in)  y,
integer  unit 
)

Definition at line 1138 of file CALLING_cpHTO.f.

1138  USE hto_acmplx_pro
1139  USE hto_full_ln
1140  USE hto_riemann
1141  USE hto_units
1142 *
1143  IMPLICIT NONE
1144 *
1145  INTEGER unit
1146  real*8 xms,atheta,theta,sr,si
1147  real*8, dimension (:) :: x,y
1148  real*8, dimension(3,2) :: aux
1149  real*8, dimension(2) :: lnc,lncs,add,value
1150  INTENT(IN) x,y
1151 *
1152  IF(abs(y(2)).ne.1.d0) THEN
1153  IF(unit.eq.1) THEN
1154  xms= x(1)*x(1)+x(2)*x(2)
1155  IF(abs(1.d0-sqrt(xms)).gt.1.d-12) THEN
1156  print*,' apparent inconsistency '
1157  ENDIF
1158  atheta= atan(abs(x(2)/x(1)))
1159  sr= x(1)/abs(x(1))
1160  si= x(2)/abs(x(2))
1161  IF(sr > 0.d0.and.si > 0.d0) THEN
1162  theta= atheta
1163  ELSEIF(sr > 0.d0.and.si < 0.d0) THEN
1164  theta= 2.d0*pi-atheta
1165  ELSEIF(sr < 0.d0.and.si > 0.d0) THEN
1166  theta= pi-atheta
1167  ELSEIF(sr < 0.d0.and.si < 0.d0) THEN
1168  theta= pi+atheta
1169  ENDIF
1170  aux= hto_poly_unit(theta)
1171  value(1:2)= aux(2,1:2)
1172  ELSE
1173  value= hto_li3(x)
1174  ENDIF
1175  RETURN
1176  ENDIF
1177 *
1178  IF(x(1) > 1.d0) THEN
1179  IF(y(1) < 1.d0) THEN
1180  print*,' anomaly'
1181  print*,x
1182  print*,y
1183  stop
1184  ENDIF
1185  lnc= x(1).fln.x(2)
1186  lncs= lnc.cp.lnc
1187  add(1)= -lncs(2)
1188  add(2)= lncs(1)
1189  IF(y(2) < 0.d0.and.x(2) > 0.d0) THEN
1190  value= hto_li3(x)-pi*add
1191  ELSEIF(y(2) > 0.d0.and.x(2) < 0.d0) THEN
1192  value= hto_li3(x)+pi*add
1193  ELSE
1194  value= hto_li3(x)
1195  ENDIF
1196  ELSE
1197  value= hto_li3(x)
1198  ENDIF
1199 *
1200  RETURN
1201 *

◆ hto_poly_unit()

real*8 function, dimension(3,2) hto_sp_fun::hto_poly_unit ( real*8  theta)

Definition at line 1208 of file CALLING_cpHTO.f.

1208  USE hto_riemann
1209  USE hto_full_ln
1210  USE hto_real_lnz
1211  USE hto_imag_lnz
1212  USE hto_linear_comb_c
1213  USE hto_acmplx_pro
1214  USE hto_bernoulli
1215  USE hto_units
1216 *
1217  IMPLICIT NONE
1218 *
1219  INTEGER n,i
1220  real*8 theta
1221  real*8, dimension(2) :: arg,ln,z,sum2,sum3,sum4,pw2,pw3,pw4,
1222  # cb2,cb3,cb4
1223  real*8, dimension(4) :: psi
1224  real*8, dimension(0:15) :: rzf
1225  real*8, dimension(3,2) :: value
1226  real*8, dimension(0:17) :: c2,c3,c4,fac
1227 *
1228  data rzf(0:15) /-0.5d0,0.d0,1.644934066848226d0,
1229  # 1.202056903159594d0,1.082323233711138d0,
1230  # 1.036927755143370d0,1.017343061984449d0,
1231  # 1.008349277381923d0,1.004077356197944d0,
1232  # 1.002008392826082d0,1.000994575127818d0,
1233  # 1.000494188604119d0,1.000246086553308d0,
1234  # 1.000122713346578d0,1.000061248135059d0,
1235  # 1.000030588236307d0/
1236 *
1237  psi(1)= eg
1238  psi(2)= psi(1)+1.d0
1239  psi(3)= psi(2)+0.5d0
1240  psi(4)= psi(3)+1.d0/3.d0
1241 *
1242  fac(0)= 1.d0
1243  DO n=1,17
1244  fac(n)= n*fac(n-1)
1245  ENDDO
1246 *
1247  arg(1)= 0.d0
1248  arg(2)= -theta
1249 *
1250  ln= arg(1).fln.arg(2)
1251 *
1252  DO n=0,17
1253  IF((2-n) >= 0) THEN
1254  c2(n)= rzf(2-n)/fac(n)
1255  ELSE
1256  c2(n)= -b_num(n-1)/(n-1)/fac(n)
1257  ENDIF
1258  IF((3-n) >= 0) THEN
1259  c3(n)= rzf(3-n)/fac(n)
1260  ELSE
1261  c3(n)= -b_num(n-2)/(n-2)/fac(n)
1262  ENDIF
1263  IF((4-n) >= 0) THEN
1264  c4(n)= rzf(4-n)/fac(n)
1265  ELSE
1266  c4(n)= -b_num(n-3)/(n-3)/fac(n)
1267  ENDIF
1268  ENDDO
1269 *
1270  z(1)= 0.d0
1271  z(2)= theta
1272 *
1273  sum2= c2.lcc.z
1274  sum3= c3.lcc.z
1275  sum4= c4.lcc.z
1276 *
1277  pw2= z
1278  pw3= z.cp.pw2
1279  pw4= z.cp.pw3
1280 *
1281  cb2(1:2)= co(1:2)*(psi(2)-psi(1))-ln(1:2)
1282  cb3(1:2)= co(1:2)*(psi(3)-psi(1))-ln(1:2)
1283  cb4(1:2)= co(1:2)*(psi(4)-psi(1))-ln(1:2)
1284 *
1285  cb2= pw2.cp.cb2
1286  cb3= pw3.cp.cb3
1287  cb4= pw4.cp.cb4
1288 *
1289  value(1,1:2)= sum2(1:2)+cb2(1:2)
1290 *
1291  value(2,1:2)= sum3(1:2)+0.5d0*cb3(1:2)
1292 *
1293  value(3,1:2)= 1.d0/6.d0*(sum4(1:2)+1.d0/6.d0*cb4(1:2))
1294 *
1295  RETURN
1296 *

◆ hto_s_niels_up4()

recursive real*8 function, dimension(6,2) hto_sp_fun::hto_s_niels_up4 ( real*8, dimension(2), intent(in)  x)

Definition at line 766 of file CALLING_cpHTO.f.

766  USE hto_riemann
767  USE hto_acmplx_pro
768  USE hto_cmplx_ln
769  USE hto_real_ln
770  USE hto_common_niels
771  USE hto_linear_comb
772  USE hto_units
773 *
774  IMPLICIT NONE
775 *
776  real*8, dimension(6,2) :: res
777  real*8, dimension(2) :: x
778  INTENT(IN) x
779 *
780  INTEGER i,j,l
781  real*8 ym,ym2,ymod,s1,s2,s3,l2,xmo,z_exp,sum,s13m1,s22m1
782  real*8, dimension(6) :: sign
783  real*8, dimension(0:15) :: sumr
784  real*8, dimension(6,2) :: aniels,bniels,add
785  real*8, dimension(2) :: ln_omx,ln2_omx,ln3_omx,ln_omy,ln2_omy,
786  # ln3_omy,ln_y,ln2_y,ln3_y,ln_yi,ln_myi,ln2_myi,ln3_myi,ln4_y,
787  # y,omy,yi,yt,omx,myi,prod,li2,li3,s12,prod1,prod2,prodl,
788  # prodl1,prodl2,add1,add2,add3,add4,add5,add6,add7,add8,add9,
789  # add10,ln_my,ln2_my,ln3_my,ln4_my,add11,add12,add13,add14,
790  # add15,add16,add17,add18,add19,add20,add21,ln4_omx,ln_ymo,
791  # ln2_ymo,ln3_ymo,ln4_ymo
792  real*8, parameter :: li4h=6.72510831970049127d-2
793 *
794  IF(abs(x(2)).ne.1.d0) THEN
795  print 1,x
796  stop
797  ENDIF
798  1 format(' wrong input for niels ',2e20.5)
799 *
800  l2= log(2.d0)
801  s13m1= li4h+l2*(7.d0/8.d0*rz3-rz2*l2/4.d0+l2*l2*l2/24.d0)-rz4
802  s22m1= 2.d0*li4h+l2*(7.d0/4.d0*rz3-0.5d0*rz2*l2+l2*l2*l2/12.d0)-
803  # 15.d0/8.d0*rz4
804 *
805  IF(x(1).eq.1.d0) THEN
806  res(1,1)= rz2
807  res(2,1)= rz3
808  res(3,1)= rz3
809  res(4,1)= rz4
810  res(5,1)= -1/2*(rz2**2-3.d0*rz4)
811  res(6,1)= rz4
812  FORALL (i=1:6) res(i,2)= 0.d0
813  RETURN
814  ENDIF
815  IF(x(1).eq.0.5d0) THEN
816  res(1,1)= 0.5d0*(rz2-l2*l2)
817  res(2,1)= 7.d0/8.d0*rz3-l2*(0.5d0*rz2-l2*l2/6.d0)
818  res(3,1)= rz3/8.d0-l2*l2*l2/6.d0
819  res(4,1)= li4h
820  res(5,1)= -l2*res(3,1)-l2**4/8.d0-rz2*rz2/4.d0+3.d0/4.d0*rz4
821  res(6,1)= -res(4,1)+l2*res(3,1)+l2**2/4.d0*(1.d0/3.d0+rz2)-
822  # l2*rz3-rz4
823  FORALL (i=1:6) res(i,2)= 0.d0
824  RETURN
825  ENDIF
826 *
827 * gets rid of arguments with negative real part
828 *
829  IF(x(1) < 0.d0) THEN
830  xmo= x(1)-1.d0
831  ymod= xmo*xmo
832  y(1)= x(1)*xmo/ymod
833  y(2)= -x(2)
834  omx= co-x
835  ln_omx= omx(1).cln.omx(2)
836  ln2_omx= ln_omx.cp.ln_omx
837  ln3_omx= ln2_omx.cp.ln_omx
838  ln4_omx= ln3_omx.cp.ln_omx
839  aniels= hto_s_niels_up4(y)
840  li2(1:2)= aniels(1,1:2)
841  li3(1:2)= aniels(2,1:2)
842  s12(1:2)= aniels(3,1:2)
843  add1= ln2_omx.cp.li2
844  add2= ln_omx.cp.li3
845  add3= ln_omx.cp.s12
846  add(1,1:2)= -0.5d0*ln2_omx(1:2)
847  sign(1)= -1.d0
848  prod= ln_omx.cp.li2
849  add(2,1:2)= aniels(3,1:2)-ln3_omx(1:2)/6.d0-prod(1:2)
850  sign(2)= -1.d0
851  add(3,1:2)= ln3_omx(1:2)/6.d0
852  sign(3)= 1.d0
853  add(4,1:2)= -0.5d0*add1(1:2)-add2(1:2)+add3(1:2)
854  # -1.d0/24.d0*ln4_omx(1:2)-aniels(6,1:2)+aniels(5,1:2)
855  # +co(1:2)*(2.d0*s13m1-s22m1+1.d0/8.d0*rz4)
856  sign(4)= -1.d0
857  add(5,1:2)= add3(1:2)+1.d0/24.d0*ln4_omx(1:2)
858  # -2.d0*aniels(6,1:2)+aniels(5,1:2)+co(1:2)*(4.d0*s13m1
859  # -2.d0*s22m1+0.25d0*rz4)
860  sign(5)= 1.d0
861  add(6,1:2)= -1.d0/24.d0*ln4_omx(1:2)-aniels(6,1:2)+
862  # co(1:2)*(4.d0*s13m1-2.d0*s22m1+0.25d0*rz4)
863  sign(6)= -1.d0
864  ELSE
865  y= x
866  add= 0.d0
867  sign= 1.d0
868  ENDIF
869 *
870  ym2= y(1)*y(1)
871  ym= abs(y(1))
872 *
873 * |y| < 1 & Re(y) < 1/2
874 *
875  IF(ym < 1.d0.and.y(1) < 0.5d0) THEN
876  omy= co-y
877  z_exp= -(omy(1).rln.omy(2))
878  DO l=1,3
879  sumr(0:15)= plr(l,0:15)
880  sum= sumr.lc.z_exp
881  res(l,1)= sum
882  res(l,2)= 0.d0
883  ENDDO
884  DO l=4,6
885  sumr(0:15)= plr_4(l-3,0:15)
886  sum= sumr.lc.z_exp
887  res(l,1)= sum
888  res(l,2)= 0.d0
889  ENDDO
890  FORALL (i=1:6,j=1:2) res(i,j)= sign(i)*res(i,j)+add(i,j)
891  RETURN
892 *
893 * |y| < 1 & Re(y) > 1/2
894 *
895  ELSEIF(ym < 1.d0.and.y(1) > 0.5d0) THEN
896  omy= co-y
897  ln_omy= omy(1).cln.omy(2)
898  ln_y= y(1).cln.y(2)
899  ln2_omy= ln_omy.cp.ln_omy
900  ln2_y= ln_y.cp.ln_y
901  ln3_y= ln2_y.cp.ln_y
902  aniels= hto_s_niels_up4(omy)
903  s12(1:2)= aniels(3,1:2)
904  li3(1:2)= aniels(2,1:2)
905  li2(1:2)= aniels(1,1:2)
906  prodl= ln_y.cp.ln_omy
907  res(1,1:2)= -li2(1:2)-prodl(1:2)+rz2*co(1:2)
908  prodl= ln_omy.cp.ln2_y
909  prod= ln_y.cp.li2
910  res(2,1:2)= -s12(1:2)-prod(1:2)-0.5d0*prodl(1:2)+
911  # rz2*ln_y(1:2)+rz3*co(1:2)
912  prodl= ln_y.cp.ln2_omy
913  prod= ln_omy.cp.li2
914  res(3,1:2)= -li3(1:2)+prod(1:2)+0.5d0*prodl(1:2)+rz3*co(1:2)
915  add1= ln2_y.cp.li2
916  add2= ln_y.cp.s12
917  add3= ln_omy.cp.ln3_y
918  add4= ln_omy.cp.s12
919  add5= ln_y.cp.ln_omy
920  add5= add5.cp.li2
921  add6= ln_y.cp.li3
922  add7= ln2_y.cp.ln2_omy
923  add8= ln_omy.cp.li3
924  add9= ln_y.cp.ln3_omy
925  add10= ln2_omy.cp.li2
926  res(4,1:2)= -aniels(6,1:2)-0.5d0*add1(1:2)-add2(1:2)-
927  # add3(1:2)/6.d0+rz3*ln_y(1:2)+0.5d0*rz2*ln2_y+rz4*co(1:2)
928  res(5,1:2)= aniels(5,1:2)-add4(1:2)-add5(1:2)+add6(1:2)
929  # -add2(1:2)-0.25d0*add7(1:2)+0.5d0*co(1:2)*(rz2*rz2-3.d0*rz4)
930  res(6,1:2)= aniels(4,1:2)-add8(1:2)+add9(1:2)/6.d0+
931  # 0.5d0*add10(1:2)-co(1:2)*s13m1
932  FORALL (i=1:6,j=1:2) res(i,j)= sign(i)*res(i,j)+add(i,j)
933  RETURN
934 *
935 * |y| > 1 & Re(y^-1) < 1/2
936 *
937  ELSEIF(ym > 1.d0.and.y(1) >= 2.d0) THEN
938  yi(1)= y(1)/ym2
939  yi(2)= -y(2)
940  ln_yi= yi(1).cln.yi(2)
941  ln_my= (-y(1)).cln.(-y(2))
942  ln2_my= ln_my.cp.ln_my
943  ln3_my= ln3_my.cp.ln_my
944  ln4_my= ln3_my.cp.ln_my
945  myi(1)= -y(1)/ym2
946  myi(2)= y(2)
947  ln_myi= myi(1).cln.myi(2)
948  ln2_myi= ln_myi.cp.ln_myi
949  ln3_myi= ln2_myi.cp.ln_myi
950  aniels= hto_s_niels_up4(yi)
951  li2(1:2)= aniels(1,1:2)
952  li3(1:2)= aniels(2,1:2)
953  s12(1:2)= aniels(3,1:2)
954  res(1,1:2)= -aniels(1,1:2)-0.5d0*ln2_myi(1:2)-rz2*co(1:2)
955  res(2,1:2)= aniels(2,1:2)+ln3_myi(1:2)/6.d0+rz2*ln_myi(1:2)
956  li2(1:2)= aniels(1,1:2)
957  prod= ln_myi.cp.li2
958  res(3,1:2)= aniels(2,1:2)-aniels(3,1:2)-ln3_myi(1:2)/6.d0-
959  # prod(1:2)+rz3*co(1:2)
960  res(4,1:2)= aniels(4,1:2)+0.5d0*rz2*ln2_my(1:2)+
961  # 1.d0/24.d0*ln4_my+7.d0/4.d0*rz4*co(1:2)
962  add1= ln_my.cp.li3
963  add1= ln_my.cp.s12
964  add3= ln2_my.cp.li2
965  res(5,1:2)= -aniels(5,1:2)+2.d0*aniels(4,1:2)+add1(1:2)-
966  # rz3*ln_my-ln4_my/24.d0+7.d0/4.d0*rz4*co(1:2)
967  res(6,1:2)= aniels(4,1:2)+aniels(6,1:2)-aniels(5,1:2)+
968  # add1(1:2)-add2(1:2)+0.5d0*add3(1:2)+ln4_my/24.d0+
969  # co(1:2)*(-2.d0*s13m1+s22m1+7.d0/8.d0*rz4)
970  FORALL (i=1:6,j=1:2) res(i,j)= sign(i)*res(i,j)+add(i,j)
971  RETURN
972 *
973 * |y| > 1 & Re(y^-1) > 1/2
974 *
975  ELSEIF(ym > 1.d0.and.y(1) < 2.d0) THEN
976  yt(1)= (y(1)-1.d0)/y(1)
977  yt(2)= y(2)
978  omy= co-y
979  ln_y= y(1).cln.y(2)
980  ln_omy= omy(1).cln.omy(2)
981  ln2_y= ln_y.cp.ln_y
982  ln3_y= ln2_y.cp.ln_y
983  ln4_y= ln3_y.cp.ln_y
984  ln2_omy= ln_omy.cp.ln_omy
985  ln_my= (-y(1)).cln.(-y(2))
986  ln2_my= ln_my.cp.ln_my
987  ln3_my= ln2_my.cp.ln_my
988  ln4_my= ln3_my.cp.ln_my
989  aniels= hto_s_niels_up4(yt)
990  bniels= 0.d0
991  prod= ln_y.cp.ln_omy
992  bniels(1,1:2)= aniels(1,1:2)-prod(1:2)+0.5d0*ln2_y(1:2)+
993  # rz2*co(1:2)
994  res(1,1:2)= aniels(1,1:2)-prod(1:2)+0.5d0*ln2_y(1:2)+
995  # rz2*co(1:2)
996  prodl= ln2_y.cp.ln_omy
997  li2(1:2)= bniels(1,1:2)
998  prod= ln_y.cp.li2
999  res(2,1:2)= -aniels(3,1:2)+prod(1:2)+0.5d0*prodl(1:2)-
1000  # ln3_y(1:2)/6.d0+rz3*co(1:2)
1001  li2(1:2)= aniels(1,1:2)
1002  prodl1= ln_omy.cp.ln2_y
1003  prodl2= ln_y.cp.ln2_omy
1004  prod1= ln_omy.cp.li2
1005  prod2= ln_y.cp.li2
1006  res(3,1:2)= rz3*co(1:2)-0.5d0*prodl1(1:2)-prod1(1:2)+
1007  # 0.5d0*prodl2(1:2)+ln3_y(1:2)/6.d0+prod2(1:2)+
1008  # aniels(2,1:2)-aniels(3,1:2)
1009  add1= ln_ymo.cp.ln3_y
1010  add2= ln_y.cp.s12
1011  add3= ln2_y.cp.li2
1012  add4= ln_ymo.cp.s12
1013  add5= ln_my.cp.ln_ymo
1014  add5= add5.cp.ln2_y
1015  add6= ln_my.cp.s12
1016  add7= ln_y.cp.ln_ymo
1017  add7= add7.cp.li2
1018  add8= ln_my.cp.ln3_y
1019  add9= ln_y.cp.ln_my
1020  add9= add9.cp.li2
1021  add10= ln_y.cp.ln_my
1022  add11= ln_y.cp.li3
1023  add12= ln2_y.cp.ln2_ymo
1024  add13= ln_ymo.cp.li3
1025  add14= ln_my.cp.ln_ymo
1026  add14= add14.cp.li2
1027  add15= ln_my.cp.li3
1028  add16= ln_y.cp.ln_ymo
1029  add16= add16.cp.ln2_my
1030  add17= ln_y.cp.ln_my
1031  add17= add17.cp.ln2_ymo
1032  add18= ln_y.cp.ln3_ymo
1033  add19= ln2_ymo.cp.li2
1034  add20= ln2_my.cp.li2
1035  add21= ln2_y.cp.ln2_my
1036  res(4,1:2)= -1.d0/6.d0*add1(1:2)-add2(1:2)+ln_y(1:2)*rz3
1037  # -0.5d0*ln2_my*rz2+0.5d0*add3(1:2)-0.5d0*ln2_y(1:2)*rz2
1038  # -1.d0/24.d0*ln4_my(1:2)+1.d0/6.d0*ln4_y(1:2)+aniels(6,1:2)
1039  # -11.d0/4.d0*rz4*co(1:2)
1040  res(5,1:2)= -5.d0/6.d0*add1(1:2)+add4(1:2)+0.5d0*add5(1:2)
1041  # -0.5d0*add8(1:2)+add6(1:2)-add7(1:2)-add9(1:2)
1042  # +add10(1:2)*rz2+add11(1:2)-3.d0*add2(1:2)+ln_y(1:2)*rz3
1043  # +1.d0/4.d0*add12(1:2)+2.d0*add3(1:2)-ln2_y(1:2)*rz2
1044  # +1.d0/24.d0*ln4_my(1:2)+7.d0/12.d0*ln4_y(1:2)
1045  # +2.d0*aniels(6,1:2)-aniels(5,1:2)-7.d0/2.d0*rz4*co(1:2)
1046  res(6,1:2)= -7.d0/6.d0*add1(1:2)-add13(1:2)+add4(1:2)
1047  # +3.d0/2.d0*add5(1:2)+add14(1:2)-add8(1:2)-add15(1:2)
1048  # +add6(1:2)-0.5d0*add16(1:2)-2.d0*add7(1:2)-0.5d0*add17(1:2)
1049  # -2.d0*add9(1:2)+add10(1:2)*rz2-1.d0/6.d0*add18(1:2)
1050  # +2.d0*add11(1:2)-2.d0*add2(1:2)+0.5d0*add19(1:2)
1051  # +0.5d0*add20(1:2)-0.5d0*ln2_my*rz2+3.d0/4.d0*add12(1:2)
1052  # +0.5d0*add21(1:2)+2.d0*add3(1:2)-0.5d0*ln2_y(1:2)*rz2
1053  # -1.d0/24.d0*ln4_my(1:2)+ 7.d0/12.d0*ln4_y(1:2)+aniels(4,1:2)
1054  # +aniels(6,1:2)-aniels(5,1:2)+co(1:2)*(2.d0*s13m1-s22m1
1055  # -21.d0/8.d0*rz4)
1056  FORALL (i=1:6,j=1:2) res(i,j)= sign(i)*res(i,j)+add(i,j)
1057  RETURN
1058  ENDIF
1059 *
hto_cmplx_ln
Definition: CALLING_cpHTO.f:439
hto_bernoulli::b_num
real *8, dimension(0:18) b_num
Definition: CALLING_cpHTO.f:210
value
pymela::gHIGGS_KAPPA value("gHIGGS_KAPPA_TILDE", pymela::gHIGGS_KAPPA_TILDE) .value("SIZE_HQQ"
hto_acmplx_pro
Definition: CALLING_cpHTO.f:257
hto_linear_comb_c::lcc
real *8 function, dimension(2) lcc(c, x)
Definition: CALLING_cpHTO.f:346
hto_real_ln
Definition: CALLING_cpHTO.f:424
hto_units::one
real *8, parameter one
Definition: CALLING_cpHTO.f:184
hto_riemann::rz3
real *8, parameter rz3
Definition: CALLING_cpHTO.f:200
hto_riemann
Definition: CALLING_cpHTO.f:195
hto_linear_comb_c
Definition: CALLING_cpHTO.f:340
hto_units::co
real *8, dimension(1:2) co
Definition: CALLING_cpHTO.f:188
hto_kountac::kp
integer kp
Definition: CALLING_cpHTO.f:725
hto_acmplx_rat
Definition: CALLING_cpHTO.f:292
hto_cmplx_ln::cln
real *8 function, dimension(2) cln(x, ep)
Definition: CALLING_cpHTO.f:445
hto_kountac
Definition: CALLING_cpHTO.f:724
hto_units
Definition: CALLING_cpHTO.f:179
hto_imag_lnz
Definition: CALLING_cpHTO.f:379
hto_full_ln
Definition: CALLING_cpHTO.f:464
hto_linear_comb::lc
real *8 function lc(c, x)
Definition: CALLING_cpHTO.f:736
hto_units::eps
real *8, parameter eps
Definition: CALLING_cpHTO.f:183
hto_common_niels::plr
real *8, dimension(3, 0:15) plr
Definition: CALLING_cpHTO.f:718
hto_real_lnz
Definition: CALLING_cpHTO.f:364
hto_full_ln::fln
real *8 function, dimension(2) fln(x, y)
Definition: CALLING_cpHTO.f:470
hto_bernoulli
Definition: CALLING_cpHTO.f:209
hto_riemann::rz4
real *8, parameter rz4
Definition: CALLING_cpHTO.f:201
hto_common_niels::plr_4
real *8, dimension(3, 0:15) plr_4
Definition: CALLING_cpHTO.f:719
hto_riemann::pis
real *8, parameter pis
Definition: CALLING_cpHTO.f:197
hto_riemann::eg
real *8, parameter eg
Definition: CALLING_cpHTO.f:203
hto_real_ln::rln
real *8 function rln(x, ep)
Definition: CALLING_cpHTO.f:430
hto_acmplx_pro::cp
real *8 function, dimension(2) cp(x, y)
Definition: CALLING_cpHTO.f:263
hto_kountac::km
integer km
Definition: CALLING_cpHTO.f:725
hto_common_niels
Definition: CALLING_cpHTO.f:717
hto_linear_comb
Definition: CALLING_cpHTO.f:730
hto_riemann::rz2
real *8, parameter rz2
Definition: CALLING_cpHTO.f:199
Generated by   doxygen 1.8.17