include/Geant4/G4PolynomialSolver.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4PolynomialSolver inline methods implementation
0027 //
0028 // Author: E.Medernach, 19.12.2000 - First implementation
0029 // --------------------------------------------------------------------
0030
0031 #define POLEPSILON 1e-12
0032 #define POLINFINITY 9.0E99
0033 #define ITERATION 12  // 20 But 8 is really enough for Newton with a good guess
0034
0035 template <class T, class F>
0036 G4PolynomialSolver<T, F>::G4PolynomialSolver(T* typeF, F func, F deriv,
0037                                              G4double precision)
0038 {
0039   Precision     = precision;
0040   FunctionClass = typeF;
0041   Function      = func;
0042   Derivative    = deriv;
0043 }
0044
0045 template <class T, class F>
0046 G4PolynomialSolver<T, F>::~G4PolynomialSolver()
0047 {}
0048
0049 template <class T, class F>
0050 G4double G4PolynomialSolver<T, F>::solve(G4double IntervalMin,
0051                                          G4double IntervalMax)
0052 {
0053   return Newton(IntervalMin, IntervalMax);
0054 }
0055
0056 /* If we want to be general this could work for any
0057    polynomial of order more that 4 if we find the (ORDER + 1)
0058    control points
0059 */
0060 #define NBBEZIER 5
0061
0062 template <class T, class F>
0063 G4int G4PolynomialSolver<T, F>::BezierClipping(/*T* typeF,F func,F deriv,*/
0064                                                G4double* IntervalMin,
0065                                                G4double* IntervalMax)
0066 {
0067   /** BezierClipping is a clipping interval Newton method **/
0068   /** It works by clipping the area where the polynomial is **/
0069
0070   G4double P[NBBEZIER][2], D[2];
0071   G4double NewMin, NewMax;
0072
0073   G4int IntervalIsVoid = 1;
0074
0075   /*** Calculating Control Points  ***/
0076   /* We see the polynomial as a Bezier curve for some control points to find */
0077
0078   /*
0079     For 5 control points (polynomial of degree 4) this is:
0080
0081     0     p0 = F((*IntervalMin))
0082     1/4   p1 = F((*IntervalMin)) + ((*IntervalMax) - (*IntervalMin))/4
0083                  * F'((*IntervalMin))
0084     2/4   p2 = 1/6 * (16*F(((*IntervalMax) + (*IntervalMin))/2)
0085                       - (p0 + 4*p1 + 4*p3 + p4))
0086     3/4   p3 = F((*IntervalMax)) - ((*IntervalMax) - (*IntervalMin))/4
0087                  * F'((*IntervalMax))
0088     1     p4 = F((*IntervalMax))
0089   */
0090
0091   /* x,y,z,dx,dy,dz are constant during searching */
0092
0093   D[0] = (FunctionClass->*Derivative)(*IntervalMin);
0094
0095   P[0][0] = (*IntervalMin);
0096   P[0][1] = (FunctionClass->*Function)(*IntervalMin);
0097
0098   if(std::fabs(P[0][1]) < Precision)
0099   {
0100     return 1;
0101   }
0102
0103   if(((*IntervalMax) - (*IntervalMin)) < POLEPSILON)
0104   {
0105     return 1;
0106   }
0107
0108   P[1][0] = (*IntervalMin) + ((*IntervalMax) - (*IntervalMin)) / 4;
0109   P[1][1] = P[0][1] + (((*IntervalMax) - (*IntervalMin)) / 4.0) * D[0];
0110
0111   D[1] = (FunctionClass->*Derivative)(*IntervalMax);
0112
0113   P[4][0] = (*IntervalMax);
0114   P[4][1] = (FunctionClass->*Function)(*IntervalMax);
0115
0116   P[3][0] = (*IntervalMax) - ((*IntervalMax) - (*IntervalMin)) / 4;
0117   P[3][1] = P[4][1] - ((*IntervalMax) - (*IntervalMin)) / 4 * D[1];
0118
0119   P[2][0] = ((*IntervalMax) + (*IntervalMin)) / 2;
0120   P[2][1] =
0121     (16 * (FunctionClass->*Function)(((*IntervalMax) + (*IntervalMin)) / 2) -
0122      (P[0][1] + 4 * P[1][1] + 4 * P[3][1] + P[4][1])) /
0123     6;
0124
0125   {
0126     G4double Intersection;
0127     G4int i, j;
0128
0129     NewMin = (*IntervalMax);
0130     NewMax = (*IntervalMin);
0131
0132     for(i = 0; i < 5; ++i)
0133       for(j = i + 1; j < 5; ++j)
0134       {
0135         /* there is an intersection only if each have different signs */
0136         if(((P[j][1] > -Precision) && (P[i][1] < Precision)) ||
0137            ((P[j][1] < Precision) && (P[i][1] > -Precision)))
0138         {
0139           IntervalIsVoid = 0;
0140           Intersection =
0141             P[j][0] - P[j][1] * ((P[i][0] - P[j][0]) / (P[i][1] - P[j][1]));
0142           if(Intersection < NewMin)
0143           {
0144             NewMin = Intersection;
0145           }
0146           if(Intersection > NewMax)
0147           {
0148             NewMax = Intersection;
0149           }
0150         }
0151       }
0152
0153     if(IntervalIsVoid != 1)
0154     {
0155       (*IntervalMax) = NewMax;
0156       (*IntervalMin) = NewMin;
0157     }
0158   }
0159
0160   if(IntervalIsVoid == 1)
0161   {
0162     return -1;
0163   }
0164
0165   return 0;
0166 }
0167
0168 template <class T, class F>
0169 G4double G4PolynomialSolver<T, F>::Newton(G4double IntervalMin,
0170                                           G4double IntervalMax)
0171 {
0172   /* So now we have a good guess and an interval where
0173      if there are an intersection the root must be */
0174
0175   G4double Value    = 0;
0176   G4double Gradient = 0;
0177   G4double Lambda;
0178
0179   G4int i = 0;
0180   G4int j = 0;
0181
0182   /* Reduce interval before applying Newton Method */
0183   {
0184     G4int NewtonIsSafe;
0185
0186     while((NewtonIsSafe = BezierClipping(&IntervalMin, &IntervalMax)) == 0)
0187       ;
0188
0189     if(NewtonIsSafe == -1)
0190     {
0191       return POLINFINITY;
0192     }
0193   }
0194
0195   Lambda = IntervalMin;
0196   Value  = (FunctionClass->*Function)(Lambda);
0197
0198   //  while ((std::fabs(Value) > Precision)) {
0199   while(j != -1)
0200   {
0201     Value = (FunctionClass->*Function)(Lambda);
0202
0203     Gradient = (FunctionClass->*Derivative)(Lambda);
0204
0205     Lambda = Lambda - Value / Gradient;
0206
0207     if(std::fabs(Value) <= Precision)
0208     {
0209       ++j;
0210       if(j == 2)
0211       {
0212         j = -1;
0213       }
0214     }
0215     else
0216     {
0217       ++i;
0218
0219       if(i > ITERATION)
0220         return POLINFINITY;
0221     }
0222   }
0223   return Lambda;
0224 }