Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/parameterisations/gflash/src/Gamma.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 /parameterisations/gflash/src/Gamma.cc (Version 11.3.0) and /parameterisations/gflash/src/Gamma.cc (Version 9.4.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: Gamma.cc,v 1.6 2006-06-29 19:14:28 gunter Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-09-04-patch-01 $
 26 //                                                 28 //
 27 //                                                 29 //
 28 // -------------------------------------------     30 // ------------------------------------------------------------
 29 // GEANT 4 class implementation                    31 // GEANT 4 class implementation
 30 // -------------------------------------------     32 // ------------------------------------------------------------
 31                                                    33 
 32 #include <cmath>                                   34 #include <cmath>
 33 #include <string.h>                                35 #include <string.h>
 34 #include "Gamma.hh"                                36 #include "Gamma.hh"
 35                                                    37 
 36 MyGamma::MyGamma(){}                               38 MyGamma::MyGamma(){}
 37                                                    39 
 38 MyGamma::~MyGamma(){}                              40 MyGamma::~MyGamma(){}
 39                                                    41 
 40 //____________________________________________     42 //____________________________________________________________________________
 41 double MyGamma::Gamma(double z)                    43 double MyGamma::Gamma(double z)
 42 {                                                  44 {
 43   if (z <= 0)                                  <<  45   // Computation of gamma(z) for all z>0.
 44       return 0;                                <<  46   //
 45                                                <<  47   // The algorithm is based on the article by C.Lanczos [1] as denoted in
 46   return std::tgamma(z);                       <<  48   // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
                                                   >>  49   //
                                                   >>  50   // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
                                                   >>  51   //
                                                   >>  52   //--- Nve 14-nov-1998 UU-SAP Utrecht
                                                   >>  53   
                                                   >>  54   if (z<=0) return 0;
                                                   >>  55   
                                                   >>  56   double v = LnGamma(z);
                                                   >>  57   return std::exp(v);
 47 }                                                  58 }
 48                                                    59 
 49 //____________________________________________     60 //____________________________________________________________________________
 50 double MyGamma::Gamma(double a, double x)      <<  61 double MyGamma::Gamma(double a,double x)
 51 {                                                  62 {
 52   // Computation of the incomplete gamma funct     63   // Computation of the incomplete gamma function P(a,x)
 53   //                                               64   //
 54   // The algorithm is based on the formulas an     65   // The algorithm is based on the formulas and code as denoted in
 55   // Numerical Recipes 2nd ed. on p. 210-212 (     66   // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
 56   //                                               67   //
 57   //--- Nve 14-nov-1998 UU-SAP Utrecht             68   //--- Nve 14-nov-1998 UU-SAP Utrecht
 58                                                <<  69   
 59   if (a <= 0 || x <= 0) return 0;                  70   if (a <= 0 || x <= 0) return 0;
 60                                                <<  71   
 61   if (x < (a + 1))                             <<  72   if (x < (a+1)) return GamSer(a,x);
 62     return GamSer(a, x);                       <<  73   else           return GamCf(a,x);
 63   else                                         << 
 64     return GamCf(a, x);                        << 
 65 }                                                  74 }
 66                                                    75 
 67 //____________________________________________     76 //____________________________________________________________________________
 68 double MyGamma::GamCf(double a, double x)      <<  77 double MyGamma::GamCf(double a,double x)
 69 {                                                  78 {
 70   // Computation of the incomplete gamma funct     79   // Computation of the incomplete gamma function P(a,x)
 71   // via its continued fraction representation     80   // via its continued fraction representation.
 72   //                                               81   //
 73   // The algorithm is based on the formulas an     82   // The algorithm is based on the formulas and code as denoted in
 74   // Numerical Recipes 2nd ed. on p. 210-212 (     83   // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
 75   //                                               84   //
 76   //--- Nve 14-nov-1998 UU-SAP Utrecht             85   //--- Nve 14-nov-1998 UU-SAP Utrecht
 77                                                <<  86   
 78   int itmax = 100;  // Maximum number of itera <<  87   int itmax    = 100;      // Maximum number of iterations
 79   double eps = 3.e-7;  // Relative accuracy    <<  88   double eps   = 3.e-7;    // Relative accuracy
 80   double fpmin = 1.e-30;  // Smallest double v <<  89   double fpmin = 1.e-30;   // Smallest double value allowed here
 81                                                <<  90   
 82   if (a <= 0 || x <= 0) return 0;                  91   if (a <= 0 || x <= 0) return 0;
 83                                                <<  92   
 84   double gln = LnGamma(a);                         93   double gln = LnGamma(a);
 85   double b = x + 1 - a;                        <<  94   double b   = x+1-a;
 86   double c = 1 / fpmin;                        <<  95   double c   = 1/fpmin;
 87   double d = 1 / b;                            <<  96   double d   = 1/b;
 88   double h = d;                                <<  97   double h   = d;
 89   double an, del;                              <<  98   double an,del;
 90   for (int i = 1; i <= itmax; i++) {           <<  99   for (int i=1; i<=itmax; i++) {
 91     an = double(-i) * (double(i) - a);         << 100     an = double(-i)*(double(i)-a);
 92     b += 2;                                       101     b += 2;
 93     d = an * d + b;                            << 102     d  = an*d+b;
 94     if (Abs(d) < fpmin) d = fpmin;                103     if (Abs(d) < fpmin) d = fpmin;
 95     c = b + an / c;                            << 104     c = b+an/c;
 96     if (Abs(c) < fpmin) c = fpmin;                105     if (Abs(c) < fpmin) c = fpmin;
 97     d = 1 / d;                                 << 106     d   = 1/d;
 98     del = d * c;                               << 107     del = d*c;
 99     h = h * del;                               << 108     h   = h*del;
100     if (Abs(del - 1) < eps) break;             << 109     if (Abs(del-1) < eps) break;
101     // if (i==itmax) cout << "*GamCf(a,x)* a t << 110     //if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
102   }                                               111   }
103   double v = Exp(-x + a * Log(x) - gln) * h;   << 112   double v = Exp(-x+a*Log(x)-gln)*h;
104   return (1 - v);                              << 113   return (1-v);
105 }                                                 114 }
106                                                   115 
107 //____________________________________________    116 //____________________________________________________________________________
108 double MyGamma::GamSer(double a, double x)     << 117 double MyGamma::GamSer(double a,double x)
109 {                                                 118 {
110   // Computation of the incomplete gamma funct    119   // Computation of the incomplete gamma function P(a,x)
111   // via its series representation.               120   // via its series representation.
112   //                                              121   //
113   // The algorithm is based on the formulas an    122   // The algorithm is based on the formulas and code as denoted in
114   // Numerical Recipes 2nd ed. on p. 210-212 (    123   // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
115   //                                              124   //
116   //--- Nve 14-nov-1998 UU-SAP Utrecht            125   //--- Nve 14-nov-1998 UU-SAP Utrecht
117                                                << 126   
118   int itmax = 100;  // Maximum number of itera << 127   int itmax  = 100;   // Maximum number of iterations
119   double eps = 3.e-7;  // Relative accuracy    << 128   double eps = 3.e-7; // Relative accuracy
120                                                << 129   
121   if (a <= 0 || x <= 0) return 0;                 130   if (a <= 0 || x <= 0) return 0;
122                                                << 131   
123   double gln = LnGamma(a);                        132   double gln = LnGamma(a);
124   double ap = a;                               << 133   double ap  = a;
125   double sum = 1 / a;                          << 134   double sum = 1/a;
126   double del = sum;                               135   double del = sum;
127   for (int n = 1; n <= itmax; n++) {           << 136   for (int n=1; n<=itmax; n++) {
128     ap += 1;                                   << 137     ap  += 1;
129     del = del * x / ap;                        << 138     del  = del*x/ap;
130     sum += del;                                   139     sum += del;
131     if (MyGamma::Abs(del) < Abs(sum * eps)) br << 140     if (MyGamma::Abs(del) < Abs(sum*eps)) break;
132     // if (n==itmax) cout << "*GamSer(a,x)* a  << 141     //if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
133   }                                               142   }
134   double v = sum * Exp(-x + a * Log(x) - gln); << 143   double v = sum*Exp(-x+a*Log(x)-gln);
135   return v;                                       144   return v;
136 }                                                 145 }
137                                                   146 
                                                   >> 147 
138 double MyGamma::LnGamma(double z)                 148 double MyGamma::LnGamma(double z)
139 {                                                 149 {
140   if (z <= 0)                                  << 150   // Computation of ln[gamma(z)] for all z>0.
141       return 0;                                << 151   //
142                                                << 152   // The algorithm is based on the article by C.Lanczos [1] as denoted in
143   return std::lgamma(z);                       << 153   // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
                                                   >> 154   //
                                                   >> 155   // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
                                                   >> 156   //
                                                   >> 157   // The accuracy of the result is better than 2e-10.
                                                   >> 158   //
                                                   >> 159   //--- Nve 14-nov-1998 UU-SAP Utrecht
                                                   >> 160   
                                                   >> 161   if (z<=0) return 0;
                                                   >> 162   
                                                   >> 163   // Coefficients for the series expansion
                                                   >> 164   double c[7] = { 2.5066282746310005, 76.18009172947146, -86.50532032941677
                                                   >> 165     ,24.01409824083091,  -1.231739572450155, 0.1208650973866179e-2
                                                   >> 166     ,-0.5395239384953e-5};
                                                   >> 167   
                                                   >> 168   double x   = z;
                                                   >> 169   double y   = x;
                                                   >> 170   double tmp = x+5.5;
                                                   >> 171   tmp = (x+0.5)*Log(tmp)-tmp;
                                                   >> 172   double ser = 1.000000000190015;
                                                   >> 173   for (int i=1; i<7; i++) {
                                                   >> 174     y   += 1;
                                                   >> 175     ser += c[i]/y;
                                                   >> 176   }
                                                   >> 177   double v = tmp+Log(c[0]*ser/x);
                                                   >> 178   return v;
144 }                                                 179 }
145                                                   180