JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
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