40{
42 const G4int maxNumber = 12;
43
44 G4int i = 1, j = 1, k = 1;
46 G4double newton1 = 0.0, temp1 = 0.0, temp2 = 0.0, temp3 = 0.0, temp = 0.0;
47 G4double piInMinusQ = std::pow(pi, -0.25);
48
52
54 {
55 if(i == 1)
56 {
57 newton0 =
58 std::sqrt((
G4double)(2 * nHermite + 1)) -
59 1.85575001 * std::pow((
G4double)(2 * nHermite + 1), -0.16666999);
60 }
61 else if(i == 2)
62 {
63 newton0 -= 1.14001 * std::pow((
G4double) nHermite, 0.425999) / newton0;
64 }
65 else if(i == 3)
66 {
67 newton0 = 1.86002 * newton0 - 0.86002 *
fAbscissa[0];
68 }
69 else if(i == 4)
70 {
71 newton0 = 1.91001 * newton0 - 0.91001 *
fAbscissa[1];
72 }
73 else
74 {
75 newton0 = 2.0 * newton0 -
fAbscissa[i - 3];
76 }
77 for(k = 1; k <= maxNumber; ++k)
78 {
79 temp1 = piInMinusQ;
80 temp2 = 0.0;
81 for(j = 1; j <= nHermite; ++j)
82 {
83 temp3 = temp2;
84 temp2 = temp1;
85 temp1 = newton0 * std::sqrt(2.0 / j) * temp2 -
86 std::sqrt(((
G4double)(j - 1)) / j) * temp3;
87 }
88 temp = std::sqrt((
G4double) 2 * nHermite) * temp2;
89 newton1 = newton0;
90 newton0 = newton1 - temp1 / temp;
91 if(std::fabs(newton0 - newton1) <= tolerance)
92 {
93 break;
94 }
95 }
96 if(k > maxNumber)
97 {
98 G4Exception(
"G4GaussHermiteQ::G4GaussHermiteQ()",
"OutOfRange",
100 "Too many iterations in Gauss-Hermite constructor.");
101 }
103 fWeight[i - 1] = 2.0 / (temp * temp);
104 }
105}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4VGaussianQuadrature(function pFunction)