Geant4: G4Bessel.cc Source File

00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // *                                                                  *
00021 // * Parts of this code which have been  developed by QinetiQ Ltd     *
00022 // * under contract to the European Space Agency (ESA) are the        *
00023 // * intellectual property of ESA. Rights to use, copy, modify and    *
00024 // * redistribute this software for general public use are granted    *
00025 // * in compliance with any licensing, distribution and development   *
00026 // * policy adopted by the Geant4 Collaboration. This code has been   *
00027 // * written by QinetiQ Ltd for the European Space Agency, under ESA  *
00028 // * contract 17191/03/NL/LvH (Aurora Programme).                     *
00029 // *                                                                  *
00030 // * By using,  copying,  modifying or  distributing the software (or *
00031 // * any work based  on the software)  you  agree  to acknowledge its *
00032 // * use  in  resulting  scientific  publications,  and indicate your *
00033 // * acceptance of all terms of the Geant4 Software license.          *
00034 // ********************************************************************
00035 //
00036 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00037 //
00038 // MODULE:              G4Bessel.cc
00039 //
00040 // Version:             B.1
00041 // Date:                15/04/04
00042 // Author:              P R Truscott
00043 // Organisation:        QinetiQ Ltd, UK
00044 // Customer:            ESA/ESTEC, NOORDWIJK
00045 // Contract:            17191/03/NL/LvH
00046 //
00047 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00048 //
00049 // CHANGE HISTORY
00050 // --------------
00051 //
00052 // 18 Noevmber 2003, P R Truscott, QinetiQ Ltd, UK
00053 // Created.
00054 //
00055 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
00056 // Beta release
00057 //
00058 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00060 //
00061 #include "G4Bessel.hh"
00062 #include "G4PhysicalConstants.hh"
00064 //
00065 G4Bessel::G4Bessel ()
00066 {;}
00068 //
00069 G4Bessel::~G4Bessel ()
00070 {;}
00072 //
00073 G4double G4Bessel::I0 (G4double x)
00074 {
00075   const G4double P1 = 1.0,
00076                  P2 = 3.5156229,
00077                  P3 = 3.0899424,
00078                  P4 = 1.2067492,
00079                  P5 = 0.2659732,
00080                  P6 = 0.0360768,
00081                  P7 = 0.0045813;
00082   const G4double Q1 = 0.39894228,
00083                  Q2 = 0.01328592,
00084                  Q3 = 0.00225319,
00085                  Q4 =-0.00157565,
00086                  Q5 = 0.00916281,
00087                  Q6 =-0.02057706,
00088                  Q7 = 0.02635537,
00089                  Q8 =-0.01647633,
00090                  Q9 = 0.00392377;
00091   
00092   G4double I = 0.0;
00093   if (std::fabs(x) < 3.75)
00094   {
00095     G4double y = std::pow(x/3.75, 2.0);
00096     I = P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
00097   }
00098   else
00099   {
00100     G4double ax = std::fabs(x);
00101     G4double y  = 3.75/ax;
00102     I  = std::exp(ax) / std::sqrt(ax) *
00103       (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
00104   }
00105   return I;
00106 }
00108 //
00109 G4double G4Bessel::K0 (G4double x)
00110 {
00111   const G4double P1 =-0.57721566,
00112                  P2 = 0.42278420,
00113                  P3 = 0.23069756,
00114                  P4 = 0.03488590,
00115                  P5 = 0.00262698,
00116                  P6 = 0.00010750,
00117                  P7 = 0.00000740;
00118   const G4double Q1 = 1.25331414,
00119                  Q2 =-0.07832358,
00120                  Q3 = 0.02189568,
00121                  Q4 =-0.01062446,
00122                  Q5 = 0.00587872,
00123                  Q6 =-0.00251540,
00124                  Q7 = 0.00053208;
00125 
00126   G4double K = 0.0;
00127   if (x <= 2.0)
00128   {
00129     G4double y = x * x / 4.0;
00130     K = (-std::log(x/2.0)) * I0(x) +
00131       P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
00132   }
00133   else
00134   {
00135     G4double y = 2.0 / x;
00136     K = std::exp(-x)  / std::sqrt(x) *
00137       (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
00138   }
00139   return K;
00140 }
00142 //
00143 G4double G4Bessel::I1 (G4double x)
00144 {
00145   const G4double P1 = 0.5,
00146                  P2 = 0.87890594,
00147                  P3 = 0.51498869,
00148                  P4 = 0.15084934,
00149                  P5 = 0.02658733,
00150                  P6 = 0.00301532,
00151                  P7 = 0.00032411;
00152   const G4double Q1 = 0.39894228,
00153                  Q2 =-0.03988024,
00154                  Q3 =-0.00362018,
00155                  Q4 = 0.00163801,
00156                  Q5 =-0.01031555,
00157                  Q6 = 0.02282967,
00158                  Q7 =-0.02895312,
00159                  Q8 = 0.01787654,
00160                  Q9 =-0.00420059;
00161 
00162   G4double I = 0.0;
00163   if (std::fabs(x) < 3.75)
00164   {
00165     G4double ax = std::fabs(x);
00166     G4double y = std::pow(x/3.75, 2.0);
00167     I = ax*(P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
00168   }
00169   else
00170   {
00171     G4double ax = std::fabs(x);
00172     G4double y  = 3.75/ax;
00173     I  = std::exp(ax) / std::sqrt(ax) *
00174       (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
00175   }
00176   if (x < 0.0) I = -I;
00177   return I;
00178 }
00180 //
00181 G4double G4Bessel::K1 (G4double x)
00182 {
00183   const G4double P1 = 1.0,
00184                  P2 = 0.15443144,
00185                  P3 =-0.67278579,
00186                  P4 =-0.18156897,
00187                  P5 =-0.01919402,
00188                  P6 =-0.00110404,
00189                  P7 =-0.00004686;
00190   const G4double Q1 = 1.25331414,
00191                  Q2 = 0.23498619,
00192                  Q3 =-0.03655620,
00193                  Q4 = 0.01504268,
00194                  Q5 =-0.00780353,
00195                  Q6 = 0.00325614,
00196                  Q7 =-0.00068245;
00197 
00198   G4double K = 0.0;
00199   if (x <= 2.0)
00200   {
00201     G4double y = x * x / 4.0;
00202     K = std::log(x/2.0)*I1(x) + 1.0/x *
00203       (P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
00204   }
00205   else
00206   {
00207     G4double y = 2.0 / x;
00208     K = std::exp(-x) / std::sqrt(x) *
00209       (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
00210   }
00211   return K;
00212 }
00214 //
00215 G4double G4Bessel::pI0 (G4double x)
00216 {
00217   const G4double A0  = 0.1250000000000E+00,
00218                  A1  = 7.0312500000000E-02,
00219                  A2  = 7.3242187500000E-02,
00220                  A3  = 1.1215209960938E-01,
00221                  A4  = 2.2710800170898E-01,
00222                  A5  = 5.7250142097473E-01,
00223                  A6  = 1.7277275025845E+00,
00224                  A7  = 6.0740420012735E+00,
00225                  A8  = 2.4380529699556E+01,
00226                  A9  = 1.1001714026925E+02,
00227                  A10 = 5.5133589612202E+02,
00228                  A11 = 3.0380905109224E+03;
00229 
00230   G4double I = 0.0;
00231   if (x == 0.0)
00232   {
00233     I = 1.0;
00234   }
00235   else if (x < 18.0)
00236   {
00237     I          = 1.0;
00238     G4double y = x * x;
00239     G4double q = 1.0;
00240     for (G4int i=1; i<101; i++)
00241     {
00242       q *= 0.25 * y / i / i;
00243       I += q;
00244       if (std::abs(q/I) < 1.0E-15) break;
00245     }
00246   }
00247   else
00248   {
00249     G4double y = 1.0 / x;
00250     I = std::exp(x) / std::sqrt(twopi*x) *
00251       (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*(A7+y*(A8+y*(A9+y*(A10+y*A11))))))))))));
00252   }
00253 
00254   return I;
00255 }
00257 //
00258 G4double G4Bessel::pI1 (G4double x)
00259 {
00260   const G4double A0  = -0.3750000000000E+00,
00261                  A1  = -1.1718750000000E-01,
00262                  A2  = -1.0253906250000E-01,
00263                  A3  = -1.4419555664063E-01,
00264                  A4  = -2.775764465332E-01,
00265                  A5  = -6.7659258842468E-01,
00266                  A6  = -1.9935317337513E+00,
00267                  A7  = -6.8839142681099E+00,
00268                  A8  = -2.7248827311269E+01,
00269                  A9  = -1.2159789187654E+02,
00270                  A10 = -6.0384407670507E+02,
00271                  A11 = -3.3022722944809E+03;
00272 
00273   G4double I = 0.0;
00274   if (x == 0.0)
00275   {
00276     I = 0.0;
00277   }
00278   else if (x < 18.0)
00279   {
00280     I          = 1.0;
00281     G4double y = x * x;
00282     G4double q = 1.0;
00283     for (G4int i=1; i<101; i++)
00284     {
00285       q *= 0.25 * y / i / (i+1.0);
00286       I += q;
00287       if (std::abs(q/I) < 1.0E-15) break;
00288     }
00289     I *= 0.5 * x;
00290     
00291   }
00292   else
00293   {
00294     G4double y = 1.0 / x;
00295     I = std::exp(x) / std::sqrt(twopi*x) *
00296       (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*(A7+y*(A8+y*(A9+y*(A10+y*A11))))))))))));
00297   }
00298 
00299   return I;
00300 }
00302 //
00303 G4double G4Bessel::pK0 (G4double x)
00304 {
00305   const G4double A0 = 0.1250000000000E+00,
00306                  A1 = 0.2109375000000E+00,
00307                  A2 = 1.0986328125000E+00,
00308                  A3 = 1.1775970458984E+01,
00309                  A4 = 2.1461706161499E+02,
00310                  A5 = 5.9511522710323E+03,
00311                  A6 = 2.3347645606175E+05,
00312                  A7 = 1.2312234987631E+07;
00313 
00314   G4double K = 0.0;
00315   if (x == 0.0)
00316   {
00317     K = 1.0E+307;
00318   }
00319   else if (x < 9.0)
00320   {
00321     G4double y = x * x;
00322     G4double C = -std::log(x/2.0) - 0.5772156649015329;
00323     G4double q = 1.0;
00324     G4double t = 0.0;
00325     for (G4int i=1; i<51; i++)
00326     {
00327       q *= 0.25 * y / i / i;
00328       t += 1.0 / i ;
00329       K += q * (t+C);
00330     }
00331     K += C;
00332   }
00333   else
00334   {
00335     G4double y = 1.0 / x / x;
00336     K = 0.5 / x / pI0(x) *
00337       (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*A7))))))));
00338   }
00339   
00340   return K;
00341 }
00343 //
00344 G4double G4Bessel::pK1 (G4double x)
00345 {
00346   G4double K = 0.0;
00347   if (x == 0.0)
00348     K = 1.0E+307;
00349   else
00350     K = (1.0/x - pI1(x)*pK0(x)) / pI0(x);
00351   return K;
00352 }
00354 //