Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/management/include/G4Pow.hh

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 ]

  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 // G4Pow
 27 //
 28 // Class description:
 29 //
 30 // Utility singleton class for the fast computation of log and pow
 31 // functions. Integer argument should be in the interval 0-512, no
 32 // check is performed inside these methods for performance reasons.
 33 // For factorial integer argument should be in the interval 0-170
 34 // Computations with double arguments are fast for the interval
 35 // 0.002-511.5 for all functions except exponent, which is computed
 36 // for the interval 0-84.4, standard library is used in the opposite case
 37 
 38 // Author: Vladimir Ivanchenko, 23.05.2009
 39 // --------------------------------------------------------------------
 40 #ifndef G4Pow_hh
 41 #define G4Pow_hh 1
 42 
 43 #include "G4DataVector.hh"
 44 #include "G4Exp.hh"
 45 #include "G4Log.hh"
 46 #include "globals.hh"
 47 
 48 class G4Pow
 49 {
 50  public:
 51   static G4Pow* GetInstance();
 52   ~G4Pow() = default;
 53 
 54   // Fast computation of Z^1/3
 55   //
 56   inline G4double Z13(G4int Z) const;
 57   G4double A13(G4double A) const;
 58 
 59   // Fast computation of Z^2/3
 60   //
 61   inline G4double Z23(G4int Z) const;
 62   inline G4double A23(G4double A) const;
 63 
 64   // Fast computation of log(Z)
 65   //
 66   inline G4double logZ(G4int Z) const;
 67   inline G4double logA(G4double A) const;
 68   inline G4double logX(G4double x) const;
 69 
 70   // Fast computation of log10(Z)
 71   //
 72   inline G4double log10Z(G4int Z) const;
 73   inline G4double log10A(G4double A) const;
 74 
 75   // Fast computation of exp(X)
 76   //
 77   inline G4double expA(G4double A) const;
 78 
 79   // Fast computation of pow(Z,X)
 80   //
 81   inline G4double powZ(G4int Z, G4double y) const;
 82   inline G4double powA(G4double A, G4double y) const;
 83   G4double powN(G4double x, G4int n) const;
 84 
 85   // Fast factorial
 86   //
 87   inline G4double factorial(G4int Z) const;
 88   inline G4double logfactorial(G4int Z) const;
 89 
 90  private:
 91   G4Pow();
 92 
 93   G4double A13Low(const G4double, const G4bool) const;
 94   G4double A13High(const G4double, const G4bool) const;
 95 
 96   inline G4double logBase(G4double x) const;
 97 
 98   static G4Pow* fpInstance;
 99 
100   const G4double onethird = 1.0 / 3.0;
101   const G4int max2        = 5;
102 
103   G4double maxA;
104   G4double maxLowA;
105   G4double maxA2;
106   G4double maxAexp;
107 
108   G4DataVector ener;
109   G4DataVector logen;
110   G4DataVector pz13;
111   G4DataVector lowa13;
112   G4DataVector lz;
113   G4DataVector lz2;
114   G4DataVector fexp;
115   G4DataVector fact;
116   G4DataVector logfact;
117 };
118 
119 // -----------------------------
120 // Inline methods implementation
121 // -----------------------------
122 
123 inline G4double G4Pow::Z13(G4int Z) const { return pz13[Z]; }
124 
125 inline G4double G4Pow::Z23(G4int Z) const
126 {
127   G4double x = Z13(Z);
128   return x * x;
129 }
130 
131 inline G4double G4Pow::A23(G4double A) const
132 {
133   G4double x = A13(A);
134   return x * x;
135 }
136 
137 inline G4double G4Pow::logZ(G4int Z) const { return lz[Z]; }
138 
139 inline G4double G4Pow::logBase(G4double a) const
140 {
141   G4double res;
142   if(a <= maxA2)
143   {
144     G4int i = G4int(max2 * (a - 1) + 0.5);
145     if(i > max2)
146     {
147       i = max2;
148     }
149     G4double x = a / (G4double(i) / max2 + 1) - 1;
150     res        = lz2[i] + x * (1.0 - (0.5 - onethird * x) * x);
151   }
152   else if(a <= maxA)
153   {
154     G4int i    = G4int(a + 0.5);
155     G4double x = a / G4double(i) - 1;
156     res        = lz[i] + x * (1.0 - (0.5 - onethird * x) * x);
157   }
158   else
159   {
160     res = G4Log(a);
161   }
162   return res;
163 }
164 
165 inline G4double G4Pow::logA(G4double A) const
166 {
167   return (1.0 <= A ? logBase(A) : -logBase(1. / A));
168 }
169 
170 inline G4double G4Pow::logX(G4double x) const
171 {
172   G4double res = 0.0;
173   G4double a   = (1.0 <= x) ? x : 1.0 / x;
174 
175   if(a <= maxA)
176   {
177     res = logBase(a);
178   }
179   else if(a <= ener[2])
180   {
181     res = logen[1] + logBase(a / ener[1]);
182   }
183   else if(a <= ener[3])
184   {
185     res = logen[2] + logBase(a / ener[2]);
186   }
187   else
188   {
189     res = G4Log(a);
190   }
191 
192   if(1.0 > x)
193   {
194     res = -res;
195   }
196   return res;
197 }
198 
199 inline G4double G4Pow::log10Z(G4int Z) const { return lz[Z] / lz[10]; }
200 
201 inline G4double G4Pow::log10A(G4double A) const { return logX(A) / lz[10]; }
202 
203 inline G4double G4Pow::expA(G4double A) const
204 {
205   G4double res;
206   G4double a = (0.0 <= A) ? A : -A;
207 
208   if(a <= maxAexp)
209   {
210     G4int i    = G4int(2 * a + 0.5);
211     G4double x = a - i * 0.5;
212     res        = fexp[i] * (1.0 + x * (1.0 + 0.5 * (1.0 + onethird * x) * x));
213   }
214   else
215   {
216     res = G4Exp(a);
217   }
218   if(0.0 > A)
219   {
220     res = 1.0 / res;
221   }
222   return res;
223 }
224 
225 inline G4double G4Pow::powZ(G4int Z, G4double y) const
226 {
227   return expA(y * lz[Z]);
228 }
229 
230 inline G4double G4Pow::powA(G4double A, G4double y) const
231 {
232   return (0.0 == A ? 0.0 : expA(y * logX(A)));
233 }
234 
235 inline G4double G4Pow::factorial(G4int Z) const { return fact[Z]; }
236 
237 inline G4double G4Pow::logfactorial(G4int Z) const { return logfact[Z]; }
238 
239 #endif
240