Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/src/G4SimpleIntegration.cc

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/src/G4SimpleIntegration.cc (Version 11.3.0) and /global/HEPNumerics/src/G4SimpleIntegration.cc (Version 6.2.p1)


  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 // G4SimpleIntegration class implementation    << 
 27 //                                                 23 //
 28 // Author: V.Grichine, 26.03.1997              <<  24 // $Id: G4SimpleIntegration.cc,v 1.3 2001/07/11 10:00:42 gunter Exp $
 29 // ------------------------------------------- <<  25 // GEANT4 tag $Name: geant4-05-02-patch-01 $
                                                   >>  26 //
                                                   >>  27 // Implementation file for simple integration methods
                                                   >>  28 //
 30                                                    29 
 31 #include "G4SimpleIntegration.hh"                  30 #include "G4SimpleIntegration.hh"
 32 #include "globals.hh"                          << 
 33                                                    31 
 34 G4SimpleIntegration::G4SimpleIntegration(funct <<  32 
 35   : fFunction(pFunction)                       <<  33 G4int G4SimpleIntegration::fMaxDepth = 100 ;
 36 {}                                             <<  34 
 37                                                <<  35 
 38 G4SimpleIntegration::G4SimpleIntegration(funct <<  36 G4SimpleIntegration::G4SimpleIntegration( function pFunction )
 39                                          G4dou <<  37 {
 40   : fFunction(pFunction)                       <<  38    fFunction = pFunction ;
 41   , fTolerance(pTolerance)                     <<  39 }
 42 {}                                             <<  40 
 43                                                <<  41 G4SimpleIntegration::G4SimpleIntegration( function pFunction,
 44 // Simple integration methods                  <<  42             G4double pTolerance)
 45                                                <<  43 {
 46 G4double G4SimpleIntegration::Trapezoidal(G4do <<  44    fFunction = pFunction ;
 47                                           G4in <<  45    fTolerance = pTolerance ;
 48 {                                              <<  46 }
 49   G4double Step = (xFinal - xInitial) / iterat <<  47 
 50   G4double mean = (fFunction(xInitial) + fFunc <<  48 
 51   G4double x    = xInitial;                    <<  49 G4SimpleIntegration::~G4SimpleIntegration() 
 52   for(G4int i = 1; i < iterationNumber; ++i)   <<  50 {
 53   {                                            <<  51    ;
 54     x += Step;                                 <<  52 }
 55     mean += fFunction(x);                      <<  53        
 56   }                                            <<  54        // Simple integration methods
 57   return mean * Step;                          <<  55        
 58 }                                              <<  56 G4double
 59                                                <<  57 G4SimpleIntegration::Trapezoidal(G4double xInitial,
 60 G4double G4SimpleIntegration::MidPoint(G4doubl <<  58                                  G4double xFinal,
 61                                        G4int i <<  59                G4int iterationNumber ) 
 62 {                                              <<  60 {
 63   G4double Step = (xFinal - xInitial) / iterat <<  61    G4int i ;
 64   G4double x    = xInitial + 0.5 * Step;       <<  62    G4double Step = (xFinal - xInitial)/iterationNumber ;
 65   G4double mean = fFunction(x);                <<  63    G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ;
 66   for(G4int i = 1; i < iterationNumber; ++i)   <<  64    G4double x = xInitial ;
 67   {                                            <<  65    for(i=1;i<iterationNumber;i++)
 68     x += Step;                                 <<  66    {
 69     mean += fFunction(x);                      <<  67       x += Step ;
 70   }                                            <<  68       mean += fFunction(x) ;
 71   return mean * Step;                          <<  69    }
 72 }                                              <<  70    return mean*Step ;
 73                                                <<  71 }
 74 G4double G4SimpleIntegration::Gauss(G4double x <<  72 
 75                                     G4int iter <<  73 G4double 
 76 {                                              <<  74 G4SimpleIntegration::MidPoint(G4double xInitial,
 77   G4double x                 = 0.;             <<  75                               G4double xFinal,
 78   static const G4double root = 1.0 / std::sqrt <<  76             G4int iterationNumber ) 
 79   G4double Step              = (xFinal - xInit <<  77 {
 80   G4double delta             = Step * root;    <<  78    G4int i ;
 81   G4double mean              = 0.0;            <<  79    G4double Step = (xFinal - xInitial)/iterationNumber ;
 82   for(G4int i = 0; i < iterationNumber; ++i)   <<  80    G4double x = xInitial + 0.5*Step;
 83   {                                            <<  81    G4double mean = fFunction(x) ;
 84     x = (2 * i + 1) * Step;                    <<  82    for(i=1;i<iterationNumber;i++)
 85     mean += (fFunction(x + delta) + fFunction( <<  83    {
 86   }                                            <<  84       x += Step ;
 87   return mean * Step;                          <<  85       mean += fFunction(x) ;
 88 }                                              <<  86    }
 89                                                <<  87    return mean*Step ;   
 90 G4double G4SimpleIntegration::Simpson(G4double <<  88 }
 91                                       G4int it <<  89 
 92 {                                              <<  90 G4double      
 93   G4double Step  = (xFinal - xInitial) / itera <<  91 G4SimpleIntegration::Gauss(G4double xInitial,
 94   G4double x     = xInitial;                   <<  92                            G4double xFinal,
 95   G4double xPlus = xInitial + 0.5 * Step;      <<  93          G4int iterationNumber ) 
 96   G4double mean  = (fFunction(xInitial) + fFun <<  94 {
 97   G4double sum   = fFunction(xPlus);           <<  95    G4int i ;
 98   for(G4int i = 1; i < iterationNumber; ++i)   <<  96    G4double x ;
 99   {                                            <<  97    static G4double root = 1.0/sqrt(3.0) ;
100     x += Step;                                 <<  98    G4double Step = (xFinal - xInitial)/(2.0*iterationNumber) ;
101     xPlus += Step;                             <<  99    G4double delta = Step*root ;
102     mean += fFunction(x);                      << 100    G4double mean = 0.0 ;
103     sum += fFunction(xPlus);                   << 101    for(i=0;i<iterationNumber;i++)
104   }                                            << 102    {
105   mean += 2.0 * sum;                           << 103       x = (2*i + 1)*Step ;
106   return mean * Step / 3.0;                    << 104       mean += (fFunction(x+delta) + fFunction(x-delta)) ;
107 }                                              << 105    }
108                                                << 106    return mean*Step ;   
109 // Adaptive Gauss integration                  << 107 }
110                                                << 108 
111 G4double G4SimpleIntegration::AdaptGaussIntegr << 109 G4double    
112                                                << 110 G4SimpleIntegration::Simpson(G4double xInitial,
113 {                                              << 111                              G4double xFinal,
114   G4int depth  = 0;                            << 112            G4int iterationNumber ) 
115   G4double sum = 0.0;                          << 113 {
116   AdaptGauss(xInitial, xFinal, sum, depth);    << 114    G4int i ;
117   return sum;                                  << 115    G4double Step = (xFinal - xInitial)/iterationNumber ;
118 }                                              << 116    G4double x = xInitial ;
119                                                << 117    G4double xPlus = xInitial + 0.5*Step ;
120 G4double G4SimpleIntegration::Gauss(G4double x << 118    G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ;
121 {                                              << 119    G4double sum = fFunction(xPlus) ;
122   static const G4double root = 1.0 / std::sqrt << 120    for(i=1;i<iterationNumber;i++)
123                                                << 121    {
124   G4double xMean = (xInitial + xFinal) / 2.0;  << 122       x     += Step ;
125   G4double Step  = (xFinal - xInitial) / 2.0;  << 123       xPlus += Step ;
126   G4double delta = Step * root;                << 124       mean  += fFunction(x) ;
127   G4double sum   = (fFunction(xMean + delta) + << 125       sum   += fFunction(xPlus) ;
128                                                << 126    }
129   return sum * Step;                           << 127    mean += 2.0*sum ;
130 }                                              << 128    return mean*Step/3.0 ;   
131                                                << 129 }
132 void G4SimpleIntegration::AdaptGauss(G4double  << 130 
133                                      G4double& << 131 
134 {                                              << 132 
135   if(depth > fMaxDepth)                        << 133        // Adaptive Gauss integration
136   {                                            << 134        
137     G4Exception("G4SimpleIntegration::AdaptGau << 135 G4double       
138                 "Function varies too rapidly ! << 136 G4SimpleIntegration::AdaptGaussIntegration( G4double xInitial,
139   }                                            << 137                                             G4double xFinal   ) 
140   G4double xMean     = (xInitial + xFinal) / 2 << 138 {
141   G4double leftHalf  = Gauss(xInitial, xMean); << 139    G4int depth = 0 ;
142   G4double rightHalf = Gauss(xMean, xFinal);   << 140    G4double sum = 0.0 ;
143   G4double full      = Gauss(xInitial, xFinal) << 141    AdaptGauss(xInitial,xFinal,sum,depth) ;
144   if(std::fabs(leftHalf + rightHalf - full) <  << 142    return sum ;
145   {                                            << 143 }
146     sum += full;                               << 144           
147   }                                            << 145 
148   else                                         << 146 G4double
149   {                                            << 147 G4SimpleIntegration::Gauss( G4double xInitial,
150     ++depth;                                   << 148                             G4double xFinal   ) 
151     AdaptGauss(xInitial, xMean, sum, depth);   << 149 {
152     AdaptGauss(xMean, xFinal, sum, depth);     << 150    static G4double root = 1.0/sqrt(3.0) ;
153   }                                            << 151    
                                                   >> 152    G4double xMean = (xInitial + xFinal)/2.0 ;
                                                   >> 153    G4double Step = (xFinal - xInitial)/2.0 ;
                                                   >> 154    G4double delta = Step*root ;
                                                   >> 155    G4double sum = (fFunction(xMean + delta) + fFunction(xMean - delta)) ;
                                                   >> 156    
                                                   >> 157    return sum*Step ;   
                                                   >> 158 }
                                                   >> 159 
                                                   >> 160 
                                                   >> 161 void    
                                                   >> 162 G4SimpleIntegration::AdaptGauss( G4double xInitial,
                                                   >> 163                                  G4double xFinal,
                                                   >> 164          G4double& sum,
                                                   >> 165                G4int& depth      ) 
                                                   >> 166 {
                                                   >> 167    if(depth >fMaxDepth)
                                                   >> 168    {
                                                   >> 169       G4Exception("Function varies too rapidly in G4SimpleIntegration::AdaptGauss") ;
                                                   >> 170    }
                                                   >> 171    G4double xMean = (xInitial + xFinal)/2.0 ;
                                                   >> 172    G4double leftHalf = Gauss(xInitial,xMean) ;
                                                   >> 173    G4double rightHalf = Gauss(xMean,xFinal) ;
                                                   >> 174    G4double full = Gauss(xInitial,xFinal) ;
                                                   >> 175    if(fabs(leftHalf+rightHalf-full) < fTolerance)
                                                   >> 176    {
                                                   >> 177       sum += full ;
                                                   >> 178    }
                                                   >> 179    else
                                                   >> 180    {
                                                   >> 181       depth++ ;
                                                   >> 182       AdaptGauss(xInitial,xMean,sum,depth) ;
                                                   >> 183       AdaptGauss(xMean,xFinal,sum,depth) ;
                                                   >> 184    }
154 }                                                 185 }
155                                                   186