Geant4  10.00.p03
G4Integrator.icc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4Integrator.icc 81649 2014-06-04 09:56:53Z gcosmo $
28 //
29 // Implementation of G4Integrator methods.
30 //
31 //
32 
33 /////////////////////////////////////////////////////////////////////
34 //
35 // Sympson integration method
36 //
37 /////////////////////////////////////////////////////////////////////
38 //
39 // Integration of class member functions T::f by Simpson method.
40 
41 template <class T, class F>
42 G4double G4Integrator<T,F>::Simpson( T& typeT,
43  F f,
44  G4double xInitial,
45  G4double xFinal,
46  G4int iterationNumber )
47 {
48  G4int i ;
49  G4double step = (xFinal - xInitial)/iterationNumber ;
50  G4double x = xInitial ;
51  G4double xPlus = xInitial + 0.5*step ;
52  G4double mean = ( (typeT.*f)(xInitial) + (typeT.*f)(xFinal) )*0.5 ;
53  G4double sum = (typeT.*f)(xPlus) ;
54 
55  for(i=1;i<iterationNumber;i++)
56  {
57  x += step ;
58  xPlus += step ;
59  mean += (typeT.*f)(x) ;
60  sum += (typeT.*f)(xPlus) ;
61  }
62  mean += 2.0*sum ;
63 
64  return mean*step/3.0 ;
65 }
66 
67 /////////////////////////////////////////////////////////////////////
68 //
69 // Integration of class member functions T::f by Simpson method.
70 // Convenient to use with 'this' pointer
71 
72 template <class T, class F>
73 G4double G4Integrator<T,F>::Simpson( T* ptrT,
74  F f,
75  G4double xInitial,
76  G4double xFinal,
77  G4int iterationNumber )
78 {
79  G4int i ;
80  G4double step = (xFinal - xInitial)/iterationNumber ;
81  G4double x = xInitial ;
82  G4double xPlus = xInitial + 0.5*step ;
83  G4double mean = ( (ptrT->*f)(xInitial) + (ptrT->*f)(xFinal) )*0.5 ;
84  G4double sum = (ptrT->*f)(xPlus) ;
85 
86  for(i=1;i<iterationNumber;i++)
87  {
88  x += step ;
89  xPlus += step ;
90  mean += (ptrT->*f)(x) ;
91  sum += (ptrT->*f)(xPlus) ;
92  }
93  mean += 2.0*sum ;
94 
95  return mean*step/3.0 ;
96 }
97 
98 /////////////////////////////////////////////////////////////////////
99 //
100 // Integration of class member functions T::f by Simpson method.
101 // Convenient to use, when function f is defined in global scope, i.e. in main()
102 // program
103 
104 template <class T, class F>
105 G4double G4Integrator<T,F>::Simpson( G4double (*f)(G4double),
106  G4double xInitial,
107  G4double xFinal,
108  G4int iterationNumber )
109 {
110  G4int i ;
111  G4double step = (xFinal - xInitial)/iterationNumber ;
112  G4double x = xInitial ;
113  G4double xPlus = xInitial + 0.5*step ;
114  G4double mean = ( (*f)(xInitial) + (*f)(xFinal) )*0.5 ;
115  G4double sum = (*f)(xPlus) ;
116 
117  for(i=1;i<iterationNumber;i++)
118  {
119  x += step ;
120  xPlus += step ;
121  mean += (*f)(x) ;
122  sum += (*f)(xPlus) ;
123  }
124  mean += 2.0*sum ;
125 
126  return mean*step/3.0 ;
127 }
128 
129 //////////////////////////////////////////////////////////////////////////
130 //
131 // Adaptive Gauss method
132 //
133 //////////////////////////////////////////////////////////////////////////
134 //
135 //
136 
137 template <class T, class F>
138 G4double G4Integrator<T,F>::Gauss( T& typeT, F f,
139  G4double xInitial, G4double xFinal )
140 {
141  static const G4double root = 1.0/std::sqrt(3.0) ;
142 
143  G4double xMean = (xInitial + xFinal)/2.0 ;
144  G4double Step = (xFinal - xInitial)/2.0 ;
145  G4double delta = Step*root ;
146  G4double sum = ((typeT.*f)(xMean + delta) +
147  (typeT.*f)(xMean - delta)) ;
148 
149  return sum*Step ;
150 }
151 
152 //////////////////////////////////////////////////////////////////////
153 //
154 //
155 
156 template <class T, class F> G4double
157 G4Integrator<T,F>::Gauss( T* ptrT, F f, G4double a, G4double b )
158 {
159  return Gauss(*ptrT,f,a,b) ;
160 }
161 
162 ///////////////////////////////////////////////////////////////////////
163 //
164 //
165 
166 template <class T, class F>
167 G4double G4Integrator<T,F>::Gauss( G4double (*f)(G4double),
168  G4double xInitial, G4double xFinal)
169 {
170  static const G4double root = 1.0/std::sqrt(3.0) ;
171 
172  G4double xMean = (xInitial + xFinal)/2.0 ;
173  G4double Step = (xFinal - xInitial)/2.0 ;
174  G4double delta = Step*root ;
175  G4double sum = ( (*f)(xMean + delta) + (*f)(xMean - delta) ) ;
176 
177  return sum*Step ;
178 }
179 
180 ///////////////////////////////////////////////////////////////////////////
181 //
182 //
183 
184 template <class T, class F>
185 void G4Integrator<T,F>::AdaptGauss( T& typeT, F f, G4double xInitial,
186  G4double xFinal, G4double fTolerance,
187  G4double& sum,
188  G4int& depth )
189 {
190  if(depth > 100)
191  {
192  G4cout<<"G4Integrator<T,F>::AdaptGauss: WARNING !!!"<<G4endl ;
193  G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps "
194  <<G4endl ;
195 
196  return ;
197  }
198  G4double xMean = (xInitial + xFinal)/2.0 ;
199  G4double leftHalf = Gauss(typeT,f,xInitial,xMean) ;
200  G4double rightHalf = Gauss(typeT,f,xMean,xFinal) ;
201  G4double full = Gauss(typeT,f,xInitial,xFinal) ;
202  if(std::fabs(leftHalf+rightHalf-full) < fTolerance)
203  {
204  sum += full ;
205  }
206  else
207  {
208  depth++ ;
209  AdaptGauss(typeT,f,xInitial,xMean,fTolerance,sum,depth) ;
210  AdaptGauss(typeT,f,xMean,xFinal,fTolerance,sum,depth) ;
211  }
212 }
213 
214 template <class T, class F>
215 void G4Integrator<T,F>::AdaptGauss( T* ptrT, F f, G4double xInitial,
216  G4double xFinal, G4double fTolerance,
217  G4double& sum,
218  G4int& depth )
219 {
220  AdaptGauss(*ptrT,f,xInitial,xFinal,fTolerance,sum,depth) ;
221 }
222 
223 /////////////////////////////////////////////////////////////////////////
224 //
225 //
226 template <class T, class F>
227 void G4Integrator<T,F>::AdaptGauss( G4double (*f)(G4double),
228  G4double xInitial, G4double xFinal,
229  G4double fTolerance, G4double& sum,
230  G4int& depth )
231 {
232  if(depth > 100)
233  {
234  G4cout<<"G4SimpleIntegration::AdaptGauss: WARNING !!!"<<G4endl ;
235  G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps "
236  <<G4endl ;
237 
238  return ;
239  }
240  G4double xMean = (xInitial + xFinal)/2.0 ;
241  G4double leftHalf = Gauss(f,xInitial,xMean) ;
242  G4double rightHalf = Gauss(f,xMean,xFinal) ;
243  G4double full = Gauss(f,xInitial,xFinal) ;
244  if(std::fabs(leftHalf+rightHalf-full) < fTolerance)
245  {
246  sum += full ;
247  }
248  else
249  {
250  depth++ ;
251  AdaptGauss(f,xInitial,xMean,fTolerance,sum,depth) ;
252  AdaptGauss(f,xMean,xFinal,fTolerance,sum,depth) ;
253  }
254 }
255 
256 ////////////////////////////////////////////////////////////////////////
257 //
258 // Adaptive Gauss integration with accuracy 'e'
259 // Convenient for using with class object typeT
260 
261 template<class T, class F>
262 G4double G4Integrator<T,F>::AdaptiveGauss( T& typeT, F f, G4double xInitial,
263  G4double xFinal, G4double e )
264 {
265  G4int depth = 0 ;
266  G4double sum = 0.0 ;
267  AdaptGauss(typeT,f,xInitial,xFinal,e,sum,depth) ;
268  return sum ;
269 }
270 
271 ////////////////////////////////////////////////////////////////////////
272 //
273 // Adaptive Gauss integration with accuracy 'e'
274 // Convenient for using with 'this' pointer
275 
276 template<class T, class F>
277 G4double G4Integrator<T,F>::AdaptiveGauss( T* ptrT, F f, G4double xInitial,
278  G4double xFinal, G4double e )
279 {
280  return AdaptiveGauss(*ptrT,f,xInitial,xFinal,e) ;
281 }
282 
283 ////////////////////////////////////////////////////////////////////////
284 //
285 // Adaptive Gauss integration with accuracy 'e'
286 // Convenient for using with global scope function f
287 
288 template <class T, class F>
289 G4double G4Integrator<T,F>::AdaptiveGauss( G4double (*f)(G4double),
290  G4double xInitial, G4double xFinal, G4double e )
291 {
292  G4int depth = 0 ;
293  G4double sum = 0.0 ;
294  AdaptGauss(f,xInitial,xFinal,e,sum,depth) ;
295  return sum ;
296 }
297 
298 ////////////////////////////////////////////////////////////////////////////
299 // Gauss integration methods involving ortogonal polynomials
300 ////////////////////////////////////////////////////////////////////////////
301 //
302 // Methods involving Legendre polynomials
303 //
304 /////////////////////////////////////////////////////////////////////////
305 //
306 // The value nLegendre set the accuracy required, i.e the number of points
307 // where the function pFunction will be evaluated during integration.
308 // The function creates the arrays for abscissas and weights that used
309 // in Gauss-Legendre quadrature method.
310 // The values a and b are the limits of integration of the function f .
311 // nLegendre MUST BE EVEN !!!
312 // Returns the integral of the function f between a and b, by 2*fNumber point
313 // Gauss-Legendre integration: the function is evaluated exactly
314 // 2*fNumber times at interior points in the range of integration.
315 // Since the weights and abscissas are, in this case, symmetric around
316 // the midpoint of the range of integration, there are actually only
317 // fNumber distinct values of each.
318 // Convenient for using with some class object dataT
319 
320 template <class T, class F>
321 G4double G4Integrator<T,F>::Legendre( T& typeT, F f, G4double a, G4double b,
322  G4int nLegendre )
323 {
324  G4double nwt, nwt1, temp1, temp2, temp3, temp ;
325  G4double xDiff, xMean, dx, integral ;
326 
327  const G4double tolerance = 1.6e-10 ;
328  G4int i, j, k = nLegendre ;
329  G4int fNumber = (nLegendre + 1)/2 ;
330 
331  if(2*fNumber != k)
332  {
333  G4Exception("G4Integrator<T,F>::Legendre(T&,F, ...)", "InvalidCall",
334  FatalException, "Invalid (odd) nLegendre in constructor.");
335  }
336 
337  G4double* fAbscissa = new G4double[fNumber] ;
338  G4double* fWeight = new G4double[fNumber] ;
339 
340  for(i=1;i<=fNumber;i++) // Loop over the desired roots
341  {
342  nwt = std::cos(CLHEP::pi*(i - 0.25)/(k + 0.5)) ; // Initial root approximation
343 
344  do // loop of Newton's method
345  {
346  temp1 = 1.0 ;
347  temp2 = 0.0 ;
348  for(j=1;j<=k;j++)
349  {
350  temp3 = temp2 ;
351  temp2 = temp1 ;
352  temp1 = ((2.0*j - 1.0)*nwt*temp2 - (j - 1.0)*temp3)/j ;
353  }
354  temp = k*(nwt*temp1 - temp2)/(nwt*nwt - 1.0) ;
355  nwt1 = nwt ;
356  nwt = nwt1 - temp1/temp ; // Newton's method
357  }
358  while(std::fabs(nwt - nwt1) > tolerance) ;
359 
360  fAbscissa[fNumber-i] = nwt ;
361  fWeight[fNumber-i] = 2.0/((1.0 - nwt*nwt)*temp*temp) ;
362  }
363 
364  //
365  // Now we ready to get integral
366  //
367 
368  xMean = 0.5*(a + b) ;
369  xDiff = 0.5*(b - a) ;
370  integral = 0.0 ;
371  for(i=0;i<fNumber;i++)
372  {
373  dx = xDiff*fAbscissa[i] ;
374  integral += fWeight[i]*( (typeT.*f)(xMean + dx) +
375  (typeT.*f)(xMean - dx) ) ;
376  }
377  delete[] fAbscissa;
378  delete[] fWeight;
379  return integral *= xDiff ;
380 }
381 
382 ///////////////////////////////////////////////////////////////////////
383 //
384 // Convenient for using with the pointer 'this'
385 
386 template <class T, class F>
387 G4double G4Integrator<T,F>::Legendre( T* ptrT, F f, G4double a,
388  G4double b, G4int nLegendre )
389 {
390  return Legendre(*ptrT,f,a,b,nLegendre) ;
391 }
392 
393 ///////////////////////////////////////////////////////////////////////
394 //
395 // Convenient for using with global scope function f
396 
397 template <class T, class F>
398 G4double G4Integrator<T,F>::Legendre( G4double (*f)(G4double),
399  G4double a, G4double b, G4int nLegendre)
400 {
401  G4double nwt, nwt1, temp1, temp2, temp3, temp ;
402  G4double xDiff, xMean, dx, integral ;
403 
404  const G4double tolerance = 1.6e-10 ;
405  G4int i, j, k = nLegendre ;
406  G4int fNumber = (nLegendre + 1)/2 ;
407 
408  if(2*fNumber != k)
409  {
410  G4Exception("G4Integrator<T,F>::Legendre(...)", "InvalidCall",
411  FatalException, "Invalid (odd) nLegendre in constructor.");
412  }
413 
414  G4double* fAbscissa = new G4double[fNumber] ;
415  G4double* fWeight = new G4double[fNumber] ;
416 
417  for(i=1;i<=fNumber;i++) // Loop over the desired roots
418  {
419  nwt = std::cos(CLHEP::pi*(i - 0.25)/(k + 0.5)) ; // Initial root approximation
420 
421  do // loop of Newton's method
422  {
423  temp1 = 1.0 ;
424  temp2 = 0.0 ;
425  for(j=1;j<=k;j++)
426  {
427  temp3 = temp2 ;
428  temp2 = temp1 ;
429  temp1 = ((2.0*j - 1.0)*nwt*temp2 - (j - 1.0)*temp3)/j ;
430  }
431  temp = k*(nwt*temp1 - temp2)/(nwt*nwt - 1.0) ;
432  nwt1 = nwt ;
433  nwt = nwt1 - temp1/temp ; // Newton's method
434  }
435  while(std::fabs(nwt - nwt1) > tolerance) ;
436 
437  fAbscissa[fNumber-i] = nwt ;
438  fWeight[fNumber-i] = 2.0/((1.0 - nwt*nwt)*temp*temp) ;
439  }
440 
441  //
442  // Now we ready to get integral
443  //
444 
445  xMean = 0.5*(a + b) ;
446  xDiff = 0.5*(b - a) ;
447  integral = 0.0 ;
448  for(i=0;i<fNumber;i++)
449  {
450  dx = xDiff*fAbscissa[i] ;
451  integral += fWeight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx) ) ;
452  }
453  delete[] fAbscissa;
454  delete[] fWeight;
455 
456  return integral *= xDiff ;
457 }
458 
459 ////////////////////////////////////////////////////////////////////////////
460 //
461 // Returns the integral of the function to be pointed by T::f between a and b,
462 // by ten point Gauss-Legendre integration: the function is evaluated exactly
463 // ten times at interior points in the range of integration. Since the weights
464 // and abscissas are, in this case, symmetric around the midpoint of the
465 // range of integration, there are actually only five distinct values of each
466 // Convenient for using with class object typeT
467 
468 template <class T, class F>
469 G4double G4Integrator<T,F>::Legendre10( T& typeT, F f,G4double a, G4double b)
470 {
471  G4int i ;
472  G4double xDiff, xMean, dx, integral ;
473 
474  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
475 
476  static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
477  0.679409568299024, 0.865063366688985,
478  0.973906528517172 } ;
479 
480  static const G4double weight[] = { 0.295524224714753, 0.269266719309996,
481  0.219086362515982, 0.149451349150581,
482  0.066671344308688 } ;
483  xMean = 0.5*(a + b) ;
484  xDiff = 0.5*(b - a) ;
485  integral = 0.0 ;
486  for(i=0;i<5;i++)
487  {
488  dx = xDiff*abscissa[i] ;
489  integral += weight[i]*( (typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ;
490  }
491  return integral *= xDiff ;
492 }
493 
494 ///////////////////////////////////////////////////////////////////////////
495 //
496 // Convenient for using with the pointer 'this'
497 
498 template <class T, class F>
499 G4double G4Integrator<T,F>::Legendre10( T* ptrT, F f,G4double a, G4double b)
500 {
501  return Legendre10(*ptrT,f,a,b) ;
502 }
503 
504 //////////////////////////////////////////////////////////////////////////
505 //
506 // Convenient for using with global scope functions
507 
508 template <class T, class F>
509 G4double G4Integrator<T,F>::Legendre10( G4double (*f)(G4double),
510  G4double a, G4double b )
511 {
512  G4int i ;
513  G4double xDiff, xMean, dx, integral ;
514 
515  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
516 
517  static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
518  0.679409568299024, 0.865063366688985,
519  0.973906528517172 } ;
520 
521  static const G4double weight[] = { 0.295524224714753, 0.269266719309996,
522  0.219086362515982, 0.149451349150581,
523  0.066671344308688 } ;
524  xMean = 0.5*(a + b) ;
525  xDiff = 0.5*(b - a) ;
526  integral = 0.0 ;
527  for(i=0;i<5;i++)
528  {
529  dx = xDiff*abscissa[i] ;
530  integral += weight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx)) ;
531  }
532  return integral *= xDiff ;
533 }
534 
535 ///////////////////////////////////////////////////////////////////////
536 //
537 // Returns the integral of the function to be pointed by T::f between a and b,
538 // by 96 point Gauss-Legendre integration: the function is evaluated exactly
539 // ten Times at interior points in the range of integration. Since the weights
540 // and abscissas are, in this case, symmetric around the midpoint of the
541 // range of integration, there are actually only five distinct values of each
542 // Convenient for using with some class object typeT
543 
544 template <class T, class F>
545 G4double G4Integrator<T,F>::Legendre96( T& typeT, F f,G4double a, G4double b)
546 {
547  G4int i ;
548  G4double xDiff, xMean, dx, integral ;
549 
550  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
551 
552  static const G4double
553  abscissa[] = {
554  0.016276744849602969579, 0.048812985136049731112,
555  0.081297495464425558994, 0.113695850110665920911,
556  0.145973714654896941989, 0.178096882367618602759, // 6
557 
558  0.210031310460567203603, 0.241743156163840012328,
559  0.273198812591049141487, 0.304364944354496353024,
560  0.335208522892625422616, 0.365696861472313635031, // 12
561 
562  0.395797649828908603285, 0.425478988407300545365,
563  0.454709422167743008636, 0.483457973920596359768,
564  0.511694177154667673586, 0.539388108324357436227, // 18
565 
566  0.566510418561397168404, 0.593032364777572080684,
567  0.618925840125468570386, 0.644163403784967106798,
568  0.668718310043916153953, 0.692564536642171561344, // 24
569 
570  0.715676812348967626225, 0.738030643744400132851,
571  0.759602341176647498703, 0.780369043867433217604,
572  0.800308744139140817229, 0.819400310737931675539, // 30
573 
574  0.837623511228187121494, 0.854959033434601455463,
575  0.871388505909296502874, 0.886894517402420416057,
576  0.901460635315852341319, 0.915071423120898074206, // 36
577 
578  0.927712456722308690965, 0.939370339752755216932,
579  0.950032717784437635756, 0.959688291448742539300,
580  0.968326828463264212174, 0.975939174585136466453, // 42
581 
582  0.982517263563014677447, 0.988054126329623799481,
583  0.992543900323762624572, 0.995981842987209290650,
584  0.998364375863181677724, 0.999689503883230766828 // 48
585  } ;
586 
587  static const G4double
588  weight[] = {
589  0.032550614492363166242, 0.032516118713868835987,
590  0.032447163714064269364, 0.032343822568575928429,
591  0.032206204794030250669, 0.032034456231992663218, // 6
592 
593  0.031828758894411006535, 0.031589330770727168558,
594  0.031316425596862355813, 0.031010332586313837423,
595  0.030671376123669149014, 0.030299915420827593794, // 12
596 
597  0.029896344136328385984, 0.029461089958167905970,
598  0.028994614150555236543, 0.028497411065085385646,
599  0.027970007616848334440, 0.027412962726029242823, // 18
600 
601  0.026826866725591762198, 0.026212340735672413913,
602  0.025570036005349361499, 0.024900633222483610288,
603  0.024204841792364691282, 0.023483399085926219842, // 24
604 
605  0.022737069658329374001, 0.021966644438744349195,
606  0.021172939892191298988, 0.020356797154333324595,
607  0.019519081140145022410, 0.018660679627411467385, // 30
608 
609  0.017782502316045260838, 0.016885479864245172450,
610  0.015970562902562291381, 0.015038721026994938006,
611  0.014090941772314860916, 0.013128229566961572637, // 36
612 
613  0.012151604671088319635, 0.011162102099838498591,
614  0.010160770535008415758, 0.009148671230783386633,
615  0.008126876925698759217, 0.007096470791153865269, // 42
616 
617  0.006058545504235961683, 0.005014202742927517693,
618  0.003964554338444686674, 0.002910731817934946408,
619  0.001853960788946921732, 0.000796792065552012429 // 48
620  } ;
621  xMean = 0.5*(a + b) ;
622  xDiff = 0.5*(b - a) ;
623  integral = 0.0 ;
624  for(i=0;i<48;i++)
625  {
626  dx = xDiff*abscissa[i] ;
627  integral += weight[i]*((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ;
628  }
629  return integral *= xDiff ;
630 }
631 
632 ///////////////////////////////////////////////////////////////////////
633 //
634 // Convenient for using with the pointer 'this'
635 
636 template <class T, class F>
637 G4double G4Integrator<T,F>::Legendre96( T* ptrT, F f,G4double a, G4double b)
638 {
639  return Legendre96(*ptrT,f,a,b) ;
640 }
641 
642 ///////////////////////////////////////////////////////////////////////
643 //
644 // Convenient for using with global scope function f
645 
646 template <class T, class F>
647 G4double G4Integrator<T,F>::Legendre96( G4double (*f)(G4double),
648  G4double a, G4double b )
649 {
650  G4int i ;
651  G4double xDiff, xMean, dx, integral ;
652 
653  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
654 
655  static const G4double
656  abscissa[] = {
657  0.016276744849602969579, 0.048812985136049731112,
658  0.081297495464425558994, 0.113695850110665920911,
659  0.145973714654896941989, 0.178096882367618602759, // 6
660 
661  0.210031310460567203603, 0.241743156163840012328,
662  0.273198812591049141487, 0.304364944354496353024,
663  0.335208522892625422616, 0.365696861472313635031, // 12
664 
665  0.395797649828908603285, 0.425478988407300545365,
666  0.454709422167743008636, 0.483457973920596359768,
667  0.511694177154667673586, 0.539388108324357436227, // 18
668 
669  0.566510418561397168404, 0.593032364777572080684,
670  0.618925840125468570386, 0.644163403784967106798,
671  0.668718310043916153953, 0.692564536642171561344, // 24
672 
673  0.715676812348967626225, 0.738030643744400132851,
674  0.759602341176647498703, 0.780369043867433217604,
675  0.800308744139140817229, 0.819400310737931675539, // 30
676 
677  0.837623511228187121494, 0.854959033434601455463,
678  0.871388505909296502874, 0.886894517402420416057,
679  0.901460635315852341319, 0.915071423120898074206, // 36
680 
681  0.927712456722308690965, 0.939370339752755216932,
682  0.950032717784437635756, 0.959688291448742539300,
683  0.968326828463264212174, 0.975939174585136466453, // 42
684 
685  0.982517263563014677447, 0.988054126329623799481,
686  0.992543900323762624572, 0.995981842987209290650,
687  0.998364375863181677724, 0.999689503883230766828 // 48
688  } ;
689 
690  static const G4double
691  weight[] = {
692  0.032550614492363166242, 0.032516118713868835987,
693  0.032447163714064269364, 0.032343822568575928429,
694  0.032206204794030250669, 0.032034456231992663218, // 6
695 
696  0.031828758894411006535, 0.031589330770727168558,
697  0.031316425596862355813, 0.031010332586313837423,
698  0.030671376123669149014, 0.030299915420827593794, // 12
699 
700  0.029896344136328385984, 0.029461089958167905970,
701  0.028994614150555236543, 0.028497411065085385646,
702  0.027970007616848334440, 0.027412962726029242823, // 18
703 
704  0.026826866725591762198, 0.026212340735672413913,
705  0.025570036005349361499, 0.024900633222483610288,
706  0.024204841792364691282, 0.023483399085926219842, // 24
707 
708  0.022737069658329374001, 0.021966644438744349195,
709  0.021172939892191298988, 0.020356797154333324595,
710  0.019519081140145022410, 0.018660679627411467385, // 30
711 
712  0.017782502316045260838, 0.016885479864245172450,
713  0.015970562902562291381, 0.015038721026994938006,
714  0.014090941772314860916, 0.013128229566961572637, // 36
715 
716  0.012151604671088319635, 0.011162102099838498591,
717  0.010160770535008415758, 0.009148671230783386633,
718  0.008126876925698759217, 0.007096470791153865269, // 42
719 
720  0.006058545504235961683, 0.005014202742927517693,
721  0.003964554338444686674, 0.002910731817934946408,
722  0.001853960788946921732, 0.000796792065552012429 // 48
723  } ;
724  xMean = 0.5*(a + b) ;
725  xDiff = 0.5*(b - a) ;
726  integral = 0.0 ;
727  for(i=0;i<48;i++)
728  {
729  dx = xDiff*abscissa[i] ;
730  integral += weight[i]*((*f)(xMean + dx) + (*f)(xMean - dx)) ;
731  }
732  return integral *= xDiff ;
733 }
734 
735 //////////////////////////////////////////////////////////////////////////////
736 //
737 // Methods involving Chebyshev polynomials
738 //
739 ///////////////////////////////////////////////////////////////////////////
740 //
741 // Integrates function pointed by T::f from a to b by Gauss-Chebyshev
742 // quadrature method.
743 // Convenient for using with class object typeT
744 
745 template <class T, class F>
746 G4double G4Integrator<T,F>::Chebyshev( T& typeT, F f, G4double a,
747  G4double b, G4int nChebyshev )
748 {
749  G4int i ;
750  G4double xDiff, xMean, dx, integral = 0.0 ;
751 
752  G4int fNumber = nChebyshev ; // Try to reduce fNumber twice ??
753  G4double cof = CLHEP::pi/fNumber ;
754  G4double* fAbscissa = new G4double[fNumber] ;
755  G4double* fWeight = new G4double[fNumber] ;
756  for(i=0;i<fNumber;i++)
757  {
758  fAbscissa[i] = std::cos(cof*(i + 0.5)) ;
759  fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ;
760  }
761 
762  //
763  // Now we ready to estimate the integral
764  //
765 
766  xMean = 0.5*(a + b) ;
767  xDiff = 0.5*(b - a) ;
768  for(i=0;i<fNumber;i++)
769  {
770  dx = xDiff*fAbscissa[i] ;
771  integral += fWeight[i]*(typeT.*f)(xMean + dx) ;
772  }
773  delete[] fAbscissa;
774  delete[] fWeight;
775  return integral *= xDiff ;
776 }
777 
778 ///////////////////////////////////////////////////////////////////////
779 //
780 // Convenient for using with 'this' pointer
781 
782 template <class T, class F>
783 G4double G4Integrator<T,F>::Chebyshev( T* ptrT, F f, G4double a,
784  G4double b, G4int n )
785 {
786  return Chebyshev(*ptrT,f,a,b,n) ;
787 }
788 
789 ////////////////////////////////////////////////////////////////////////
790 //
791 // For use with global scope functions f
792 
793 template <class T, class F>
794 G4double G4Integrator<T,F>::Chebyshev( G4double (*f)(G4double),
795  G4double a, G4double b, G4int nChebyshev )
796 {
797  G4int i ;
798  G4double xDiff, xMean, dx, integral = 0.0 ;
799 
800  G4int fNumber = nChebyshev ; // Try to reduce fNumber twice ??
801  G4double cof = CLHEP::pi/fNumber ;
802  G4double* fAbscissa = new G4double[fNumber] ;
803  G4double* fWeight = new G4double[fNumber] ;
804  for(i=0;i<fNumber;i++)
805  {
806  fAbscissa[i] = std::cos(cof*(i + 0.5)) ;
807  fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ;
808  }
809 
810  //
811  // Now we ready to estimate the integral
812  //
813 
814  xMean = 0.5*(a + b) ;
815  xDiff = 0.5*(b - a) ;
816  for(i=0;i<fNumber;i++)
817  {
818  dx = xDiff*fAbscissa[i] ;
819  integral += fWeight[i]*(*f)(xMean + dx) ;
820  }
821  delete[] fAbscissa;
822  delete[] fWeight;
823  return integral *= xDiff ;
824 }
825 
826 //////////////////////////////////////////////////////////////////////
827 //
828 // Method involving Laguerre polynomials
829 //
830 //////////////////////////////////////////////////////////////////////
831 //
832 // Integral from zero to infinity of std::pow(x,alpha)*std::exp(-x)*f(x).
833 // The value of nLaguerre sets the accuracy.
834 // The function creates arrays fAbscissa[0,..,nLaguerre-1] and
835 // fWeight[0,..,nLaguerre-1] .
836 // Convenient for using with class object 'typeT' and (typeT.*f) function
837 // (T::f)
838 
839 template <class T, class F>
840 G4double G4Integrator<T,F>::Laguerre( T& typeT, F f, G4double alpha,
841  G4int nLaguerre )
842 {
843  const G4double tolerance = 1.0e-10 ;
844  const G4int maxNumber = 12 ;
845  G4int i, j, k ;
846  G4double nwt=0., nwt1, temp1, temp2, temp3, temp, cofi ;
847  G4double integral = 0.0 ;
848 
849  G4int fNumber = nLaguerre ;
850  G4double* fAbscissa = new G4double[fNumber] ;
851  G4double* fWeight = new G4double[fNumber] ;
852 
853  for(i=1;i<=fNumber;i++) // Loop over the desired roots
854  {
855  if(i == 1)
856  {
857  nwt = (1.0 + alpha)*(3.0 + 0.92*alpha)
858  / (1.0 + 2.4*fNumber + 1.8*alpha) ;
859  }
860  else if(i == 2)
861  {
862  nwt += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
863  }
864  else
865  {
866  cofi = i - 2 ;
867  nwt += ((1.0+2.55*cofi)/(1.9*cofi)
868  + 1.26*cofi*alpha/(1.0+3.5*cofi))
869  * (nwt - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
870  }
871  for(k=1;k<=maxNumber;k++)
872  {
873  temp1 = 1.0 ;
874  temp2 = 0.0 ;
875 
876  for(j=1;j<=fNumber;j++)
877  {
878  temp3 = temp2 ;
879  temp2 = temp1 ;
880  temp1 = ((2*j - 1 + alpha - nwt)*temp2 - (j - 1 + alpha)*temp3)/j ;
881  }
882  temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/nwt ;
883  nwt1 = nwt ;
884  nwt = nwt1 - temp1/temp ;
885 
886  if(std::fabs(nwt - nwt1) <= tolerance)
887  {
888  break ;
889  }
890  }
891  if(k > maxNumber)
892  {
893  G4Exception("G4Integrator<T,F>::Laguerre(T,F, ...)", "Error",
894  FatalException, "Too many (>12) iterations.");
895  }
896 
897  fAbscissa[i-1] = nwt ;
898  fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) -
899  GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
900  }
901 
902  //
903  // Integral evaluation
904  //
905 
906  for(i=0;i<fNumber;i++)
907  {
908  integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ;
909  }
910  delete[] fAbscissa;
911  delete[] fWeight;
912  return integral ;
913 }
914 
915 
916 
917 //////////////////////////////////////////////////////////////////////
918 //
919 //
920 
921 template <class T, class F> G4double
922 G4Integrator<T,F>::Laguerre( T* ptrT, F f, G4double alpha, G4int nLaguerre )
923 {
924  return Laguerre(*ptrT,f,alpha,nLaguerre) ;
925 }
926 
927 ////////////////////////////////////////////////////////////////////////
928 //
929 // For use with global scope functions f
930 
931 template <class T, class F> G4double
932 G4Integrator<T,F>::Laguerre( G4double (*f)(G4double),
933  G4double alpha, G4int nLaguerre )
934 {
935  const G4double tolerance = 1.0e-10 ;
936  const G4int maxNumber = 12 ;
937  G4int i, j, k ;
938  G4double nwt=0., nwt1, temp1, temp2, temp3, temp, cofi ;
939  G4double integral = 0.0 ;
940 
941  G4int fNumber = nLaguerre ;
942  G4double* fAbscissa = new G4double[fNumber] ;
943  G4double* fWeight = new G4double[fNumber] ;
944 
945  for(i=1;i<=fNumber;i++) // Loop over the desired roots
946  {
947  if(i == 1)
948  {
949  nwt = (1.0 + alpha)*(3.0 + 0.92*alpha)
950  / (1.0 + 2.4*fNumber + 1.8*alpha) ;
951  }
952  else if(i == 2)
953  {
954  nwt += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
955  }
956  else
957  {
958  cofi = i - 2 ;
959  nwt += ((1.0+2.55*cofi)/(1.9*cofi)
960  + 1.26*cofi*alpha/(1.0+3.5*cofi))
961  * (nwt - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
962  }
963  for(k=1;k<=maxNumber;k++)
964  {
965  temp1 = 1.0 ;
966  temp2 = 0.0 ;
967 
968  for(j=1;j<=fNumber;j++)
969  {
970  temp3 = temp2 ;
971  temp2 = temp1 ;
972  temp1 = ((2*j - 1 + alpha - nwt)*temp2 - (j - 1 + alpha)*temp3)/j ;
973  }
974  temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/nwt ;
975  nwt1 = nwt ;
976  nwt = nwt1 - temp1/temp ;
977 
978  if(std::fabs(nwt - nwt1) <= tolerance)
979  {
980  break ;
981  }
982  }
983  if(k > maxNumber)
984  {
985  G4Exception("G4Integrator<T,F>::Laguerre( ...)", "Error",
986  FatalException, "Too many (>12) iterations.");
987  }
988 
989  fAbscissa[i-1] = nwt ;
990  fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) -
991  GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
992  }
993 
994  //
995  // Integral evaluation
996  //
997 
998  for(i=0;i<fNumber;i++)
999  {
1000  integral += fWeight[i]*(*f)(fAbscissa[i]) ;
1001  }
1002  delete[] fAbscissa;
1003  delete[] fWeight;
1004  return integral ;
1005 }
1006 
1007 ///////////////////////////////////////////////////////////////////////
1008 //
1009 // Auxiliary function which returns the value of std::log(gamma-function(x))
1010 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
1011 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
1012 // (Adapted from Numerical Recipes in C)
1013 //
1014 
1015 template <class T, class F>
1016 G4double G4Integrator<T,F>::GammaLogarithm(G4double xx)
1017 {
1018  static const G4double cof[6] = { 76.18009172947146, -86.50532032941677,
1019  24.01409824083091, -1.231739572450155,
1020  0.1208650973866179e-2, -0.5395239384953e-5 };
1021  G4int j;
1022  G4double x = xx - 1.0 ;
1023  G4double tmp = x + 5.5 ;
1024  tmp -= (x + 0.5) * std::log(tmp) ;
1025  G4double ser = 1.000000000190015 ;
1026 
1027  for ( j = 0; j <= 5; j++ )
1028  {
1029  x += 1.0 ;
1030  ser += cof[j]/x ;
1031  }
1032  return -tmp + std::log(2.5066282746310005*ser) ;
1033 }
1034 
1035 ///////////////////////////////////////////////////////////////////////
1036 //
1037 // Method involving Hermite polynomials
1038 //
1039 ///////////////////////////////////////////////////////////////////////
1040 //
1041 //
1042 // Gauss-Hermite method for integration of std::exp(-x*x)*f(x)
1043 // from minus infinity to plus infinity .
1044 //
1045 
1046 template <class T, class F>
1047 G4double G4Integrator<T,F>::Hermite( T& typeT, F f, G4int nHermite )
1048 {
1049  const G4double tolerance = 1.0e-12 ;
1050  const G4int maxNumber = 12 ;
1051 
1052  G4int i, j, k ;
1053  G4double integral = 0.0 ;
1054  G4double nwt=0., nwt1, temp1, temp2, temp3, temp ;
1055 
1056  G4double piInMinusQ = std::pow(CLHEP::pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ??
1057 
1058  G4int fNumber = (nHermite +1)/2 ;
1059  G4double* fAbscissa = new G4double[fNumber] ;
1060  G4double* fWeight = new G4double[fNumber] ;
1061 
1062  for(i=1;i<=fNumber;i++)
1063  {
1064  if(i == 1)
1065  {
1066  nwt = std::sqrt((G4double)(2*nHermite + 1)) -
1067  1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
1068  }
1069  else if(i == 2)
1070  {
1071  nwt -= 1.14001*std::pow((G4double)nHermite,0.425999)/nwt ;
1072  }
1073  else if(i == 3)
1074  {
1075  nwt = 1.86002*nwt - 0.86002*fAbscissa[0] ;
1076  }
1077  else if(i == 4)
1078  {
1079  nwt = 1.91001*nwt - 0.91001*fAbscissa[1] ;
1080  }
1081  else
1082  {
1083  nwt = 2.0*nwt - fAbscissa[i - 3] ;
1084  }
1085  for(k=1;k<=maxNumber;k++)
1086  {
1087  temp1 = piInMinusQ ;
1088  temp2 = 0.0 ;
1089 
1090  for(j=1;j<=nHermite;j++)
1091  {
1092  temp3 = temp2 ;
1093  temp2 = temp1 ;
1094  temp1 = nwt*std::sqrt(2.0/j)*temp2 -
1095  std::sqrt(((G4double)(j - 1))/j)*temp3 ;
1096  }
1097  temp = std::sqrt((G4double)2*nHermite)*temp2 ;
1098  nwt1 = nwt ;
1099  nwt = nwt1 - temp1/temp ;
1100 
1101  if(std::fabs(nwt - nwt1) <= tolerance)
1102  {
1103  break ;
1104  }
1105  }
1106  if(k > maxNumber)
1107  {
1108  G4Exception("G4Integrator<T,F>::Hermite(T,F, ...)", "Error",
1109  FatalException, "Too many (>12) iterations.");
1110  }
1111  fAbscissa[i-1] = nwt ;
1112  fWeight[i-1] = 2.0/(temp*temp) ;
1113  }
1114 
1115  //
1116  // Integral calculation
1117  //
1118 
1119  for(i=0;i<fNumber;i++)
1120  {
1121  integral += fWeight[i]*( (typeT.*f)(fAbscissa[i]) +
1122  (typeT.*f)(-fAbscissa[i]) ) ;
1123  }
1124  delete[] fAbscissa;
1125  delete[] fWeight;
1126  return integral ;
1127 }
1128 
1129 
1130 ////////////////////////////////////////////////////////////////////////
1131 //
1132 // For use with 'this' pointer
1133 
1134 template <class T, class F>
1135 G4double G4Integrator<T,F>::Hermite( T* ptrT, F f, G4int n )
1136 {
1137  return Hermite(*ptrT,f,n) ;
1138 }
1139 
1140 ////////////////////////////////////////////////////////////////////////
1141 //
1142 // For use with global scope f
1143 
1144 template <class T, class F>
1145 G4double G4Integrator<T,F>::Hermite( G4double (*f)(G4double), G4int nHermite)
1146 {
1147  const G4double tolerance = 1.0e-12 ;
1148  const G4int maxNumber = 12 ;
1149 
1150  G4int i, j, k ;
1151  G4double integral = 0.0 ;
1152  G4double nwt=0., nwt1, temp1, temp2, temp3, temp ;
1153 
1154  G4double piInMinusQ = std::pow(CLHEP::pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ??
1155 
1156  G4int fNumber = (nHermite +1)/2 ;
1157  G4double* fAbscissa = new G4double[fNumber] ;
1158  G4double* fWeight = new G4double[fNumber] ;
1159 
1160  for(i=1;i<=fNumber;i++)
1161  {
1162  if(i == 1)
1163  {
1164  nwt = std::sqrt((G4double)(2*nHermite + 1)) -
1165  1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
1166  }
1167  else if(i == 2)
1168  {
1169  nwt -= 1.14001*std::pow((G4double)nHermite,0.425999)/nwt ;
1170  }
1171  else if(i == 3)
1172  {
1173  nwt = 1.86002*nwt - 0.86002*fAbscissa[0] ;
1174  }
1175  else if(i == 4)
1176  {
1177  nwt = 1.91001*nwt - 0.91001*fAbscissa[1] ;
1178  }
1179  else
1180  {
1181  nwt = 2.0*nwt - fAbscissa[i - 3] ;
1182  }
1183  for(k=1;k<=maxNumber;k++)
1184  {
1185  temp1 = piInMinusQ ;
1186  temp2 = 0.0 ;
1187 
1188  for(j=1;j<=nHermite;j++)
1189  {
1190  temp3 = temp2 ;
1191  temp2 = temp1 ;
1192  temp1 = nwt*std::sqrt(2.0/j)*temp2 -
1193  std::sqrt(((G4double)(j - 1))/j)*temp3 ;
1194  }
1195  temp = std::sqrt((G4double)2*nHermite)*temp2 ;
1196  nwt1 = nwt ;
1197  nwt = nwt1 - temp1/temp ;
1198 
1199  if(std::fabs(nwt - nwt1) <= tolerance)
1200  {
1201  break ;
1202  }
1203  }
1204  if(k > maxNumber)
1205  {
1206  G4Exception("G4Integrator<T,F>::Hermite(...)", "Error",
1207  FatalException, "Too many (>12) iterations.");
1208  }
1209  fAbscissa[i-1] = nwt ;
1210  fWeight[i-1] = 2.0/(temp*temp) ;
1211  }
1212 
1213  //
1214  // Integral calculation
1215  //
1216 
1217  for(i=0;i<fNumber;i++)
1218  {
1219  integral += fWeight[i]*( (*f)(fAbscissa[i]) + (*f)(-fAbscissa[i]) ) ;
1220  }
1221  delete[] fAbscissa;
1222  delete[] fWeight;
1223  return integral ;
1224 }
1225 
1226 ////////////////////////////////////////////////////////////////////////////
1227 //
1228 // Method involving Jacobi polynomials
1229 //
1230 ////////////////////////////////////////////////////////////////////////////
1231 //
1232 // Gauss-Jacobi method for integration of ((1-x)^alpha)*((1+x)^beta)*f(x)
1233 // from minus unit to plus unit .
1234 //
1235 
1236 template <class T, class F>
1237 G4double G4Integrator<T,F>::Jacobi( T& typeT, F f, G4double alpha,
1238  G4double beta, G4int nJacobi)
1239 {
1240  const G4double tolerance = 1.0e-12 ;
1241  const G4double maxNumber = 12 ;
1242  G4int i, k, j ;
1243  G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ;
1244  G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root=0., rootTemp ;
1245 
1246  G4int fNumber = nJacobi ;
1247  G4double* fAbscissa = new G4double[fNumber] ;
1248  G4double* fWeight = new G4double[fNumber] ;
1249 
1250  for (i=1;i<=nJacobi;i++)
1251  {
1252  if (i == 1)
1253  {
1254  alphaReduced = alpha/nJacobi ;
1255  betaReduced = beta/nJacobi ;
1256  root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
1257  0.767999*alphaReduced/nJacobi) ;
1258  root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
1259  0.451998*alphaReduced*alphaReduced +
1260  0.83001*alphaReduced*betaReduced ;
1261  root = 1.0-root1/root2 ;
1262  }
1263  else if (i == 2)
1264  {
1265  root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
1266  root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
1267  root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
1268  root -= (1.0-root)*root1*root2*root3 ;
1269  }
1270  else if (i == 3)
1271  {
1272  root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
1273  root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
1274  root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
1275  root -= (fAbscissa[0]-root)*root1*root2*root3 ;
1276  }
1277  else if (i == nJacobi-1)
1278  {
1279  root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
1280  root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
1281  root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
1282  root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
1283  }
1284  else if (i == nJacobi)
1285  {
1286  root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
1287  root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
1288  root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
1289  root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
1290  }
1291  else
1292  {
1293  root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
1294  }
1295  alphaBeta = alpha + beta ;
1296  for (k=1;k<=maxNumber;k++)
1297  {
1298  temp = 2.0 + alphaBeta ;
1299  nwt1 = (alpha-beta+temp*root)/2.0 ;
1300  nwt2 = 1.0 ;
1301  for (j=2;j<=nJacobi;j++)
1302  {
1303  nwt3 = nwt2 ;
1304  nwt2 = nwt1 ;
1305  temp = 2*j+alphaBeta ;
1306  a = 2*j*(j+alphaBeta)*(temp-2.0) ;
1307  b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
1308  c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
1309  nwt1 = (b*nwt2-c*nwt3)/a ;
1310  }
1311  nwt = (nJacobi*(alpha - beta - temp*root)*nwt1 +
1312  2.0*(nJacobi + alpha)*(nJacobi + beta)*nwt2)/
1313  (temp*(1.0 - root*root)) ;
1314  rootTemp = root ;
1315  root = rootTemp - nwt1/nwt ;
1316  if (std::fabs(root-rootTemp) <= tolerance)
1317  {
1318  break ;
1319  }
1320  }
1321  if (k > maxNumber)
1322  {
1323  G4Exception("G4Integrator<T,F>::Jacobi(T,F, ...)", "Error",
1324  FatalException, "Too many (>12) iterations.");
1325  }
1326  fAbscissa[i-1] = root ;
1327  fWeight[i-1] = std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) +
1328  GammaLogarithm((G4double)(beta+nJacobi)) -
1329  GammaLogarithm((G4double)(nJacobi+1.0)) -
1330  GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
1331  *temp*std::pow(2.0,alphaBeta)/(nwt*nwt2) ;
1332  }
1333 
1334  //
1335  // Calculation of the integral
1336  //
1337 
1338  G4double integral = 0.0 ;
1339  for(i=0;i<fNumber;i++)
1340  {
1341  integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ;
1342  }
1343  delete[] fAbscissa;
1344  delete[] fWeight;
1345  return integral ;
1346 }
1347 
1348 
1349 /////////////////////////////////////////////////////////////////////////
1350 //
1351 // For use with 'this' pointer
1352 
1353 template <class T, class F>
1354 G4double G4Integrator<T,F>::Jacobi( T* ptrT, F f, G4double alpha,
1355  G4double beta, G4int n)
1356 {
1357  return Jacobi(*ptrT,f,alpha,beta,n) ;
1358 }
1359 
1360 /////////////////////////////////////////////////////////////////////////
1361 //
1362 // For use with global scope f
1363 
1364 template <class T, class F>
1365 G4double G4Integrator<T,F>::Jacobi( G4double (*f)(G4double), G4double alpha,
1366  G4double beta, G4int nJacobi)
1367 {
1368  const G4double tolerance = 1.0e-12 ;
1369  const G4double maxNumber = 12 ;
1370  G4int i, k, j ;
1371  G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ;
1372  G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root=0., rootTemp ;
1373 
1374  G4int fNumber = nJacobi ;
1375  G4double* fAbscissa = new G4double[fNumber] ;
1376  G4double* fWeight = new G4double[fNumber] ;
1377 
1378  for (i=1;i<=nJacobi;i++)
1379  {
1380  if (i == 1)
1381  {
1382  alphaReduced = alpha/nJacobi ;
1383  betaReduced = beta/nJacobi ;
1384  root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
1385  0.767999*alphaReduced/nJacobi) ;
1386  root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
1387  0.451998*alphaReduced*alphaReduced +
1388  0.83001*alphaReduced*betaReduced ;
1389  root = 1.0-root1/root2 ;
1390  }
1391  else if (i == 2)
1392  {
1393  root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
1394  root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
1395  root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
1396  root -= (1.0-root)*root1*root2*root3 ;
1397  }
1398  else if (i == 3)
1399  {
1400  root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
1401  root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
1402  root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
1403  root -= (fAbscissa[0]-root)*root1*root2*root3 ;
1404  }
1405  else if (i == nJacobi-1)
1406  {
1407  root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
1408  root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
1409  root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
1410  root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
1411  }
1412  else if (i == nJacobi)
1413  {
1414  root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
1415  root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
1416  root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
1417  root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
1418  }
1419  else
1420  {
1421  root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
1422  }
1423  alphaBeta = alpha + beta ;
1424  for (k=1;k<=maxNumber;k++)
1425  {
1426  temp = 2.0 + alphaBeta ;
1427  nwt1 = (alpha-beta+temp*root)/2.0 ;
1428  nwt2 = 1.0 ;
1429  for (j=2;j<=nJacobi;j++)
1430  {
1431  nwt3 = nwt2 ;
1432  nwt2 = nwt1 ;
1433  temp = 2*j+alphaBeta ;
1434  a = 2*j*(j+alphaBeta)*(temp-2.0) ;
1435  b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
1436  c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
1437  nwt1 = (b*nwt2-c*nwt3)/a ;
1438  }
1439  nwt = (nJacobi*(alpha - beta - temp*root)*nwt1 +
1440  2.0*(nJacobi + alpha)*(nJacobi + beta)*nwt2) /
1441  (temp*(1.0 - root*root)) ;
1442  rootTemp = root ;
1443  root = rootTemp - nwt1/nwt ;
1444  if (std::fabs(root-rootTemp) <= tolerance)
1445  {
1446  break ;
1447  }
1448  }
1449  if (k > maxNumber)
1450  {
1451  G4Exception("G4Integrator<T,F>::Jacobi(...)", "Error",
1452  FatalException, "Too many (>12) iterations.");
1453  }
1454  fAbscissa[i-1] = root ;
1455  fWeight[i-1] =
1456  std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) +
1457  GammaLogarithm((G4double)(beta+nJacobi)) -
1458  GammaLogarithm((G4double)(nJacobi+1.0)) -
1459  GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
1460  *temp*std::pow(2.0,alphaBeta)/(nwt*nwt2);
1461  }
1462 
1463  //
1464  // Calculation of the integral
1465  //
1466 
1467  G4double integral = 0.0 ;
1468  for(i=0;i<fNumber;i++)
1469  {
1470  integral += fWeight[i]*(*f)(fAbscissa[i]) ;
1471  }
1472  delete[] fAbscissa;
1473  delete[] fWeight;
1474  return integral ;
1475 }
1476 
1477 //
1478 //
1479 ///////////////////////////////////////////////////////////////////