Geant4 10.7.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 60 of file G4AnalyticalPolSolver.hh.

Constructor & Destructor Documentation

◆ G4AnalyticalPolSolver()

G4AnalyticalPolSolver::G4AnalyticalPolSolver ( )

Definition at line 38 of file G4AnalyticalPolSolver.cc.

38{ ; }

◆ ~G4AnalyticalPolSolver()

G4AnalyticalPolSolver::~G4AnalyticalPolSolver ( )

Definition at line 42 of file G4AnalyticalPolSolver.cc.

42{ ; }

Member Function Documentation

◆ BiquadRoots()

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

Definition at line 203 of file G4AnalyticalPolSolver.cc.

204{
205 G4double a, b, c, d, e;
206 G4int i, k, j;
207
208 if(p[0] != 1.0)
209 {
210 for(k = 1; k < 5; k++)
211 {
212 p[k] = p[k] / p[0];
213 }
214 p[0] = 1.;
215 }
216 e = 0.25 * p[1];
217 b = 2 * e;
218 c = b * b;
219 d = 0.75 * c;
220 b = p[3] + b * (c - p[2]);
221 a = p[2] - d;
222 c = p[4] + e * (e * a - p[3]);
223 a = a - d;
224
225 p[1] = 0.5 * a;
226 p[2] = (p[1] * p[1] - c) * 0.25;
227 p[3] = b * b / (-64.0);
228
229 if(p[3] < 0.)
230 {
231 CubicRoots(p, r);
232
233 for(k = 1; k < 4; k++)
234 {
235 if(r[2][k] == 0. && r[1][k] > 0)
236 {
237 d = r[1][k] * 4;
238 a = a + d;
239
240 if(a >= 0. && b >= 0.)
241 {
242 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
253 b = 0.5 * (a + b / p[1]);
254
255 p[2] = c / b;
256 QuadRoots(p, r);
257
258 for(i = 1; i < 3; i++)
259 {
260 for(j = 1; j < 3; j++)
261 {
262 r[j][i + 2] = r[j][i];
263 }
264 }
265 p[1] = -p[1];
266 p[2] = b;
267 QuadRoots(p, r);
268
269 for(i = 1; i < 5; i++)
270 {
271 r[1][i] = r[1][i] - e;
272 }
273
274 return 4;
275 }
276 }
277 }
278 if(p[2] < 0.)
279 {
280 b = std::sqrt(c);
281 d = b + b - a;
282 p[1] = 0.;
283
284 if(d > 0.)
285 {
286 p[1] = std::sqrt(d);
287 }
288 }
289 else
290 {
291 if(p[1] > 0.)
292 {
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
300 if(b != 0.)
301 {
302 p[1] = 0;
303 }
304 else
305 {
306 for(k = 1; k < 5; k++)
307 {
308 r[1][k] = -e;
309 r[2][k] = 0;
310 }
311 return 0;
312 }
313 }
314
315 p[2] = c / b;
316 QuadRoots(p, r);
317
318 for(k = 1; k < 3; k++)
319 {
320 for(j = 1; j < 3; j++)
321 {
322 r[j][k + 2] = r[j][k];
323 }
324 }
325 p[1] = -p[1];
326 p[2] = b;
327 QuadRoots(p, r);
328
329 for(k = 1; k < 5; k++)
330 {
331 r[1][k] = r[1][k] - e;
332 }
333
334 return 4;
335}
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4AnalyticalPolSolver::CubicRoots
G4int CubicRoots(G4double p[5], G4double r[3][5])
Definition: G4AnalyticalPolSolver.cc:86
G4AnalyticalPolSolver::QuadRoots
G4int QuadRoots(G4double p[5], G4double r[3][5])
Definition: G4AnalyticalPolSolver.cc:51

◆ CubicRoots()

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

Definition at line 86 of file G4AnalyticalPolSolver.cc.

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

Referenced by BiquadRoots(), and QuarticRoots().

◆ QuadRoots()

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

Definition at line 51 of file G4AnalyticalPolSolver.cc.

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

Referenced by BiquadRoots().

◆ QuarticRoots()

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

Definition at line 339 of file G4AnalyticalPolSolver.cc.

340{
341 G4double a0, a1, a2, a3, y1;
342 G4double R2, D2, E2, D, E, R = 0.;
343 G4double a, b, c, d, ds;
344
345 G4double reRoot[4];
346 G4int k, noReRoots = 0;
347
348 for(k = 0; k < 4; k++)
349 {
350 reRoot[k] = DBL_MAX;
351 }
352
353 if(p[0] != 1.0)
354 {
355 for(k = 1; k < 5; k++)
356 {
357 p[k] = p[k] / p[0];
358 }
359 p[0] = 1.;
360 }
361 a3 = p[1];
362 a2 = p[2];
363 a1 = p[3];
364 a0 = p[4];
365
366 // resolvent cubic equation cofs:
367
368 p[1] = -a2;
369 p[2] = a1 * a3 - 4 * a0;
370 p[3] = 4 * a2 * a0 - a1 * a1 - a3 * a3 * a0;
371
372 CubicRoots(p, r);
373
374 for(k = 1; k < 4; k++)
375 {
376 if(r[2][k] == 0.) // find a real root
377 {
378 noReRoots++;
379 reRoot[k] = r[1][k];
380 }
381 else
382 reRoot[k] = DBL_MAX; // kInfinity;
383 }
384 y1 = DBL_MAX; // kInfinity;
385 for(k = 1; k < 4; k++)
386 {
387 if(reRoot[k] < y1)
388 {
389 y1 = reRoot[k];
390 }
391 }
392
393 R2 = 0.25 * a3 * a3 - a2 + y1;
394 b = 0.25 * (4 * a3 * a2 - 8 * a1 - a3 * a3 * a3);
395 c = 0.75 * a3 * a3 - 2 * a2;
396 a = c - R2;
397 d = 4 * y1 * y1 - 16 * a0;
398
399 if(R2 > 0.)
400 {
401 R = std::sqrt(R2);
402 D2 = a + b / R;
403 E2 = a - b / R;
404
405 if(D2 >= 0.)
406 {
407 D = std::sqrt(D2);
408 r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * D;
409 r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * D;
410 r[2][1] = 0.;
411 r[2][2] = 0.;
412 }
413 else
414 {
415 D = std::sqrt(-D2);
416 r[1][1] = -0.25 * a3 + 0.5 * R;
417 r[1][2] = -0.25 * a3 + 0.5 * R;
418 r[2][1] = 0.5 * D;
419 r[2][2] = -0.5 * D;
420 }
421 if(E2 >= 0.)
422 {
423 E = std::sqrt(E2);
424 r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * E;
425 r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * E;
426 r[2][3] = 0.;
427 r[2][4] = 0.;
428 }
429 else
430 {
431 E = std::sqrt(-E2);
432 r[1][3] = -0.25 * a3 - 0.5 * R;
433 r[1][4] = -0.25 * a3 - 0.5 * R;
434 r[2][3] = 0.5 * E;
435 r[2][4] = -0.5 * E;
436 }
437 }
438 else if(R2 < 0.)
439 {
440 R = std::sqrt(-R2);
441 G4complex CD2(a, -b / R);
442 G4complex CD = std::sqrt(CD2);
443
444 r[1][1] = -0.25 * a3 + 0.5 * real(CD);
445 r[1][2] = -0.25 * a3 - 0.5 * real(CD);
446 r[2][1] = 0.5 * R + 0.5 * imag(CD);
447 r[2][2] = 0.5 * R - 0.5 * imag(CD);
448 G4complex CE2(a, b / R);
449 G4complex CE = std::sqrt(CE2);
450
451 r[1][3] = -0.25 * a3 + 0.5 * real(CE);
452 r[1][4] = -0.25 * a3 - 0.5 * real(CE);
453 r[2][3] = -0.5 * R + 0.5 * imag(CE);
454 r[2][4] = -0.5 * R - 0.5 * imag(CE);
455 }
456 else // R2=0 case
457 {
458 if(d >= 0.)
459 {
460 D2 = c + std::sqrt(d);
461 E2 = c - std::sqrt(d);
462
463 if(D2 >= 0.)
464 {
465 D = std::sqrt(D2);
466 r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * D;
467 r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * D;
468 r[2][1] = 0.;
469 r[2][2] = 0.;
470 }
471 else
472 {
473 D = std::sqrt(-D2);
474 r[1][1] = -0.25 * a3 + 0.5 * R;
475 r[1][2] = -0.25 * a3 + 0.5 * R;
476 r[2][1] = 0.5 * D;
477 r[2][2] = -0.5 * D;
478 }
479 if(E2 >= 0.)
480 {
481 E = std::sqrt(E2);
482 r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * E;
483 r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * E;
484 r[2][3] = 0.;
485 r[2][4] = 0.;
486 }
487 else
488 {
489 E = std::sqrt(-E2);
490 r[1][3] = -0.25 * a3 - 0.5 * R;
491 r[1][4] = -0.25 * a3 - 0.5 * R;
492 r[2][3] = 0.5 * E;
493 r[2][4] = -0.5 * E;
494 }
495 }
496 else
497 {
498 ds = std::sqrt(-d);
499 G4complex CD2(c, ds);
500 G4complex CD = std::sqrt(CD2);
501
502 r[1][1] = -0.25 * a3 + 0.5 * real(CD);
503 r[1][2] = -0.25 * a3 - 0.5 * real(CD);
504 r[2][1] = 0.5 * R + 0.5 * imag(CD);
505 r[2][2] = 0.5 * R - 0.5 * imag(CD);
506
507 G4complex CE2(c, -ds);
508 G4complex CE = std::sqrt(CE2);
509
510 r[1][3] = -0.25 * a3 + 0.5 * real(CE);
511 r[1][4] = -0.25 * a3 - 0.5 * real(CE);
512 r[2][3] = -0.5 * R + 0.5 * imag(CE);
513 r[2][4] = -0.5 * R - 0.5 * imag(CE);
514 }
515 }
516 return 4;
517}
D
double D(double temp)
Definition: G4DNAElectronHoleRecombination.cc:91
a0
const G4double a0
Definition: G4IonCoulombCrossSection.cc:68
G4complex
std::complex< G4double > G4complex
Definition: G4Types.hh:88
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62
The documentation for this class was generated from the following files: