Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/src/G4GaussJacobiQ.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/G4GaussJacobiQ.cc (Version 11.3.0) and /global/HEPNumerics/src/G4GaussJacobiQ.cc (Version 4.1.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 // G4GaussJacobiQ class implementation         << 
 27 //                                                 23 //
 28 // Author: V.Grichine, 13.05.1997              <<  24 // $Id: G4GaussJacobiQ.cc,v 1.4 2001/07/11 10:00:41 gunter Exp $
 29 // ------------------------------------------- <<  25 // GEANT4 tag $Name: geant4-04-01 $
 30                                                <<  26 //
 31 #include "G4GaussJacobiQ.hh"                       27 #include "G4GaussJacobiQ.hh"
 32                                                    28 
                                                   >>  29 
 33 // -------------------------------------------     30 // -------------------------------------------------------------
 34 //                                                 31 //
 35 // Constructor for Gauss-Jacobi integration me <<  32 // Constructor for Gauss-Jacobi integration method. 
 36 //                                                 33 //
 37                                                    34 
 38 G4GaussJacobiQ::G4GaussJacobiQ(function pFunct <<  35 G4GaussJacobiQ::G4GaussJacobiQ(       function pFunction,
 39                                G4double beta,  <<  36                     G4double alpha,
 40   : G4VGaussianQuadrature(pFunction)           <<  37                                       G4double beta, 
                                                   >>  38                     G4int nJacobi           ) 
                                                   >>  39    : G4VGaussianQuadrature(pFunction)
 41                                                    40 
 42 {                                                  41 {
 43   const G4double tolerance = 1.0e-12;          <<  42   const G4double tolerance = 1.0e-12 ;
 44   const G4double maxNumber = 12;               <<  43   const G4double maxNumber = 12 ;
 45   G4int i = 1, k = 1;                          <<  44   G4int i, k, j ;
 46   G4double root      = 0.;                     <<  45   G4double root=0.;
 47   G4double alphaBeta = 0.0, alphaReduced = 0.0 <<  46   G4double alphaBeta, alphaReduced, betaReduced, root1, root2, root3 ;
 48            root2 = 0.0, root3 = 0.0;           <<  47   G4double a, b, c, newton1, newton2, newton3, newton, temp, rootTemp ;
 49   G4double a = 0.0, b = 0.0, c = 0.0, newton1  <<  48 
 50            newton3 = 0.0, newton0 = 0.0, temp  <<  49   fNumber   = nJacobi ;
 51                                                <<  50   fAbscissa = new G4double[fNumber] ;
 52   fNumber   = nJacobi;                         <<  51   fWeight   = new G4double[fNumber] ;
 53   fAbscissa = new G4double[fNumber];           << 
 54   fWeight   = new G4double[fNumber];           << 
 55                                                    52 
 56   for(i = 1; i <= nJacobi; ++i)                <<  53   for (i=1;i<=nJacobi;i++)
 57   {                                                54   {
 58     if(i == 1)                                 <<  55      if (i == 1)
 59     {                                          <<  56      {
 60       alphaReduced = alpha / nJacobi;          <<  57   alphaReduced = alpha/nJacobi ;
 61       betaReduced  = beta / nJacobi;           <<  58   betaReduced = beta/nJacobi ;
 62       root1        = (1.0 + alpha) * (2.78002  <<  59   root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
 63                                0.767999 * alph <<  60         0.767999*alphaReduced/nJacobi) ;
 64       root2        = 1.0 + 1.48 * alphaReduced <<  61   root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
 65               0.451998 * alphaReduced * alphaR <<  62        0.451998*alphaReduced*alphaReduced+0.83001*alphaReduced*betaReduced ;
 66               0.83001 * alphaReduced * betaRed <<  63   root  = 1.0-root1/root2 ;
 67       root = 1.0 - root1 / root2;              <<  64      } 
 68     }                                          <<  65      else if (i == 2)
 69     else if(i == 2)                            <<  66      {
 70     {                                          <<  67   root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
 71       root1 = (4.1002 + alpha) / ((1.0 + alpha <<  68   root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
 72       root2 = 1.0 + 0.06 * (nJacobi - 8.0) * ( <<  69   root3=1.0+0.012002*beta*(1.0+0.24997*fabs(alpha))/nJacobi ;
 73       root3 =                                  <<  70   root -= (1.0-root)*root1*root2*root3 ;
 74         1.0 + 0.012002 * beta * (1.0 + 0.24997 <<  71      } 
 75       root -= (1.0 - root) * root1 * root2 * r <<  72      else if (i == 3) 
 76     }                                          <<  73      {
 77     else if(i == 3)                            <<  74   root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
 78     {                                          <<  75   root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
 79       root1 = (1.67001 + 0.27998 * alpha) / (1 <<  76   root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
 80       root2 = 1.0 + 0.22 * (nJacobi - 8.0) / n <<  77   root -= (fAbscissa[0]-root)*root1*root2*root3 ;
 81       root3 = 1.0 + 8.0 * beta / ((6.28001 + b <<  78      }
 82       root -= (fAbscissa[0] - root) * root1 *  <<  79      else if (i == nJacobi-1)
 83     }                                          <<  80      {
 84     else if(i == nJacobi - 1)                  <<  81   root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
 85     {                                          <<  82   root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
 86       root1 = (1.0 + 0.235002 * beta) / (0.766 <<  83   root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
 87       root2 = 1.0 / (1.0 + 0.639002 * (nJacobi <<  84   root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
 88                              (1.0 + 0.71001 *  <<  85      } 
 89       root3 = 1.0 / (1.0 + 20.0 * alpha / ((7. <<  86      else if (i == nJacobi) 
 90       root += (root - fAbscissa[nJacobi - 4])  <<  87      {
 91     }                                          <<  88   root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
 92     else if(i == nJacobi)                      <<  89   root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
 93     {                                          <<  90   root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
 94       root1 = (1.0 + 0.37002 * beta) / (1.6700 <<  91   root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
 95       root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8 <<  92      } 
 96       root3 =                                  <<  93      else
 97         1.0 / (1.0 + 8.0 * alpha / ((6.28002 + <<  94      {
 98       root += (root - fAbscissa[nJacobi - 3])  <<  95   root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
 99     }                                          <<  96      }
100     else                                       <<  97      alphaBeta = alpha + beta ;
101     {                                          <<  98      for (k=1;k<=maxNumber;k++)
102       root = 3.0 * fAbscissa[i - 2] - 3.0 * fA <<  99      {
103     }                                          << 100   temp = 2.0 + alphaBeta ;
104     alphaBeta = alpha + beta;                  << 101   newton1 = (alpha-beta+temp*root)/2.0 ;
105     for(k = 1; k <= maxNumber; ++k)            << 102   newton2 = 1.0 ;
106     {                                          << 103   for (j=2;j<=nJacobi;j++)
107       temp    = 2.0 + alphaBeta;               << 104   {
108       newton1 = (alpha - beta + temp * root) / << 105      newton3 = newton2 ;
109       newton2 = 1.0;                           << 106      newton2 = newton1 ;
110       for(G4int j = 2; j <= nJacobi; ++j)      << 107      temp = 2*j+alphaBeta ;
111       {                                        << 108      a = 2*j*(j+alphaBeta)*(temp-2.0) ;
112         newton3 = newton2;                     << 109      b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
113         newton2 = newton1;                     << 110      c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
114         temp    = 2 * j + alphaBeta;           << 111      newton1 = (b*newton2-c*newton3)/a ;
115         a       = 2 * j * (j + alphaBeta) * (t << 112   }
116         b       = (temp - 1.0) *               << 113   newton = (nJacobi*(alpha - beta - temp*root)*newton1 +
117             (alpha * alpha - beta * beta + tem << 114         2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/
118         c       = 2.0 * (j - 1 + alpha) * (j - << 115        (temp*(1.0 - root*root)) ;
119         newton1 = (b * newton2 - c * newton3)  << 116   rootTemp = root ;
120       }                                        << 117   root = rootTemp - newton1/newton ;
121       newton0 = (nJacobi * (alpha - beta - tem << 118   if (fabs(root-rootTemp) <= tolerance)
122                  2.0 * (nJacobi + alpha) * (nJ << 119   {
123                 (temp * (1.0 - root * root));  << 120      break ;
124       rootTemp = root;                         << 121   }
125       root     = rootTemp - newton1 / newton0; << 122      }
126       if(std::fabs(root - rootTemp) <= toleran << 123      if (k > maxNumber) 
127       {                                        << 124      {
128         break;                                 << 125         G4Exception("Too many iterations in G4GaussJacobiQ::G4GaussJacobiQ") ;
129       }                                        << 126      }
130     }                                          << 127      fAbscissa[i-1] = root ;
131     if(k > maxNumber)                          << 128      fWeight[i-1] = exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 
132     {                                          << 129             GammaLogarithm((G4double)(beta+nJacobi)) - 
133       G4Exception("G4GaussJacobiQ::G4GaussJaco << 130             GammaLogarithm((G4double)(nJacobi+1.0)) -
134                   FatalException, "Too many it << 131             GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
135     }                                          << 132             *temp*pow(2.0,alphaBeta)/(newton*newton2)                          ;
136     fAbscissa[i - 1] = root;                   << 
137     fWeight[i - 1] =                           << 
138       std::exp(GammaLogarithm((G4double)(alpha << 
139                GammaLogarithm((G4double)(beta  << 
140                GammaLogarithm((G4double)(nJaco << 
141                GammaLogarithm((G4double)(nJaco << 
142       temp * std::pow(2.0, alphaBeta) / (newto << 
143   }                                               133   }
144 }                                                 134 }
145                                                   135 
                                                   >> 136 
146 // -------------------------------------------    137 // ----------------------------------------------------------
147 //                                                138 //
148 // Gauss-Jacobi method for integration of      << 139 // Gauss-Jacobi method for integration of ((1-x)^alpha)*((1+x)^beta)*pFunction(x)
149 // ((1-x)^alpha)*((1+x)^beta)*pFunction(x)     << 
150 // from minus unit to plus unit .                 140 // from minus unit to plus unit .
151                                                   141 
152 G4double G4GaussJacobiQ::Integral() const      << 142 
                                                   >> 143 G4double 
                                                   >> 144 G4GaussJacobiQ::Integral() const 
153 {                                                 145 {
154   G4double integral = 0.0;                     << 146    G4int i ;
155   for(G4int i = 0; i < fNumber; ++i)           << 147    G4double integral = 0.0 ;
156   {                                            << 148    for(i=0;i<fNumber;i++)
157     integral += fWeight[i] * fFunction(fAbscis << 149    {
158   }                                            << 150       integral += fWeight[i]*fFunction(fAbscissa[i]) ;
159   return integral;                             << 151    }
                                                   >> 152    return integral ;
160 }                                                 153 }
                                                   >> 154 
161                                                   155