47 G4double alphaBeta=0.0, alphaReduced=0.0, betaReduced=0.0,
48 root1=0.0, root2=0.0, root3=0.0 ;
50 newton1=0.0, newton2=0.0, newton3=0.0, newton0=0.0,
51 temp=0.0, rootTemp=0.0 ;
57 for (i=1;i<=nJacobi;i++)
61 alphaReduced = alpha/nJacobi ;
62 betaReduced = beta/nJacobi ;
63 root1 = (1.0+
alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
64 0.767999*alphaReduced/nJacobi) ;
65 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced
66 + 0.451998*alphaReduced*alphaReduced
67 + 0.83001*alphaReduced*betaReduced ;
68 root = 1.0-root1/root2 ;
72 root1=(4.1002+
alpha)/((1.0+alpha)*(1.0+0.155998*
alpha)) ;
73 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*
alpha)/nJacobi ;
74 root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
75 root -= (1.0-root)*root1*root2*root3 ;
79 root1=(1.67001+0.27998*
alpha)/(1.0+0.37002*alpha) ;
80 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
81 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
82 root -= (
fAbscissa[0]-root)*root1*root2*root3 ;
84 else if (i == nJacobi-1)
86 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
87 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
88 root3=1.0/(1.0+20.0*alpha/((7.5+
alpha)*nJacobi*nJacobi)) ;
89 root += (root-
fAbscissa[nJacobi-4])*root1*root2*root3 ;
91 else if (i == nJacobi)
93 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
94 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
95 root3 = 1.0/(1.0+8.0*alpha/((6.28002+
alpha)*nJacobi*nJacobi)) ;
96 root += (root-
fAbscissa[nJacobi-3])*root1*root2*root3 ;
102 alphaBeta = alpha + beta ;
103 for (
k=1;
k<=maxNumber;
k++)
105 temp = 2.0 + alphaBeta ;
106 newton1 = (alpha-beta+temp*root)/2.0 ;
108 for (
G4int j=2;j<=nJacobi;j++)
112 temp = 2*j+alphaBeta ;
113 a = 2*j*(j+alphaBeta)*(temp-2.0) ;
114 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
115 c = 2.0*(j-1+
alpha)*(j-1+beta)*temp ;
116 newton1 = (
b*newton2-
c*newton3)/a ;
118 newton0 = (nJacobi*(alpha - beta - temp*root)*newton1 +
119 2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/
120 (temp*(1.0 - root*root)) ;
122 root = rootTemp - newton1/newton0 ;
123 if (std::fabs(root-rootTemp) <= tolerance)
130 G4Exception(
"G4GaussJacobiQ::G4GaussJacobiQ()",
"OutOfRange",
138 *temp*std::pow(2.0,alphaBeta)/(newton0*newton2) ;