Geant4 Cross Reference |
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 #include "G4ios.hh" 27 #include "G4LegendrePolynomial.hh" 28 #include "G4Pow.hh" 29 #include "G4Exp.hh" 30 #include "G4Log.hh" 31 32 using namespace std; 33 34 G4double G4LegendrePolynomial::GetCoefficient(std::size_t i, std::size_t order) 35 { 36 if(order >= fCoefficients.size()) BuildUpToOrder(order); 37 if(order >= fCoefficients.size() || 38 i/2 >= fCoefficients[order].size() || 39 (i%2) != order %2) return 0; 40 return fCoefficients[order][i/2]; 41 } 42 43 G4double G4LegendrePolynomial::EvalLegendrePoly(G4int order, G4double x) 44 { 45 // Call EvalAssocLegendrePoly with m=0 46 return (EvalAssocLegendrePoly(order,0,x)); 47 } 48 49 G4double G4LegendrePolynomial::EvalAssocLegendrePoly(G4int l, G4int m, G4double x) 50 { 51 // Invalid calls --> 0. (Keeping for backward compatibility; should we throw instead?) 52 if (l<0 || m<-l || m>l) return 0.0; 53 54 // New check: g4pow doesn't handle integer arguments over 512. 55 // For us that means: 56 if ((l+m) > 512 || (l-m) > 512 || (2*m) > 512) return 0.0; 57 58 G4Pow* g4pow = G4Pow::GetInstance(); 59 G4double x2 = x*x; 60 61 // hard-code the first few orders for speed 62 switch (l) { 63 case 0 : 64 return 1; 65 case 1 : 66 switch (m) { 67 case -1 : return 0.5 * sqrt(1.-x2); 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 - x2); 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 - x2) * sqrt(1.0 - x2); 84 case -2 : return 0.125 * x * (1.0 - x2); 85 case -1 : return 0.125 * (5.0 * x2 - 1.0) * sqrt(1.0 - x2); 86 case 0 : return 0.5*(5.*x*x2 - 3.*x); 87 case 1 : return -1.5*(5.*x2-1.)*sqrt(1.-x2); 88 case 2 : return 15.*x*(1.-x2); 89 case 3 : return -15.*(1.-x2)*sqrt(1.-x2); 90 }; 91 break; 92 case 4 : 93 switch (m) { 94 case -4 : return (105.0/40320.0)*(1. - 2.*x2 + x2*x2); 95 case -3 : return (105.0/5040.0)*x*(1.-x2)*sqrt(1.-x2); 96 case -2 : return (15.0/720.0)*(7.*x2-1.)*(1.-x2); 97 case -1 : return 0.125*(7.*x*x2-3.*x)*sqrt(1.-x2); 98 case 0 : return 0.125*(35.*x2*x2 - 30.*x2 + 3.); 99 case 1 : return -2.5*(7.*x*x2-3.*x)*sqrt(1.-x2); 100 case 2 : return 7.5*(7.*x2-1.)*(1.-x2); 101 case 3 : return -105.*x*(1.-x2)*sqrt(1.-x2); 102 case 4 : return 105.*(1. - 2.*x2 + 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 = EvalAssocLegendrePoly(l, -m, x); 111 return complementary_value * (m%2==0 ? 1.0 : -1.0) * g4pow->factorial(l+m)/g4pow->factorial(l-m); 112 } 113 114 // Iteratively walk up from P[m,m,x] to P[l,m,x] 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->logfactorial(m)) * 126 G4Exp(G4Log((1.0-(x2))*0.25)*0.5*G4double(m)); 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))*previous + x*G4double(2*i-1)*current )/G4double(i-m); 133 previous = current; 134 current = next; 135 } 136 137 return current; 138 } 139 140 void G4LegendrePolynomial::BuildUpToOrder(std::size_t orderMax) 141 { 142 if(orderMax > 30) { 143 G4cout << "G4LegendrePolynomial::GetCoefficient(): " 144 << "I refuse to make a Legendre Polynomial of order " 145 << orderMax << G4endl; 146 return; 147 } 148 while(fCoefficients.size() < orderMax+1) { /* Loop checking, 30-Oct-2015, G.Folger */ 149 std::size_t order = fCoefficients.size(); 150 fCoefficients.resize(order+1); 151 if(order <= 1) fCoefficients[order].push_back(1.); 152 else { 153 for(std::size_t iCoeff = 0; iCoeff < order+1; ++iCoeff) { 154 if((order % 2) == (iCoeff % 2)) { 155 G4double coeff = 0; 156 if(iCoeff <= order-2) coeff -= fCoefficients[order-2][iCoeff/2]*G4double(order-1); 157 if(iCoeff > 0) coeff += fCoefficients[order-1][(iCoeff-1)/2]*G4double(2*order-1); 158 coeff /= G4double(order); 159 fCoefficients[order].push_back(coeff); 160 } 161 } 162 } 163 } 164 } 165 166