Geant4  10.01.p01
G4AnalyticalPolSolver.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4AnalyticalPolSolver.cc 67970 2013-03-13 10:10:06Z gcosmo $
28 //
29 
30 #include "globals.hh"
31 #include <complex>
32 
33 #include "G4AnalyticalPolSolver.hh"
34 
36 
38 
40 
42 
44 //
45 // Array r[3][5] p[5]
46 // Roots of poly p[0] x^2 + p[1] x+p[2]=0
47 //
48 // x = r[1][k] + i r[2][k]; k = 1, 2
49 
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 }
77 
79 //
80 // Array r[3][5] p[5]
81 // Roots of poly p[0] x^3 + p[1] x^2...+p[3]=0
82 // x=r[1][k] + i r[2][k] k=1,...,3
83 // Assumes 0<arctan(x)<pi/2 for x>0
84 
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 }
165 
167 //
168 // Array r[3][5] p[5]
169 // Roots of poly p[0] x^4 + p[1] x^3...+p[4]=0
170 // x=r[1][k] + i r[2][k] k=1,...,4
171 
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 }
269 
271 
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 }
441 
442 //
443 //
const G4double a0
static const G4double a1
G4double a
Definition: TRTMaterials.hh:39
G4int QuadRoots(G4double p[5], G4double r[3][5])
int G4int
Definition: G4Types.hh:78
std::complex< G4double > G4complex
Definition: G4Types.hh:81
static const G4double a3
G4int QuarticRoots(G4double p[5], G4double r[3][5])
double G4double
Definition: G4Types.hh:76
static const G4double d2
G4int BiquadRoots(G4double p[5], G4double r[3][5])
#define DBL_MAX
Definition: templates.hh:83
G4int CubicRoots(G4double p[5], G4double r[3][5])
static const G4double a2