00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00046 #include "G4INCLRootFinder.hh"
00047 #include "G4INCLGlobals.hh"
00048 #include "G4INCLLogger.hh"
00049 #include <utility>
00050 #include <cmath>
00051
00052 namespace G4INCL {
00053
00054 std::pair<G4double,G4double> RootFinder::solution;
00055
00056 const G4double RootFinder::toleranceY = 1.e-4;
00057
00058 G4bool RootFinder::solve(RootFunctor const * const f, const G4double x0) {
00059
00060 const G4double y0 = (*f)(x0);
00061 if( std::abs(y0) < toleranceY ) {
00062 solution = std::make_pair(x0,y0);
00063 return true;
00064 }
00065
00066
00067 std::pair<G4double,G4double> bracket = bracketRoot(f,x0);
00068 G4double x1 = bracket.first;
00069 G4double x2 = bracket.second;
00070
00071 if(x1>x2) {
00072
00073 G4double y_at_zero = (*f)(0.);
00074 if(std::abs(y_at_zero)<=toleranceY) {
00075 f->cleanUp(true);
00076 solution = std::make_pair(0.,y_at_zero);
00077 return true;
00078 } else {
00079 WARN("Root-finding algorithm could not bracket the root." << std::endl);
00080 f->cleanUp(false);
00081 return false;
00082 }
00083 }
00084
00085 G4double y1 = (*f)(x1);
00086 G4double y2 = (*f)(x2);
00087 G4double x = x1;
00088 G4double y = y1;
00089
00090
00091
00092
00093
00094
00095 G4int lastUpdated = 0;
00096
00097 for(G4int iterations=0; std::abs(y) > toleranceY; iterations++) {
00098
00099 if(iterations > maxIterations) {
00100 WARN("Root-finding algorithm did not converge." << std::endl);
00101 f->cleanUp(false);
00102 return false;
00103 }
00104
00105
00106 x = (y1*x2-y2*x1)/(y1-y2);
00107
00108
00109 y = (*f)(x);
00110
00111
00112 if(Math::sign(y) == Math::sign(y1)) {
00113 x1=x;
00114 y1=y;
00115 if(lastUpdated==-1) y2 *= 0.5;
00116 lastUpdated = -1;
00117 } else {
00118 x2=x;
00119 y2=y;
00120 if(lastUpdated==1) y1 *= 0.5;
00121 lastUpdated = 1;
00122 }
00123 }
00124
00125
00126
00127
00128
00129 solution = std::make_pair(x,y);
00130 f->cleanUp(true);
00131 return true;
00132 }
00133
00134 std::pair<G4double,G4double> RootFinder::bracketRoot(RootFunctor const * const f, G4double x0) {
00135 G4double y0 = (*f)(x0);
00136
00137 const G4double scaleFactor = 1.5;
00138
00139 G4double x1;
00140 if(x0!=0.)
00141 x1=scaleFactor*x0;
00142 else
00143 x1=1.;
00144 G4double y1 = (*f)(x1);
00145
00146 if(Math::sign(y0)!=Math::sign(y1))
00147 return std::make_pair(x0,x1);
00148
00149 const G4double scaleFactorMinus1 = 1./scaleFactor;
00150 G4double oldx0, oldx1, oldy1;
00151 G4int iterations=0;
00152 do {
00153 if(iterations > maxIterations) {
00154 DEBUG("Could not bracket the root." << std::endl);
00155 return std::make_pair((G4double) 1.,(G4double) -1.);
00156 }
00157
00158 oldx0=x0;
00159 oldx1=x1;
00160 oldy1=y1;
00161
00162 x0 *= scaleFactorMinus1;
00163 x1 *= scaleFactor;
00164 y0 = (*f)(x0);
00165 y1 = (*f)(x1);
00166 iterations++;
00167 } while(Math::sign(y0)==Math::sign(y1));
00168
00169 if(Math::sign(y1)==Math::sign(oldy1))
00170 return std::make_pair(x0,oldx0);
00171 else
00172 return std::make_pair(oldx1,x1);
00173 }
00174
00175 }