Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/include/G4Integrator.icc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /global/HEPNumerics/include/G4Integrator.icc (Version 11.3.0) and /global/HEPNumerics/include/G4Integrator.icc (Version 7.1)


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