Geant4  10.01.p03
G4SimpleIntegration.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4SimpleIntegration.cc 69546 2013-05-08 09:50:34Z gcosmo $
28 //
29 // Implementation file for simple integration methods
30 //
31 
32 #include "globals.hh"
33 #include "G4SimpleIntegration.hh"
34 
35 
37  : fFunction(pFunction),
38  fTolerance(.0001),
39  fMaxDepth(100)
40 {
41 }
42 
44  G4double pTolerance)
45  : fFunction(pFunction),
46  fTolerance(pTolerance),
47  fMaxDepth(100)
48 {
49 }
50 
51 
53 {
54 }
55 
56  // Simple integration methods
57 
60  G4double xFinal,
61  G4int iterationNumber )
62 {
63  G4double Step = (xFinal - xInitial)/iterationNumber ;
64  G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ;
65  G4double x = xInitial ;
66  for(G4int i=1;i<iterationNumber;i++)
67  {
68  x += Step ;
69  mean += fFunction(x) ;
70  }
71  return mean*Step ;
72 }
73 
74 G4double
76  G4double xFinal,
77  G4int iterationNumber )
78 {
79  G4double Step = (xFinal - xInitial)/iterationNumber ;
80  G4double x = xInitial + 0.5*Step;
81  G4double mean = fFunction(x) ;
82  for(G4int i=1;i<iterationNumber;i++)
83  {
84  x += Step ;
85  mean += fFunction(x) ;
86  }
87  return mean*Step ;
88 }
89 
90 G4double
92  G4double xFinal,
93  G4int iterationNumber )
94 {
95  G4double x=0.;
96  static const G4double root = 1.0/std::sqrt(3.0) ;
97  G4double Step = (xFinal - xInitial)/(2.0*iterationNumber) ;
98  G4double delta = Step*root ;
99  G4double mean = 0.0 ;
100  for(G4int i=0;i<iterationNumber;i++)
101  {
102  x = (2*i + 1)*Step ;
103  mean += (fFunction(x+delta) + fFunction(x-delta)) ;
104  }
105  return mean*Step ;
106 }
107 
108 G4double
110  G4double xFinal,
111  G4int iterationNumber )
112 {
113  G4double Step = (xFinal - xInitial)/iterationNumber ;
114  G4double x = xInitial ;
115  G4double xPlus = xInitial + 0.5*Step ;
116  G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ;
117  G4double sum = fFunction(xPlus) ;
118  for(G4int i=1;i<iterationNumber;i++)
119  {
120  x += Step ;
121  xPlus += Step ;
122  mean += fFunction(x) ;
123  sum += fFunction(xPlus) ;
124  }
125  mean += 2.0*sum ;
126  return mean*Step/3.0 ;
127 }
128 
129 
130 
131  // Adaptive Gauss integration
132 
133 G4double
135  G4double xFinal )
136 {
137  G4int depth = 0 ;
138  G4double sum = 0.0 ;
139  AdaptGauss(xInitial,xFinal,sum,depth) ;
140  return sum ;
141 }
142 
143 
144 G4double
146  G4double xFinal )
147 {
148  static const G4double root = 1.0/std::sqrt(3.0) ;
149 
150  G4double xMean = (xInitial + xFinal)/2.0 ;
151  G4double Step = (xFinal - xInitial)/2.0 ;
152  G4double delta = Step*root ;
153  G4double sum = (fFunction(xMean + delta) + fFunction(xMean - delta)) ;
154 
155  return sum*Step ;
156 }
157 
158 
159 void
161  G4double xFinal,
162  G4double& sum,
163  G4int& depth )
164 {
165  if(depth >fMaxDepth)
166  {
167  G4Exception("G4SimpleIntegration::AdaptGauss()", "Error",
168  FatalException, "Function varies too rapidly !") ;
169  }
170  G4double xMean = (xInitial + xFinal)/2.0 ;
171  G4double leftHalf = Gauss(xInitial,xMean) ;
172  G4double rightHalf = Gauss(xMean,xFinal) ;
173  G4double full = Gauss(xInitial,xFinal) ;
174  if(std::fabs(leftHalf+rightHalf-full) < fTolerance)
175  {
176  sum += full ;
177  }
178  else
179  {
180  depth++ ;
181  AdaptGauss(xInitial,xMean,sum,depth) ;
182  AdaptGauss(xMean,xFinal,sum,depth) ;
183  }
184 }
G4double MidPoint(G4double xInitial, G4double xFinal, G4int iterationNumber)
G4SimpleIntegration(function pFunction)
G4double Gauss(G4double xInitial, G4double xFinal, G4int iterationNumber)
G4double AdaptGaussIntegration(G4double xInitial, G4double xFinal)
int G4int
Definition: G4Types.hh:78
void AdaptGauss(G4double xInitial, G4double xFinal, G4double &sum, G4int &depth)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double Simpson(G4double xInitial, G4double xFinal, G4int iterationNumber)
Definition: Step.hh:41
G4double Trapezoidal(G4double xInitial, G4double xFinal, G4int iterationNumber)
double G4double
Definition: G4Types.hh:76