64 #
clxe,
clxmu,
clxtau,
clxu,
clxd,
clxc,
clxs,
clxt,
clxb,
muhcp,
65 #
scalec,
cxqs,
cxq,
cxls,
cxl,
cxtmb,
cxtmbs,
cxbi,
cxtmbi,
xq,
74 real*8,
parameter ::
m0= 0.d0
75 real*8,
parameter ::
mw= 80.398d0
76 real*8,
parameter ::
mz= 91.1876d0
77 real*8,
parameter ::
me= 0.51099907d-3
78 real*8,
parameter ::
mm= 0.105658389d0
79 real*8,
parameter ::
mtl= 1.77684d0
80 real*8,
parameter ::
muq= 0.190d0
81 real*8,
parameter ::
mdq= 0.190d0
82 real*8,
parameter ::
mcq= 1.55d0
83 real*8,
parameter ::
msq= 0.190d0
84 real*8,
parameter ::
mbq= 4.69d0
85 real*8,
parameter ::
mb= 4.69d0
86 real*8,
parameter ::
mtiny= 1.d-10
87 real*8,
parameter ::
imw= 2.08872d0
88 real*8,
parameter ::
imz= 2.4952d0
124 CHARACTER (len=8) ampm
129 CHARACTER (len= 9),
parameter,
dimension(12) :: month= (/
130 #
'January ',
'February ',
'March ',
'April ',
131 #
'May ',
'June ',
'July ',
'August ',
132 #
'September',
'October ',
'November ',
'December ' /)
138 CALL date_and_time (values= values)
150 ELSE IF (h== 12)
THEN
151 IF(n== 0 .and. s== 0)
THEN
160 ELSE IF (h== 12)
THEN
161 IF(n== 0 .and. s== 0)
THEN
169 WRITE (*,
'(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)')
170 # d,trim(month(m)),y,h,
':',n,
':',s,
'.',mm,trim(ampm)
183 real*8,
parameter ::
eps= -1.d0
184 real*8,
parameter ::
one= 1.d0
185 real*8,
parameter ::
zero= 0.d0
186 real*8,
parameter ::
qeps= 1.d-25
187 real*8,
parameter ::
ez= 1.d-15
188 real*8 ::
co(1:2)=(/1.d0,0.d0/)
189 real*8 ::
ci(1:2)=(/0.d0,1.d0/)
190 real*8 ::
cz(1:2)=(/0.d0,0.d0/)
196 real*8,
parameter :: pi= 3.141592653589793238462643d0
197 real*8,
parameter ::
pis= pi*pi
199 real*8,
parameter ::
rz2= 1.64493406684823d0
200 real*8,
parameter ::
rz3= 1.20205690315959d0
201 real*8,
parameter ::
rz4= 1.08232323371114d0
202 real*8,
parameter ::
rz5= 1.03692775514337d0
203 real*8,
parameter ::
eg= 0.5772156649d0
204 real*8,
parameter ::
ln_pi= 0.114472988584940016d1
210 real*8 ::
b_num(0:18)=(/1.d0,-5.d-1,1.66666666666666d-1,
211 # 0.d0,-3.33333333333333d-2,0.d0,2.38095238095238d-2,
212 # 0.d0,-3.33333333333333d-2,0.d0,7.57575757575757d-2,
213 # 0.d0,-2.53113553113553d-1,0.d0,1.16666666666666d0,
214 # 0.d0,-7.09215686274509d0,0.d0,5.49711779448621d1/)
258 INTERFACE OPERATOR ( .
cp. )
266 real*8,
dimension(2) ::
cp,z1,z2
267 real*8,
INTENT(IN),
dimension(2) :: x,y
268 IF(abs(x(2)).eq.1.d0.and.abs(y(2)).eq.1.d0)
THEN
272 IF(abs(x(2)).eq.1.d0)
THEN
278 IF(abs(y(2)).eq.1.d0)
THEN
284 cp(1)= z1(1)*z2(1)-z1(2)*z2(2)
285 cp(2)= z1(1)*z2(2)+z1(2)*z2(1)
293 INTERFACE OPERATOR ( .
cq. )
299 real*8,
dimension(2) ::
cq,z1,z2
300 real*8,
INTENT(IN),
dimension(2) :: x,y
302 IF(x(2).eq.0.d0.and.y(2).eq.0.d0)
THEN
307 IF(abs(x(2)).eq.1.d0.and.abs(y(2)).eq.1.d0)
THEN
311 IF(abs(x(2)).eq.1.d0)
then
316 ELSEIF(abs(y(2)).eq.1.d0)
then
325 IF(abs(z2(1)) > abs(z2(2)))
THEN
327 cq(1)= (z1(1)+theta*z1(2))/(theta*z2(2)+z2(1))
328 cq(2)= (z1(2)-theta*z1(1))/(theta*z2(2)+z2(1))
331 cq(1)= (theta*z1(1)+z1(2))/(theta*z2(1)+z2(2))
332 cq(2)= (theta*z1(2)-z1(1))/(theta*z2(1)+z2(2))
341 INTERFACE OPERATOR ( .
lcc. )
349 real*8,
dimension(:) :: c
350 real*8,
dimension(2) :: x,
lcc,sum
356 sum= (sum.
cp.x)+c(i)*
co
366 MODULE PROCEDURE rlnz
372 real*8,
INTENT(IN) :: x,y
373 rlnz= 0.5d0*log(x*x+y*y)
381 MODULE PROCEDURE ilnz
387 real*8,
INTENT(IN) :: x,y
388 real*8 teta,tnteta,sr,si,ax,ay,
ilnz
399 ELSE IF(y.eq.0.d0)
THEN
425 INTERFACE OPERATOR ( .
rln. )
432 real*8,
INTENT(IN) :: x,ep
440 INTERFACE OPERATOR ( .
cln. )
447 real*8,
dimension(2) ::
cln
448 real*8,
INTENT(IN) :: x,ep
449 IF(abs(ep).ne.1.d0)
THEN
450 print*,
' Wrong argument for CLN '
465 INTERFACE OPERATOR ( .
fln. )
472 real*8 ax,ay,teta,tnteta,sr,si
473 real*8,
INTENT(IN) :: x,y
474 real*8,
dimension(2) ::
fln
475 IF(abs(y).eq.1.d0)
THEN
483 fln(1)= 0.5d0*log(x*x+y*y)
493 ELSE IF(y.eq.0.d0)
THEN
519 MODULE PROCEDURE lnsrs
527 real*8,
dimension(2) ::
lnsrs,x,y,olnsrs
528 real*8 ax,ay,teta,tnteta,ilnx,sr,si,sx,sy,xms
531 IF(abs(y(2)).ne.1.d0)
THEN
533 xms= x(1)*x(1)+x(2)*x(2)
534 IF(abs(1.d0-sqrt(xms)).lt.1.d-12)
THEN
535 teta= atan(abs(x(2)/x(1)))
542 lnsrs(2)= si*(pi-teta)
547 olnsrs= y(1).
fln.y(2)
552 lnsrs(1)= 0.5d0*log(x(1)*x(1)+x(2)*x(2))
555 IF(x(1).eq.0.d0)
THEN
562 ELSE IF(x(2).eq.0.d0)
THEN
583 IF((x(2) > 0.d0.and.y(2) > 0).or.
584 # (x(2) < 0.d0.and.y(2) < 0))
THEN
586 ELSEIF(x(2) > 0.d0.and.y(2) < 0)
THEN
589 print*,
'+++++++++++++++++++++++++++++++++++'
590 print*,
' anomaly ln '
593 print*,
'+++++++++++++++++++++++++++++++++++'
599 IF(abs(y(2)).ne.1.d0)
THEN
600 print*,
' Wrong argument for LNSRS '
604 IF((x(2)*y(2)) < 0.d0)
THEN
605 lnsrs(2)= ilnx+2.d0*y(2)*pi
620 MODULE PROCEDURE crsrs
628 real*8,
dimension(2) ::
crsrs
629 real*8,
INTENT(IN),
dimension(2):: x,y
632 rs= x(1)*x(1)+x(2)*x(2)
640 rs= x(1)*x(1)+x(2)*x(2)
645 IF((x(2)*y(2)) > 0.d0)
THEN
648 ELSEIF((x(2)*y(2)) < 0.d0)
THEN
660 INTERFACE OPERATOR ( .
cr. )
666 real*8,
dimension(2) ::
cr
667 real*8,
INTENT(IN) :: x,ep
673 cr(2)= ep*sqrt(abs(x))
681 INTERFACE OPERATOR ( .
crz. )
688 real*8,
dimension(2) ::
crz
689 real*8,
INTENT(IN) :: x,y
691 IF(abs(y/x).lt.1.d-8)
THEN
693 crz(1)= sqrt(x)*(1.d0+0.5d0*y*y/(x*x))
694 crz(2)= 0.5d0*y/sqrt(x)
695 ELSEIF(x.lt.0.d0)
THEN
697 crz(1)= 0.5d0*abs(y)/sax
698 crz(2)= sign(
one,y)*sax*(1.d0+0.5d0*y*y/(x*x))
705 IF(x > 0.d0.and.y.eq.0.d0)
THEN
709 crz(2)= sign(
one,y)*sqrt(v)
718 real*8,
dimension(3,0:15) ::
plr
731 INTERFACE OPERATOR ( .
lc. )
737 real*8,
dimension(:) :: c
776 real*8,
dimension(6,2) :: res
777 real*8,
dimension(2) :: x
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
794 IF(abs(x(2)).ne.1.d0)
THEN
798 1
format(
' wrong input for niels ',2e20.5)
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)-
805 IF(x(1).eq.1.d0)
THEN
810 res(5,1)= -1/2*(
rz2**2-3.d0*
rz4)
812 FORALL (i=1:6) res(i,2)= 0.d0
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
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)-
823 FORALL (i=1:6) res(i,2)= 0.d0
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
840 li2(1:2)= aniels(1,1:2)
841 li3(1:2)= aniels(2,1:2)
842 s12(1:2)= aniels(3,1:2)
846 add(1,1:2)= -0.5d0*ln2_omx(1:2)
849 add(2,1:2)= aniels(3,1:2)-ln3_omx(1:2)/6.d0-prod(1:2)
851 add(3,1:2)= ln3_omx(1:2)/6.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)
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)
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)
875 IF(ym < 1.d0.and.y(1) < 0.5d0)
THEN
877 z_exp= -(omy(1).
rln.omy(2))
879 sumr(0:15)=
plr(l,0:15)
885 sumr(0:15)=
plr_4(l-3,0:15)
890 FORALL (i=1:6,j=1:2) res(i,j)= sign(i)*res(i,j)+add(i,j)
895 ELSEIF(ym < 1.d0.and.y(1) > 0.5d0)
THEN
897 ln_omy= omy(1).
cln.omy(2)
899 ln2_omy= ln_omy.
cp.ln_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
910 res(2,1:2)= -s12(1:2)-prod(1:2)-0.5d0*prodl(1:2)+
912 prodl= ln_y.
cp.ln2_omy
914 res(3,1:2)= -li3(1:2)+prod(1:2)+0.5d0*prodl(1:2)+
rz3*
co(1:2)
917 add3= ln_omy.
cp.ln3_y
922 add7= ln2_y.
cp.ln2_omy
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)
937 ELSEIF(ym > 1.d0.and.y(1) >= 2.d0)
THEN
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
947 ln_myi= myi(1).
cln.myi(2)
948 ln2_myi= ln_myi.
cp.ln_myi
949 ln3_myi= ln2_myi.
cp.ln_myi
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)
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)
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)
975 ELSEIF(ym > 1.d0.and.y(1) < 2.d0)
THEN
976 yt(1)= (y(1)-1.d0)/y(1)
980 ln_omy= omy(1).
cln.omy(2)
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
992 bniels(1,1:2)= aniels(1,1:2)-prod(1:2)+0.5d0*ln2_y(1:2)+
994 res(1,1:2)= aniels(1,1:2)-prod(1:2)+0.5d0*ln2_y(1:2)+
996 prodl= ln2_y.
cp.ln_omy
997 li2(1:2)= bniels(1,1:2)
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
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
1013 add5= ln_my.
cp.ln_ymo
1016 add7= ln_y.
cp.ln_ymo
1018 add8= ln_my.
cp.ln3_y
1021 add10= ln_y.
cp.ln_my
1023 add12= ln2_y.
cp.ln2_ymo
1024 add13= ln_ymo.
cp.li3
1025 add14= ln_my.
cp.ln_ymo
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
1056 FORALL (i=1:6,j=1:2) res(i,j)= sign(i)*res(i,j)+add(i,j)
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
1080 IF(abs(y(2)).ne.1.d0)
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 '
1086 atheta= atan(abs(x(2)/x(1)))
1089 IF(sr > 0.d0.and.si > 0.d0)
THEN
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
1095 ELSEIF(sr < 0.d0.and.si < 0.d0)
THEN
1099 value(1:2)= aux(1,1:2)
1106 IF(x(1) > 1.d0)
THEN
1107 IF(y(1) < 1.d0)
THEN
1109 print*,
'+++++++++++++++++++++++++++++++++++'
1110 print*,
' anomaly Li2 '
1113 print*,
'+++++++++++++++++++++++++++++++++++'
1118 IF(y(2) < 0.d0.and.x(2) > 0.d0)
THEN
1121 ELSEIF(y(2) > 0.d0.and.x(2) < 0.d0)
THEN
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
1152 IF(abs(y(2)).ne.1.d0)
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 '
1158 atheta= atan(abs(x(2)/x(1)))
1161 IF(sr > 0.d0.and.si > 0.d0)
THEN
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
1167 ELSEIF(sr < 0.d0.and.si < 0.d0)
THEN
1171 value(1:2)= aux(2,1:2)
1178 IF(x(1) > 1.d0)
THEN
1179 IF(y(1) < 1.d0)
THEN
1189 IF(y(2) < 0.d0.and.x(2) > 0.d0)
THEN
1191 ELSEIF(y(2) > 0.d0.and.x(2) < 0.d0)
THEN
1221 real*8,
dimension(2) :: arg,ln,z,sum2,sum3,sum4,pw2,pw3,pw4,
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
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/
1239 psi(3)= psi(2)+0.5d0
1240 psi(4)= psi(3)+1.d0/3.d0
1250 ln= arg(1).
fln.arg(2)
1254 c2(n)= rzf(2-n)/fac(n)
1256 c2(n)= -
b_num(n-1)/(n-1)/fac(n)
1259 c3(n)= rzf(3-n)/fac(n)
1261 c3(n)= -
b_num(n-2)/(n-2)/fac(n)
1264 c4(n)= rzf(4-n)/fac(n)
1266 c4(n)= -
b_num(n-3)/(n-3)/fac(n)
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)
1289 value(1,1:2)= sum2(1:2)+cb2(1:2)
1291 value(2,1:2)= sum3(1:2)+0.5d0*cb3(1:2)
1293 value(3,1:2)= 1.d0/6.d0*(sum4(1:2)+1.d0/6.d0*cb4(1:2))
1301 FUNCTION hto_li2(x)
RESULT(value)
1311 real*8,
dimension(2) ::x,
value
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
1323 IF(x(1) < 0.d0)
THEN
1327 clnomx= omx(1).
fln.omx(2)
1328 add1=
pis/6.d0*
co-(clnx.
cp.clnomx)
1335 ym2= y(1)*y(1)+y(2)*y(2)
1342 clnoy= oy(1).
fln.oy(2)
1343 add2= -
pis/6.d0*
co-0.5d0*(clnoy.
cp.clnoy)
1356 clnomz= omz(1).
fln.omz(2)
1357 add3=
pis/6.d0*
co-(clnz.
cp.clnomz)
1364 par= omt(1).
fln.omt(2)
1367 bf(n)=
b_num(n)/fact
1372 parm2= parr*parr+pari*pari
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)
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)+
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)+
1392 value= sign1*(sign2*(sign3*res+add3)+add2)+add1
1406 real*8,
dimension (:) :: x
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/)
1432 FORALL (n=0:14) bf(n)= b(n)/(n+1.d0)
1434 xm2= x(1)*x(1)+x(2)*x(2)
1440 IF(xtst <= 1.d-20)
THEN
1443 ELSE IF(xm > 1.d-20)
THEN
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))
1457 y2r= y(1)*y(1)-y(2)*y(2)
1458 IF(y(1) >= 0.d0.or.y2r < 0.d0)
THEN
1460 IF(ytst <= 1.d-20)
THEN
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)
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))))))))))))))))
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
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)
1517 p12= pr1*pr1+pj1*pj1
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)
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))))))))))))))))
1548 p22= pr2*pr2+pj2*pj2
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)
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))))))))))))))))
1579 ELSE IF(y(1) < 0.d0)
THEN
1584 IF(t(1) <= 0.5d0)
THEN
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)
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
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)
1639 p12= pr1*pr1+pj1*pj1
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)
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))))))))))))))))
1670 p22= pr2*pr2+pj2*pj2
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)
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
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
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)
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
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)
1756 p12= pr1*pr1+pj1*pj1
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)
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))))))))))))))))
1787 p22= pr2*pr2+pj2*pj2
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)
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
1813 hto_li3= -resa+0.25d0*resb+addx
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
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
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
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
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
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
1916 real*8 teta,zm,zm2,tnteta,sr,si
1917 real*8,
dimension(2) :: arg,aarg,res
1920 IF((abs(arg(2))-
eps).eq.0.d0)
THEN
1921 res(1)= log(abs(arg(1)))
1922 IF(arg(1) > 0.d0)
THEN
1925 res(2)= pi*sign(
one,arg(2))
1931 zm2= (arg(1))**2+(arg(2))**2
1934 IF(arg(1).eq.0.d0)
THEN
1935 IF(arg(2) > 0.d0)
THEN
1942 ELSE IF(arg(2).eq.0.d0)
THEN
1943 IF(arg(1) > 0.d0)
THEN
1951 tnteta= aarg(2)/aarg(1)
1958 res(2)= si*(pi-teta)
1975 real*8,
dimension(2) :: arg,omarg,res,ares
1976 real*8,
dimension(10) :: ct,sn
1977 INTENT(IN) arg,omarg
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)
1990 ares(1)= ct(10)/10.d0
1991 ares(2)= sn(10)/10.d0
1993 ares(1)= ares(1)*zm+ct(k)/k
1994 ares(2)= ares(2)*zm+sn(k)/k
1996 ares(1)= -ares(1)*zm
1997 ares(2)= -ares(2)*zm
2101 real*8 fr2,mur,asmur,mc,mb,mt,a,b,
hto_dzero,
2103 real*8,
parameter :: eps=1.d-10
2104 INTEGER,
parameter :: MAXF=10000
2105 INTEGER,
parameter :: MODE=1
2108 IF(mur*sqrt(fr2).LE.mc)
THEN
2115 r0c = 1.d0/sqrt(fr2)
2177 real*8 fr2,r0,asi,mc,mb,mt,r20,mc2,mb2,mt2
2178 real*8,
parameter :: pi= 3.14159265358979d0
2188 zeta(1) = 0.5772 15664 90153 d0
2189 zeta(2) = 1.64493 40668 48226 d0
2190 zeta(3) = 1.20205 69031 59594 d0
2191 zeta(4) = 1.08232 32337 11138 d0
2192 zeta(5) = 1.03692 77551 43370 d0
2193 zeta(6) = 1.01734 30619 84449 d0
2210 IF( (
ivfns .EQ. 0) .AND. (
nff .LT. 3) )
THEN
2211 print*,
'Wrong flavour number for FFNS evolution. STOP'
2214 IF( (
ivfns .EQ. 0) .AND. (
nff .GT. 5) )
THEN
2215 print*,
'Wrong flavour number for FFNS evolution. STOP'
2219 IF(
naord .GT. 3 )
THEN
2220 print*,
'Specified order in a_s too high. STOP'
2224 IF( (
ivfns .NE. 0) .AND. (fr2 .GT. 4.001d0) )
THEN
2225 print*,
'Too low mu_r for VFNS evolution. STOP'
2229 IF( (
ivfns .EQ. 1) .AND. (
m20 .GT. mc2) )
THEN
2230 print*,
'Too high mu_0 for VFNS evolution. STOP'
2234 IF( (asi .GT. 2.d0) .OR. (asi .LT. 2.d-2) )
THEN
2235 print*,
'alpha_s out of range. STOP'
2239 IF( (
ivfns .EQ. 1) .AND. (mc2 .GT. mb2) )
THEN
2240 print*,
'Wrong charm-bottom mass hierarchy. STOP'
2243 IF( (
ivfns .EQ. 1) .AND. (mb2 .GT. mt2) )
THEN
2244 print*,
'Wrong bottom-top mass hierarchy. STOP'
2255 as0 = asi / (4.d0* pi)
2259 IF(
ivfns .NE. 0)
THEN
2278 real*8 m2,mur,r2,asi,asf,r20,r2t,r2b,r2c,
hto_alphas
2279 real*8,
parameter :: pi= 3.1415 9265358979d0
2282 FUNCTION hto_as(R2,R20,AS0,NF)
2295 m2 = r2 * exp(+
logfr)
2296 IF(
ivfns .EQ. 0)
THEN
2309 IF(m2 .GT.
m2t)
THEN
2315 ELSE IF(m2 .GT.
m2b)
THEN
2321 ELSE IF(m2 .GT.
m2c)
THEN
2372 INTEGER nf,prvcll,k1,k2
2373 real*8 asnf,logrh,cmci30,cmcf30,cmcf31,
2375 real*8,
dimension(3,0:3) :: cmc
2383 SAVE cmc,cmci30,cmcf30,cmcf31,cmci31,prvcll
2390 IF(prvcll .NE. 1)
THEN
2393 cmc(1,1) = 2.d0/3.d0
2395 cmc(2,0) = 14.d0/3.d0
2396 cmc(2,1) = 38.d0/3.d0
2397 cmc(2,2) = 4.d0/9.d0
2399 cmci30 = + 80507.d0/432.d0 *
zeta(3) + 58933.d0/1944.d0
2400 # + 128.d0/3.d0 *
zeta(2) * (1.+ dlog(2.d0)/3.d0)
2401 cmcf30 = - 64.d0/9.d0 * (
zeta(2) + 2479.d0/3456.d0)
2402 cmci31 = 8941.d0/27.d0
2403 cmcf31 = - 409.d0/27.d0
2404 cmc(3,2) = 511.d0/9.d0
2405 cmc(3,3) = 8.d0/27.d0
2415 cmc(3,0) = cmci30 +
nf * cmcf30
2416 cmc(3,1) = cmci31 +
nf * cmcf31
2419 IF(
naord .EQ. 0)
GO TO 1
2457 # asc3,asb4,ast5,sc,sb,st
2478 sc = log(
as0 / asc3)
2486 sb = log(
asc / asb4)
2494 st = log(
asb / ast5)
2520 FUNCTION hto_as(R2, R20,AS0,NF)
2526 real*8 r2,r20,as0,as1,as2,as3,as4,fbeta1,fbeta2,fbeta3,fbeta4,a,
2527 # lrrat,dlr,xk0,xk1,xk2,xk3,
hto_as
2528 INTEGER,
parameter :: nfmin =3
2529 INTEGER,
parameter :: nfmax =6
2530 real*8,
parameter :: sxth =0.166666666666666d0
2551 IF(
naord .EQ. 0)
THEN
2555 ELSE IF(
naord .EQ. 1)
THEN
2561 as2=
hto_as + 0.5d0 * xk0
2564 as3=
hto_as + 0.5d0 * xk1
2570 hto_as =
hto_as + sxth * (xk0 + 2.d0* xk1 + 2.d0* xk2 + xk3)
2573 ELSE IF(
naord .EQ. 2)
THEN
2580 as2=
hto_as + 0.5d0 * xk0
2584 as3=
hto_as + 0.5d0 * xk1
2592 hto_as =
hto_as + sxth * (xk0 + 2.d0* xk1 + 2.d0* xk2 + xk3)
2595 ELSE IF(
naord .EQ. 3)
THEN
2602 as2=
hto_as + 0.5d0 * xk0
2606 as3=
hto_as + 0.5d0 * xk1
2614 hto_as =
hto_as + sxth * (xk0 + 2.d0* xk1 + 2.d0* xk2 + xk3)
2649 real*8 b00,b01,b10,b11
2650 INTEGER,
parameter :: NFMIN=3
2651 INTEGER,
parameter :: NFMAX=6
2657 b00 = 11.d0/3.d0 *
ca
2658 b01 = -4.d0/3.d0 *
tr
2659 b10 = 34.d0/3.d0 *
ca**2
2660 b11 = -20.d0/3.d0 *
ca*
tr - 4.*
cf*
tr
2666 beta0(nf) = b00 + b01 * nf
2667 beta1(nf) = b10 + b11 * nf
2669 beta2(nf) = 1428.50d0 - 279.611d0 * nf + 6.01852d0 * nf**2
2670 beta3(nf) = 29243.0d0 - 6946.30d0 * nf + 405.089d0 * nf**2
2671 # + 1.49931d0 * nf**3
2682 FUNCTION hto_dzero(A0,B0,EPS,MAXF,F,MODE)
2683 IMPLICIT REAL*8 (a-h,o-z)
2695 dimension im1(2),im2(2),lmt(2)
2696 parameter(z1 = 1, half = z1/2)
2697 DATA im1 /2,3/, im2 /-1,3/
2700 IF(mode .NE. 1 .AND. mode .NE. 2)
THEN
2709 IF(fa*fb .GT. 0)
THEN
2724 3
IF(abs(fc) .LT. abs(fb))
THEN
2739 IF(abs(hb) .GT. tol)
THEN
2740 IF(ie .GT. im1(mode))
THEN
2759 IF(ie .EQ. im2(mode)) p=p+p
2760 IF(p .EQ. 0 .OR. p .LE. q*tol)
THEN
2762 ELSEIF(p .LT. hb*q)
THEN
2774 IF(mf .GT. maxf)
THEN
2776 print*,
"Error in HTO_DZERO: TOO MANY FUNCTION CALLS"
2780 IF(fb .EQ. 0 .OR. sign(z1,fc) .EQ. sign(z1,fb))
GO TO 1
2781 IF(w .EQ. hb)
GO TO 2
2788 101
FORMAT(
'Error in DZERO: MODE = ',i3,
' ILLEGAL')
2789 102
FORMAT(
'Error in DZERO: F(A) AND F(B) HAVE THE SAME SIGN, A = ',
2790 1 1p,d15.8,
', B = ',d15.8)
2799 FUNCTION hto_b0af_em(scal,psi,ps0i,xmsi,xm0i)
RESULT(value)
2810 real*8 scal,ps0i,xm0i,scals,ps0,xm0
2811 real*8,
dimension(2) ::
value,psi,ps,xmsi,xms,betasc,betas,
2812 # betac,beta,argc,arg,lbet
2820 betasc=
co-4.d0*(xms.
cq.ps)
2821 IF(betasc(2).eq.0.d0)
THEN
2823 betac= (betasc(1)).
cr.(betasc(2))
2825 betac= (betasc(1)).
crz.(betasc(2))
2827 argc= (betac+
co).
cq.(betac-
co)
2828 IF(argc(2).eq.0.d0)
THEN
2830 lbet= argc(1).
fln.argc(2)
2832 betas(1)= 1.d0-4.d0*(xm0*xm0)/ps0
2834 beta= (betas(1)).
cr.(betas(2))
2835 arg= (beta+
co).
cq.(beta-
co)
2836 IF(arg(2).eq.0.d0) arg(2)=
eps
2837 lbet= argc.
lnsrs.arg
2840 value= 2.d0*
co-(betac.
cp.lbet)
2850 FUNCTION hto_quarkqcd(scal,psi,ps0i,xmsi,xm0i,type)
RESULT(value)
2864 INTEGER it,
type,unit
2865 real*8 scal,ps0i,xm0i,scals,ps0,xm0,sgn
2866 real*8,
dimension(2) ::
value,psi,ps,xmsi,xms,betasc,betas,
2867 # betac,beta,argc,arg,lq,lm,cx,comx,cxs,x,omx,xs,clx,clomx,
2868 # clxs,li2cx,li3cx,li2cxs,li3cxs,copx,clopx,lqs,qcd,lms,
2869 # opx,tau,taus,clxx,clxxs
2870 real*8,
dimension(6,2) :: aux,auxs
2878 IF(psi(2).eq.0.d0.and.xmsi(2).eq.0.d0.
2879 # and.psi(1).le.4.d0*xmsi(1))
THEN
2885 IF(abs(ps(2)/ps(1)).lt.1.d-10.and.xms(2).eq.0.d0)
THEN
2886 betasc(1)= 1.d0-4.d0*xms(1)/ps(1)
2887 betasc(2)= 4.d0/(ps(1)*ps(1))*xms(1)*ps(2)
2889 betasc=
co-4.d0*(xms.
cq.ps)
2891 IF(betasc(2).eq.0.d0)
THEN
2893 betac= (betasc(1)).
cr.(betasc(2))
2895 betac= (betasc(1)).
crz.(betasc(2))
2897 argc= (betac+
co).
cq.(betac-
co)
2899 betas(1)= 1.d0-4.d0*xm0*xm0/ps0
2901 beta= (betas(1)).
cr.(betas(2))
2902 arg= (beta+
co).
cq.(beta-
co)
2904 IF(arg(2).eq.0.d0)
THEN
2911 x= (beta-
co).
cq.(beta+
co)
2917 IF(arg(2).eq.0.d0) arg(2)=
eps
2918 IF(argc(2).eq.0.d0)
THEN
2921 lq= argc(1).
fln.argc(2)
2933 sgn= sign(
one,cx(1))
2936 clx= cx(1).
fln.cx(2)
2938 clxx= cxs(1).
fln.cxs(2)
2940 clomx= comx(1).
fln.comx(2)
2941 clopx= copx(1).
fln.copx(2)
2943 li2cx(1:2)= aux(1,1:2)
2944 li3cx(1:2)= aux(2,1:2)
2946 li2cxs(1:2)= auxs(1,1:2)
2947 li3cxs(1:2)= auxs(2,1:2)
2948 ELSEIF(it.eq.1.d0)
THEN
2949 cx= (betac-
co).
cq.(betac+
co)
2957 clomx= comx.
lnsrs.omx
2958 clopx= copx.
lnsrs.opx
2965 IF(xms(2).eq.0.d0)
THEN
2969 lm= xms(1).
fln.xms(2)
2978 qcd= -3.d0/4.d0*(
co-12.d0*tau)*
rz2
2979 # +1.d0/16.d0*(3.d0*
co+344.d0*tau)
2980 # -3.d0/2.d0*((
co-6.d0*tau).
cp.lms)
2981 # +3.d0/4.d0*((3.d0*
co-14.d0*tau).
cp.(betac.
cp.clx))
2982 # -3.d0*((li3cxs-2.d0*li3cx+4.d0/3.d0*(li2cx.
cp.clx)
2983 # +1.d0/3.d0*(clxs.
cp.clomx)+
rz3*
co
2984 # -2.d0/3.d0*(clxx.
cp.li2cxs)-1.d0/6.d0*(clxxs.
cp.clomx)
2986 # (clxxs.
cp.clopx)).
cp.((
co-4.d0*tau).
cp.(
co-2.d0*tau)))
2987 # -1.d0/2.d0*((4.d0*li2cxs-4.d0*li2cx-4.d0*(clomx.
cp.clx)
2988 # -2.d0*(cx.
cp.(clxs.
cq.comx))+4.d0*(clxx.
cp.clomx)
2989 # +4.d0*(clxx.
cp.clopx)+(clxxs.
cp.(cxs.
cq.comx))
2990 # +(clxxs.
cp.(cxs.
cq.copx))).
cp.(betac.
cp.(
co-4.d0*tau)))
2991 # +1.d0/4.d0*(lm.
cp.(11.d0*
co-108.d0*tau))
2992 # +1.d0/4.d0*(clxs.
cp.(3.d0*
co+58.d0*taus-28.d0*tau))
2994 ELSEIF(type.eq.1)
THEN
2996 qcd= -3.d0/4.d0*(
co-12.d0*tau)*
rz2
2997 # +1.d0/16.d0*(67.d0*
co-40.d0*tau)
2998 # +3.d0/4.d0*((3.d0*
co-14.d0*tau).
cp.(betac.
cp.clx))
2999 # -3.d0*((li3cxs-2.d0*li3cx+4.d0/3.d0*(li2cx.
cp.clx)
3000 # +1.d0/3.d0*(clxs.
cp.clomx)+
rz3*
co
3001 # -2.d0/3.d0*(clxx.
cp.li2cxs)-1.d0/6.d0*(clxxs.
cp.clomx)
3003 # (clxxs.
cp.clopx)).
cp.((
co-4.d0*tau).
cp.(
co-2.d0*tau)))
3004 # -1.d0/2.d0*((4.d0*li2cxs-4.d0*li2cx-4.d0*(clomx.
cp.clx)
3005 # -2.d0*(cx.
cp.(clxs.
cq.comx))+4.d0*(clxx.
cp.clomx)
3006 # +4.d0*(clxx.
cp.clopx)+(clxxs.
cp.(cxs.
cq.comx))
3007 # +(clxxs.
cp.(cxs.
cq.copx))).
cp.(betac.
cp.(
co-4.d0*tau)))
3009 # (((
co-8.d0*tau).
cp.(betac.
cp.lq)).
cp.(
co-3.d0/4.d0*lm))
3010 # -3.d0*((betac.
cp.lq).
cp.(lm.
cp.tau))
3011 # +4.d0*((
co.
cp.betac).
cp.(lq.
cp.tau))
3012 # -9.d0/4.d0*(lm.
cp.(
co+4.d0*tau))
3013 # +9.d0*(lms.
cp.tau)
3014 # +1.d0/4.d0*(clxs.
cp.(3.d0*
co+58.d0*taus-28.d0*tau))
3039 real*8
hto_ralphas,rs0,rs,als,x1,x2,xacc,qcdl,rat,rl,rll,rb,fac,
3040 # facp,rhs,ratl2,pqcdl4,qcdb0,qcdb1,qcdb2,qcda,pqcdl3,
3042 real*8,
dimension(5) :: b0,b1,b2
3043 INTEGER,
parameter ::
nf=5
3047 IF(als.eq.0.d0)
THEN
3056 b0(i)= (11.d0-2.d0/3.d0*i)/4.d0
3057 b1(i)= (102.d0-38.d0/3.d0*i)/16.d0
3058 b2(i)= 0.5d0*(2857.d0-i*(5033.d0/9.d0-325.d0/27.d0*i))/64.d0
3065 rb= log(b0(5)/b0(4))
3066 fac= b1(5)/b0(5)-b1(4)/b0(4)
3067 facp= b2(5)/b0(5)-b2(4)/b0(4)
3068 fac2= b1(5)*b1(5)/b0(5)/b0(5)-b1(4)*b1(4)/b0(4)/b0(4)
3069 rhs= (b0(5)-b0(4))*rl+fac*rll-b1(4)/b0(4)*rb+
3070 # b1(5)/b0(5)/b0(5)*fac*rll/rl+1.d0/b0(5)/rl*(
3071 # fac2-facp-7.d0/72.d0)
3074 qcdl= qcdl/sqrt(ratl2)
3081 rb= log(b0(4)/b0(3))
3082 fac= b1(4)/b0(4)-b1(3)/b0(3)
3083 facp= b2(4)/b0(4)-b2(3)/b0(3)
3084 fac2= b1(4)*b1(4)/b0(4)/b0(4)-b1(3)*b1(3)/b0(3)/b0(3)
3085 rhs= (b0(4)-b0(3))*rl+fac*rll-b1(3)/b0(3)*rb+
3086 # b1(4)/b0(4)/b0(4)*fac*rll/rl+1.d0/b0(4)/rl*(
3087 # fac2-facp-7.d0/72.d0)
3090 qcdl= qcdl/sqrt(ratl2)
3098 qcdb0= 11.d0-2.d0/3.d0*nfe
3099 qcdb1= 102.d0-38.d0/3.d0*nfe
3100 qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nfe+325.d0/27.d0*nfe*nfe)
3101 qcda= 2.d0*log(rs/qcdl)
3103 hto_ralphas= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*
3104 # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)-
3105 # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0))
3116 real*8
hto_qcdlam,als,rs,x1,x2,xacc,fmid,f,dx,xmid
3117 INTEGER,
parameter :: jmax=50
3121 IF(f*fmid >= 0.d0) stop
3134 IF(abs(dx) < xacc.or.fmid.eq.0.d0)
RETURN
3144 qcdb0= 11.d0-2.d0/3.d0*
nf
3145 qcdb1= 102.d0-38.d0/3.d0*
nf
3146 qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*
nf+325.d0/27.d0*
nf*
nf)
3147 qcda= 2.d0*log(rs/x)
3148 hto_qcdscale= als-(4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*
3149 # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)-
3150 # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0)))
3163 real*8 scal,scal2,alsr,aexps,scaling1g,als1g,alsc4,alsc42,
3164 # als1g2,alsb,alsb2,alsz,alsz2,fn,b03,b13,b23,g03,g13,
3165 # g23,rex3,b03s,b03c,b13s,sm1g,cfm13,cfm23,rmsc,
3166 # b04,b04c,b14,b24,g04,g14,g24,rex4,b04s,b14s,asldf,cfm14,
3167 # cfm24,rsmb,zero,x1,x2,xacc,cmm4,cmb,b05,b15,b25,g05,g15,
3168 # g25,rex5,b05s,b05c,b15s,cfm15,cfm25,rcqm,bmm5,rbqm,scal1g,
3170 real*8,
dimension(2) ::
value
3192 b03= (11.d0-2.d0/3.d0*fn)/4.d0
3193 b13= (102.d0-38.d0/3.d0*fn)/16.d0
3194 b23= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0*fn*fn)/64.d0
3196 g13= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0
3197 g23= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*
rz3)*fn-
3198 # 140.d0/81.d0*fn*fn)/64.d0
3205 cfm13= g13/b03-b13*g03/b03s
3206 cfm23= g23/b03-b13*g13/b03s-b23*g03/b03s+b13s*g03/b03c
3207 rsmc= sm1g*(alsc4/als1g)**rex3*(1.d0+cfm13*
3208 # asldf+0.5d0*cfm13*cfm13*asldf*asldf+0.5d0*cfm23*
3212 b04= (11.d0-2.d0/3.d0*fn)/4.d0
3213 b14= (102.d0-38.d0/3.d0*fn)/16.d0
3214 b24= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0*fn*fn)/64.d0
3216 g14= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0
3217 g24= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*
rz3)*fn-
3218 # 140.d0/81.d0*fn*fn)/64.d0
3224 cfm14= g14/b04-b14*g04/b04s
3225 cfm24= g24/b04-b14*g14/b04s-b24*g04/b04s+
3227 rsmb= rsmc*(alsb/alsc4)**rex4*(1.d0+cfm14*
3228 # asldf+0.5d0*cfm14*cfm14*asldf*asldf+0.5d0*cfm24*
3242 cmb= cmm4*(alsb/alsc4)**rex4*(1.d0+cfm14*
3243 # (alsb-alsc4)+0.5d0*cfm14*cfm14*
3244 # (alsb-alsc4)**2+0.5d0*cfm24*(alsb2-alsc42))
3249 b05= (11.d0-2.d0/3.d0*fn)/4.d0
3250 b15= (102.d0-38.d0/3.d0*fn)/16.d0
3251 b25= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0*fn*fn)/64.d0
3253 g15= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0
3254 g25= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*
rz3)*fn-
3255 # 140.d0/81.d0*fn*fn)/64.d0
3260 cfm15= g15/b05-b15*g05/b05s
3261 cfm25= g25/b05-b15*g15/b05s-b25*g05/b05s+b15s*g05/b05c
3262 rcqm= cmb*(alsz/alsb)**rex5*(1.d0+cfm15*
3263 # (alsz-alsb)+0.5d0*cfm15*cfm15*
3264 # (alsz-alsb)**2+0.5d0*cfm25*(alsz2-alsb2))
3277 rbqm= bmm5*(alsz/alsb)**rex5*(1.d0+cfm15*
3278 # asldf+0.5d0*cfm15*cfm15*asldf*asldf+0.5d0*cfm25*
3291 FUNCTION hto_rrunm(x1,x2,xacc,qm,als,rm1,rm2,fn)
3294 INTEGER,
parameter :: jmax=50
3296 real*8
hto_rrunm,x1,x2,xacc,qm,als,rm1,rm2,fn,fmid,f,dx,xmid
3300 IF(f*fmid.ge.0.d0)
then
3301 print*,
'root must be bracketed for bisection'
3303 1
format(/
' error detected by HTO_rrunm ',/
3304 #
' current value of quark mass = ',e20.5)
3319 IF(abs(dx) < xacc.or.fmid.eq.0.d0)
RETURN
3331 real*8
hto_qcdmass,qm,als,rm1,rm2,fn,x,rln,rlns,r1,r2,delta0,
3338 delta0= 3.d0/4.d0*
rz2-3.d0/8.d0
3339 delta1= r1*(
pis/8.d0+r1*(-0.597d0+0.230d0*r1))
3340 delta2= r2*(
pis/8.d0+r2*(-0.597d0+0.230d0*r2))
3342 rhs= 1.d0+als*(4.d0/3.d0+rln+als*(3817.d0/288.d0-8.d0/3.d0+
3343 # 2.d0/3.d0*(2.d0+log(2.d0))*
rz2-
rz3/6.d0-fn/3.d0*(
rz2+
3344 # 71.d0/48.d0)+4.d0/3.d0*(delta0+delta1+delta2)+(173.d0/
3345 # 24.d0-13.d0/36.d0*fn)*rln+(15.d0/8.d0-fn/12.d0)*rlns))
3372 real*8 scal,scals,ps0i,ps0,xm1,xm2
3373 real*8,
intent(in),
dimension(2) :: xm0i
3374 real*8,
intent(in),
dimension(2,2) :: xmsi
3375 real*8,
dimension(2) ::
value,xm1c,xm2c
3376 real*8,
dimension(2,2) :: xms
3377 real*8,
dimension(2) :: xm0
3378 real*8,
dimension(2) :: psi,ps,aroot,root,lambdasc,lambdas,
3379 # lambdac,lambda,argc,arg,llam,l1,l2
3386 xm1c(1:2)= xms(1,1:2)
3387 xm2c(1:2)= xms(2,1:2)
3390 xm1= xm0(1)*xm0(1)/ps0
3391 xm2= xm0(2)*xm0(2)/ps0
3394 root= (aroot(1).
crz.aroot(2))
3396 lambdasc=
co+(xm1c.
cp.xm1c)+(xm2c.
cp.xm2c)-2.d0*(
3397 # xm1c+xm2c+(xm1c.
cp.xm2c))
3398 lambdas(1)= 1.d0+xm1*xm1+xm2*xm2-2.d0*(
3401 lambdac= (lambdasc(1)).
crz.(lambdasc(2))
3402 lambda= (lambdas(1)).
cr.(lambdas(2))
3403 IF(lambda(2).eq.0.d0) lambda(2)= -
eps
3405 argc= 0.5d0*((-
co+xm1c+xm2c-lambdac).
cq.root)
3407 arg(1)= 0.5d0*(-1.d0+xm1+xm2-lambda(1))/sqrt(xm1*xm2)
3410 llam= argc.
lnsrs.arg
3412 l1= xm1c(1).
fln.xm1c(2)
3413 l2= xm2c(1).
fln.xm2c(2)
3415 value= (((xm1c-xm2c-
co).
cq.lambdac).
cp.llam)+0.5d0*(l1-l2)
3427 FUNCTION hto_zeroin(f, ax, bx, aerr, rerr)
RESULT(fn_val)
3461 real*8,
INTENT(IN) :: ax
3462 real*8,
INTENT(IN) :: bx
3463 real*8,
INTENT(IN) :: aerr
3464 real*8,
INTENT(IN) :: rerr
3470 FUNCTION f(x)
RESULT(fn_val)
3472 real*8,
INTENT(IN) :: x
3477 real*8 :: a,b,c,d,e,eps,fa,fb,fc,tol,xm,p,q,r,s,atol,rtol
3489 IF(fa*fb.gt.0.d0)
THEN
3493 rtol= max(0.5d0*rerr, 2.0d0*eps)
3501 20
IF(abs(fc) < abs(fb))
THEN
3512 tol= rtol*max(abs(b),abs(c))+atol
3514 IF(abs(xm) > tol)
THEN
3515 IF(fb /= 0.0d0)
THEN
3519 IF(abs(e) >= tol)
THEN
3520 IF(abs(fa) > abs(fb))
THEN
3538 p= s*((c-b)*q*(q-r)-(b-a)*(r-1.0d0))
3539 q= (q-1.0d0)*(r-1.0d0)*(s-1.0d0)
3549 IF(2.0*p < (3.0*xm*q-abs(tol*q)))
THEN
3550 IF(p < abs(0.5d0*e*q))
THEN
3568 IF(abs(d) > tol) b= b+d
3569 IF(abs(d) <= tol) b= b+sign(tol,xm)
3571 IF(fb*(fc/abs(fc)) > 0.0d0)
GO TO 10
3597 real*8,
INTENT(IN) :: x
3610 FUNCTION hto_shh(muhr,scal,rgh)
3632 real*8
hto_shh,muhr,rgh,muhs,scal,scals,p2,xm0,str,sti,ewc,
3633 # sconv,asmur,emc,emb,emt,as_nlo,rmbs,rmcs,rmss,
3634 # crmbs,crmcs,lcxb,lcxc,lcxbs,lcxcs,lclxb,lclxc,qcdtop,neg
3635 real*8,
dimension(2) :: axm0
3636 real*8,
dimension(2,2) :: axms
3637 real*8,
dimension(3,2,2) :: bmcb
3638 real*8,
dimension(8,2,2) :: bmcf
3639 real*8,
dimension(1,2,2) :: bmct
3640 real*8,
dimension(3,2) :: bm0b,coefbb
3641 real*8,
dimension(8,2) :: bm0f,coefbf
3642 real*8,
dimension(1,2) :: bm0t,coefbt
3643 real*8,
dimension(2) :: sh,shs,clh,b0sumb,b0sumf,cxp,ksumb,
3644 # ksumf,coefb1,coefb2,coefb3,coefb4,coefb5,coefb6,coefb7,
3645 # coefb8,coefb9,coefb10,coefb11,coefb12,ab2,ab3,
3646 # totalf,totalt,totalb,total,b0sumt,ksumt,xms,
3647 # b0part,cpt,ccxt,deltag,sww0,coefw1,coefw2,
3648 # ksumw,ccxts,cctsi,shi,clt,cxhw,cstsi,csmcts,cltmw,
3649 # b0sumw,sww,dw,b0sumw1,b0sumw2,cmxw,clmw,cxtw,
3650 # kfmsb,ktmsb,kbmsb,nloqcd,runbc,ttqcd,wccxt,wclt,wccxts,
3651 # wcxtw,wcltmw,ascal
3690 real*8 scal,ps0i,xm0i
3691 real*8,
dimension(2) ::
value,psi,xmsi
3708 crmbs= runbc(2)*runbc(2)
3709 crmcs= runbc(1)*runbc(1)
3722 ELSEIF(
gtop == 0)
THEN
3733 wccxts= wccxt.
cp.wccxt
3742 clmw= cmxw(1).
fln.cmxw(2)
3743 clmw(2)= clmw(2)-2.d0*pi
3746 sh(2)= -muhr*rgh/scals
3752 clh= sh(1).
fln.sh(2)
3753 clt= ccxt(1).
fln.ccxt(2)
3754 wclt= wccxt(1).
fln.wccxt(2)
3757 cltmw= cxtw(1).
fln.cxtw(2)
3758 wcltmw= wcxtw(1).
fln.wcxtw(2)
3761 coefb1= (12.d0*
cxw-4.d0*sh+(shs.
cq.
cxw))/64.d0
3768 coefb3= 9.d0/128.d0*(shs.
cq.
cxw)
3771 coefb4= -3.d0/32.d0*((4.d0*ccxt-sh).
cp.(ccxt.
cq.
cxw))
3775 coefb5= -3.d0/32.d0*((4.d0*
co*lcxb-sh).
cq.
cxw)*lcxb
3779 coefb7= -3.d0/32.d0*((4.d0*
co*lcxc-sh).
cq.
cxw)*lcxc
3801 b0sumb= (coefb1.
cp.b0part)
3806 b0sumb= b0sumb+(coefb2.
cp.b0part)
3811 b0sumb= b0sumb+(coefb3.
cp.b0part)
3813 xms(1:2)= crmbs*
co(1:2)
3816 b0sumf= (coefb5.
cp.b0part)
3820 b0sumf= b0sumf+(coefb6.
cp.b0part)
3821 xms(1:2)= crmcs*
co(1:2)
3824 b0sumf= b0sumf+(coefb7.
cp.b0part)
3828 b0sumf= b0sumf+(coefb8.
cp.b0part)
3832 b0sumf= b0sumf+(coefb9.
cp.b0part)
3836 b0sumf= b0sumf+(coefb10.
cp.b0part)
3840 b0sumf= b0sumf+(coefb11.
cp.b0part)
3844 b0sumf= b0sumf+(coefb12.
cp.b0part)
3849 b0sumt= (coefb4.
cp.b0part)
3852 xms(1:2)= crmbs*
co(1:2)
3855 xms(1:2)= crmcs*
co(1:2)
3857 nloqcd= nloqcd+crmcs*
hto_quarkqcd(scal,cxp,p2,xms,xm0,iz)
3858 IF(
qcdc==0) nloqcd= 0.d0
3865 IF(
qcdc==0) ttqcd= 0.d0
3867 ELSEIF(
ifb.eq.1)
THEN
3869 xms(1:2)= crmbs*
co(1:2)
3872 b0sumf= (coefb5.
cp.b0part)
3876 b0sumf= b0sumf+(coefb6.
cp.b0part)
3877 xms(1:2)= crmcs*
co(1:2)
3880 b0sumf= b0sumf+(coefb7.
cp.b0part)
3884 b0sumf= b0sumf+(coefb8.
cp.b0part)
3888 b0sumf= b0sumf+(coefb9.
cp.b0part)
3892 b0sumf= b0sumf+(coefb10.
cp.b0part)
3896 b0sumf= b0sumf+(coefb11.
cp.b0part)
3900 b0sumf= b0sumf+(coefb12.
cp.b0part)
3903 xms(1:2)= crmbs*
co(1:2)
3906 xms(1:2)= crmcs*
co(1:2)
3908 nloqcd= nloqcd+crmcs*
hto_quarkqcd(scal,cxp,p2,xms,xm0,iz)
3909 IF(
qcdc==0) nloqcd= 0.d0
3911 ELSEIF(
ifb.eq.2)
THEN
3916 b0sumt= (coefb4.
cp.b0part)
3923 IF(
qcdc==0) ttqcd= 0.d0
3925 ELSEIF(
ifb.eq.3)
THEN
3931 b0sumb= (coefb1.
cp.b0part)
3936 b0sumb= b0sumb+(coefb2.
cp.b0part)
3941 b0sumb= b0sumb+(coefb3.
cp.b0part)
3945 IF(
ifb.eq.0.or.
ifb.eq.1)
THEN
3948 #
cxu*
clxu+lcxc*lclxc))/32.d0+
3952 IF(
ifb.eq.0.or.
ifb.eq.2)
THEN
3953 ksumt= 3.d0/32.d0*(4.d0*((ccxt.
cq.
cxw).
cp.ccxt)
3957 IF(
ifb.eq.0.or.
ifb.eq.3)
THEN
3959 ksumb= (-2.d0*((2.d0*
co+cctsi+3.d0*(sh.
cq.
cxw)).
cp.sh)
3962 # -3.d0*(clh.
cp.(shs.
cq.
cxw)))/128.d0
3966 ewc= 4.d0*sqrt(2.d0)*
g_f/
pis
3969 IF(
ifb.eq.0.or.
ifb.eq.1)
THEN
3970 totalf= b0sumf+ksumf
3973 IF(
ifb.eq.0.or.
ifb.eq.2)
THEN
3974 totalt= b0sumt+ksumt
3976 neg= ewc*(
swr*totalt(2)+
swi*totalt(1))/sconv
3977 # -as_nlo*
g_f/(sqrt(2.d0)*
pis)*
3978 # (sh(1)*ttqcd(2)+sh(2)*ttqcd(1))/sconv
3985 qcdtop= as_nlo*
g_f/(sqrt(2.d0)*
pis)*
3986 # (sh(1)*ttqcd(2)+sh(2)*ttqcd(1))/sconv
3992 IF(
ifb.eq.0.or.
ifb.eq.3)
THEN
3993 totalb= b0sumb+ksumb
3994 IF((
swr*totalb(2)+
swi*totalb(1)) < 0.d0)
THEN
4002 total= totalf+totalt+totalb
4008 sww0= -(38.d0*
cxw+6.d0*wccxt+7.d0*sh
4009 # -48.d0*(((wccxt.
cq.sh).
cq.
cxw).
cp.wccxt)+8.d0*(
cxw.
cq.sh))/128.d0
4011 # +3.d0/32.d0*(((
co-4.d0*((wccxt.
cq.sh).
cq.
cxw).
cp.wccxt)).
cp.wclt)
4012 # +((((8.d0*
co-17.d0*cstsi+3.d0*cctsi).
cp.
cxw)
4016 coefw1= -(((8.d0*
co-(sh.
cq.
cxw)).
cp.sh)*sh
4017 # -4.d0*((-12.d0*
cxw+7.d0*sh).
cp.
cxw))/192.d0
4030 axms(2,2)= -muhr*rgh
4033 b0sumw1= (coefw1.
cp.b0part)
4034 b0sumw= (coefw1.
cp.b0part)
4043 b0sumw2= (coefw2.
cp.b0part)
4044 b0sumw= b0sumw+(coefw2.
cp.b0part)
4047 # -0.5d0*((3.d0*
co-(wccxts.
cq.
cxws)).
cp.wccxt)).
cp.wcltmw)
4049 # +((36.d0*
cxw-14.d0*sh-18.d0*((
co-4.d0*(wccxt.
cq.sh)).
cp.wccxt)
4051 # -6.d0*(((2.d0*
co+((
cxwi-12.d0*shi).
cp.wccxt)).
cp.wccxt)
4052 # +1.d0/6.d0*((15.d0*
co-(sh.
cq.
cxw)).
cp.sh)
4056 # -48.d0*
csts-5.d0*cctsi))
4058 # -72.d0*((wclt.
cp.wccxts).
cp.(shi-1.d0/12.d0*(wccxt.
cq.
cxws)))
4061 ksumw= ksumw/192.d0+3.d0/16.d0*(
cxw.
cp.(
clw-clmw))
4065 dw= (-sww+sww0)/sconv+deltag/16.d0
4070 # (1.d0+ewc*(
swr*dw(1)-
swi*dw(2)))
4071 # -ewc*(
swr*total(2)+
swi*total(1))/sconv
4072 # +as_nlo*
g_f/(sqrt(2.d0)*
pis)*
4073 # (sh(1)*nloqcd(2)+sh(2)*nloqcd(1))/sconv
4076 ELSEIF(
ifb.eq.1)
THEN
4078 # (1.d0+ewc*(
swr*dw(1)-
swi*dw(2)))
4079 # -ewc*(
swr*totalf(2)+
swi*totalf(1))/sconv
4080 # +as_nlo*
g_f/(sqrt(2.d0)*
pis)*
4081 # (sh(1)*nloqcd(2)+sh(2)*nloqcd(1))/sconv
4082 ELSEIF(
ifb.eq.2)
THEN
4084 # (1.d0+ewc*(
swr*dw(1)-
swi*dw(2)))
4085 # -ewc*(
swr*totalt(2)+
swi*totalt(1))/sconv
4087 ELSEIF(
ifb.eq.3)
THEN
4089 # (1.d0+ewc*(
swr*dw(1)-
swi*dw(2)))
4090 # -ewc*(
swr*totalb(2)+
swi*totalb(1))/sconv
4110 real*8 u,
value,evalue,mass
4111 real*8,
dimension(103) :: bc,cc,dc
4134 INTEGER k,n,i,top,gdim,l
4136 real*8,
dimension(103) :: xc,yc
4137 real*8,
dimension(103) :: x,y
4139 real*8,
DIMENSION(gdim) :: b
4142 real*8,
DIMENSION(gdim) :: c
4145 real*8,
DIMENSION(gdim) :: d
4149 real*8,
PARAMETER:: zero=0.0, two=2.0, three=3.0
4154 DATA (xc(i),i=1,103)/
4155 # 90.d0,95.d0,100.d0,105.d0,110.d0,115.d0,120.d0,
4156 # 125.d0,130.d0,135.d0,140.d0,145.d0,150.d0,155.d0,160.d0,165.d0,
4157 # 170.d0,175.d0,180.d0,185.d0,190.d0,195.d0,200.d0,210.d0,220.d0,
4158 # 230.d0,240.d0,250.d0,260.d0,270.d0,280.d0,290.d0,300.d0,310.d0,
4159 # 320.d0,330.d0,340.d0,350.d0,360.d0,370.d0,380.d0,390.d0,400.d0,
4160 # 410.d0,420.d0,430.d0,440.d0,450.d0,460.d0,470.d0,480.d0,490.d0,
4161 # 500.d0,510.d0,520.d0,530.d0,540.d0,550.d0,560.d0,570.d0,580.d0,
4162 # 590.d0,600.d0,610.d0,620.d0,630.d0,640.d0,650.d0,660.d0,670.d0,
4163 # 680.d0,690.d0,700.d0,710.d0,720.d0,730.d0,740.d0,750.d0,760.d0,
4164 # 770.d0,780.d0,790.d0,800.d0,810.d0,820.d0,830.d0,840.d0,850.d0,
4165 # 860.d0,870.d0,880.d0,890.d0,900.d0,910.d0,920.d0,930.d0,940.d0,
4166 # 950.d0,960.d0,970.d0,980.d0,990.d0,1000.d0/
4168 DATA (yc(i),i=1,103)/
4169 # 2.20d-3,2.32d-3,2.46d-3,2.62d-3,2.82d-3,3.09d-3,3.47d-3,4.03d-3,
4170 # 4.87d-3,6.14d-3,8.12d-3,1.14d-2,1.73d-2,3.02d-2,8.29d-2,2.46d-1,
4171 # 3.80d-1,5.00d-1,6.31d-1,8.32d-1,1.04d0,1.24d0,1.43d0,1.85d0,
4172 # 2.31d0,2.82d0,3.40d0,4.04d0,4.76d0,5.55d0,6.43d0,7.39d0,8.43d0,
4173 # 9.57d0,10.8d0,12.1d0,13.5d0,15.2d0,17.6d0,20.2d0,23.1d0,26.1d0,
4174 # 29.2d0,32.5d0,35.9d0,39.4d0,43.1d0,46.9d0,50.8d0,54.9d0,59.1d0,
4175 # 63.5d0,68.0d0,72.7d0,77.6d0,82.6d0,87.7d0,93.1d0,98.7d0,104.d0,
4176 # 110.d0,116.d0,123.d0,129.d0,136.d0,143.d0,150.d0,158.d0,166.d0,
4177 # 174.d0,182.d0,190.d0,199.d0,208.d0,218.d0,227.d0,237.d0,248.d0,
4178 # 258.d0,269.d0,281.d0,292.d0,304.d0,317.d0,330.d0,343.d0,357.d0,
4179 # 371.d0,386.d0,401.d0,416.d0,432.d0,449.d0,466.d0,484.d0,502.d0,
4180 # 521.d0,540.d0,560.d0,581.d0,602.d0,624.d0,647.d0/
4192 c(2)= (y(2)-y(1))/d(1)
4195 b(k)= two*(d(k-1)+d(k))
4196 c(k+1)= (y(k+1)-y(k))/d(k)
4208 c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
4209 c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
4210 c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
4211 c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
4225 c(k)= (c(k)-d(k)*c(k+1))/b(k)
4230 b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
4232 b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
4233 d(k)= (c(k+1)-c(k))/d(k)
4249 real*8,
INTENT(IN) :: u
4252 INTEGER j,k,n,l,top,gdim
4254 real*8,
dimension(103) :: xc,yc
4255 real*8,
dimension(103) :: x,y
4256 real*8,
DIMENSION(gdim) :: b,c,d
4259 real*8,
INTENT(OUT),
OPTIONAL:: f,fp,fpp,fppp
4262 INTEGER,
SAVE :: i=1
4264 real*8,
PARAMETER:: two=2.0, three=3.0, six=6.0
4268 DATA (xc(l),l=1,103)/
4269 # 90.d0,95.d0,100.d0,105.d0,110.d0,115.d0,120.d0,
4270 # 125.d0,130.d0,135.d0,140.d0,145.d0,150.d0,155.d0,160.d0,165.d0,
4271 # 170.d0,175.d0,180.d0,185.d0,190.d0,195.d0,200.d0,210.d0,220.d0,
4272 # 230.d0,240.d0,250.d0,260.d0,270.d0,280.d0,290.d0,300.d0,310.d0,
4273 # 320.d0,330.d0,340.d0,350.d0,360.d0,370.d0,380.d0,390.d0,400.d0,
4274 # 410.d0,420.d0,430.d0,440.d0,450.d0,460.d0,470.d0,480.d0,490.d0,
4275 # 500.d0,510.d0,520.d0,530.d0,540.d0,550.d0,560.d0,570.d0,580.d0,
4276 # 590.d0,600.d0,610.d0,620.d0,630.d0,640.d0,650.d0,660.d0,670.d0,
4277 # 680.d0,690.d0,700.d0,710.d0,720.d0,730.d0,740.d0,750.d0,760.d0,
4278 # 770.d0,780.d0,790.d0,800.d0,810.d0,820.d0,830.d0,840.d0,850.d0,
4279 # 860.d0,870.d0,880.d0,890.d0,900.d0,910.d0,920.d0,930.d0,940.d0,
4280 # 950.d0,960.d0,970.d0,980.d0,990.d0,1000.d0/
4282 DATA (yc(l),l=1,103)/
4283 # 2.20d-3,2.32d-3,2.46d-3,2.62d-3,2.82d-3,3.09d-3,3.47d-3,4.03d-3,
4284 # 4.87d-3,6.14d-3,8.12d-3,1.14d-2,1.73d-2,3.02d-2,8.29d-2,2.46d-1,
4285 # 3.80d-1,5.00d-1,6.31d-1,8.32d-1,1.04d0,1.24d0,1.43d0,1.85d0,
4286 # 2.31d0,2.82d0,3.40d0,4.04d0,4.76d0,5.55d0,6.43d0,7.39d0,8.43d0,
4287 # 9.57d0,10.8d0,12.1d0,13.5d0,15.2d0,17.6d0,20.2d0,23.1d0,26.1d0,
4288 # 29.2d0,32.5d0,35.9d0,39.4d0,43.1d0,46.9d0,50.8d0,54.9d0,59.1d0,
4289 # 63.5d0,68.0d0,72.7d0,77.6d0,82.6d0,87.7d0,93.1d0,98.7d0,104.d0,
4290 # 110.d0,116.d0,123.d0,129.d0,136.d0,143.d0,150.d0,158.d0,166.d0,
4291 # 174.d0,182.d0,190.d0,199.d0,208.d0,218.d0,227.d0,237.d0,248.d0,
4292 # 258.d0,269.d0,281.d0,292.d0,304.d0,317.d0,330.d0,343.d0,357.d0,
4293 # 371.d0,386.d0,401.d0,416.d0,432.d0,449.d0,466.d0,484.d0,502.d0,
4294 # 521.d0,540.d0,560.d0,581.d0,602.d0,624.d0,647.d0/
4305 IF ( (i<1) .OR. (i >= n) ) i=1
4306 IF ( (u < x(i)) .OR. (u >= x(i+1)) )
THEN
4327 IF (
Present(f)) f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
4328 IF (
Present(fp)) fp= b(i)+dx*(two*c(i) + dx*three*d(i))
4329 IF (
Present(fpp)) fpp= two*c(i) + dx*six*d(i)
4330 IF (
Present(fppp)) fppp= six*d(i)
4347 real*8 u,
value,evalue,mass
4348 real*8,
dimension(22) :: bc,cc,dc
4372 INTEGER k,n,i,top,gdim,l
4374 real*8,
dimension(22) :: xc,yc
4375 real*8,
dimension(22) :: x,y
4377 real*8,
DIMENSION(gdim) :: b
4380 real*8,
DIMENSION(gdim) :: c
4383 real*8,
DIMENSION(gdim) :: d
4387 real*8,
PARAMETER:: zero=0.0, two=2.0, three=3.0
4392 DATA (xc(i),i=1,22)/
4393 # 190.d0,191.d0,192.d0,193.d0,194.d0,195.d0,196.d0,197.d0,198.d0,
4394 # 199.d0,200.d0,240.d0,241.d0,242.d0,243.d0,244.d0,245.d0,246.d0,
4395 # 247.d0,248.d0,249.d0,250.d0/
4397 DATA (yc(i),i=1,22)/
4398 # 0.10346d1,0.10750d1,0.11151d1,0.11549d1,0.11942d1,0.12329d1,
4399 # 0.12708d1,0.13083d1,0.13456d1,0.13832d1,0.14212d1,
4400 # 0.32401d1,0.32994d1,0.33593d1,0.34200d1,0.34813d1,0.35433d1,
4401 # 0.36061d1,0.36695d1,0.37336d1,0.37984d1,0.38640d1/
4413 c(2)= (y(2)-y(1))/d(1)
4416 b(k)= two*(d(k-1)+d(k))
4417 c(k+1)= (y(k+1)-y(k))/d(k)
4429 c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
4430 c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
4431 c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
4432 c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
4446 c(k)= (c(k)-d(k)*c(k+1))/b(k)
4451 b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
4453 b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
4454 d(k)= (c(k+1)-c(k))/d(k)
4470 real*8,
INTENT(IN) :: u
4473 INTEGER j,k,n,l,top,gdim
4475 real*8,
dimension(22) :: xc,yc
4476 real*8,
dimension(22) :: x,y
4477 real*8,
DIMENSION(gdim) :: b,c,d
4480 real*8,
INTENT(OUT),
OPTIONAL:: f,fp,fpp,fppp
4483 INTEGER,
SAVE :: i=1
4485 real*8,
PARAMETER:: two=2.0, three=3.0, six=6.0
4489 DATA (xc(l),l=1,22)/
4490 # 190.d0,191.d0,192.d0,193.d0,194.d0,195.d0,196.d0,197.d0,198.d0,
4491 # 199.d0,200.d0,240.d0,241.d0,242.d0,243.d0,244.d0,245.d0,246.d0,
4492 # 247.d0,248.d0,249.d0,250.d0/
4494 DATA (yc(l),l=1,22)/
4495 # 0.10346d1,0.10750d1,0.11151d1,0.11549d1,0.11942d1,0.12329d1,
4496 # 0.12708d1,0.13083d1,0.13456d1,0.13832d1,0.14212d1,
4497 # 0.32401d1,0.32994d1,0.33593d1,0.34200d1,0.34813d1,0.35433d1,
4498 # 0.36061d1,0.36695d1,0.37336d1,0.37984d1,0.38640d1/
4509 IF ( (i<1) .OR. (i >= n) ) i=1
4510 IF ( (u < x(i)) .OR. (u >= x(i+1)) )
THEN
4531 IF (
Present(f)) f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
4532 IF (
Present(fp)) fp= b(i)+dx*(two*c(i) + dx*three*d(i))
4533 IF (
Present(fpp)) fpp= two*c(i) + dx*six*d(i)
4534 IF (
Present(fppp)) fppp= six*d(i)
4544 FUNCTION hto_deriv(scal,rhm)
RESULT(value)
4557 real*8 scal,scals,rhm,ps,xw,xz,lw,lz,ls,xe,lxe,xmu,lxmu,xtau,
4558 # lxtau,xu,lxu,xd,lxd,xc,lxc,xs,lxs,xt,lxt,xb,lxb,rmbs,rmcs,
4559 # xes,xmus,xtaus,xus,xds,xcs,xss,xts,xbs,
4560 # xec,xmuc,xtauc,xuc,xdc,xcc,xsc,xtc,xbc,
4562 # bxeis,bxmuis,bxtauis,bxuis,bxdis,bxcis,bxsis,bxtis,bxbis,
4563 # bxes,bxmus,bxtaus,bxus,bxds,bxcs,bxss,bxts,bxbs,
4565 real*8,
dimension(2) :: bs,bw,bz,arg,lbw,lbz,dmz,bzi,bh,lbh,bwi,
4566 # dmw,dmh,sh,bxe,lbxe,bxmu,lbxmu,bxtau,lbxtau,bxu,lbxu,bxd,
4567 # lbxd,bxc,lbxc,bxs,lbxs,bxt,lbxt,bxb,lbxb,runbc,bwic,bzic,
4568 # d2mh,dmwh,dmzh,dmzw,dmww,dmzz,
4569 # bxei,bxmui,bxtaui,bxui,bxdi,bxci,bxsi,bxti,bxbi
4570 real*8,
dimension(10,2) ::
value
4602 fw= 1.d0-4.d0*xw*(1.d0-3.d0*xw)
4603 fz= 1.d0-4.d0*xz*(1.d0-3.d0*xz)
4613 rmbs= runbc(2)*runbc(2)
4614 rmcs= runbc(1)*runbc(1)
4657 bh= (bs(1)).
cr.(bs(2))
4659 IF(arg(2).eq.0.d0) arg(2)=
eps
4660 lbh= arg(1).
fln.arg(2)
4665 bs(1)= 1.d0-4.d0*(
mw*
mw)/ps
4667 bw= (bs(1)).
cr.(bs(2))
4668 IF(bw(2).eq.0.d0) bw(2)= -
eps
4670 IF(bwi(2).eq.0.d0) bwi(2)=
eps
4671 bwis= 1.d0/(1.d0-4.d0*(
mw*
mw)/ps)
4674 IF(arg(2).eq.0.d0) arg(2)=
eps
4675 lbw= arg(1).
fln.arg(2)
4679 bs(1)= 1.d0-4.d0*(
mz*
mz)/ps
4681 bz= (bs(1)).
cr.(bs(2))
4682 IF(bz(2).eq.0.d0) bz(2)= -
eps
4684 IF(bzi(2).eq.0.d0) bzi(2)=
eps
4685 bzis= 1.d0/(1.d0-4.d0*(
mz*
mz)/ps)
4688 IF(arg(2).eq.0.d0) arg(2)=
eps
4689 lbz= arg(1).
fln.arg(2)
4693 bs(1)= 1.d0-4.d0*(
me*
me)/ps
4695 bxe= (bs(1)).
cr.(bs(2))
4696 arg= (bxe+
co).
cq.(bxe-
co)
4697 IF(arg(2).eq.0.d0) arg(2)=
eps
4698 lbxe= arg(1).
fln.arg(2)
4700 bs(1)= 1.d0-4.d0*(
mm*
mm)/ps
4702 bxmu= (bs(1)).
cr.(bs(2))
4703 arg= (bxmu+
co).
cq.(bxmu-
co)
4704 IF(arg(2).eq.0.d0) arg(2)=
eps
4705 lbxmu= arg(1).
fln.arg(2)
4707 bs(1)= 1.d0-4.d0*(
mtl*
mtl)/ps
4709 bxtau= (bs(1)).
cr.(bs(2))
4710 arg= (bxtau+
co).
cq.(bxtau-
co)
4711 IF(arg(2).eq.0.d0) arg(2)=
eps
4712 lbxtau= arg(1).
fln.arg(2)
4714 bs(1)= 1.d0-4.d0*(
muq*
muq)/ps
4716 bxu= (bs(1)).
cr.(bs(2))
4717 arg= (bxu+
co).
cq.(bxu-
co)
4718 IF(arg(2).eq.0.d0) arg(2)=
eps
4719 lbxu= arg(1).
fln.arg(2)
4721 bs(1)= 1.d0-4.d0*(
mdq*
mdq)/ps
4723 bxd= (bs(1)).
cr.(bs(2))
4724 arg= (bxd+
co).
cq.(bxd-
co)
4725 IF(arg(2).eq.0.d0) arg(2)=
eps
4726 lbxd= arg(1).
fln.arg(2)
4728 bs(1)= 1.d0-4.d0*rmcs/ps
4730 bxc= (bs(1)).
cr.(bs(2))
4731 arg= (bxc+
co).
cq.(bxc-
co)
4732 IF(arg(2).eq.0.d0) arg(2)=
eps
4733 lbxc= arg(1).
fln.arg(2)
4735 bs(1)= 1.d0-4.d0*(
msq*
msq)/ps
4737 bxs= (bs(1)).
cr.(bs(2))
4738 arg= (bxs+
co).
cq.(bxs-
co)
4739 IF(arg(2).eq.0.d0) arg(2)=
eps
4740 lbxs= arg(1).
fln.arg(2)
4742 bs(1)= 1.d0-4.d0*(
mt*
mt)/ps
4744 bxt= (bs(1)).
cr.(bs(2))
4745 arg= (bxt+
co).
cq.(bxt-
co)
4746 IF(arg(2).eq.0.d0) arg(2)=
eps
4747 lbxt= arg(1).
fln.arg(2)
4749 bs(1)= 1.d0-4.d0*rmbs/ps
4751 bxb= (bs(1)).
cr.(bs(2))
4752 arg= (bxb+
co).
cq.(bxb-
co)
4753 IF(arg(2).eq.0.d0) arg(2)=
eps
4754 lbxb= arg(1).
fln.arg(2)
4765 bxeis= 1.d0/(1.d0-4.d0*(
me*
me)/ps)
4766 bxmuis= 1.d0/(1.d0-4.d0*(
mm*
mm)/ps)
4767 bxtauis= 1.d0/(1.d0-4.d0*(
mtl*
mtl)/ps)
4768 bxuis= 1.d0/(1.d0-4.d0*(
muq*
muq)/ps)
4769 bxdis= 1.d0/(1.d0-4.d0*(
mdq*
mdq)/ps)
4770 bxcis= 1.d0/(1.d0-4.d0*rmcs/ps)
4771 bxsis= 1.d0/(1.d0-4.d0*(
msq*
msq)/ps)
4772 bxtis= 1.d0/(1.d0-4.d0*(
mt*
mt)/ps)
4773 bxbis= 1.d0/(1.d0-4.d0*rmbs/ps)
4775 bxes= 1.d0-4.d0*(
me*
me)/ps
4776 bxmus= 1.d0-4.d0*(
mm*
mm)/ps
4777 bxtaus= 1.d0-4.d0*(
mtl*
mtl)/ps
4778 bxus= 1.d0-4.d0*(
muq*
muq)/ps
4779 bxds= 1.d0-4.d0*(
mdq*
mdq)/ps
4780 bxcs= 1.d0-4.d0*rmcs/ps
4781 bxss= 1.d0-4.d0*(
msq*
msq)/ps
4782 bxts= 1.d0-4.d0*(
mt*
mt)/ps
4783 bxbs= 1.d0-4.d0*rmbs/ps
4785 sh= - 9.d0*(lbh.
cp.bh)
4787 # - 2.d0*fw*(lbw.
cp.bw)
4789 # - (1.d0-6.d0*xz)*lz
4790 # - 2.d0*(1.d0-6.d0*xw)*lw
4791 # - 6.d0*(1.d0-2.d0*xw)*ls
4792 # +2.d0*(9.d0+xw*(-10.d0+24.d0*xw))
4793 # - 2.d0*xz*(5.d0-12.d0*xz)
4796 sh= sh/(128.d0*xw)*ps+ps/(32.d0*xw)*(
4797 # - 3.d0*xt*bxts*(bxt.
cp.lbxt)-3.d0*xt*lxt*
co
4798 # - 3.d0*xb*bxbs*(bxb.
cp.lbxb)-3.d0*xb*lxb*
co
4799 # - 3.d0*xs*bxss*(bxs.
cp.lbxs)-3.d0*xs*lxs*
co
4800 # - 3.d0*xc*bxcs*(bxc.
cp.lbxc)-3.d0*xc*lxc*
co
4801 # - 3.d0*xd*bxds*(bxd.
cp.lbxd)-3.d0*xd*lxd*
co
4802 # - 3.d0*xu*bxus*(bxu.
cp.lbxu)-3.d0*xd*lxu*
co
4803 # - xtau*bxtaus*(bxtau.
cp.lbxtau)-xtau*lxtau*
co
4804 # - xmu*bxmus*(bxmu.
cp.lbxmu)-xmu*lxmu*
co
4805 # - xe*bxes*(bxe.
cp.lbxe)-xe*lxe*
co+
co*(
4806 # - ls*(3.d0*(xt+xb+xs+xc+xd+xu)+xtau+xmu+xe)
4807 # +6.d0*xt*(1.d0-2.d0*xt)
4808 # +6.d0*xb*(1.d0-2.d0*xb)
4809 # +6.d0*xs*(1.d0-2.d0*xs)
4810 # +6.d0*xc*(1.d0-2.d0*xc)
4811 # +6.d0*xd*(1.d0-2.d0*xd)
4812 # +6.d0*xu*(1.d0-2.d0*xu)
4813 # +2.d0*xtau*(1.d0-2.d0*xtau)
4814 # +2.d0*xmu*(1.d0-2.d0*xmu)
4815 # +2.d0*xe*(1.d0-2.d0*xe)))
4817 dmh= - 18.d0*(bh.
cp.lbh)
4818 # - 2.d0*xz*fz*(lbz.
cq.bz)
4819 # - 2.d0*(1.d0-2.d0*xz)*(lbz.
cp.bz)
4820 # - 4.d0*xw*fw*(lbw.
cq.bw)
4821 # - 4.d0*(1.d0-2.d0*xw)*(lbw.
cp.bw)
4823 # - 2.d0*(1.d0-3.d0*xz)*lz
4824 # - 4.d0*(1.d0-3.d0*xw)*lw
4825 # - 12.d0*(1.d0-xw)*ls
4826 # +6.d0*(5.d0-2.d0*xw*(1.d0+2.d0*xw))
4827 # - 6.d0*xz*(1.d0+2.d0*xz)
4830 dmh= dmh/(128.d0*xw)
4831 # +1.d0/(32.d0*xw)*(
4832 # - 2.d0*xes*bxes*(lbxe.
cq.bxe)
4833 # - 2.d0*xmus*bxmus*(lbxmu.
cq.bxmu)
4834 # - 2.d0*xtaus*bxtaus*(lbxtau.
cq.bxtau)
4835 # - 6.d0*xus*bxus*(lbxu.
cq.bxu)
4836 # - 6.d0*xds*bxds*(lbxd.
cq.bxd)
4837 # - 6.d0*xcs*bxcs*(lbxc.
cq.bxc)
4838 # - 6.d0*xss*bxss*(lbxs.
cq.bxs)
4839 # - 6.d0*xts*bxts*(lbxt.
cq.bxt)
4840 # - 6.d0*xbs*bxbs*(lbxb.
cq.bxb)
4841 # - 3.d0*xt*((lbxt.
cp.bxt)+lxt*
co)
4842 # - 3.d0*xb*((lbxb.
cp.bxb)+lxb*
co)
4843 # - 3.d0*xc*((lbxc.
cp.bxc)+lxc*
co)
4844 # - 3.d0*xs*((lbxs.
cp.bxs)+lxs*
co)
4845 # - 3.d0*xu*((lbxu.
cp.bxu)+lxu*
co)
4846 # - 3.d0*xd*((lbxd.
cp.bxd)+lxd*
co)
4847 # - xtau*((lbxtau.
cp.bxtau)+lxtau*
co)
4848 # - xmu*((lbxmu.
cp.bxmu)+lxmu*
co)
4849 # - xe*((lbxe.
cp.bxe)+lxe*
co)+
co*(
4850 # - ls*(3.d0*(xt+xb+xs+xc+xd+xu)+xtau+xmu+xe)
4851 # +3.d0*xt*(2.d0+xts*bxts)
4852 # +3.d0*xb*(2.d0+xbs*bxbs)
4853 # +3.d0*xc*(2.d0+xcs*bxcs)
4854 # +3.d0*xs*(2.d0+xss*bxss)
4855 # +3.d0*xu*(2.d0+xus*bxus)
4856 # +3.d0*xd*(2.d0+xds*bxds)
4857 # +xtau*(2.d0+xtaus*bxtaus)
4858 # +xmu*(2.d0+xmus*bxmus)
4859 # +xe*(2.d0+xes*bxes)))
4861 dmw= (9.d0*(bh.
cp.lbh)
4863 # - xw*fw*(lbw.
cq.bw)
4864 # +2.d0*(1.d0-12.d0*xw*xw)*(bw.
cp.lbw)
4866 # 2.d0*lw+(1.d0-6.d0*xz)*lz+6.d0*(1.d0-xz)*ls
4867 # - 2.d0*(9.d0-xw*(2.d0+36.d0*xw))
4868 # +2.d0*xz*(5.d0-12.d0*xz)))/(128.d0*xw*xw)
4869 dmw= dmw-1.d0/(32.d0*xw*xw)*(
4870 # - 3.d0*xt*bxts*(bxt.
cp.lbxt)-3.d0*xt*lxt*
co
4871 # - 3.d0*xb*bxbs*(bxb.
cp.lbxb)-3.d0*xb*lxb*
co
4872 # - 3.d0*xs*bxss*(bxs.
cp.lbxs)-3.d0*xs*lxs*
co
4873 # - 3.d0*xc*bxcs*(bxc.
cp.lbxc)-3.d0*xc*lxc*
co
4874 # - 3.d0*xd*bxds*(bxd.
cp.lbxd)-3.d0*xd*lxd*
co
4875 # - 3.d0*xu*bxus*(bxu.
cp.lbxu)-3.d0*xd*lxu*
co
4876 # - xtau*bxtaus*(bxtau.
cp.lbxtau)-xtau*lxtau*
co
4877 # - xmu*bxmus*(bxmu.
cp.lbxmu)-xmu*lxmu*
co
4878 # - xe*bxes*(bxe.
cp.lbxe)-xe*lxe*
co
4880 # - ls*(3.d0*(xt+xb+xs+xc+xd+xu)+xtau+xmu+xe)
4881 # +6.d0*xt*(1.d0-2.d0*xt)
4882 # +6.d0*xb*(1.d0-2.d0*xb)
4883 # +6.d0*xs*(1.d0-2.d0*xs)
4884 # +6.d0*xc*(1.d0-2.d0*xc)
4885 # +6.d0*xd*(1.d0-2.d0*xd)
4886 # +6.d0*xu*(1.d0-2.d0*xu)
4887 # +2.d0*xtau*(1.d0-2.d0*xtau)
4888 # +2.d0*xmu*(1.d0-2.d0*xmu)
4889 # +2.d0*xe*(1.d0-2.d0*xe)))
4891 dmz= (-0.25d0+(xz-3.d0*xz))*(lbz.
cq.bz)
4892 # +2.d0*(1-6.d0*xz)*(bz.
cp.lbz)+
4893 #
co*(-4.d0+30.d0*xz+3.d0*(lz+ls))
4897 d2mh= (2.d0*fw*bwis*
co
4900 # +2.d0*xz*xz*fz*(lbz.
cp.bzic)
4901 # -2.d0*xz*(1.d0-12.d0*xz)*(lbz.
cp.bzi)
4903 # +4.d0*xw*xw*fz*(lbw.
cp.bwic)
4904 # -4.d0*xw*(1.d0-12.d0*xw)*(lbw.
cp.bwi)
4907 # -6.d0*ls-lz-2.d0*lw+9.d0+4.d0*xw*(1+3.d0*xw)+
4908 # 2.d0*xz*(1.d0+3.d0*xz)))/(64.d0*xw*ps)
4910 d2mh= d2mh+1.d0/(32.d0*xw*ps)*(
4911 # 4.d0*xec*(lbxe*bxeis)
4913 # -16.d0*xec*(lbxe.
cp.bxei)
4914 # +4.d0*xmuc*(lbxmu*bxmuis)
4916 # -16.d0*xmuc*(lbxmu.
cp.bxmui)
4917 # +4.d0*xtauc*(lbxtau*bxtauis)
4919 # -16.d0*xtauc*(lbxtau.
cp.bxtaui)
4920 # +12.d0*xuc*(lbxu*bxuis)
4922 # -48.d0*xuc*(lbxu.
cp.bxui)
4923 # +12.d0*xdc*(lbxd*bxdis)
4925 # -48.d0*xdc*(lbxd.
cp.bxdi)
4926 # +12.d0*xcc*(lbxc*bxcis)
4928 # -48.d0*xcc*(lbxc.
cp.bxci)
4929 # +12.d0*xsc*(lbxs*bxsis)
4931 # -48.d0*xsc*(lbxs.
cp.bxsi)
4932 # +12.d0*xtc*(lbxt*bxtis)
4934 # -48.d0*xtc*(lbxu.
cp.bxti)
4935 # +12.d0*xbc*(lbxb*bxbis)
4937 # -48.d0*xbc*(lbxb.
cp.bxbi)+
4949 dmwh= (2.d0*xw*fw*bwis*
co
4951 # +xz*fz*(lbz.
cq.bz)
4952 # +(1.d0-2.d0*xz)*(lbz.
cp.bz)
4953 # - 4.d0*xw*xw*fw*(lbw.
cp.bwic)
4954 # - xw*(1.d0+xw*xw*(-10.d0+48.d0*xw))*(lbw.
cq.bw)
4958 # +(1.d0-3.d0*xz)*lz
4959 # +3.d0*(2.d0-xz)*ls
4960 # - (15.d0+xw*(-2.d0+12.d0*xw))
4961 # +3.d0*xz*(1.d0+2.d0*xz)))/(64.d0*xw*xw*ps)
4963 dmzh= (2.d0*fz*bzis*
co
4964 # - 4.d0*xs*fz*(lbz.
cp.bzic)
4965 # - 3.d0*(1.d0+xz*(-6.d0+24.d0*xz))*(lbz.
cq.bz)
4970 # - 4.d0*(1.d0+6.d0*xz)))/(128.d0*xw*ps)
4972 dmzw= (fz*(lbz.
cq.bz)
4973 # - 8.d0*(1.d0-6.d0*xz)*(lbz.
cp.bz)
4977 # +8.d0*(2.d0-15.d0*xz)))/(256.d0*xw*xw*ps)
4979 dmww= (xw*(1.d0-4.d0*xw*(1.d0-2.d0*xw))*bwis*
co
4980 # - 18.d0*(lbh.
cp.bh)
4981 # - 2.d0*fz*(lbz.
cp.bz)
4982 # - 2.d0*xw*xw*fw*(lbw.
cp.bwic)
4983 # +2.d0*xw*(1.d0-12.d0*xw*xw)*(lbw.
cq.bw)
4984 # - 4.d0*(lbw.
cp.bw)
4987 # - 2.d0*(1.d0-6.d0*xz)*lz
4988 # - 12.d0*(1.d0-xz)*ls
4989 # +4.d0*(9.d0-xw*(1.d0-6.d0*xw))
4990 # - 4.d0*xz*(5.d0-12.d0*xz)))/(128.d0*xw*xw*xw*ps)
4991 dmww= dmww+1.d0/(162.d0*xw*xw*xw*ps)*(
4992 # - 3.d0*xt*bxts*(bxt.
cp.lbxt)-3.d0*xt*lxt*
co
4993 # - 3.d0*xb*bxbs*(bxb.
cp.lbxb)-3.d0*xb*lxb*
co
4994 # - 3.d0*xs*bxss*(bxs.
cp.lbxs)-3.d0*xs*lxs*
co
4995 # - 3.d0*xc*bxcs*(bxc.
cp.lbxc)-3.d0*xc*lxc*
co
4996 # - 3.d0*xd*bxds*(bxd.
cp.lbxd)-3.d0*xd*lxd*
co
4997 # - 3.d0*xu*bxus*(bxu.
cp.lbxu)-3.d0*xd*lxu*
co
4998 # - xtau*bxtaus*(bxtau.
cp.lbxtau)-xtau*lxtau*
co
4999 # - xmu*bxmus*(bxmu.
cp.lbxmu)-xmu*lxmu*
co
5000 # - xe*bxes*(bxe.
cp.lbxe)-xe*lxe*
co+
co*(
5001 # - ls*(3.d0*(xt+xb+xs+xc+xd+xu)+xtau+xmu+xe)
5002 # +6.d0*xt*(1.d0-2.d0*xt)
5003 # +6.d0*xb*(1.d0-2.d0*xb)
5004 # +6.d0*xs*(1.d0-2.d0*xs)
5005 # +6.d0*xc*(1.d0-2.d0*xc)
5006 # +6.d0*xd*(1.d0-2.d0*xd)
5007 # +6.d0*xu*(1.d0-2.d0*xu)
5008 # +2.d0*xtau*(1.d0-2.d0*xtau)
5009 # +2.d0*xmu*(1.d0-2.d0*xmu)
5010 # +2.d0*xe*(1.d0-2.d0*xe)))
5013 # - 2.d0*xz*fz*(lbz.
cp.bzic)
5014 # +8.d0*xz*(1.d0-6.d0*xz)*(lbz.
cq.bz)
5015 # - 48.d0*xz*(lbz.
cp.bz)
5016 # +4.d0*(1.d0+42.d0*xz)*
co)/(256.d0*xw*xz*ps)
5018 value(1,1:2)= sh(1:2)
5019 value(2,1:2)= dmh(1:2)
5020 value(3,1:2)= dmw(1:2)
5021 value(4,1:2)= dmz(1:2)
5022 value(5,1:2)= d2mh(1:2)
5023 value(6,1:2)= dmwh(1:2)
5024 value(7,1:2)= dmzh(1:2)
5025 value(8,1:2)= dmzw(1:2)
5026 value(9,1:2)= dmww(1:2)
5027 value(10,1:2)= dmzz(1:2)
5046 real*8 muh,cpgh,gos,mhb,ghb,m,sclaec,expghi,ewc,ghi
5047 real*8,
dimension(10,2) :: expc
5050 FUNCTION hto_deriv(scal,rhm)
RESULT(value)
5060 real*8,
dimension(10,2) ::
value
5070 IF(muh.le.200.d0)
THEN
5072 ewc= 4.d0*sqrt(2.d0)*1.16637d-5*(
mw*
mw)/
pis
5074 # -ewc*ewc*expc(2,2)*expc(2,2)*gos
5075 # +0.5d0*ewc*(expc(2,2)/muh-expc(5,2)*muh)*gos*gos
5077 ELSEIF((muh > 200.d0).and.(muh < 240.d0))
THEN
5084 mhb= sqrt(muh*muh+cpgh*cpgh)
5093 SUBROUTINE hto_gh(muh,cpgh)
5111 real*8 muh,scals,x1,x2,xacc,tgh,itgh,ghf,ghb,ght,cpgh
5112 real*8,
dimension(2) :: cmw,cmz
5194 IF(
muhcp.ge.900.d0)
THEN
5196 ELSEIF((
muhcp.gt.600.d0).and.(
muhcp.lt.900.d0))
THEN
5198 ELSEIF(
muhcp.lt.160.d0)
THEN
5233 print*,
'---- Warning ------------------------ '
5235 print*,
' inconsistent interval for tot '
5237 print*,
' inconsistent interval for lf '
5239 print*,
' inconsistent interval for top '
5241 print*,
' inconsistent interval for bos '
5243 print*,
'------------------------------------ '
5276 print*,
'---- Warning ------------------------ '
5278 print*,
' inconsistent interval for tot '
5280 print*,
' inconsistent interval for lf '
5282 print*,
' inconsistent interval for top '
5284 print*,
' inconsistent interval for bos '
5292 116
format(
' gH Total [GeV] =',e20.7)
5293 117
format(
' gH LightF [GeV] =',e20.5)
5294 1118
format(
' gH top [GeV] =',e20.5)
5295 118
format(
' gH Bos [GeV] =',e20.5)
5310 SUBROUTINE hto_hbrd(HTO_FCN,n,x,fvec,epsfcn,tol,info,diag)
5315 INTEGER,
INTENT(IN) :: n
5316 real*8,
INTENT(IN OUT) :: x(n)
5317 real*8,
INTENT(IN OUT) :: fvec(n)
5318 real*8,
INTENT(IN) :: epsfcn
5319 real*8,
INTENT(IN) :: tol
5320 INTEGER,
INTENT(OUT) :: info
5321 real*8,
INTENT(OUT) :: diag(n)
5326 SUBROUTINE hto_fcn(N,X,FVEC,IFLAG)
5328 INTEGER,
INTENT(IN) :: n
5329 real*8,
INTENT(IN) :: x(n)
5330 real*8,
INTENT(OUT) :: fvec(n)
5331 INTEGER,
INTENT(IN OUT) :: iflag
5332 END SUBROUTINE hto_fcn
5419 INTEGER :: maxfev, ml,mode,mu,nfev,nprint
5421 real*8,
PARAMETER :: factor= 100.0d0,zero= 0.0d0
5427 IF(n <= 0.or.epsfcn < zero.or.tol < zero)
GO TO 20
5437 CALL hto_hybrd(hto_fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag,
5438 # mode,factor,nprint,info,nfev)
5439 IF(info == 5) info= 4
5448 SUBROUTINE hto_hybrd(HTO_FCN,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,
5449 # diag,mode,factor,nprint,info,nfev)
5451 INTEGER,
INTENT(IN) :: n
5452 real*8,
INTENT(IN OUT) :: x(n)
5453 real*8,
INTENT(IN OUT) :: fvec(n)
5454 real*8,
INTENT(IN) :: xtol
5455 INTEGER,
INTENT(IN OUT) :: maxfev
5456 INTEGER,
INTENT(IN OUT) :: ml
5457 INTEGER,
INTENT(IN) :: mu
5458 real*8,
INTENT(IN) :: epsfcn
5459 real*8,
INTENT(OUT) :: diag(n)
5460 INTEGER,
INTENT(IN) :: mode
5461 real*8,
INTENT(IN) :: factor
5462 INTEGER,
INTENT(IN OUT) :: nprint
5463 INTEGER,
INTENT(OUT) :: info
5464 INTEGER,
INTENT(OUT) :: nfev
5469 SUBROUTINE hto_fcn(N,X,FVEC,IFLAG)
5471 INTEGER,
INTENT(IN) :: n
5472 real*8,
INTENT(IN) :: x(n)
5473 real*8,
INTENT(OUT) :: fvec(n)
5474 INTEGER,
INTENT(IN OUT) :: iflag
5475 END SUBROUTINE hto_fcn
5623 INTEGER :: i,iflag,iter,j,jm1,l,lr,msum,ncfail,ncsuc,
5626 LOGICAL :: jeval,sing
5627 real*8 :: actred,delta,epsmch,fnorm,fnorm1,pnorm,prered,
5628 # ratio,sum,temp,xnorm
5629 real*8,
PARAMETER :: one= 1.0d0,p1= 0.1d0,p5= 0.5d0,
5630 # p001= 0.001d0,p0001= 0.0001d0,zero= 0.0d0
5631 real*8 :: fjac(n,n),r(n*(n+1)/2),qtf(n),wa1(n),wa2(n),
5634 epsmch= epsilon(1.0d0)
5639 IF(n > 0.and.xtol >= zero.and.maxfev > 0.and.ml >= 0
5640 # .and.mu >= 0.and.factor > zero )
THEN
5645 CALL hto_fcn(n,x,fvec,iflag)
5649 msum= min(ml+mu+1,n)
5657 CALL hto_fdjac1(hto_fcn,n,x,fvec,fjac,n,iflag,ml,mu,epsfcn,
5661 CALL hto_qrfac(n,n,fjac,n,.false.,iwa,1,wa1,wa2,wa3)
5666 IF(wa2(j) == zero) diag(j)= one
5669 wa3(1:n)= diag(1:n)*x(1:n)
5672 IF(delta == zero) delta= factor
5676 IF(fjac(j,j) /= zero)
THEN
5679 sum= sum+fjac(i,j)*qtf(i)
5681 temp= -sum/fjac(j,j)
5683 qtf(i)= qtf(i)+fjac(i,j)*temp
5698 IF(wa1(j) == zero) sing= .true.
5703 diag(j)= max(diag(j),wa2(j))
5706 120
IF(nprint > 0)
THEN
5708 IF(mod(iter-1,nprint) == 0)
CALL hto_fcn(n,x,fvec,iflag)
5709 IF(iflag < 0)
GO TO 190
5711 CALL hto_dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
5715 wa3(j)= diag(j)*wa1(j)
5718 IF(iter == 1) delta= min(delta,pnorm)
5720 CALL hto_fcn(n,wa2,wa4,iflag)
5725 IF(fnorm1 < fnorm) actred= one-(fnorm1/fnorm)**2
5730 sum= sum+r(l)*wa1(j)
5737 IF(temp < fnorm) prered= one-(temp/fnorm)**2
5739 IF(prered > zero) ratio= actred/prered
5747 IF(ratio >= p5.or.ncsuc > 1) delta= max(delta,pnorm/p5)
5748 IF(abs(ratio-one) <= p1) delta= pnorm/p5
5750 IF(ratio >= p0001)
THEN
5753 wa2(j)= diag(j)*x(j)
5761 IF(actred >= p001) nslow1= 0
5762 IF(jeval) nslow2= nslow2+1
5763 IF(actred >= p1) nslow2= 0
5764 IF(delta <= xtol*xnorm.or.fnorm == zero) info= 1
5766 IF(nfev >= maxfev) info= 2
5767 IF(p1*max(p1*delta,pnorm) <= epsmch*xnorm) info= 3
5768 IF(nslow2 == 5) info= 4
5769 IF(nslow1 == 10) info= 5
5771 IF(ncfail /= 2)
THEN
5775 sum= sum+fjac(i,j)*wa4(i)
5777 wa2(j)= (sum-wa3(j))/pnorm
5778 wa1(j)= diag(j)*((diag(j)*wa1(j))/pnorm)
5779 IF(ratio >= p0001) qtf(j)= sum
5795 190
IF(iflag < 0) info= iflag
5797 IF(nprint > 0)
CALL hto_fcn(n,x,fvec,iflag)
5804 SUBROUTINE hto_dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
5806 INTEGER,
INTENT(IN) :: n
5807 INTEGER,
INTENT(IN) :: lr
5808 real*8,
INTENT(IN) :: r(lr)
5809 real*8,
INTENT(IN) :: diag(n)
5810 real*8,
INTENT(IN) :: qtb(n)
5811 real*8,
INTENT(IN) :: delta
5812 real*8,
INTENT(IN OUT) :: x(n)
5813 real*8,
INTENT(OUT) :: wa1(n)
5814 real*8,
INTENT(OUT) :: wa2(n)
5874 INTEGER :: i,j,jj,jp1,k,l
5875 real*8 :: alpha,bnorm,epsmch,gnorm,qnorm,sgnorm,sum,temp
5877 epsmch= epsilon(1.0d0)
5892 IF(temp == 0.0)
THEN
5895 temp= max(temp,abs(r(l)))
5899 IF(temp == 0.0) temp= epsmch
5901 x(j)= (qtb(j)-sum)/temp
5906 wa2(j)= diag(j)*x(j)
5910 IF(qnorm > delta)
THEN
5915 wa1(i)= wa1(i)+r(l)*temp
5918 wa1(j)= wa1(j)/diag(j)
5923 IF(gnorm /= 0.0)
THEN
5925 wa1(j)= (wa1(j)/gnorm)/diag(j)
5931 sum= sum+r(l)*wa1(i)
5937 sgnorm= (gnorm/temp)/temp
5939 IF(sgnorm < delta)
THEN
5941 temp= (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
5942 temp= temp-(delta/qnorm)*(sgnorm/delta)**2+
5943 # sqrt((temp-(delta/qnorm))**2+(1.0d0-(delta/qnorm)**2)*
5944 # (1.0d0-( sgnorm/delta)**2))
5945 alpha= ((delta/qnorm)*(1.0d0-(sgnorm/delta)**2))/temp
5948 temp= (1.0d0-alpha)*min(sgnorm,delta)
5950 x(j)= temp*wa1(j)+alpha*x(j)
5960 SUBROUTINE hto_fdjac1(HTO_FCN,n,x,fvec,fjac,ldfjac,iflag,ml,mu,
5963 INTEGER,
INTENT(IN) :: n
5964 REAL*8,
INTENT(IN OUT) :: x(n)
5965 real*8,
INTENT(IN) :: fvec(n)
5966 INTEGER,
INTENT(IN) :: ldfjac
5967 real*8,
INTENT(OUT) :: fjac(ldfjac,n)
5968 INTEGER,
INTENT(IN OUT) :: iflag
5969 INTEGER,
INTENT(IN) :: ml
5970 INTEGER,
INTENT(IN) :: mu
5971 real*8,
INTENT(IN) :: epsfcn
5972 real*8,
INTENT(IN OUT) :: wa1(n)
5973 real*8,
INTENT(OUT) :: wa2(n)
5978 SUBROUTINE hto_fcn(N,X,FVEC,IFLAG)
5980 INTEGER,
PARAMETER :: dp= selected_real_kind(14,60)
5981 INTEGER,
INTENT(IN) :: n
5982 real*8,
INTENT(IN) :: x(n)
5983 real*8,
INTENT(OUT) :: fvec(n)
5984 INTEGER,
INTENT(IN OUT) :: iflag
5985 END SUBROUTINE hto_fcn
6072 INTEGER :: i,j,k,msum
6073 REAL*8 :: eps,epsmch,h,temp
6074 real*8,
PARAMETER :: zero= 0.0d0
6076 epsmch= epsilon(1.0d0)
6077 eps= sqrt(max(epsfcn,epsmch))
6083 IF(h == zero) h= eps
6085 CALL hto_fcn(n,x,wa1,iflag)
6089 fjac(i,j)= (wa1(i)-fvec(i))/h
6097 IF(h == zero) h= eps
6100 CALL hto_fcn(n,x,wa1,iflag)
6105 IF(h == zero) h= eps
6108 IF(i >= j-mu.and.i <= j+ml) fjac(i,j)= (wa1(i)-fvec(i))/h
6122 INTEGER,
INTENT(IN) :: m
6123 INTEGER,
INTENT(IN) :: n
6124 INTEGER,
INTENT(IN) :: ldq
6125 real*8,
INTENT(OUT) :: q(ldq,m)
6126 real*8,
INTENT(OUT) :: wa(m)
6164 INTEGER :: i,j,jm1,k,l,minmn,np1
6166 real*8,
PARAMETER :: one= 1.0d0,zero= 0.0d0
6193 IF(wa(k) /= zero)
THEN
6197 sum= sum+q(i,j)*wa(i)
6201 q(i,j)= q(i,j)-temp*wa(i)
6213 SUBROUTINE hto_qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
6215 INTEGER,
INTENT(IN) :: m
6216 INTEGER,
INTENT(IN) :: n
6217 INTEGER,
INTENT(IN) :: lda
6218 real*8,
INTENT(IN OUT) :: a(lda,n)
6219 LOGICAL,
INTENT(IN) :: pivot
6220 INTEGER,
INTENT(IN) :: lipvt
6221 INTEGER,
INTENT(OUT) :: ipvt(lipvt)
6222 real*8,
INTENT(OUT) :: rdiag(n)
6223 real*8,
INTENT(OUT) :: acnorm(n)
6224 real*8,
INTENT(OUT) :: wa(n)
6297 INTEGER :: i,j,jp1,k,kmax,minmn
6298 REAL*8 :: ajnorm,epsmch,sum,temp
6299 real*8,
PARAMETER :: one= 1.0d0,p05= 0.05d0,zero= 0.0d0
6301 epsmch= epsilon(1.0d0)
6306 IF(pivot) ipvt(j)= j
6313 IF(rdiag(k) > rdiag(kmax)) kmax= k
6321 rdiag(kmax)= rdiag(j)
6329 IF(ajnorm /= zero)
THEN
6330 IF(a(j,j) < zero) ajnorm= -ajnorm
6332 a(i,j)= a(i,j)/ajnorm
6340 sum= sum+a(i,j)*a(i,k)
6344 a(i,k)= a(i,k)-temp*a(i,j)
6346 IF(.NOT.(.NOT.pivot.OR.rdiag(k) == zero))
THEN
6347 temp= a(j,k)/rdiag(k)
6348 rdiag(k)= rdiag(k)*sqrt(max(zero,one-temp**2))
6349 IF(p05*(rdiag(k)/wa(k))**2 <= epsmch)
THEN
6368 INTEGER,
INTENT(IN) :: m
6369 INTEGER,
INTENT(IN) :: n
6370 INTEGER,
INTENT(IN) :: lda
6371 real*8,
INTENT(IN OUT) :: a(lda,n)
6372 real*8,
INTENT(IN) :: v(n)
6373 real*8,
INTENT(IN) :: w(n)
6421 INTEGER :: i,j,nmj,nm1
6422 REAL*8 :: COS,SIN,temp
6423 real*8,
PARAMETER :: one= 1.0d0
6429 IF(abs(v(j)) > one) cos= one/v(j)
6430 IF(abs(v(j)) > one) sin= sqrt(one-cos**2)
6431 IF(abs(v(j)) <= one) sin= v(j)
6432 IF(abs(v(j)) <= one) cos= sqrt(one-sin**2)
6434 temp= cos*a(i,j)-sin*a(i,n)
6435 a(i,n)= sin*a(i,j)+cos*a(i,n)
6440 IF(abs(w(j)) > one) cos= one/w(j)
6441 IF(abs(w(j)) > one) sin= sqrt(one-cos**2)
6442 IF(abs(w(j)) <= one) sin= w(j)
6443 IF(abs(w(j)) <= one) cos= sqrt(one-sin**2)
6445 temp= cos*a(i,j)+sin*a(i,n)
6446 a(i,n)= -sin*a(i,j)+cos*a(i,n)
6460 INTEGER,
INTENT(IN) :: m
6461 INTEGER,
INTENT(IN) :: n
6462 INTEGER,
INTENT(IN) :: ls
6463 real*8,
INTENT(IN OUT) :: s(ls)
6464 real*8,
INTENT(IN) :: u(m)
6465 real*8,
INTENT(IN OUT) :: v(n)
6466 real*8,
INTENT(OUT) :: w(m)
6467 LOGICAL,
INTENT(OUT) :: sing
6533 INTEGER :: i,j,jj,l,nmj,nm1
6534 REAL*8 :: COS,cotan,giant,SIN,TAN,tau,temp
6535 real*8,
PARAMETER :: one= 1.0d0,p5= 0.5d0,p25= 0.25d0,zero= 0.0d0
6538 jj= (n*(2*m-n+1))/2-(m-n)
6550 IF(v(j) /= zero)
THEN
6551 IF(abs(v(n)) < abs(v(j)))
THEN
6553 sin= p5/sqrt(p25+p25*cotan**2)
6556 IF(abs(cos)*giant > one) tau= one/cos
6559 cos= p5/sqrt(p25+p25*tan**2)
6563 v(n)= sin*v(j)+cos*v(n)
6567 temp= cos*s(l)-sin*w(i)
6568 w(i)= sin*s(l)+cos*w(i)
6576 w(i)= w(i)+v(n)*u(i)
6581 IF(w(j) /= zero)
THEN
6582 IF(abs(s(jj)) < abs(w(j)))
THEN
6584 sin= p5/sqrt(p25+p25*cotan**2)
6587 IF(abs(cos)*giant > one) tau= one/cos
6590 cos= p5/sqrt(p25+p25*tan**2)
6596 temp= cos*s(l)+sin*w(i)
6597 w(i)= -sin*s(l)+cos*w(i)
6603 IF(s(jj) == zero) sing= .true.
6612 IF(s(jj) == zero) sing= .true.
6622 INTEGER,
INTENT(IN) :: n
6623 REAL*8,
INTENT(IN) :: x(n)
6662 REAL*8 :: agiant,floatn,s1,s2,s3,xabs,x1max,x3max
6663 real*8,
PARAMETER :: rdwarf= 1.0d-100,rgiant= 1.0d+100
6671 agiant= rgiant/floatn
6674 IF(xabs <= rdwarf.or.xabs >= agiant)
THEN
6675 IF(xabs > rdwarf)
THEN
6676 IF(xabs > x1max)
THEN
6677 s1= 1.0d0+s1*(x1max/xabs)**2
6680 s1= s1+(xabs/x1max)**2
6683 IF(xabs > x3max)
THEN
6684 s3= 1.0d0+s3*(x3max/xabs)**2
6687 IF(xabs /= 0.0d0) s3= s3+(xabs/x3max)**2
6694 IF(s1 /= 0.0d0)
THEN
6695 fn_val= x1max*sqrt(s1+(s2/x1max)/x1max)
6697 IF(s2 /= 0.0d0)
THEN
6698 IF(s2 >= x3max) fn_val= sqrt(s2*(1.0d0+(x3max/s2)*(x3max*s3)))
6699 IF(s2 < x3max) fn_val= sqrt(x3max*((s2/x3max)+(x3max*s3)))
6701 fn_val= x3max*sqrt(s3)
6723 real*8,
dimension(3) :: xcp,fvcp,diag
6729 xcp(1)= sqrt(20.d0/
mw)
6736 print*,
'ALGORITHM ESTIMATES THAT THE RELATIVE ERROR'
6737 print*,
'BETWEEN X AND THE SOLUTION IS AT MOST TOL'
6738 ELSEIF(info == 2)
THEN
6739 print*,
'NUMBER OF CALLS TO FCN HAS REACHED'
6740 print*,
'OR EXCEEDED 200*(N+1)'
6741 ELSEIF(info == 3)
THEN
6742 print*,
'TOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN'
6743 print*,
'THE APPROXIMATE SOLUTION X IS POSSIBLE'
6744 ELSEIF(info == 4)
THEN
6745 print*,
'ITERATION IS NOT MAKING GOOD PROGRESS'
6748 print 2005,1.d-2*xcp(2)*xcp(2)*m
6749 print 2004,1.d-1*xcp(1)*xcp(1)*
mw,
6752 print 2006,1.d-2*xcp(3)*xcp(3)*
mt,
imt
6755 2003
format(
' info =',i2)
6757 2004
format(
' gammaW = ',2e20.5)
6758 2005
format(
' gammaH = ',e20.5)
6759 2006
format(
' gammat = ',2e20.5)
6784 INTEGER nc,n,iflag,iz
6785 REAL*8 muh,rgh,muhs,scal,scals,p2,xm0,str,sti,rgw,
6786 # lswr,lswi,lmw,bxm0,emc,emb,emt,asmur,rgt,as_nlo,ewc,
6787 # crmbs,crmcs,lcxb,lcxc,lcxbs,lcxcs,lclxb,lclxc,as_lo
6788 real*8,
dimension(n) :: xcp,fvcp
6789 real*8,
dimension(2,2) :: axms
6790 real*8,
dimension(2) :: axm0,bxms,runbc
6791 real*8,
dimension(2) :: sh,shs,clh,b0sumb,b0sumf,cxp,ksumb,
6792 # ksumf,coefb1,coefb2,coefb3,coefb4,coefb5,coefb6,coefb7,
6793 # coefb8,coefb9,coefb10,coefb11,coefb12,
6794 # shhf,shht,shhb,shh,b0sumt,ksumt,xms,
6795 # b0part,cpt,cpts,ccxt,deltag,sww0,coefw1,
6796 # coefw2,ksumw,ccxts,cctsi,shi,clt,cxhw,cstsi,csmcts,cltmw,
6797 # b0sumw,sww,dw,b0sumw1,b0sumw2,cxw,cxz,ccts,clw,csts,
6798 # cctq,cxws,cmxw,clmw,cxtw,kfmsb,ktmsb,kbmsb,kwmsb,
6799 # coeft1,coeft2,coeft3,coeft4,b0sumtop,ksumtop,stt,
6800 # coeft1s,b0sumtops,ksumtops,stts,cctqi,ucpt,nloqcd,ttqcd,
6801 # dww,sww0w,swww,ksumww,lcxwi,lclcts
6840 real*8 scal,ps0i,xm0i
6841 real*8,
dimension(2) ::
value,psi,xmsi
6846 FUNCTION hto_lb0af_em(scal,psi,ps0i,xmsi,xm0i)
RESULT(value)
6856 real*8 scal,ps0i,xm0i
6857 real*8,
dimension(2) ::
value,psi,xmsi
6862 FUNCTION hto_lb0af_dm(scal,psi,ps0i,xmsi,xm0i)
RESULT(value)
6873 real*8,
dimension(2,2) :: xmsi
6874 real*8,
dimension(2) ::
value,psi,xm0i
6890 real*8,
intent(in),
dimension(2) :: xm0i
6891 real*8,
intent(in),
dimension(2,2) :: xmsi
6892 real*8,
dimension(2) ::
value
6893 real*8,
dimension(2) :: psi
6902 cxw(1)=
mw*
mw/scals*(1.d0-0.5d-2*(xcp(1)*xcp(1))**2)
6903 cxw(2)= -1.d-1*
mw*
mw/scals*xcp(1)*xcp(1)
6909 rgw= 1.d-1*xcp(1)*xcp(1)*
mw
6910 rgh= 1.d-2*xcp(2)*xcp(2)*muh
6912 rgt= 1.d-2*xcp(3)*xcp(3)*
mt
6917 lswr=
mw*
mw*(1.d0-0.5d-2*(xcp(1)*xcp(1))**2)
6920 clw= cxw(1).
fln.cxw(2)
6921 clmw= cmxw(1).
fln.cmxw(2)
6922 clmw(2)= clmw(2)-2.d0*pi
6930 lclcts= ccts(1).
fln.ccts(2)
6940 crmbs= runbc(2)*runbc(2)
6941 crmcs= runbc(1)*runbc(1)
6964 sh(2)= -muh*rgh/scals
6971 clh= sh(1).
fln.sh(2)
6972 clt= ccxt(1).
fln.ccxt(2)
6973 cltmw= cxtw(1).
fln.cxtw(2)
6976 coefb1= (12.d0*cxw-4.d0*sh+(shs.
cq.cxw))/64.d0
6979 coefb2= (-4.d0*(sh.
cq.ccts)+(shs.
cq.cxw)+
6980 # 12.d0*(cxw.
cq.cctq))/128.d0
6983 coefb3= 9.d0/128.d0*(shs.
cq.cxw)
6986 coefb4= -3.d0/32.d0*((4.d0*ccxt-sh).
cp.(ccxt.
cq.cxw))
6990 coefb5= -3.d0/32.d0*((4.d0*
co*lcxb-sh).
cq.cxw)*lcxb
6994 coefb7= -3.d0/32.d0*((4.d0*
co*lcxc-sh).
cq.cxw)*lcxc
6998 coefb9= -3.d0/32.d0*((4.d0*
co*
cxs-sh).
cq.cxw)*
cxs
7000 coefb10= -3.d0/32.d0*((4.d0*
co*
cxd-sh).
cq.cxw)*
cxd
7002 coefb11= -3.d0/32.d0*((4.d0*
co*
cxu-sh).
cq.cxw)*
cxu
7004 coefb12= -1.d0/32.d0*((4.d0*
co*
cxe-sh).
cq.cxw)*
cxe
7014 b0sumb= coefb1.
cp.b0part
7020 b0sumb= b0sumb+(coefb2.
cp.b0part)
7026 b0sumb= b0sumb+(coefb3.
cp.b0part)
7028 xms(1:2)= crmbs*
co(1:2)
7031 b0sumf= coefb5.
cp.b0part
7036 b0sumf= b0sumf+(coefb6.
cp.b0part)
7038 xms(1:2)= crmcs*
co(1:2)
7041 b0sumf= b0sumf+(coefb7.
cp.b0part)
7046 b0sumf= b0sumf+(coefb8.
cp.b0part)
7051 b0sumf= b0sumf+(coefb9.
cp.b0part)
7056 b0sumf= b0sumf+(coefb10.
cp.b0part)
7061 b0sumf= b0sumf+(coefb11.
cp.b0part)
7066 b0sumf= b0sumf+(coefb12.
cp.b0part)
7071 b0sumt= coefb4.
cp.b0part
7074 xms(1:2)= crmbs*
co(1:2)
7077 xms(1:2)= crmcs*
co(1:2)
7080 IF(
qcdc==0) nloqcd= 0.d0
7087 IF(
qcdc==0) ttqcd= 0.d0
7094 ksumt= 3.d0/8.d0*((ccxt.
cq.cxw).
cp.ccxt)
7095 # -3.d0/32.d0*(clt.
cp.((sh.
cq.cxw).
cp.ccxt))
7097 ksumb= -1.d0/64.d0*((2.d0*
co+cctsi+3.d0*(sh.
cq.cxw)).
cp.sh)
7098 # -1.d0/128.d0*(lclcts.
cp.((6.d0*cctsi-(sh.
cq.cxw)).
cp.sh))
7099 # +3.d0/128.d0*(clw.
cp.((4.d0*
co+2.d0*cctsi-(sh.
cq.cxw)).
cp.sh))
7100 # -3.d0/128.d0*(clh.
cp.(shs.
cq.cxw))
7104 IF((lswr*shht(2)+lswi*shht(1)) < 0.d0) shht= 0.d0
7107 IF((lswr*shhb(2)+lswi*shhb(1)) < 0.d0) shhb= 0.d0
7112 deltag= 6.d0*
co+0.5d0*(((7.d0*
co-4.d0*csts).
cq.csts).
cp.lclcts)
7114 sww0= -(38.d0*cxw+6.d0*ccxt+7.d0*sh
7115 # -48.d0*(((ccxt.
cq.sh).
cq.cxw).
cp.ccxt)+8.d0*(cxw.
cq.sh))/128.d0
7116 # -3.d0/64.d0*((cxw-sh+(cxws.
cq.cxhw)).
cp.clh)
7117 # +3.d0/32.d0*(((
co-4.d0*((ccxt.
cq.sh).
cq.cxw).
cp.ccxt)).
cp.clt)
7118 # +((((8.d0*
co-17.d0*cstsi+3.d0*cctsi).
cp.cxw)
7119 # -6.d0*((cxw.
cq.sh).
cq.cctq)).
cp.lclcts)/64.d0
7120 # -((cxw.
cq.sh).
cq.cctq)/32.d0+5.d0/128.d0*(cxw.
cq.ccts)
7122 sww0w= -(38.d0*cxw+6.d0*ccxt+7.d0*sh
7123 # -48.d0*(((ccxt.
cq.sh).
cq.cxw).
cp.ccxt)+8.d0*(cxw.
cq.sh))/128.d0
7124 # -3.d0/64.d0*((cxw-sh+(cxws.
cq.cxhw)).
cp.(clh-clw))
7126 # (((
co-4.d0*((ccxt.
cq.sh).
cq.cxw).
cp.ccxt)).
cp.(clt-clw))
7127 # +((((8.d0*
co-17.d0*cstsi+3.d0*cctsi).
cp.cxw)
7128 # -6.d0*((cxw.
cq.sh).
cq.cctq)).
cp.lclcts)/64.d0
7129 # -((cxw.
cq.sh).
cq.cctq)/32.d0+5.d0/128.d0*(cxw.
cq.ccts)
7131 coefw1= -(((8.d0*
co-(sh.
cq.cxw)).
cp.sh)*sh
7132 # -4.d0*((-12.d0*cxw+7.d0*sh).
cp.cxw))/192.d0
7134 coefw2= -((cxws.
cq.csmcts).
cp.(416.d0*
co-192.d0*csts
7135 # -((132.d0*
co-((12.d0*
co+cctsi).
cq.ccts)).
cq.ccts)))/192.d0
7148 b0sumw1= (coefw1.
cp.b0part)
7149 b0sumw= (coefw1.
cp.b0part)
7158 b0sumw2= (coefw2.
cp.b0part)
7159 b0sumw= b0sumw+(coefw2.
cp.b0part)
7162 # -0.5d0*((3.d0*
co-(ccxts.
cq.cxws)).
cp.ccxt)).
cp.cltmw)
7163 # -((24.d0*cxw-((14.d0*
co-(sh.
cq.cxw)).
cp.sh)).
cp.clh)
7164 # +((36.d0*cxw-14.d0*sh-18.d0*((
co-4.d0*(ccxt.
cq.sh)).
cp.ccxt)
7165 # +(shs.
cq.cxw)).
cp.clw)
7166 # -6.d0*(((2.d0*
co+((lcxwi-12.d0*shi).
cp.ccxt)).
cp.ccxt)
7167 # +1.d0/6.d0*((15.d0*
co-(sh.
cq.cxw)).
cp.sh)
7168 # +2.d0/9.d0*((97.d0*
co+9.d0*(cxw.
cq.sh)).
cp.cxw))
7169 # +(((cxw.
cq.ccts).
cp.(
co-6.d0*(cxw.
cq.sh))).
cq.ccts)
7170 # -2.d0*(((cxw.
cq.csmcts).
cp.lclcts).
cp.(62.d0*
co
7171 # -48.d0*csts-5.d0*cctsi))
7172 # -18.d0*(((cxws.
cq.sh).
cq.cctq).
cp.lclcts)
7173 # -72.d0*((clt.
cp.ccxts).
cp.(shi-1.d0/12.d0*(ccxt.
cq.cxws)))
7174 # +23.d0*(cxw.
cq.ccts)
7175 ksumw= ksumw/192.d0+3.d0/16.d0*(cxw.
cp.(clw-clmw))
7177 ksumww= -12.d0*((cxw
7178 # -0.5d0*((3.d0*
co-(ccxts.
cq.cxws)).
cp.ccxt)).
cp.(cltmw-clw))
7179 # -((24.d0*cxw-((14.d0*
co-(sh.
cq.cxw)).
cp.sh)).
cp.(clh-clw))
7180 # -6.d0*(((2.d0*
co+((lcxwi-12.d0*shi).
cp.ccxt)).
cp.ccxt)
7181 # +1.d0/6.d0*((15.d0*
co-(sh.
cq.cxw)).
cp.sh)
7182 # +2.d0/9.d0*((97.d0*
co+9.d0*(cxw.
cq.sh)).
cp.cxw))
7183 # +(((cxw.
cq.ccts).
cp.(
co-6.d0*(cxw.
cq.sh))).
cq.ccts)
7184 # -2.d0*(((cxw.
cq.csmcts).
cp.lclcts).
cp.(62.d0*
co
7185 # -48.d0*csts-5.d0*cctsi))
7186 # -18.d0*(((cxws.
cq.sh).
cq.cctq).
cp.lclcts)
7188 # (((clt-clw).
cp.ccxts).
cp.(shi-1.d0/12.d0*(ccxt.
cq.cxws)))
7189 # +23.d0*(cxw.
cq.ccts)
7190 ksumww= ksumww/192.d0+3.d0/16.d0*(cxw.
cp.(clw-clmw))
7195 dw= -sww+sww0+deltag/16.d0
7196 dww= -swww+sww0w+deltag/16.d0
7199 ksumtop= 1.d0/128.d0*(-(2.d0*cxw+3.d0*
cxb*
co+
cxbs*lcxwi)+48.d0*
7202 # +1.d0/576.d0*(-(50.d0*cxw-9.d0*
cxb*
co-9.d0*
cxbs*lcxwi
7203 # +17.d0*(cxw.
cq.cctq)-40.d0*(cxw.
cq.ccts))+18.d0*(cpts.
cq.cxw)
7204 # +(cpt.
cp.(
co+17.d0*cctsi+9.d0*
cxb*lcxwi-216.d0*((
7205 # cpts.
cq.cxw).
cq.sh)+54.d0*((cxw.
cq.sh).
cq.cctq)
7206 # +18.d0*(sh.
cq.cxw)-216.d0*(shi.
cp.cxw)*
cxbs+108.d0*(cxw.
cq.sh))))
7207 # +1.d0/1152.d0*(clt.
cp.(432.d0*((cpts.
cq.cxw).
cp.(cpt.
cq.sh))
7208 # -(cxw.
cp.(32.d0*
co-40.d0*cctsi+17.d0*cctqi))
7209 # +(cpt.
cp.(32.d0-64.d0*csts-41.d0*cctsi-9.d0*(sh.
cq.cxw)))))
7210 # -1.d0/128.d0*(clh.
cp.((-4.d0*cpt+5.d0*sh).
cp.(cpt.
cq.cxw)))
7211 # +1.d0/1152.d0*(clw.
cp.((50.d0*cxw+27.d0*
cxb*
co+9.d0*
cxbs*
7212 # lcxwi+17.d0*(cxw.
cq.cctq)-40.d0*(cxw.
cq.ccts))
7213 # +9.d0*(cpts.
cq.cxw)+(cpt.
cp.(7.d0*
co+64.d0*csts-7.d0*cctsi
7214 # -18.d0*
cxb*lcxwi-108.d0*((cxw.
cq.sh).
cq.cctq)
7215 # -216.d0*(cxw.
cq.sh)))))
7216 # -1.d0/1152.d0*(lclcts.
cp.((cxw.
cp.(32.d0*
co-40.d0*
7217 # cctsi+17.d0*cctqi))+(cpt.
cp.(16.d0*
co+64.d0*csts-7.d0*
7218 # cctsi-108.d0*((cxw.
cq.sh).
cq.cctq)))))
7220 coeft1= -1.d0/576.d0*(-(cxw.
cp.(32.d0*
co-40.d0*cctsi+17.d0*
7221 # cctqi))+(cpt.
cp.(16.d0*
co+64.d0*csts-7.d0*cctsi)))
7223 coeft2= 1.d0/64.d0*((-4.d0*cpt+sh).
cp.(cpt.
cq.cxw))
7225 coeft3= 1.d0/9.d0*(csts.
cp.cpt)
7227 coeft4= 1.d0/64.d0*((-2.d0*cxw+
cxb*
co+
cxbs*lcxwi)
7228 # +(cpts.
cq.cxw)+(cpt.
cp.(
co-2.d0*
cxb*lcxwi)))
7241 b0sumtop= (coeft1.
cp.b0part)
7250 b0sumtop= b0sumtop+(coeft2.
cp.b0part)
7253 b0sumtop= b0sumtop+(coeft3.
cp.b0part)
7262 b0sumtop= b0sumtop+(coeft4.
cp.b0part)
7264 ksumtops= cpt/6.d0-(clt.
cp.cpt)/2.d0
7266 coeft1s= 1.d0/3.d0*cpt
7269 b0sumtops= coeft1s.
cp.b0part
7271 stt= b0sumtop+ksumtop
7272 stts= b0sumtops+ksumtops
7274 ewc= 4.d0*sqrt(2.d0)*
g_f/
pis
7284 fvcp(1)= rgw/
mw*(1.d0+ewc*(lswr*dww(1)-lswi*dww(2)))
7285 # -ewc*(scal/
mw)**2*(lswr*swww(2)+lswi*swww(1))
7287 fvcp(2)= rgh/muh*(1.d0+ewc*(lswr*dw(1)-lswi*dw(2)))
7288 # -ewc*(lswr*shh(2)+lswi*shh(1))*scals/muhs
7289 # +as_nlo*
g_f/(sqrt(2.d0)*
pis)*
7290 # (sh(1)*nloqcd(2)+sh(2)*nloqcd(1))*scals/muhs
7291 # +as_nlo*
g_f/(sqrt(2.d0)*
pis)*
7292 # (sh(1)*ttqcd(2)+sh(2)*ttqcd(1))*scals/muhs
7295 fvcp(3)= rgt/
mt*(1.d0+2.d0*ewc*(lswr*dw(1)-lswi*dw(2)))
7296 # -scals/(
mt*
mt)*(ewc*(lswr*stt(2)+lswi*stt(1))+
7297 # 4.d0*(scals/(
mt*
mt))*as_lo*stts(2))
7306 FUNCTION hto_lquarkqcd(scal,psi,ps0i,xmsi,xm0i,type)
RESULT(value)
7320 INTEGER it,
type,unit
7321 REAL*8 scal,ps0i,xm0i,scals,ps0,xm0,sgn
7322 real*8,
dimension(2) ::
value,psi,ps,xmsi,xms,betasc,betas,
7323 # betac,beta,argc,arg,lq,lm,cx,comx,cxs,x,omx,xs,clx,clomx,
7324 # clxs,li2cx,li3cx,li2cxs,li3cxs,copx,clopx,lqs,qcd,lms,
7325 # opx,tau,taus,clxx,clxxs
7326 real*8,
dimension(6,2) :: aux,auxs
7334 IF(psi(2).eq.0.d0.and.xmsi(2).eq.0.d0.
7335 # and.psi(1).le.4.d0*xmsi(1))
THEN
7341 IF(abs(ps(2)/ps(1)).lt.1.d-10.and.xms(2).eq.0.d0)
THEN
7342 betasc(1)= 1.d0-4.d0*xms(1)/ps(1)
7343 betasc(2)= 4.d0/(ps(1)*ps(1))*xms(1)*ps(2)
7345 betasc=
co-4.d0*(xms.
cq.ps)
7347 IF(betasc(2).eq.0.d0)
THEN
7349 betac= (betasc(1)).
cr.(betasc(2))
7351 betac= (betasc(1)).
crz.(betasc(2))
7353 argc= (betac+
co).
cq.(betac-
co)
7355 betas(1)= 1.d0-4.d0*xm0*xm0/ps0
7357 beta= (betas(1)).
cr.(betas(2))
7358 arg= (beta+
co).
cq.(beta-
co)
7360 IF(arg(2).eq.0.d0)
THEN
7367 x= (beta-
co).
cq.(beta+
co)
7373 IF(arg(2).eq.0.d0) arg(2)=
eps
7374 IF(argc(2).eq.0.d0)
THEN
7377 lq= argc(1).
fln.argc(2)
7389 sgn= sign(
one,cx(1))
7392 clx= cx(1).
fln.cx(2)
7394 clxx= cxs(1).
fln.cxs(2)
7396 clomx= comx(1).
fln.comx(2)
7397 clopx= copx(1).
fln.copx(2)
7399 li2cx(1:2)= aux(1,1:2)
7400 li3cx(1:2)= aux(2,1:2)
7402 li2cxs(1:2)= auxs(1,1:2)
7403 li3cxs(1:2)= auxs(2,1:2)
7404 ELSEIF(it.eq.1.d0)
THEN
7405 cx= (betac-
co).
cq.(betac+
co)
7413 clomx= comx.
lnsrs.omx
7414 clopx= copx.
lnsrs.opx
7421 IF(xms(2).eq.0.d0)
THEN
7425 lm= xms(1).
fln.xms(2)
7434 qcd= -3.d0/4.d0*(
co-12.d0*tau)*
rz2
7435 # +1.d0/16.d0*(3.d0*
co+344.d0*tau)
7436 # -3.d0/2.d0*((
co-6.d0*tau).
cp.lms)
7437 # +3.d0/4.d0*((3.d0*
co-14.d0*tau).
cp.(betac.
cp.clx))
7438 # -3.d0*((li3cxs-2.d0*li3cx+4.d0/3.d0*(li2cx.
cp.clx)
7439 # +1.d0/3.d0*(clxs.
cp.clomx)+
rz3*
co
7440 # -2.d0/3.d0*(clxx.
cp.li2cxs)-1.d0/6.d0*(clxxs.
cp.clomx)
7442 # (clxxs.
cp.clopx)).
cp.((
co-4.d0*tau).
cp.(
co-2.d0*tau)))
7443 # -1.d0/2.d0*((4.d0*li2cxs-4.d0*li2cx-4.d0*(clomx.
cp.clx)
7444 # -2.d0*(cx.
cp.(clxs.
cq.comx))+4.d0*(clxx.
cp.clomx)
7445 # +4.d0*(clxx.
cp.clopx)+(clxxs.
cp.(cxs.
cq.comx))
7446 # +(clxxs.
cp.(cxs.
cq.copx))).
cp.(betac.
cp.(
co-4.d0*tau)))
7447 # +1.d0/4.d0*(lm.
cp.(11.d0*
co-108.d0*tau))
7448 # +1.d0/4.d0*(clxs.
cp.(3.d0*
co+58.d0*taus-28.d0*tau))
7450 ELSEIF(type.eq.1)
THEN
7452 qcd= -3.d0/4.d0*(
co-12.d0*tau)*
rz2
7453 # +1.d0/16.d0*(67.d0*
co-40.d0*tau)
7454 # +3.d0/4.d0*((3.d0*
co-14.d0*tau).
cp.(betac.
cp.clx))
7455 # -3.d0*((li3cxs-2.d0*li3cx+4.d0/3.d0*(li2cx.
cp.clx)
7456 # +1.d0/3.d0*(clxs.
cp.clomx)+
rz3*
co
7457 # -2.d0/3.d0*(clxx.
cp.li2cxs)-1.d0/6.d0*(clxxs.
cp.clomx)
7459 # (clxxs.
cp.clopx)).
cp.((
co-4.d0*tau).
cp.(
co-2.d0*tau)))
7460 # -1.d0/2.d0*((4.d0*li2cxs-4.d0*li2cx-4.d0*(clomx.
cp.clx)
7461 # -2.d0*(cx.
cp.(clxs.
cq.comx))+4.d0*(clxx.
cp.clomx)
7462 # +4.d0*(clxx.
cp.clopx)+(clxxs.
cp.(cxs.
cq.comx))
7463 # +(clxxs.
cp.(cxs.
cq.copx))).
cp.(betac.
cp.(
co-4.d0*tau)))
7465 # (((
co-8.d0*tau).
cp.(betac.
cp.lq)).
cp.(
co-3.d0/4.d0*lm))
7466 # -3.d0*((betac.
cp.lq).
cp.(lm.
cp.tau))
7467 # +4.d0*((
co.
cp.betac).
cp.(lq.
cp.tau))
7468 # -9.d0/4.d0*(lm.
cp.(
co+4.d0*tau))
7469 # +9.d0*(lms.
cp.tau)
7470 # +1.d0/4.d0*(clxs.
cp.(3.d0*
co+58.d0*taus-28.d0*tau))
7482 FUNCTION hto_lb0af_dm(scal,psi,ps0i,xmsi,xm0i)
RESULT(value)
7494 real*8 scal,ps0i,scals,ps0,xm1,xm2
7495 real*8,
dimension(2,2) :: xmsi,xms
7496 real*8,
dimension(2) ::
value,psi,ps,lambdasc,lambdas,
7497 # lambdac,lambda,argc,arg,llam,lm,xm0,xm0i,aroot,
7498 # root,rat,lnr,xm1c,xm2c
7505 xm1c(1:2)= xms(1,1:2)
7506 xm2c(1:2)= xms(2,1:2)
7510 root= (aroot(1).
crz.aroot(2))
7512 lambdasc= (ps.
cp.ps)+(xm1c.
cp.xm1c)+(xm2c.
cp.xm2c)-2.d0*(
7513 # (ps.
cp.xm1c)+(ps.
cp.xm2c)+(xm1c.
cp.xm2c))
7514 lambdas(1)= ps0*ps0+xm1*xm1+xm2*xm2-2.d0*(
7515 # ps0*xm1+ps0*xm2+xm1*xm2)
7517 lambdac= (lambdasc(1)).
crz.(lambdasc(2))
7518 lambda= (lambdas(1)).
cr.(lambdas(2))
7519 IF(lambda(2).eq.0.d0) lambda(2)= -
eps
7521 argc= 0.5d0*((-ps+xm1c+xm2c-lambdac).
cq.root)
7523 arg(1)= 0.5d0*(-ps0+xm1+xm2-lambda(1))/sqrt(xm1*xm2)
7526 llam= argc.
lnsrs.arg
7529 lnr= rat(1).
fln.rat(2)
7531 value= 2.d0*
co-0.5d0*(((xm1c-xm2c).
cq.ps).
cp.lnr)-
7532 # ((lambdac.
cq.ps).
cp.llam)
7540 FUNCTION hto_lb0af_em(scal,psi,ps0i,xmsi,xm0i)
RESULT(value)
7552 real*8 scal,ps0i,xm0i,scals,ps0,xm0
7553 real*8,
dimension(2) ::
value,psi,ps,xmsi,xms,betasc,betas,
7554 # betac,beta,argc,arg,lbet
7562 betasc=
co-4.d0*(xms.
cq.ps)
7563 IF(betasc(2).eq.0.d0)
THEN
7565 betac= (betasc(1)).
cr.(betasc(2))
7567 betac= (betasc(1)).
crz.(betasc(2))
7569 argc= (betac+
co).
cq.(betac-
co)
7570 IF(argc(2).eq.0.d0)
THEN
7572 lbet= argc(1).
fln.argc(2)
7574 betas(1)= 1.d0-4.d0*(xm0*xm0)/ps0
7576 beta= (betas(1)).
cr.(betas(2))
7577 arg= (beta+
co).
cq.(beta-
co)
7578 IF(arg(2).eq.0.d0) arg(2)=
eps
7579 lbet= argc.
lnsrs.arg
7582 value= 2.d0*
co-(betac.
cp.lbet)
7603 REAL*8 scal,scals,ps0i,ps0,xm1,xm2
7604 real*8,
intent(in),
dimension(2) :: xm0i
7605 real*8,
intent(in),
dimension(2,2) :: xmsi
7606 real*8,
dimension(2) ::
value,xm1c,xm2c
7607 real*8,
dimension(2,2) :: xms
7608 real*8,
dimension(2) :: xm0
7609 real*8,
dimension(2) :: psi,ps,aroot,root,lambdasc,lambdas,
7610 # lambdac,lambda,argc,arg,llam,l1,l2
7617 xm1c(1:2)= xms(1,1:2)
7618 xm2c(1:2)= xms(2,1:2)
7621 xm1= xm0(1)*xm0(1)/ps0
7622 xm2= xm0(2)*xm0(2)/ps0
7625 root= (aroot(1).
crz.aroot(2))
7627 lambdasc=
co+(xm1c.
cp.xm1c)+(xm2c.
cp.xm2c)-2.d0*(
7628 # xm1c+xm2c+(xm1c.
cp.xm2c))
7629 lambdas(1)= 1.d0+xm1*xm1+xm2*xm2-2.d0*(
7632 lambdac= (lambdasc(1)).
crz.(lambdasc(2))
7633 lambda= (lambdas(1)).
cr.(lambdas(2))
7634 IF(lambda(2).eq.0.d0) lambda(2)= -
eps
7636 argc= 0.5d0*((-
co+xm1c+xm2c-lambdac).
cq.root)
7638 arg(1)= 0.5d0*(-1.d0+xm1+xm2-lambda(1))/sqrt(xm1*xm2)
7641 llam= argc.
lnsrs.arg
7643 l1= xm1c(1).
fln.xm1c(2)
7644 l2= xm2c(1).
fln.xm2c(2)
7646 value= (((xm1c-xm2c-
co).
cq.lambdac).
cp.llam)+0.5d0*(l1-l2)
7655 SUBROUTINE call_hto(mhiggs,mtop,mhb,ghb)
7660 real*8 mhiggs,mtop,mhb,ghb