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