Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/src/G4AnalyticalPolSolver.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /global/HEPNumerics/src/G4AnalyticalPolSolver.cc (Version 11.3.0) and /global/HEPNumerics/src/G4AnalyticalPolSolver.cc (Version 9.1)


  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 // G4AnalyticalPolSolver class implementation  << 
 27 //                                                 26 //
 28 // Author: V.Grichine, 24.04.97                <<  27 // $Id: G4AnalyticalPolSolver.cc,v 1.7 2007/11/13 17:35:06 gcosmo Exp $
 29 // ------------------------------------------- <<  28 // GEANT4 tag $Name: geant4-09-01 $
                                                   >>  29 //
 30                                                    30 
 31 #include "globals.hh"                          <<  31 #include  "globals.hh"
 32 #include <complex>                             <<  32 #include  <complex>
 33                                                    33 
 34 #include "G4AnalyticalPolSolver.hh"                34 #include "G4AnalyticalPolSolver.hh"
 35                                                    35 
 36 //////////////////////////////////////////////     36 //////////////////////////////////////////////////////////////////////////////
 37                                                    37 
 38 G4AnalyticalPolSolver::G4AnalyticalPolSolver() <<  38 G4AnalyticalPolSolver::G4AnalyticalPolSolver() {;}
 39                                                    39 
 40 //////////////////////////////////////////////     40 //////////////////////////////////////////////////////////////////////////////
 41                                                    41 
 42 G4AnalyticalPolSolver::~G4AnalyticalPolSolver( <<  42 G4AnalyticalPolSolver::~G4AnalyticalPolSolver() {;}
 43                                                    43 
 44 //////////////////////////////////////////////     44 //////////////////////////////////////////////////////////////////////////////
 45 //                                                 45 //
 46 // Array r[3][5]  p[5]                             46 // Array r[3][5]  p[5]
 47 // Roots of poly p[0] x^2 + p[1] x+p[2]=0          47 // Roots of poly p[0] x^2 + p[1] x+p[2]=0
 48 //                                                 48 //
 49 // x = r[1][k] + i r[2][k];  k = 1, 2              49 // x = r[1][k] + i r[2][k];  k = 1, 2
 50                                                    50 
 51 G4int G4AnalyticalPolSolver::QuadRoots(G4doubl <<  51 G4int G4AnalyticalPolSolver::QuadRoots( G4double p[5], G4double r[3][5] )
 52 {                                                  52 {
 53   G4double b, c, d2, d;                            53   G4double b, c, d2, d;
 54                                                    54 
 55   b  = -p[1] / p[0] / 2.;                      <<  55   b  = -p[1]/p[0]/2.;
 56   c  = p[2] / p[0];                            <<  56   c  =  p[2]/p[0];
 57   d2 = b * b - c;                              <<  57   d2 =  b*b - c;
 58                                                <<  58   
 59   if(d2 >= 0.)                                 <<  59   if( d2 >= 0. )
 60   {                                                60   {
 61     d       = std::sqrt(d2);                       61     d       = std::sqrt(d2);
 62     r[1][1] = b - d;                           <<  62     r[1][1] = b - d;   
 63     r[1][2] = b + d;                           <<  63     r[1][2] = b + d;   
 64     r[2][1] = 0.;                              <<  64     r[2][1] = 0.; 
 65     r[2][2] = 0.;                                  65     r[2][2] = 0.;
 66   }                                                66   }
 67   else                                             67   else
 68   {                                                68   {
 69     d       = std::sqrt(-d2);                      69     d       = std::sqrt(-d2);
 70     r[2][1] = d;                               <<  70     r[2][1] =  d; 
 71     r[2][2] = -d;                                  71     r[2][2] = -d;
 72     r[1][1] = b;                               <<  72     r[1][1] =  b; 
 73     r[1][2] = b;                               <<  73     r[1][2] =  b;
 74   }                                                74   }
 75                                                    75 
 76   return 2;                                        76   return 2;
 77 }                                                  77 }
 78                                                    78 
 79 //////////////////////////////////////////////     79 //////////////////////////////////////////////////////////////////////////////
 80 //                                                 80 //
 81 // Array r[3][5]  p[5]                             81 // Array r[3][5]  p[5]
 82 // Roots of poly p[0] x^3 + p[1] x^2...+p[3]=0     82 // Roots of poly p[0] x^3 + p[1] x^2...+p[3]=0
 83 // x=r[1][k] + i r[2][k]  k=1,...,3                83 // x=r[1][k] + i r[2][k]  k=1,...,3
 84 // Assumes 0<arctan(x)<pi/2 for x>0                84 // Assumes 0<arctan(x)<pi/2 for x>0
 85                                                    85 
 86 G4int G4AnalyticalPolSolver::CubicRoots(G4doub <<  86 G4int G4AnalyticalPolSolver::CubicRoots( G4double p[5], G4double r[3][5] )
 87 {                                                  87 {
 88   G4double x, t, b, c, d;                      <<  88   G4double x,t,b,c,d;
 89   G4int k;                                         89   G4int k;
 90                                                    90 
 91   if(p[0] != 1.)                               <<  91   if( p[0] != 1. )
 92   {                                                92   {
 93     for(k = 1; k < 4; k++)                     <<  93     for(k = 1; k < 4; k++ ) { p[k] = p[k]/p[0]; }
 94     {                                          << 
 95       p[k] = p[k] / p[0];                      << 
 96     }                                          << 
 97     p[0] = 1.;                                     94     p[0] = 1.;
 98   }                                                95   }
 99   x = p[1] / 3.0;                              <<  96   x = p[1]/3.0; 
100   t = x * p[1];                                <<  97   t = x*p[1];
101   b = 0.5 * (x * (t / 1.5 - p[2]) + p[3]);     <<  98   b = 0.5*( x*( t/1.5 - p[2] ) + p[3] ); 
102   t = (t - p[2]) / 3.0;                        <<  99   t = ( t - p[2] )/3.0;
103   c = t * t * t;                               << 100   c = t*t*t; 
104   d = b * b - c;                               << 101   d = b*b - c;
105                                                   102 
106   if(d >= 0.)                                  << 103   if( d >= 0. )
107   {                                               104   {
108     d = std::pow((std::sqrt(d) + std::fabs(b)) << 105     d = std::pow( (std::sqrt(d) + std::fabs(b) ), 1.0/3.0 );
109                                                << 106     
110     if(d != 0.)                                << 107     if( d != 0. )
111     {                                             108     {
112       if(b > 0.)                               << 109        if( b > 0. ) { b = -d; }
113       {                                        << 110        else         { b =  d; }
114         b = -d;                                << 111        c =  t/b;
115       }                                        << 
116       else                                     << 
117       {                                        << 
118         b = d;                                 << 
119       }                                        << 
120       c = t / b;                               << 
121     }                                             112     }
122     d       = std::sqrt(0.75) * (b - c);       << 113     d       =  std::sqrt(0.75)*(b - c); 
123     r[2][2] = d;                               << 114     r[2][2] =  d; 
124     b       = b + c;                           << 115     b       =  b + c;
125     c       = -0.5 * b - x;                    << 116     c       = -0.5*b-x;
126     r[1][2] = c;                               << 117     r[1][2] =  c;
127                                                << 118 
128     if((b > 0. && x <= 0.) || (b < 0. && x > 0 << 119     if( ( b > 0. &&  x <= 0. ) || ( b < 0. && x > 0. ) )
129     {                                          << 120     {
130       r[1][1] = c;                             << 121        r[1][1] =  c; 
131       r[2][1] = -d;                            << 122        r[2][1] = -d; 
132       r[1][3] = b - x;                         << 123        r[1][3] =  b - x;
133       r[2][3] = 0;                             << 124        r[2][3] =  0;
134     }                                             125     }
135     else                                          126     else
136     {                                             127     {
137       r[1][1] = b - x;                         << 128        r[1][1] =  b - x; 
138       r[2][1] = 0.;                            << 129        r[2][1] =  0.; 
139       r[1][3] = c;                             << 130        r[1][3] =  c;
140       r[2][3] = -d;                            << 131        r[2][3] = -d;
141     }                                             132     }
142   }  // end of 2 equal or complex roots        << 133   }              // end of 2 equal or complex roots 
143   else                                            134   else
144   {                                               135   {
145     if(b == 0.)                                << 136     if( b == 0. ) { d =  std::atan(1.0)/1.5; }
146     {                                          << 137     else          { d =  std::atan( std::sqrt(-d)/std::fabs(b) )/3.0; }
147       d = std::atan(1.0) / 1.5;                << 
148     }                                          << 
149     else                                       << 
150     {                                          << 
151       d = std::atan(std::sqrt(-d) / std::fabs( << 
152     }                                          << 
153                                                   138 
154     if(b < 0.)                                 << 139     if( b < 0. )  { b =  std::sqrt(t)*2.0; }
155     {                                          << 140     else          { b = -2.0*std::sqrt(t); }
156       b = std::sqrt(t) * 2.0;                  << 
157     }                                          << 
158     else                                       << 
159     {                                          << 
160       b = -2.0 * std::sqrt(t);                 << 
161     }                                          << 
162                                                   141 
163     c = std::cos(d) * b;                       << 142     c =  std::cos(d)*b; 
164     t = -std::sqrt(0.75) * std::sin(d) * b - 0 << 143     t = -std::sqrt(0.75)*std::sin(d)*b - 0.5*c;
165     d = -t - c - x;                            << 144     d = -t - c - x; 
166     c = c - x;                                 << 145     c =  c - x; 
167     t = t - x;                                 << 146     t =  t - x;
168                                                   147 
169     if(std::fabs(c) > std::fabs(t))            << 148     if( std::fabs(c) > std::fabs(t) ) { r[1][3] = c; }
170     {                                          << 
171       r[1][3] = c;                             << 
172     }                                          << 
173     else                                          149     else
174     {                                             150     {
175       r[1][3] = t;                             << 151       r[1][3] = t; 
176       t       = c;                                152       t       = c;
177     }                                             153     }
178     if(std::fabs(d) > std::fabs(t))            << 154     if( std::fabs(d) > std::fabs(t) ) { r[1][2] = d; }
179     {                                          << 
180       r[1][2] = d;                             << 
181     }                                          << 
182     else                                          155     else
183     {                                             156     {
184       r[1][2] = t;                             << 157       r[1][2] = t; 
185       t       = d;                                158       t       = d;
186     }                                             159     }
187     r[1][1] = t;                                  160     r[1][1] = t;
188                                                   161 
189     for(k = 1; k < 4; k++)                     << 162     for(k = 1; k < 4; k++ ) { r[2][k] = 0.; }
190     {                                          << 
191       r[2][k] = 0.;                            << 
192     }                                          << 
193   }                                               163   }
194   return 0;                                       164   return 0;
195 }                                                 165 }
196                                                   166 
197 //////////////////////////////////////////////    167 //////////////////////////////////////////////////////////////////////////////
198 //                                                168 //
199 // Array r[3][5]  p[5]                            169 // Array r[3][5]  p[5]
200 // Roots of poly p[0] x^4 + p[1] x^3...+p[4]=0    170 // Roots of poly p[0] x^4 + p[1] x^3...+p[4]=0
201 // x=r[1][k] + i r[2][k]  k=1,...,4               171 // x=r[1][k] + i r[2][k]  k=1,...,4
202                                                   172 
203 G4int G4AnalyticalPolSolver::BiquadRoots(G4dou << 173 G4int G4AnalyticalPolSolver::BiquadRoots( G4double p[5], G4double r[3][5] )
204 {                                                 174 {
205   G4double a, b, c, d, e;                         175   G4double a, b, c, d, e;
206   G4int i, k, j;                               << 176   G4int i, k, j, noRoots;
207                                                   177 
208   if(p[0] != 1.0)                                 178   if(p[0] != 1.0)
209   {                                               179   {
210     for(k = 1; k < 5; k++)                     << 180     for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
211     {                                          << 
212       p[k] = p[k] / p[0];                      << 
213     }                                          << 
214     p[0] = 1.;                                    181     p[0] = 1.;
215   }                                               182   }
216   e = 0.25 * p[1];                             << 183   e = 0.25*p[1];
217   b = 2 * e;                                   << 184   b = 2*e;
218   c = b * b;                                   << 185   c = b*b;
219   d = 0.75 * c;                                << 186   d = 0.75*c;
220   b = p[3] + b * (c - p[2]);                   << 187   b = p[3] + b*( c - p[2] );
221   a = p[2] - d;                                   188   a = p[2] - d;
222   c = p[4] + e * (e * a - p[3]);               << 189   c = p[4] + e*( e*a - p[3] );
223   a = a - d;                                      190   a = a - d;
224                                                   191 
225   p[1] = 0.5 * a;                              << 192   p[1] = 0.5*a;
226   p[2] = (p[1] * p[1] - c) * 0.25;             << 193   p[2] = (p[1]*p[1]-c)*0.25;
227   p[3] = b * b / (-64.0);                      << 194   p[3] = b*b/(-64.0);
228                                                   195 
229   if(p[3] < 0.)                                << 196   if( p[3] < 0. )
230   {                                               197   {
231     CubicRoots(p, r);                          << 198     noRoots = CubicRoots(p,r);
232                                                   199 
233     for(k = 1; k < 4; k++)                     << 200     for( k = 1; k < 4; k++ )
234     {                                             201     {
235       if(r[2][k] == 0. && r[1][k] > 0)         << 202       if( r[2][k] == 0. && r[1][k] > 0 )
236       {                                           203       {
237         d = r[1][k] * 4;                       << 204         d = r[1][k]*4; 
238         a = a + d;                                205         a = a + d;
239                                                   206 
240         if(a >= 0. && b >= 0.)                 << 207         if     ( a >= 0. && b >= 0.) { p[1] =  std::sqrt(d); }
241         {                                      << 208         else if( a <= 0. && b <= 0.) { p[1] =  std::sqrt(d); }
242           p[1] = std::sqrt(d);                 << 209         else                         { p[1] = -std::sqrt(d); }
243         }                                      << 
244         else if(a <= 0. && b <= 0.)            << 
245         {                                      << 
246           p[1] = std::sqrt(d);                 << 
247         }                                      << 
248         else                                   << 
249         {                                      << 
250           p[1] = -std::sqrt(d);                << 
251         }                                      << 
252                                                   210 
253         b = 0.5 * (a + b / p[1]);              << 211         b = 0.5*( a + b/p[1] );
254                                                   212 
255         p[2] = c / b;                          << 213         p[2]    = c/b; 
256         QuadRoots(p, r);                       << 214         noRoots = QuadRoots(p,r);
257                                                   215 
258         for(i = 1; i < 3; i++)                 << 216         for( i = 1; i < 3; i++ )
259         {                                         217         {
260           for(j = 1; j < 3; j++)               << 218           for( j = 1; j < 3; j++ ) { r[j][i+2] = r[j][i]; }
261           {                                    << 
262             r[j][i + 2] = r[j][i];             << 
263           }                                    << 
264         }                                         219         }
265         p[1] = -p[1];                          << 220         p[1]    = -p[1]; 
266         p[2] = b;                              << 221         p[2]    =  b; 
267         QuadRoots(p, r);                       << 222         noRoots = QuadRoots(p,r);
268                                                   223 
269         for(i = 1; i < 5; i++)                 << 224         for( i = 1; i < 5; i++ ) { r[1][i] = r[1][i] - e; }
270         {                                      << 
271           r[1][i] = r[1][i] - e;               << 
272         }                                      << 
273                                                   225 
274         return 4;                                 226         return 4;
275       }                                           227       }
276     }                                             228     }
277   }                                               229   }
278   if(p[2] < 0.)                                << 230   if( p[2] < 0. )
279   {                                               231   {
280     b    = std::sqrt(c);                       << 232     b    = std::sqrt(c); 
281     d    = b + b - a;                             233     d    = b + b - a;
282     p[1] = 0.;                                 << 234     p[1] = 0.; 
283                                                << 235     
284     if(d > 0.)                                 << 236     if( d > 0. ) { p[1] = std::sqrt(d); }
285     {                                          << 
286       p[1] = std::sqrt(d);                     << 
287     }                                          << 
288   }                                               237   }
289   else                                            238   else
290   {                                               239   {
291     if(p[1] > 0.)                              << 240     if( p[1] > 0.) { b =  std::sqrt(p[2])*2.0 + p[1]; }
292     {                                          << 241     else           { b = -std::sqrt(p[2])*2.0 + p[1]; }
293       b = std::sqrt(p[2]) * 2.0 + p[1];        << 
294     }                                          << 
295     else                                       << 
296     {                                          << 
297       b = -std::sqrt(p[2]) * 2.0 + p[1];       << 
298     }                                          << 
299                                                   242 
300     if(b != 0.)                                << 243     if( b != 0.) { p[1] = 0; }
301     {                                          << 
302       p[1] = 0;                                << 
303     }                                          << 
304     else                                          244     else
305     {                                             245     {
306       for(k = 1; k < 5; k++)                   << 246       for(k = 1; k < 5; k++ )
307       {                                           247       {
308         r[1][k] = -e;                          << 248           r[1][k] = -e;
309         r[2][k] = 0;                           << 249           r[2][k] =  0;
310       }                                           250       }
311       return 0;                                   251       return 0;
312     }                                             252     }
313   }                                               253   }
314                                                   254 
315   p[2] = c / b;                                << 255   p[2]    = c/b; 
316   QuadRoots(p, r);                             << 256   noRoots = QuadRoots(p,r);
317                                                   257 
318   for(k = 1; k < 3; k++)                       << 258   for( k = 1; k < 3; k++ )
319   {                                               259   {
320     for(j = 1; j < 3; j++)                     << 260     for( j = 1; j < 3; j++ ) { r[j][k+2] = r[j][k]; }
321     {                                          << 
322       r[j][k + 2] = r[j][k];                   << 
323     }                                          << 
324   }                                               261   }
325   p[1] = -p[1];                                << 262   p[1]    = -p[1]; 
326   p[2] = b;                                    << 263   p[2]    =  b; 
327   QuadRoots(p, r);                             << 264   noRoots = QuadRoots(p,r);
328                                                   265 
329   for(k = 1; k < 5; k++)                       << 266   for( k = 1; k < 5; k++ ) { r[1][k] = r[1][k] - e; }
330   {                                            << 
331     r[1][k] = r[1][k] - e;                     << 
332   }                                            << 
333                                                   267 
334   return 4;                                       268   return 4;
335 }                                                 269 }
336                                                   270 
337 //////////////////////////////////////////////    271 //////////////////////////////////////////////////////////////////////////////
338                                                   272 
339 G4int G4AnalyticalPolSolver::QuarticRoots(G4do << 273 G4int G4AnalyticalPolSolver::QuarticRoots( G4double p[5], G4double r[3][5])
340 {                                                 274 {
341   G4double a0, a1, a2, a3, y1;                    275   G4double a0, a1, a2, a3, y1;
342   G4double R2, D2, E2, D, E, R = 0.;              276   G4double R2, D2, E2, D, E, R = 0.;
343   G4double a, b, c, d, ds;                        277   G4double a, b, c, d, ds;
344                                                   278 
345   G4double reRoot[4];                             279   G4double reRoot[4];
346   G4int k;                                     << 280   G4int k, noRoots, noReRoots = 0;
                                                   >> 281   
                                                   >> 282   for( k = 0; k < 4; k++ ) { reRoot[k] = DBL_MAX; }
347                                                   283 
348   for(k = 0; k < 4; k++)                       << 284   if( p[0] != 1.0 )
349   {                                            << 
350     reRoot[k] = DBL_MAX;                       << 
351   }                                            << 
352                                                << 
353   if(p[0] != 1.0)                              << 
354   {                                               285   {
355     for(k = 1; k < 5; k++)                     << 286     for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
356     {                                          << 
357       p[k] = p[k] / p[0];                      << 
358     }                                          << 
359     p[0] = 1.;                                    287     p[0] = 1.;
360   }                                               288   }
361   a3 = p[1];                                      289   a3 = p[1];
362   a2 = p[2];                                      290   a2 = p[2];
363   a1 = p[3];                                      291   a1 = p[3];
364   a0 = p[4];                                      292   a0 = p[4];
365                                                   293 
366   // resolvent cubic equation cofs:               294   // resolvent cubic equation cofs:
367                                                   295 
368   p[1] = -a2;                                     296   p[1] = -a2;
369   p[2] = a1 * a3 - 4 * a0;                     << 297   p[2] = a1*a3 - 4*a0;
370   p[3] = 4 * a2 * a0 - a1 * a1 - a3 * a3 * a0; << 298   p[3] = 4*a2*a0 - a1*a1 - a3*a3*a0;
371                                                   299 
372   CubicRoots(p, r);                            << 300   noRoots = CubicRoots(p,r);
373                                                   301 
374   for(k = 1; k < 4; k++)                       << 302   for( k = 1; k < 4; k++ )
375   {                                               303   {
376     if(r[2][k] == 0.)  // find a real root     << 304     if( r[2][k] == 0. ) // find a real root
377     {                                             305     {
                                                   >> 306       noReRoots++;
378       reRoot[k] = r[1][k];                        307       reRoot[k] = r[1][k];
379     }                                             308     }
380     else                                       << 309     else reRoot[k] = DBL_MAX; // kInfinity;
381     {                                          << 
382       reRoot[k] = DBL_MAX;  // kInfinity;      << 
383     }                                          << 
384   }                                               310   }
385   y1 = DBL_MAX;  // kInfinity;                 << 311   y1 = DBL_MAX; // kInfinity;  
386   for(k = 1; k < 4; k++)                       << 312   for( k = 1; k < 4; k++ )
387   {                                               313   {
388     if(reRoot[k] < y1)                         << 314     if ( reRoot[k] < y1 ) { y1 = reRoot[k]; }
389     {                                          << 
390       y1 = reRoot[k];                          << 
391     }                                          << 
392   }                                               315   }
393                                                   316 
394   R2 = 0.25 * a3 * a3 - a2 + y1;               << 317   R2 = 0.25*a3*a3 - a2 + y1;
395   b  = 0.25 * (4 * a3 * a2 - 8 * a1 - a3 * a3  << 318   b  = 0.25*(4*a3*a2 - 8*a1 - a3*a3*a3);
396   c  = 0.75 * a3 * a3 - 2 * a2;                << 319   c  = 0.75*a3*a3 - 2*a2;
397   a  = c - R2;                                    320   a  = c - R2;
398   d  = 4 * y1 * y1 - 16 * a0;                  << 321   d  = 4*y1*y1 - 16*a0;
399                                                   322 
400   if(R2 > 0.)                                  << 323   if( R2 > 0.)
401   {                                               324   {
402     R  = std::sqrt(R2);                        << 325     R = std::sqrt(R2);
403     D2 = a + b / R;                            << 326     D2 = a + b/R;
404     E2 = a - b / R;                            << 327     E2 = a - b/R;
405                                                   328 
406     if(D2 >= 0.)                               << 329     if( D2 >= 0. )
407     {                                             330     {
408       D       = std::sqrt(D2);                    331       D       = std::sqrt(D2);
409       r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * D << 332       r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
410       r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * D << 333       r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
411       r[2][1] = 0.;                               334       r[2][1] = 0.;
412       r[2][2] = 0.;                               335       r[2][2] = 0.;
413     }                                             336     }
414     else                                          337     else
415     {                                             338     {
416       D       = std::sqrt(-D2);                   339       D       = std::sqrt(-D2);
417       r[1][1] = -0.25 * a3 + 0.5 * R;          << 340       r[1][1] = -0.25*a3 + 0.5*R;
418       r[1][2] = -0.25 * a3 + 0.5 * R;          << 341       r[1][2] = -0.25*a3 + 0.5*R;
419       r[2][1] = 0.5 * D;                       << 342       r[2][1] =  0.5*D;
420       r[2][2] = -0.5 * D;                      << 343       r[2][2] = -0.5*D;
421     }                                             344     }
422     if(E2 >= 0.)                               << 345     if( E2 >= 0. )
423     {                                             346     {
424       E       = std::sqrt(E2);                    347       E       = std::sqrt(E2);
425       r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * E << 348       r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
426       r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * E << 349       r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
427       r[2][3] = 0.;                               350       r[2][3] = 0.;
428       r[2][4] = 0.;                               351       r[2][4] = 0.;
429     }                                             352     }
430     else                                          353     else
431     {                                             354     {
432       E       = std::sqrt(-E2);                   355       E       = std::sqrt(-E2);
433       r[1][3] = -0.25 * a3 - 0.5 * R;          << 356       r[1][3] = -0.25*a3 - 0.5*R;
434       r[1][4] = -0.25 * a3 - 0.5 * R;          << 357       r[1][4] = -0.25*a3 - 0.5*R;
435       r[2][3] = 0.5 * E;                       << 358       r[2][3] =  0.5*E;
436       r[2][4] = -0.5 * E;                      << 359       r[2][4] = -0.5*E;
437     }                                             360     }
438   }                                               361   }
439   else if(R2 < 0.)                             << 362   else if( R2 < 0.)
440   {                                               363   {
441     R = std::sqrt(-R2);                           364     R = std::sqrt(-R2);
442     G4complex CD2(a, -b / R);                  << 365     G4complex CD2(a,-b/R);
443     G4complex CD = std::sqrt(CD2);                366     G4complex CD = std::sqrt(CD2);
444                                                   367 
445     r[1][1] = -0.25 * a3 + 0.5 * real(CD);     << 368     r[1][1] = -0.25*a3 + 0.5*real(CD);
446     r[1][2] = -0.25 * a3 - 0.5 * real(CD);     << 369     r[1][2] = -0.25*a3 - 0.5*real(CD);
447     r[2][1] = 0.5 * R + 0.5 * imag(CD);        << 370     r[2][1] =  0.5*R + 0.5*imag(CD);
448     r[2][2] = 0.5 * R - 0.5 * imag(CD);        << 371     r[2][2] =  0.5*R - 0.5*imag(CD);
449     G4complex CE2(a, b / R);                   << 372     G4complex CE2(a,b/R);
450     G4complex CE = std::sqrt(CE2);                373     G4complex CE = std::sqrt(CE2);
451                                                   374 
452     r[1][3] = -0.25 * a3 + 0.5 * real(CE);     << 375     r[1][3] = -0.25*a3 + 0.5*real(CE);
453     r[1][4] = -0.25 * a3 - 0.5 * real(CE);     << 376     r[1][4] = -0.25*a3 - 0.5*real(CE);
454     r[2][3] = -0.5 * R + 0.5 * imag(CE);       << 377     r[2][3] =  -0.5*R + 0.5*imag(CE);
455     r[2][4] = -0.5 * R - 0.5 * imag(CE);       << 378     r[2][4] =  -0.5*R - 0.5*imag(CE);
456   }                                               379   }
457   else  // R2=0 case                           << 380   else // R2=0 case
458   {                                               381   {
459     if(d >= 0.)                                   382     if(d >= 0.)
460     {                                             383     {
461       D2 = c + std::sqrt(d);                      384       D2 = c + std::sqrt(d);
462       E2 = c - std::sqrt(d);                      385       E2 = c - std::sqrt(d);
463                                                   386 
464       if(D2 >= 0.)                             << 387       if( D2 >= 0. )
465       {                                           388       {
466         D       = std::sqrt(D2);                  389         D       = std::sqrt(D2);
467         r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * << 390         r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
468         r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * << 391         r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
469         r[2][1] = 0.;                             392         r[2][1] = 0.;
470         r[2][2] = 0.;                             393         r[2][2] = 0.;
471       }                                           394       }
472       else                                        395       else
473       {                                           396       {
474         D       = std::sqrt(-D2);                 397         D       = std::sqrt(-D2);
475         r[1][1] = -0.25 * a3 + 0.5 * R;        << 398         r[1][1] = -0.25*a3 + 0.5*R;
476         r[1][2] = -0.25 * a3 + 0.5 * R;        << 399         r[1][2] = -0.25*a3 + 0.5*R;
477         r[2][1] = 0.5 * D;                     << 400         r[2][1] =  0.5*D;
478         r[2][2] = -0.5 * D;                    << 401         r[2][2] = -0.5*D;
479       }                                           402       }
480       if(E2 >= 0.)                             << 403       if( E2 >= 0. )
481       {                                           404       {
482         E       = std::sqrt(E2);                  405         E       = std::sqrt(E2);
483         r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * << 406         r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
484         r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * << 407         r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
485         r[2][3] = 0.;                             408         r[2][3] = 0.;
486         r[2][4] = 0.;                             409         r[2][4] = 0.;
487       }                                           410       }
488       else                                        411       else
489       {                                           412       {
490         E       = std::sqrt(-E2);                 413         E       = std::sqrt(-E2);
491         r[1][3] = -0.25 * a3 - 0.5 * R;        << 414         r[1][3] = -0.25*a3 - 0.5*R;
492         r[1][4] = -0.25 * a3 - 0.5 * R;        << 415         r[1][4] = -0.25*a3 - 0.5*R;
493         r[2][3] = 0.5 * E;                     << 416         r[2][3] =  0.5*E;
494         r[2][4] = -0.5 * E;                    << 417         r[2][4] = -0.5*E;
495       }                                           418       }
496     }                                             419     }
497     else                                          420     else
498     {                                             421     {
499       ds = std::sqrt(-d);                         422       ds = std::sqrt(-d);
500       G4complex CD2(c, ds);                    << 423       G4complex CD2(c,ds);
501       G4complex CD = std::sqrt(CD2);              424       G4complex CD = std::sqrt(CD2);
502                                                   425 
503       r[1][1] = -0.25 * a3 + 0.5 * real(CD);   << 426       r[1][1] = -0.25*a3 + 0.5*real(CD);
504       r[1][2] = -0.25 * a3 - 0.5 * real(CD);   << 427       r[1][2] = -0.25*a3 - 0.5*real(CD);
505       r[2][1] = 0.5 * R + 0.5 * imag(CD);      << 428       r[2][1] =  0.5*R + 0.5*imag(CD);
506       r[2][2] = 0.5 * R - 0.5 * imag(CD);      << 429       r[2][2] =  0.5*R - 0.5*imag(CD);
507                                                   430 
508       G4complex CE2(c, -ds);                   << 431       G4complex CE2(c,-ds);
509       G4complex CE = std::sqrt(CE2);              432       G4complex CE = std::sqrt(CE2);
510                                                   433 
511       r[1][3] = -0.25 * a3 + 0.5 * real(CE);   << 434       r[1][3] = -0.25*a3 + 0.5*real(CE);
512       r[1][4] = -0.25 * a3 - 0.5 * real(CE);   << 435       r[1][4] = -0.25*a3 - 0.5*real(CE);
513       r[2][3] = -0.5 * R + 0.5 * imag(CE);     << 436       r[2][3] =  -0.5*R + 0.5*imag(CE);
514       r[2][4] = -0.5 * R - 0.5 * imag(CE);     << 437       r[2][4] =  -0.5*R - 0.5*imag(CE);
515     }                                          << 438     }  
516   }                                               439   }
517   return 4;                                       440   return 4;
518 }                                                 441 }
519                                                   442 
520 //                                                443 //
521 //                                                444 //
522 //////////////////////////////////////////////    445 //////////////////////////////////////////////////////////////////////////////
523                                                   446