include/Geant4/G4PolynomialSolver.hh

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
0027 //
0028 // Class description:
0029 //
0030 //   G4PolynomialSolver allows the user to solve a polynomial equation
0031 //   with a great precision. This is used by Implicit Equation solver.
0032 //
0033 //   The Bezier clipping method is used to solve the polynomial.
0034 //
0035 // How to use it:
0036 //   Create a class that is the function to be solved.
0037 //   This class could have internal parameters to allow to change
0038 //   the equation to be solved without recreating a new one.
0039 //
0040 //   Define a Polynomial solver, example:
0041 //   G4PolynomialSolver<MyFunctionClass,G4double(MyFunctionClass::*)(G4double)>
0042 //     PolySolver (&MyFunction,
0043 //                 &MyFunctionClass::Function,
0044 //                 &MyFunctionClass::Derivative,
0045 //                 precision);
0046 //
0047 //   The precision is relative to the function to solve.
0048 //
0049 //   In MyFunctionClass, provide the function to solve and its derivative:
0050 //   Example of function to provide :
0051 //
0052 //   x,y,z,dx,dy,dz,Rmin,Rmax are internal variables of MyFunctionClass
0053 //
0054 //   G4double MyFunctionClass::Function(G4double value)
0055 //   {
0056 //     G4double Lx,Ly,Lz;
0057 //     G4double result;
0058 //
0059 //     Lx = x + value*dx;
0060 //     Ly = y + value*dy;
0061 //     Lz = z + value*dz;
0062 //
0063 //     result = TorusEquation(Lx,Ly,Lz,Rmax,Rmin);
0064 //
0065 //     return result ;
0066 //   }
0067 //
0068 //   G4double MyFunctionClass::Derivative(G4double value)
0069 //   {
0070 //     G4double Lx,Ly,Lz;
0071 //     G4double result;
0072 //
0073 //     Lx = x + value*dx;
0074 //     Ly = y + value*dy;
0075 //     Lz = z + value*dz;
0076 //
0077 //     result = dx*TorusDerivativeX(Lx,Ly,Lz,Rmax,Rmin);
0078 //     result += dy*TorusDerivativeY(Lx,Ly,Lz,Rmax,Rmin);
0079 //     result += dz*TorusDerivativeZ(Lx,Ly,Lz,Rmax,Rmin);
0080 //
0081 //     return result;
0082 //   }
0083 //
0084 //   Then to have a root inside an interval [IntervalMin,IntervalMax] do the
0085 //   following:
0086 //
0087 //   MyRoot = PolySolver.solve(IntervalMin,IntervalMax);
0088
0089 // Author: E.Medernach, 19.12.2000 - First implementation
0090 // --------------------------------------------------------------------
0091 #ifndef G4POL_SOLVER_HH
0092 #define G4POL_SOLVER_HH 1
0093
0094 #include "globals.hh"
0095
0096 template <class T, class F>
0097 class G4PolynomialSolver
0098 {
0099  public:
0100   G4PolynomialSolver(T* typeF, F func, F deriv, G4double precision);
0101   ~G4PolynomialSolver();
0102
0103   G4double solve(G4double IntervalMin, G4double IntervalMax);
0104
0105  private:
0106   G4double Newton(G4double IntervalMin, G4double IntervalMax);
0107   // General Newton method with Bezier Clipping
0108
0109   // Works for polynomial of order less or equal than 4.
0110   // But could be changed to work for polynomial of any order providing
0111   // that we find the bezier control points.
0112
0113   G4int BezierClipping(G4double* IntervalMin, G4double* IntervalMax);
0114   // This is just one iteration of Bezier Clipping
0115
0116   T* FunctionClass;
0117   F Function;
0118   F Derivative;
0119
0120   G4double Precision;
0121 };
0122
0123 #include "G4PolynomialSolver.icc"
0124
0125 #endif