Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/util/src/G4LegendrePolynomial.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 /processes/hadronic/util/src/G4LegendrePolynomial.cc (Version 11.3.0) and /processes/hadronic/util/src/G4LegendrePolynomial.cc (Version 11.0)


  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 #include "G4ios.hh"                                26 #include "G4ios.hh"
 27 #include "G4LegendrePolynomial.hh"                 27 #include "G4LegendrePolynomial.hh"
 28 #include "G4Pow.hh"                                28 #include "G4Pow.hh"
 29 #include "G4Exp.hh"                                29 #include "G4Exp.hh"
 30 #include "G4Log.hh"                                30 #include "G4Log.hh"
 31                                                    31 
 32 using namespace std;                               32 using namespace std;
 33                                                    33 
 34 G4double G4LegendrePolynomial::GetCoefficient( <<  34 G4double G4LegendrePolynomial::GetCoefficient(size_t i, size_t order)
 35 {                                                  35 {
 36   if(order >= fCoefficients.size()) BuildUpToO     36   if(order >= fCoefficients.size()) BuildUpToOrder(order);
 37   if(order >= fCoefficients.size() ||              37   if(order >= fCoefficients.size() ||
 38      i/2 >= fCoefficients[order].size() ||         38      i/2 >= fCoefficients[order].size() ||
 39      (i%2) != order %2) return 0;                  39      (i%2) != order %2) return 0;
 40   return fCoefficients[order][i/2];                40   return fCoefficients[order][i/2];
 41 }                                                  41 }
 42                                                    42 
 43 G4double G4LegendrePolynomial::EvalLegendrePol     43 G4double G4LegendrePolynomial::EvalLegendrePoly(G4int order, G4double x)
 44 {                                                  44 {
 45   // Call EvalAssocLegendrePoly with m=0           45   // Call EvalAssocLegendrePoly with m=0
 46   return (EvalAssocLegendrePoly(order,0,x));       46   return (EvalAssocLegendrePoly(order,0,x));
 47 }                                                  47 }
 48                                                    48 
 49 G4double G4LegendrePolynomial::EvalAssocLegend <<  49 G4double G4LegendrePolynomial::EvalAssocLegendrePoly(G4int l, G4int m, G4double x, 
                                                   >>  50                                                      map<G4int, map<G4int, G4double> >* cache)
 50 {                                                  51 {
 51   // Invalid calls --> 0. (Keeping for backwar <<  52   // Calculate P_l^m(x).
 52   if (l<0 || m<-l || m>l) return 0.0;          <<  53   // If cache ptr is non-null, use cache[l][m] if it exists, otherwise compute
                                                   >>  54   // P_l^m(x) and cache it in that position. The cache speeds up calculations
                                                   >>  55   // where many P_l^m computations are need at the same value of x.
                                                   >>  56 
                                                   >>  57   if(l<0 || m<-l || m>l) return 0;
                                                   >>  58   G4Pow* g4pow =  G4Pow::GetInstance();
                                                   >>  59 
                                                   >>  60   // Use non-log factorial for low l, m: it is more efficient until 
                                                   >>  61   // l and m get above 10 or so.
                                                   >>  62   // FIXME: G4Pow doesn't check whether the argument gets too large, 
                                                   >>  63   // which is unsafe! Max is 512; VI: It is assume that Geant4 does not
                                                   >>  64   // need higher order
                                                   >>  65   if(m<0) {
                                                   >>  66     G4double value = (m%2 ? -1. : 1.) * EvalAssocLegendrePoly(l, -m, x);
                                                   >>  67     if(l < 10) return value * g4pow->factorial(l+m)/g4pow->factorial(l-m);
                                                   >>  68     else { return value * G4Exp(g4pow->logfactorial(l+m) - g4pow->logfactorial(l-m));
                                                   >>  69     }
                                                   >>  70   }
 53                                                    71 
 54   // New check: g4pow doesn't handle integer a <<  72   // hard-code the first few orders for speed
 55   // For us that means:                        <<  73   if(l==0)   return 1;
 56   if ((l+m) > 512 || (l-m) > 512 || (2*m) > 51 <<  74   if(l==1) {
                                                   >>  75     if(m==0){return x;}
                                                   >>  76     /*m==1*/ return -sqrt(1.-x*x);
                                                   >>  77   }
                                                   >>  78   if(l<5) {
                                                   >>  79     G4double x2 = x*x;
                                                   >>  80     if(l==2) {
                                                   >>  81       if(m==0){return 0.5*(3.*x2 - 1.);}
                                                   >>  82       if(m==1){return -3.*x*sqrt(1.-x2);}
                                                   >>  83       /*m==2*/ return 3.*(1.-x2);
                                                   >>  84     }
                                                   >>  85     if(l==3) {
                                                   >>  86       if(m==0){return 0.5*(5.*x*x2 - 3.*x);}
                                                   >>  87       if(m==1){return -1.5*(5.*x2-1.)*sqrt(1.-x2);}
                                                   >>  88       if(m==2){return 15.*x*(1.-x2);}
                                                   >>  89       /*m==3*/ return -15.*(1.-x2)*sqrt(1.-x2);
                                                   >>  90     }
                                                   >>  91     if(l==4) {
                                                   >>  92       if(m==0){return 0.125*(35.*x2*x2 - 30.*x2 + 3.);}
                                                   >>  93       if(m==1){return -2.5*(7.*x*x2-3.*x)*sqrt(1.-x2);}
                                                   >>  94       if(m==2){return 7.5*(7.*x2-1.)*(1.-x2);}
                                                   >>  95       if(m==3){return -105.*x*(1.-x2)*sqrt(1.-x2);}
                                                   >>  96       /*m==4*/ return 105.*(1. - 2.*x2 + x2*x2);
                                                   >>  97     }
                                                   >>  98   }
 57                                                    99 
 58   G4Pow* g4pow = G4Pow::GetInstance();         << 100   // Easy special cases
 59   G4double x2 = x*x;                           << 101   // FIXME: G4Pow doesn't check whether the argument gets too large, which is unsafe! Max is 512. 
                                                   >> 102   if(m==l) return (l%2 ? -1. : 1.) * 
                                                   >> 103        G4Exp(g4pow->logfactorial(2*l) - g4pow->logfactorial(l)) * 
                                                   >> 104        G4Exp(G4Log((1.-x*x)*0.25)*0.5*G4double(l));
                                                   >> 105   if(m==l-1) return x*(2.*G4double(m)+1.)*EvalAssocLegendrePoly(m,m,x);
                                                   >> 106 
                                                   >> 107   // See if we have this value cached.
                                                   >> 108   if(cache != NULL && cache->count(l) > 0 && (*cache)[l].count(m) > 0) {
                                                   >> 109     return (*cache)[l][m];
                                                   >> 110   }
 60                                                   111 
 61   // hard-code the first few orders for speed  << 112   // Otherwise calculate recursively
 62   switch (l) {                                 << 113   G4double value = (x*G4double(2*l-1)*EvalAssocLegendrePoly(l-1,m,x) - 
 63     case 0 :                                   << 114              (G4double(l+m-1))*EvalAssocLegendrePoly(l-2,m,x))/G4double(l-m);
 64       return 1;                                << 115 
 65     case 1 :                                   << 116   // If we are working with a cache, cache this value.
 66       switch (m) {                             << 117   if(cache != NULL) {
 67         case -1 : return 0.5 * sqrt(1.-x2);    << 118     (*cache)[l][m] = value;
 68         case 0 : return x;                     << 
 69         case 1 : return -sqrt(1.-x2);          << 
 70       };                                       << 
 71       break;                                   << 
 72     case 2 :                                   << 
 73       switch (m) {                             << 
 74         case -2 : return 0.125 * (1.0 - x2);   << 
 75         case -1 : return 0.5 * x * sqrt(1.0 -  << 
 76         case 0 : return 0.5*(3.*x2 - 1.);      << 
 77         case 1 : return -3.*x*sqrt(1.-x2);     << 
 78         case 2 : return 3.*(1.-x2);            << 
 79       };                                       << 
 80       break;                                   << 
 81     case 3 :                                   << 
 82       switch (m) {                             << 
 83         case -3 : return (1.0/48.0) * (1.0 - x << 
 84         case -2 : return 0.125 * x * (1.0 - x2 << 
 85         case -1 : return 0.125 * (5.0 * x2 - 1 << 
 86         case 0 : return 0.5*(5.*x*x2 - 3.*x);  << 
 87         case 1 : return -1.5*(5.*x2-1.)*sqrt(1 << 
 88         case 2 : return 15.*x*(1.-x2);         << 
 89         case 3 : return -15.*(1.-x2)*sqrt(1.-x << 
 90       };                                       << 
 91       break;                                   << 
 92     case 4 :                                   << 
 93       switch (m) {                             << 
 94         case -4 : return (105.0/40320.0)*(1. - << 
 95         case -3 : return (105.0/5040.0)*x*(1.- << 
 96         case -2 : return (15.0/720.0)*(7.*x2-1 << 
 97         case -1 : return 0.125*(7.*x*x2-3.*x)* << 
 98         case 0 : return 0.125*(35.*x2*x2 - 30. << 
 99         case 1 : return -2.5*(7.*x*x2-3.*x)*sq << 
100         case 2 : return 7.5*(7.*x2-1.)*(1.-x2) << 
101         case 3 : return -105.*x*(1.-x2)*sqrt(1 << 
102         case 4 : return 105.*(1. - 2.*x2 + x2* << 
103       };                                       << 
104       break;                                   << 
105   };                                           << 
106                                                << 
107   // if m<0, compute P[l,-m,x] and correct     << 
108   if (m < 0)                                   << 
109   {                                            << 
110     G4double complementary_value = EvalAssocLe << 
111     return complementary_value * (m%2==0 ? 1.0 << 
112   }                                            << 
113                                                << 
114   // Iteratively walk up from P[m,m,x] to P[l, << 
115                                                << 
116   // prime the pump: P[l<m,m,x] = 0            << 
117   G4double previous = 0.0;                     << 
118                                                << 
119   // prime the pump: P[m,m,x]                  << 
120   G4double current;                            << 
121   if (m == 0) current = 1.0;                   << 
122   else if (m == 1) current = -sqrt((1.0 - (x2) << 
123   else {                                       << 
124     current = (m%2==0 ? 1.0 : -1.0) *          << 
125        G4Exp(g4pow->logfactorial(2*m) - g4pow- << 
126        G4Exp(G4Log((1.0-(x2))*0.25)*0.5*G4doub << 
127   }                                            << 
128                                                << 
129   // Work up to P[l,m,x]                       << 
130   for(G4int i=m+1; i<=l; i++)                  << 
131   {                                            << 
132     G4double next = (-(G4double(i+m-1))*previo << 
133     previous = current;                        << 
134     current = next;                            << 
135   }                                               119   }
136                                                << 120   return value;
137   return current;                              << 
138 }                                                 121 }
139                                                   122 
140 void G4LegendrePolynomial::BuildUpToOrder(std: << 123 void G4LegendrePolynomial::BuildUpToOrder(size_t orderMax)
141 {                                                 124 {
142   if(orderMax > 30) {                             125   if(orderMax > 30) {
143     G4cout << "G4LegendrePolynomial::GetCoeffi    126     G4cout << "G4LegendrePolynomial::GetCoefficient(): "
144            << "I refuse to make a Legendre Pol    127            << "I refuse to make a Legendre Polynomial of order " 
145            << orderMax << G4endl;                 128            << orderMax << G4endl;
146     return;                                       129     return;
147   }                                               130   }
148   while(fCoefficients.size() < orderMax+1) {      131   while(fCoefficients.size() < orderMax+1) {  /* Loop checking, 30-Oct-2015, G.Folger */
149     std::size_t order = fCoefficients.size();  << 132     size_t order = fCoefficients.size();
150     fCoefficients.resize(order+1);                133     fCoefficients.resize(order+1);
151     if(order <= 1) fCoefficients[order].push_b    134     if(order <= 1) fCoefficients[order].push_back(1.);
152     else {                                        135     else {
153       for(std::size_t iCoeff = 0; iCoeff < ord << 136       for(size_t iCoeff = 0; iCoeff < order+1; ++iCoeff) {
154         if((order % 2) == (iCoeff % 2)) {         137         if((order % 2) == (iCoeff % 2)) {
155           G4double coeff = 0;                     138           G4double coeff = 0;
156           if(iCoeff <= order-2) coeff -= fCoef    139           if(iCoeff <= order-2) coeff -= fCoefficients[order-2][iCoeff/2]*G4double(order-1);
157           if(iCoeff > 0) coeff += fCoefficient    140           if(iCoeff > 0) coeff += fCoefficients[order-1][(iCoeff-1)/2]*G4double(2*order-1);
158           coeff /= G4double(order);               141           coeff /= G4double(order);
159           fCoefficients[order].push_back(coeff    142           fCoefficients[order].push_back(coeff);
160         }                                         143         }
161       }                                           144       }
162     }                                             145     }
163   }                                               146   }
164 }                                                 147 }
165                                                   148 
166                                                   149