Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4AnalyticalPolSolver Class Reference

#include <G4AnalyticalPolSolver.hh>

Public Member Functions

 G4AnalyticalPolSolver ()
 
 ~G4AnalyticalPolSolver ()
 
G4int QuadRoots (G4double p[5], G4double r[3][5])
 
G4int CubicRoots (G4double p[5], G4double r[3][5])
 
G4int BiquadRoots (G4double p[5], G4double r[3][5])
 
G4int QuarticRoots (G4double p[5], G4double r[3][5])
 

Detailed Description

Definition at line 64 of file G4AnalyticalPolSolver.hh.

Constructor & Destructor Documentation

◆ G4AnalyticalPolSolver()

G4AnalyticalPolSolver::G4AnalyticalPolSolver ( )

Definition at line 37 of file G4AnalyticalPolSolver.cc.

37{;}

◆ ~G4AnalyticalPolSolver()

G4AnalyticalPolSolver::~G4AnalyticalPolSolver ( )

Definition at line 41 of file G4AnalyticalPolSolver.cc.

41{;}

Member Function Documentation

◆ BiquadRoots()

G4int G4AnalyticalPolSolver::BiquadRoots ( G4double  p[5],
G4double  r[3][5] 
)

Definition at line 172 of file G4AnalyticalPolSolver.cc.

173{
174 G4double a, b, c, d, e;
175 G4int i, k, j;
176
177 if(p[0] != 1.0)
178 {
179 for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
180 p[0] = 1.;
181 }
182 e = 0.25*p[1];
183 b = 2*e;
184 c = b*b;
185 d = 0.75*c;
186 b = p[3] + b*( c - p[2] );
187 a = p[2] - d;
188 c = p[4] + e*( e*a - p[3] );
189 a = a - d;
190
191 p[1] = 0.5*a;
192 p[2] = (p[1]*p[1]-c)*0.25;
193 p[3] = b*b/(-64.0);
194
195 if( p[3] < 0. )
196 {
197 CubicRoots(p,r);
198
199 for( k = 1; k < 4; k++ )
200 {
201 if( r[2][k] == 0. && r[1][k] > 0 )
202 {
203 d = r[1][k]*4;
204 a = a + d;
205
206 if ( a >= 0. && b >= 0.) { p[1] = std::sqrt(d); }
207 else if( a <= 0. && b <= 0.) { p[1] = std::sqrt(d); }
208 else { p[1] = -std::sqrt(d); }
209
210 b = 0.5*( a + b/p[1] );
211
212 p[2] = c/b;
213 QuadRoots(p,r);
214
215 for( i = 1; i < 3; i++ )
216 {
217 for( j = 1; j < 3; j++ ) { r[j][i+2] = r[j][i]; }
218 }
219 p[1] = -p[1];
220 p[2] = b;
221 QuadRoots(p,r);
222
223 for( i = 1; i < 5; i++ ) { r[1][i] = r[1][i] - e; }
224
225 return 4;
226 }
227 }
228 }
229 if( p[2] < 0. )
230 {
231 b = std::sqrt(c);
232 d = b + b - a;
233 p[1] = 0.;
234
235 if( d > 0. ) { p[1] = std::sqrt(d); }
236 }
237 else
238 {
239 if( p[1] > 0.) { b = std::sqrt(p[2])*2.0 + p[1]; }
240 else { b = -std::sqrt(p[2])*2.0 + p[1]; }
241
242 if( b != 0.) { p[1] = 0; }
243 else
244 {
245 for(k = 1; k < 5; k++ )
246 {
247 r[1][k] = -e;
248 r[2][k] = 0;
249 }
250 return 0;
251 }
252 }
253
254 p[2] = c/b;
255 QuadRoots(p,r);
256
257 for( k = 1; k < 3; k++ )
258 {
259 for( j = 1; j < 3; j++ ) { r[j][k+2] = r[j][k]; }
260 }
261 p[1] = -p[1];
262 p[2] = b;
263 QuadRoots(p,r);
264
265 for( k = 1; k < 5; k++ ) { r[1][k] = r[1][k] - e; }
266
267 return 4;
268}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4AnalyticalPolSolver::CubicRoots
G4int CubicRoots(G4double p[5], G4double r[3][5])
Definition: G4AnalyticalPolSolver.cc:85
G4AnalyticalPolSolver::QuadRoots
G4int QuadRoots(G4double p[5], G4double r[3][5])
Definition: G4AnalyticalPolSolver.cc:50

◆ CubicRoots()

G4int G4AnalyticalPolSolver::CubicRoots ( G4double  p[5],
G4double  r[3][5] 
)

Definition at line 85 of file G4AnalyticalPolSolver.cc.

86{
87 G4double x,t,b,c,d;
88 G4int k;
89
90 if( p[0] != 1. )
91 {
92 for(k = 1; k < 4; k++ ) { p[k] = p[k]/p[0]; }
93 p[0] = 1.;
94 }
95 x = p[1]/3.0;
96 t = x*p[1];
97 b = 0.5*( x*( t/1.5 - p[2] ) + p[3] );
98 t = ( t - p[2] )/3.0;
99 c = t*t*t;
100 d = b*b - c;
101
102 if( d >= 0. )
103 {
104 d = std::pow( (std::sqrt(d) + std::fabs(b) ), 1.0/3.0 );
105
106 if( d != 0. )
107 {
108 if( b > 0. ) { b = -d; }
109 else { b = d; }
110 c = t/b;
111 }
112 d = std::sqrt(0.75)*(b - c);
113 r[2][2] = d;
114 b = b + c;
115 c = -0.5*b-x;
116 r[1][2] = c;
117
118 if( ( b > 0. && x <= 0. ) || ( b < 0. && x > 0. ) )
119 {
120 r[1][1] = c;
121 r[2][1] = -d;
122 r[1][3] = b - x;
123 r[2][3] = 0;
124 }
125 else
126 {
127 r[1][1] = b - x;
128 r[2][1] = 0.;
129 r[1][3] = c;
130 r[2][3] = -d;
131 }
132 } // end of 2 equal or complex roots
133 else
134 {
135 if( b == 0. ) { d = std::atan(1.0)/1.5; }
136 else { d = std::atan( std::sqrt(-d)/std::fabs(b) )/3.0; }
137
138 if( b < 0. ) { b = std::sqrt(t)*2.0; }
139 else { b = -2.0*std::sqrt(t); }
140
141 c = std::cos(d)*b;
142 t = -std::sqrt(0.75)*std::sin(d)*b - 0.5*c;
143 d = -t - c - x;
144 c = c - x;
145 t = t - x;
146
147 if( std::fabs(c) > std::fabs(t) ) { r[1][3] = c; }
148 else
149 {
150 r[1][3] = t;
151 t = c;
152 }
153 if( std::fabs(d) > std::fabs(t) ) { r[1][2] = d; }
154 else
155 {
156 r[1][2] = t;
157 t = d;
158 }
159 r[1][1] = t;
160
161 for(k = 1; k < 4; k++ ) { r[2][k] = 0.; }
162 }
163 return 0;
164}

Referenced by BiquadRoots(), and QuarticRoots().

◆ QuadRoots()

G4int G4AnalyticalPolSolver::QuadRoots ( G4double  p[5],
G4double  r[3][5] 
)

Definition at line 50 of file G4AnalyticalPolSolver.cc.

51{
52 G4double b, c, d2, d;
53
54 b = -p[1]/p[0]/2.;
55 c = p[2]/p[0];
56 d2 = b*b - c;
57
58 if( d2 >= 0. )
59 {
60 d = std::sqrt(d2);
61 r[1][1] = b - d;
62 r[1][2] = b + d;
63 r[2][1] = 0.;
64 r[2][2] = 0.;
65 }
66 else
67 {
68 d = std::sqrt(-d2);
69 r[2][1] = d;
70 r[2][2] = -d;
71 r[1][1] = b;
72 r[1][2] = b;
73 }
74
75 return 2;
76}

Referenced by BiquadRoots().

◆ QuarticRoots()

G4int G4AnalyticalPolSolver::QuarticRoots ( G4double  p[5],
G4double  r[3][5] 
)

Definition at line 272 of file G4AnalyticalPolSolver.cc.

273{
274 G4double a0, a1, a2, a3, y1;
275 G4double R2, D2, E2, D, E, R = 0.;
276 G4double a, b, c, d, ds;
277
278 G4double reRoot[4];
279 G4int k, noReRoots = 0;
280
281 for( k = 0; k < 4; k++ ) { reRoot[k] = DBL_MAX; }
282
283 if( p[0] != 1.0 )
284 {
285 for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
286 p[0] = 1.;
287 }
288 a3 = p[1];
289 a2 = p[2];
290 a1 = p[3];
291 a0 = p[4];
292
293 // resolvent cubic equation cofs:
294
295 p[1] = -a2;
296 p[2] = a1*a3 - 4*a0;
297 p[3] = 4*a2*a0 - a1*a1 - a3*a3*a0;
298
299 CubicRoots(p,r);
300
301 for( k = 1; k < 4; k++ )
302 {
303 if( r[2][k] == 0. ) // find a real root
304 {
305 noReRoots++;
306 reRoot[k] = r[1][k];
307 }
308 else reRoot[k] = DBL_MAX; // kInfinity;
309 }
310 y1 = DBL_MAX; // kInfinity;
311 for( k = 1; k < 4; k++ )
312 {
313 if ( reRoot[k] < y1 ) { y1 = reRoot[k]; }
314 }
315
316 R2 = 0.25*a3*a3 - a2 + y1;
317 b = 0.25*(4*a3*a2 - 8*a1 - a3*a3*a3);
318 c = 0.75*a3*a3 - 2*a2;
319 a = c - R2;
320 d = 4*y1*y1 - 16*a0;
321
322 if( R2 > 0.)
323 {
324 R = std::sqrt(R2);
325 D2 = a + b/R;
326 E2 = a - b/R;
327
328 if( D2 >= 0. )
329 {
330 D = std::sqrt(D2);
331 r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
332 r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
333 r[2][1] = 0.;
334 r[2][2] = 0.;
335 }
336 else
337 {
338 D = std::sqrt(-D2);
339 r[1][1] = -0.25*a3 + 0.5*R;
340 r[1][2] = -0.25*a3 + 0.5*R;
341 r[2][1] = 0.5*D;
342 r[2][2] = -0.5*D;
343 }
344 if( E2 >= 0. )
345 {
346 E = std::sqrt(E2);
347 r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
348 r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
349 r[2][3] = 0.;
350 r[2][4] = 0.;
351 }
352 else
353 {
354 E = std::sqrt(-E2);
355 r[1][3] = -0.25*a3 - 0.5*R;
356 r[1][4] = -0.25*a3 - 0.5*R;
357 r[2][3] = 0.5*E;
358 r[2][4] = -0.5*E;
359 }
360 }
361 else if( R2 < 0.)
362 {
363 R = std::sqrt(-R2);
364 G4complex CD2(a,-b/R);
365 G4complex CD = std::sqrt(CD2);
366
367 r[1][1] = -0.25*a3 + 0.5*real(CD);
368 r[1][2] = -0.25*a3 - 0.5*real(CD);
369 r[2][1] = 0.5*R + 0.5*imag(CD);
370 r[2][2] = 0.5*R - 0.5*imag(CD);
371 G4complex CE2(a,b/R);
372 G4complex CE = std::sqrt(CE2);
373
374 r[1][3] = -0.25*a3 + 0.5*real(CE);
375 r[1][4] = -0.25*a3 - 0.5*real(CE);
376 r[2][3] = -0.5*R + 0.5*imag(CE);
377 r[2][4] = -0.5*R - 0.5*imag(CE);
378 }
379 else // R2=0 case
380 {
381 if(d >= 0.)
382 {
383 D2 = c + std::sqrt(d);
384 E2 = c - std::sqrt(d);
385
386 if( D2 >= 0. )
387 {
388 D = std::sqrt(D2);
389 r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
390 r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
391 r[2][1] = 0.;
392 r[2][2] = 0.;
393 }
394 else
395 {
396 D = std::sqrt(-D2);
397 r[1][1] = -0.25*a3 + 0.5*R;
398 r[1][2] = -0.25*a3 + 0.5*R;
399 r[2][1] = 0.5*D;
400 r[2][2] = -0.5*D;
401 }
402 if( E2 >= 0. )
403 {
404 E = std::sqrt(E2);
405 r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
406 r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
407 r[2][3] = 0.;
408 r[2][4] = 0.;
409 }
410 else
411 {
412 E = std::sqrt(-E2);
413 r[1][3] = -0.25*a3 - 0.5*R;
414 r[1][4] = -0.25*a3 - 0.5*R;
415 r[2][3] = 0.5*E;
416 r[2][4] = -0.5*E;
417 }
418 }
419 else
420 {
421 ds = std::sqrt(-d);
422 G4complex CD2(c,ds);
423 G4complex CD = std::sqrt(CD2);
424
425 r[1][1] = -0.25*a3 + 0.5*real(CD);
426 r[1][2] = -0.25*a3 - 0.5*real(CD);
427 r[2][1] = 0.5*R + 0.5*imag(CD);
428 r[2][2] = 0.5*R - 0.5*imag(CD);
429
430 G4complex CE2(c,-ds);
431 G4complex CE = std::sqrt(CE2);
432
433 r[1][3] = -0.25*a3 + 0.5*real(CE);
434 r[1][4] = -0.25*a3 - 0.5*real(CE);
435 r[2][3] = -0.5*R + 0.5*imag(CE);
436 r[2][4] = -0.5*R - 0.5*imag(CE);
437 }
438 }
439 return 4;
440}
G4complex
std::complex< G4double > G4complex
Definition: G4Types.hh:69
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
The documentation for this class was generated from the following files: