Geant4 Cross Reference |
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 // G4PolynomialSolver inline methods implement << 27 // 23 // 28 // Author: E.Medernach, 19.12.2000 - First imp << 24 // $Id: G4PolynomialSolver.icc,v 1.7 2004/11/12 17:38:29 gcosmo Exp $ 29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-07-01 $ >> 26 // >> 27 // class G4PolynomialSolver >> 28 // >> 29 // 19.12.00 E.Medernach, First implementation >> 30 // 30 31 31 #define POLEPSILON 1e-12 << 32 #define POLEPSILON 1e-12 32 #define POLINFINITY 9.0E99 << 33 #define POLINFINITY 9.0E99 33 #define ITERATION 12 // 20 But 8 is really en << 34 #define ITERATION 12 // 20 But 8 is really enough for Newton with a good guess 34 35 35 template <class T, class F> 36 template <class T, class F> 36 G4PolynomialSolver<T, F>::G4PolynomialSolver(T << 37 G4PolynomialSolver<T,F>::G4PolynomialSolver (T* typeF, F func, F deriv, 37 G 38 G4double precision) 38 { 39 { 39 Precision = precision; << 40 Precision = precision ; 40 FunctionClass = typeF; << 41 FunctionClass = typeF ; 41 Function = func; << 42 Function = func ; 42 Derivative = deriv; << 43 Derivative = deriv ; 43 } 44 } 44 45 45 template <class T, class F> 46 template <class T, class F> 46 G4PolynomialSolver<T, F>::~G4PolynomialSolver( << 47 G4PolynomialSolver<T,F>::~G4PolynomialSolver () 47 {} << 48 { >> 49 } 48 50 49 template <class T, class F> 51 template <class T, class F> 50 G4double G4PolynomialSolver<T, F>::solve(G4dou << 52 G4double G4PolynomialSolver<T,F>::solve(G4double IntervalMin, 51 G4dou << 53 G4double IntervalMax) 52 { 54 { 53 return Newton(IntervalMin, IntervalMax); << 55 return Newton(IntervalMin,IntervalMax); 54 } 56 } 55 57 >> 58 56 /* If we want to be general this could work fo 59 /* If we want to be general this could work for any 57 polynomial of order more that 4 if we find 60 polynomial of order more that 4 if we find the (ORDER + 1) 58 control points 61 control points 59 */ 62 */ 60 #define NBBEZIER 5 63 #define NBBEZIER 5 61 64 62 template <class T, class F> 65 template <class T, class F> 63 G4int G4PolynomialSolver<T, F>::BezierClipping << 66 G4int 64 << 67 G4PolynomialSolver<T,F>::BezierClipping(/*T* typeF,F func,F deriv,*/ 65 << 68 G4double *IntervalMin, >> 69 G4double *IntervalMax) 66 { 70 { 67 /** BezierClipping is a clipping interval Ne 71 /** BezierClipping is a clipping interval Newton method **/ 68 /** It works by clipping the area where the 72 /** It works by clipping the area where the polynomial is **/ 69 73 70 G4double P[NBBEZIER][2], D[2]; << 74 G4double P[NBBEZIER][2],D[2]; 71 G4double NewMin, NewMax; << 75 G4double NewMin,NewMax; 72 76 73 G4int IntervalIsVoid = 1; 77 G4int IntervalIsVoid = 1; 74 << 78 75 /*** Calculating Control Points ***/ 79 /*** Calculating Control Points ***/ 76 /* We see the polynomial as a Bezier curve f 80 /* We see the polynomial as a Bezier curve for some control points to find */ 77 81 78 /* 82 /* 79 For 5 control points (polynomial of degree 83 For 5 control points (polynomial of degree 4) this is: 80 << 84 81 0 p0 = F((*IntervalMin)) 85 0 p0 = F((*IntervalMin)) 82 1/4 p1 = F((*IntervalMin)) + ((*Interval 86 1/4 p1 = F((*IntervalMin)) + ((*IntervalMax) - (*IntervalMin))/4 83 * F'((*IntervalMin)) 87 * F'((*IntervalMin)) 84 2/4 p2 = 1/6 * (16*F(((*IntervalMax) + ( 88 2/4 p2 = 1/6 * (16*F(((*IntervalMax) + (*IntervalMin))/2) 85 - (p0 + 4*p1 + 4*p3 + p4 << 89 - (p0 + 4*p1 + 4*p3 + p4)) 86 3/4 p3 = F((*IntervalMax)) - ((*Interval 90 3/4 p3 = F((*IntervalMax)) - ((*IntervalMax) - (*IntervalMin))/4 87 * F'((*IntervalMax)) 91 * F'((*IntervalMax)) 88 1 p4 = F((*IntervalMax)) 92 1 p4 = F((*IntervalMax)) 89 */ << 93 */ 90 94 91 /* x,y,z,dx,dy,dz are constant during search 95 /* x,y,z,dx,dy,dz are constant during searching */ 92 96 93 D[0] = (FunctionClass->*Derivative)(*Interva 97 D[0] = (FunctionClass->*Derivative)(*IntervalMin); 94 << 98 95 P[0][0] = (*IntervalMin); 99 P[0][0] = (*IntervalMin); 96 P[0][1] = (FunctionClass->*Function)(*Interv 100 P[0][1] = (FunctionClass->*Function)(*IntervalMin); >> 101 97 102 98 if(std::fabs(P[0][1]) < Precision) << 103 if (std::fabs(P[0][1]) < Precision) { 99 { << 100 return 1; 104 return 1; 101 } 105 } 102 << 106 103 if(((*IntervalMax) - (*IntervalMin)) < POLEP << 107 if (((*IntervalMax) - (*IntervalMin)) < POLEPSILON) { 104 { << 105 return 1; 108 return 1; 106 } 109 } 107 110 108 P[1][0] = (*IntervalMin) + ((*IntervalMax) - << 111 P[1][0] = (*IntervalMin) + ((*IntervalMax) - (*IntervalMin))/4; 109 P[1][1] = P[0][1] + (((*IntervalMax) - (*Int << 112 P[1][1] = P[0][1] + (((*IntervalMax) - (*IntervalMin))/4.0) * D[0]; 110 113 111 D[1] = (FunctionClass->*Derivative)(*Interva 114 D[1] = (FunctionClass->*Derivative)(*IntervalMax); 112 115 113 P[4][0] = (*IntervalMax); 116 P[4][0] = (*IntervalMax); 114 P[4][1] = (FunctionClass->*Function)(*Interv 117 P[4][1] = (FunctionClass->*Function)(*IntervalMax); >> 118 >> 119 P[3][0] = (*IntervalMax) - ((*IntervalMax) - (*IntervalMin))/4; >> 120 P[3][1] = P[4][1] - ((*IntervalMax) - (*IntervalMin))/4 * D[1]; >> 121 >> 122 P[2][0] = ((*IntervalMax) + (*IntervalMin))/2; >> 123 P[2][1] = (16*(FunctionClass->*Function)(((*IntervalMax)+(*IntervalMin))/2) >> 124 - (P[0][1] + 4*P[1][1] + 4*P[3][1] + P[4][1]))/6 ; >> 125 >> 126 { >> 127 G4double Intersection ; >> 128 G4int i,j; >> 129 >> 130 NewMin = (*IntervalMax) ; >> 131 NewMax = (*IntervalMin) ; >> 132 >> 133 for (i=0;i<5;i++) >> 134 for (j=i+1;j<5;j++) >> 135 { >> 136 /* there is an intersection only if each have different signs */ >> 137 if (((P[j][1] > -Precision) && (P[i][1] < Precision)) || >> 138 ((P[j][1] < Precision) && (P[i][1] > -Precision))) { >> 139 IntervalIsVoid = 0; >> 140 Intersection = P[j][0] - P[j][1]*((P[i][0] - P[j][0])/ >> 141 (P[i][1] - P[j][1])); >> 142 if (Intersection < NewMin) { >> 143 NewMin = Intersection; >> 144 } >> 145 if (Intersection > NewMax) { >> 146 NewMax = Intersection; >> 147 } >> 148 } >> 149 } 115 150 116 P[3][0] = (*IntervalMax) - ((*IntervalMax) - << 151 117 P[3][1] = P[4][1] - ((*IntervalMax) - (*Inte << 152 if (IntervalIsVoid != 1) { 118 << 119 P[2][0] = ((*IntervalMax) + (*IntervalMin)) << 120 P[2][1] = << 121 (16 * (FunctionClass->*Function)(((*Interv << 122 (P[0][1] + 4 * P[1][1] + 4 * P[3][1] + P[ << 123 6; << 124 << 125 { << 126 G4double Intersection; << 127 G4int i, j; << 128 << 129 NewMin = (*IntervalMax); << 130 NewMax = (*IntervalMin); << 131 << 132 for(i = 0; i < 5; ++i) << 133 for(j = i + 1; j < 5; ++j) << 134 { << 135 /* there is an intersection only if ea << 136 if(((P[j][1] > -Precision) && (P[i][1] << 137 ((P[j][1] < Precision) && (P[i][1] << 138 { << 139 IntervalIsVoid = 0; << 140 Intersection = << 141 P[j][0] - P[j][1] * ((P[i][0] - P[ << 142 if(Intersection < NewMin) << 143 { << 144 NewMin = Intersection; << 145 } << 146 if(Intersection > NewMax) << 147 { << 148 NewMax = Intersection; << 149 } << 150 } << 151 } << 152 << 153 if(IntervalIsVoid != 1) << 154 { << 155 (*IntervalMax) = NewMax; 153 (*IntervalMax) = NewMax; 156 (*IntervalMin) = NewMin; 154 (*IntervalMin) = NewMin; 157 } 155 } 158 } 156 } 159 << 157 160 if(IntervalIsVoid == 1) << 158 if (IntervalIsVoid == 1) { 161 { << 162 return -1; 159 return -1; 163 } 160 } 164 << 161 165 return 0; 162 return 0; 166 } 163 } 167 164 168 template <class T, class F> 165 template <class T, class F> 169 G4double G4PolynomialSolver<T, F>::Newton(G4do << 166 G4double G4PolynomialSolver<T,F>::Newton (G4double IntervalMin, 170 G4do 167 G4double IntervalMax) 171 { 168 { 172 /* So now we have a good guess and an interv 169 /* So now we have a good guess and an interval where 173 if there are an intersection the root mus 170 if there are an intersection the root must be */ 174 171 175 G4double Value = 0; << 172 G4double Value = 0; 176 G4double Gradient = 0; 173 G4double Gradient = 0; 177 G4double Lambda; << 174 G4double Lambda ; 178 << 179 G4int i = 0; << 180 G4int j = 0; << 181 175 >> 176 G4int i=0; >> 177 G4int j=0; >> 178 >> 179 182 /* Reduce interval before applying Newton Me 180 /* Reduce interval before applying Newton Method */ 183 { 181 { 184 G4int NewtonIsSafe; << 182 G4int NewtonIsSafe ; 185 183 186 while((NewtonIsSafe = BezierClipping(&Inte << 184 while ((NewtonIsSafe = BezierClipping(&IntervalMin,&IntervalMax)) == 0) ; 187 ; << 188 185 189 if(NewtonIsSafe == -1) << 186 if (NewtonIsSafe == -1) { 190 { << 191 return POLINFINITY; 187 return POLINFINITY; 192 } 188 } 193 } 189 } 194 190 195 Lambda = IntervalMin; 191 Lambda = IntervalMin; 196 Value = (FunctionClass->*Function)(Lambda); << 192 Value = (FunctionClass->*Function)(Lambda); >> 193 197 194 198 // while ((std::fabs(Value) > Precision)) { 195 // while ((std::fabs(Value) > Precision)) { 199 while(j != -1) << 196 while (j != -1) { 200 { << 197 201 Value = (FunctionClass->*Function)(Lambda) 198 Value = (FunctionClass->*Function)(Lambda); 202 199 203 Gradient = (FunctionClass->*Derivative)(La 200 Gradient = (FunctionClass->*Derivative)(Lambda); 204 201 205 Lambda = Lambda - Value / Gradient; << 202 Lambda = Lambda - Value/Gradient ; 206 203 207 if(std::fabs(Value) <= Precision) << 204 if (std::fabs(Value) <= Precision) { 208 { << 205 j ++; 209 ++j; << 206 if (j == 2) { 210 if(j == 2) << 207 j = -1; 211 { << 208 } 212 j = -1; << 209 } else { 213 } << 210 i ++; 214 } << 211 215 else << 212 if (i > ITERATION) 216 { << 213 return POLINFINITY; 217 ++i; << 214 } 218 << 219 if(i > ITERATION) << 220 return POLINFINITY; << 221 } << 222 } 215 } 223 return Lambda; << 216 return Lambda ; 224 } 217 } 225 218