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 #include "globals.hh"
00031 #include <complex>
00032
00033 #include "G4AnalyticalPolSolver.hh"
00034
00036
00037 G4AnalyticalPolSolver::G4AnalyticalPolSolver() {;}
00038
00040
00041 G4AnalyticalPolSolver::~G4AnalyticalPolSolver() {;}
00042
00044
00045
00046
00047
00048
00049
00050 G4int G4AnalyticalPolSolver::QuadRoots( G4double p[5], G4double r[3][5] )
00051 {
00052 G4double b, c, d2, d;
00053
00054 b = -p[1]/p[0]/2.;
00055 c = p[2]/p[0];
00056 d2 = b*b - c;
00057
00058 if( d2 >= 0. )
00059 {
00060 d = std::sqrt(d2);
00061 r[1][1] = b - d;
00062 r[1][2] = b + d;
00063 r[2][1] = 0.;
00064 r[2][2] = 0.;
00065 }
00066 else
00067 {
00068 d = std::sqrt(-d2);
00069 r[2][1] = d;
00070 r[2][2] = -d;
00071 r[1][1] = b;
00072 r[1][2] = b;
00073 }
00074
00075 return 2;
00076 }
00077
00079
00080
00081
00082
00083
00084
00085 G4int G4AnalyticalPolSolver::CubicRoots( G4double p[5], G4double r[3][5] )
00086 {
00087 G4double x,t,b,c,d;
00088 G4int k;
00089
00090 if( p[0] != 1. )
00091 {
00092 for(k = 1; k < 4; k++ ) { p[k] = p[k]/p[0]; }
00093 p[0] = 1.;
00094 }
00095 x = p[1]/3.0;
00096 t = x*p[1];
00097 b = 0.5*( x*( t/1.5 - p[2] ) + p[3] );
00098 t = ( t - p[2] )/3.0;
00099 c = t*t*t;
00100 d = b*b - c;
00101
00102 if( d >= 0. )
00103 {
00104 d = std::pow( (std::sqrt(d) + std::fabs(b) ), 1.0/3.0 );
00105
00106 if( d != 0. )
00107 {
00108 if( b > 0. ) { b = -d; }
00109 else { b = d; }
00110 c = t/b;
00111 }
00112 d = std::sqrt(0.75)*(b - c);
00113 r[2][2] = d;
00114 b = b + c;
00115 c = -0.5*b-x;
00116 r[1][2] = c;
00117
00118 if( ( b > 0. && x <= 0. ) || ( b < 0. && x > 0. ) )
00119 {
00120 r[1][1] = c;
00121 r[2][1] = -d;
00122 r[1][3] = b - x;
00123 r[2][3] = 0;
00124 }
00125 else
00126 {
00127 r[1][1] = b - x;
00128 r[2][1] = 0.;
00129 r[1][3] = c;
00130 r[2][3] = -d;
00131 }
00132 }
00133 else
00134 {
00135 if( b == 0. ) { d = std::atan(1.0)/1.5; }
00136 else { d = std::atan( std::sqrt(-d)/std::fabs(b) )/3.0; }
00137
00138 if( b < 0. ) { b = std::sqrt(t)*2.0; }
00139 else { b = -2.0*std::sqrt(t); }
00140
00141 c = std::cos(d)*b;
00142 t = -std::sqrt(0.75)*std::sin(d)*b - 0.5*c;
00143 d = -t - c - x;
00144 c = c - x;
00145 t = t - x;
00146
00147 if( std::fabs(c) > std::fabs(t) ) { r[1][3] = c; }
00148 else
00149 {
00150 r[1][3] = t;
00151 t = c;
00152 }
00153 if( std::fabs(d) > std::fabs(t) ) { r[1][2] = d; }
00154 else
00155 {
00156 r[1][2] = t;
00157 t = d;
00158 }
00159 r[1][1] = t;
00160
00161 for(k = 1; k < 4; k++ ) { r[2][k] = 0.; }
00162 }
00163 return 0;
00164 }
00165
00167
00168
00169
00170
00171
00172 G4int G4AnalyticalPolSolver::BiquadRoots( G4double p[5], G4double r[3][5] )
00173 {
00174 G4double a, b, c, d, e;
00175 G4int i, k, j;
00176
00177 if(p[0] != 1.0)
00178 {
00179 for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
00180 p[0] = 1.;
00181 }
00182 e = 0.25*p[1];
00183 b = 2*e;
00184 c = b*b;
00185 d = 0.75*c;
00186 b = p[3] + b*( c - p[2] );
00187 a = p[2] - d;
00188 c = p[4] + e*( e*a - p[3] );
00189 a = a - d;
00190
00191 p[1] = 0.5*a;
00192 p[2] = (p[1]*p[1]-c)*0.25;
00193 p[3] = b*b/(-64.0);
00194
00195 if( p[3] < 0. )
00196 {
00197 CubicRoots(p,r);
00198
00199 for( k = 1; k < 4; k++ )
00200 {
00201 if( r[2][k] == 0. && r[1][k] > 0 )
00202 {
00203 d = r[1][k]*4;
00204 a = a + d;
00205
00206 if ( a >= 0. && b >= 0.) { p[1] = std::sqrt(d); }
00207 else if( a <= 0. && b <= 0.) { p[1] = std::sqrt(d); }
00208 else { p[1] = -std::sqrt(d); }
00209
00210 b = 0.5*( a + b/p[1] );
00211
00212 p[2] = c/b;
00213 QuadRoots(p,r);
00214
00215 for( i = 1; i < 3; i++ )
00216 {
00217 for( j = 1; j < 3; j++ ) { r[j][i+2] = r[j][i]; }
00218 }
00219 p[1] = -p[1];
00220 p[2] = b;
00221 QuadRoots(p,r);
00222
00223 for( i = 1; i < 5; i++ ) { r[1][i] = r[1][i] - e; }
00224
00225 return 4;
00226 }
00227 }
00228 }
00229 if( p[2] < 0. )
00230 {
00231 b = std::sqrt(c);
00232 d = b + b - a;
00233 p[1] = 0.;
00234
00235 if( d > 0. ) { p[1] = std::sqrt(d); }
00236 }
00237 else
00238 {
00239 if( p[1] > 0.) { b = std::sqrt(p[2])*2.0 + p[1]; }
00240 else { b = -std::sqrt(p[2])*2.0 + p[1]; }
00241
00242 if( b != 0.) { p[1] = 0; }
00243 else
00244 {
00245 for(k = 1; k < 5; k++ )
00246 {
00247 r[1][k] = -e;
00248 r[2][k] = 0;
00249 }
00250 return 0;
00251 }
00252 }
00253
00254 p[2] = c/b;
00255 QuadRoots(p,r);
00256
00257 for( k = 1; k < 3; k++ )
00258 {
00259 for( j = 1; j < 3; j++ ) { r[j][k+2] = r[j][k]; }
00260 }
00261 p[1] = -p[1];
00262 p[2] = b;
00263 QuadRoots(p,r);
00264
00265 for( k = 1; k < 5; k++ ) { r[1][k] = r[1][k] - e; }
00266
00267 return 4;
00268 }
00269
00271
00272 G4int G4AnalyticalPolSolver::QuarticRoots( G4double p[5], G4double r[3][5])
00273 {
00274 G4double a0, a1, a2, a3, y1;
00275 G4double R2, D2, E2, D, E, R = 0.;
00276 G4double a, b, c, d, ds;
00277
00278 G4double reRoot[4];
00279 G4int k, noReRoots = 0;
00280
00281 for( k = 0; k < 4; k++ ) { reRoot[k] = DBL_MAX; }
00282
00283 if( p[0] != 1.0 )
00284 {
00285 for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
00286 p[0] = 1.;
00287 }
00288 a3 = p[1];
00289 a2 = p[2];
00290 a1 = p[3];
00291 a0 = p[4];
00292
00293
00294
00295 p[1] = -a2;
00296 p[2] = a1*a3 - 4*a0;
00297 p[3] = 4*a2*a0 - a1*a1 - a3*a3*a0;
00298
00299 CubicRoots(p,r);
00300
00301 for( k = 1; k < 4; k++ )
00302 {
00303 if( r[2][k] == 0. )
00304 {
00305 noReRoots++;
00306 reRoot[k] = r[1][k];
00307 }
00308 else reRoot[k] = DBL_MAX;
00309 }
00310 y1 = DBL_MAX;
00311 for( k = 1; k < 4; k++ )
00312 {
00313 if ( reRoot[k] < y1 ) { y1 = reRoot[k]; }
00314 }
00315
00316 R2 = 0.25*a3*a3 - a2 + y1;
00317 b = 0.25*(4*a3*a2 - 8*a1 - a3*a3*a3);
00318 c = 0.75*a3*a3 - 2*a2;
00319 a = c - R2;
00320 d = 4*y1*y1 - 16*a0;
00321
00322 if( R2 > 0.)
00323 {
00324 R = std::sqrt(R2);
00325 D2 = a + b/R;
00326 E2 = a - b/R;
00327
00328 if( D2 >= 0. )
00329 {
00330 D = std::sqrt(D2);
00331 r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
00332 r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
00333 r[2][1] = 0.;
00334 r[2][2] = 0.;
00335 }
00336 else
00337 {
00338 D = std::sqrt(-D2);
00339 r[1][1] = -0.25*a3 + 0.5*R;
00340 r[1][2] = -0.25*a3 + 0.5*R;
00341 r[2][1] = 0.5*D;
00342 r[2][2] = -0.5*D;
00343 }
00344 if( E2 >= 0. )
00345 {
00346 E = std::sqrt(E2);
00347 r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
00348 r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
00349 r[2][3] = 0.;
00350 r[2][4] = 0.;
00351 }
00352 else
00353 {
00354 E = std::sqrt(-E2);
00355 r[1][3] = -0.25*a3 - 0.5*R;
00356 r[1][4] = -0.25*a3 - 0.5*R;
00357 r[2][3] = 0.5*E;
00358 r[2][4] = -0.5*E;
00359 }
00360 }
00361 else if( R2 < 0.)
00362 {
00363 R = std::sqrt(-R2);
00364 G4complex CD2(a,-b/R);
00365 G4complex CD = std::sqrt(CD2);
00366
00367 r[1][1] = -0.25*a3 + 0.5*real(CD);
00368 r[1][2] = -0.25*a3 - 0.5*real(CD);
00369 r[2][1] = 0.5*R + 0.5*imag(CD);
00370 r[2][2] = 0.5*R - 0.5*imag(CD);
00371 G4complex CE2(a,b/R);
00372 G4complex CE = std::sqrt(CE2);
00373
00374 r[1][3] = -0.25*a3 + 0.5*real(CE);
00375 r[1][4] = -0.25*a3 - 0.5*real(CE);
00376 r[2][3] = -0.5*R + 0.5*imag(CE);
00377 r[2][4] = -0.5*R - 0.5*imag(CE);
00378 }
00379 else
00380 {
00381 if(d >= 0.)
00382 {
00383 D2 = c + std::sqrt(d);
00384 E2 = c - std::sqrt(d);
00385
00386 if( D2 >= 0. )
00387 {
00388 D = std::sqrt(D2);
00389 r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
00390 r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
00391 r[2][1] = 0.;
00392 r[2][2] = 0.;
00393 }
00394 else
00395 {
00396 D = std::sqrt(-D2);
00397 r[1][1] = -0.25*a3 + 0.5*R;
00398 r[1][2] = -0.25*a3 + 0.5*R;
00399 r[2][1] = 0.5*D;
00400 r[2][2] = -0.5*D;
00401 }
00402 if( E2 >= 0. )
00403 {
00404 E = std::sqrt(E2);
00405 r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
00406 r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
00407 r[2][3] = 0.;
00408 r[2][4] = 0.;
00409 }
00410 else
00411 {
00412 E = std::sqrt(-E2);
00413 r[1][3] = -0.25*a3 - 0.5*R;
00414 r[1][4] = -0.25*a3 - 0.5*R;
00415 r[2][3] = 0.5*E;
00416 r[2][4] = -0.5*E;
00417 }
00418 }
00419 else
00420 {
00421 ds = std::sqrt(-d);
00422 G4complex CD2(c,ds);
00423 G4complex CD = std::sqrt(CD2);
00424
00425 r[1][1] = -0.25*a3 + 0.5*real(CD);
00426 r[1][2] = -0.25*a3 - 0.5*real(CD);
00427 r[2][1] = 0.5*R + 0.5*imag(CD);
00428 r[2][2] = 0.5*R - 0.5*imag(CD);
00429
00430 G4complex CE2(c,-ds);
00431 G4complex CE = std::sqrt(CE2);
00432
00433 r[1][3] = -0.25*a3 + 0.5*real(CE);
00434 r[1][4] = -0.25*a3 - 0.5*real(CE);
00435 r[2][3] = -0.5*R + 0.5*imag(CE);
00436 r[2][4] = -0.5*R - 0.5*imag(CE);
00437 }
00438 }
00439 return 4;
00440 }
00441
00442
00443