ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AnalyticalPolSolver.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4AnalyticalPolSolver.cc
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 //
28 
29 #include "globals.hh"
30 #include <complex>
31 
32 #include "G4AnalyticalPolSolver.hh"
33 
35 
37 
39 
41 
43 //
44 // Array r[3][5] p[5]
45 // Roots of poly p[0] x^2 + p[1] x+p[2]=0
46 //
47 // x = r[1][k] + i r[2][k]; k = 1, 2
48 
50 {
51  G4double b, c, d2, d;
52 
53  b = -p[1]/p[0]/2.;
54  c = p[2]/p[0];
55  d2 = b*b - c;
56 
57  if( d2 >= 0. )
58  {
59  d = std::sqrt(d2);
60  r[1][1] = b - d;
61  r[1][2] = b + d;
62  r[2][1] = 0.;
63  r[2][2] = 0.;
64  }
65  else
66  {
67  d = std::sqrt(-d2);
68  r[2][1] = d;
69  r[2][2] = -d;
70  r[1][1] = b;
71  r[1][2] = b;
72  }
73 
74  return 2;
75 }
76 
78 //
79 // Array r[3][5] p[5]
80 // Roots of poly p[0] x^3 + p[1] x^2...+p[3]=0
81 // x=r[1][k] + i r[2][k] k=1,...,3
82 // Assumes 0<arctan(x)<pi/2 for x>0
83 
85 {
86  G4double x,t,b,c,d;
87  G4int k;
88 
89  if( p[0] != 1. )
90  {
91  for(k = 1; k < 4; k++ ) { p[k] = p[k]/p[0]; }
92  p[0] = 1.;
93  }
94  x = p[1]/3.0;
95  t = x*p[1];
96  b = 0.5*( x*( t/1.5 - p[2] ) + p[3] );
97  t = ( t - p[2] )/3.0;
98  c = t*t*t;
99  d = b*b - c;
100 
101  if( d >= 0. )
102  {
103  d = std::pow( (std::sqrt(d) + std::fabs(b) ), 1.0/3.0 );
104 
105  if( d != 0. )
106  {
107  if( b > 0. ) { b = -d; }
108  else { b = d; }
109  c = t/b;
110  }
111  d = std::sqrt(0.75)*(b - c);
112  r[2][2] = d;
113  b = b + c;
114  c = -0.5*b-x;
115  r[1][2] = c;
116 
117  if( ( b > 0. && x <= 0. ) || ( b < 0. && x > 0. ) )
118  {
119  r[1][1] = c;
120  r[2][1] = -d;
121  r[1][3] = b - x;
122  r[2][3] = 0;
123  }
124  else
125  {
126  r[1][1] = b - x;
127  r[2][1] = 0.;
128  r[1][3] = c;
129  r[2][3] = -d;
130  }
131  } // end of 2 equal or complex roots
132  else
133  {
134  if( b == 0. ) { d = std::atan(1.0)/1.5; }
135  else { d = std::atan( std::sqrt(-d)/std::fabs(b) )/3.0; }
136 
137  if( b < 0. ) { b = std::sqrt(t)*2.0; }
138  else { b = -2.0*std::sqrt(t); }
139 
140  c = std::cos(d)*b;
141  t = -std::sqrt(0.75)*std::sin(d)*b - 0.5*c;
142  d = -t - c - x;
143  c = c - x;
144  t = t - x;
145 
146  if( std::fabs(c) > std::fabs(t) ) { r[1][3] = c; }
147  else
148  {
149  r[1][3] = t;
150  t = c;
151  }
152  if( std::fabs(d) > std::fabs(t) ) { r[1][2] = d; }
153  else
154  {
155  r[1][2] = t;
156  t = d;
157  }
158  r[1][1] = t;
159 
160  for(k = 1; k < 4; k++ ) { r[2][k] = 0.; }
161  }
162  return 0;
163 }
164 
166 //
167 // Array r[3][5] p[5]
168 // Roots of poly p[0] x^4 + p[1] x^3...+p[4]=0
169 // x=r[1][k] + i r[2][k] k=1,...,4
170 
172 {
173  G4double a, b, c, d, e;
174  G4int i, k, j;
175 
176  if(p[0] != 1.0)
177  {
178  for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
179  p[0] = 1.;
180  }
181  e = 0.25*p[1];
182  b = 2*e;
183  c = b*b;
184  d = 0.75*c;
185  b = p[3] + b*( c - p[2] );
186  a = p[2] - d;
187  c = p[4] + e*( e*a - p[3] );
188  a = a - d;
189 
190  p[1] = 0.5*a;
191  p[2] = (p[1]*p[1]-c)*0.25;
192  p[3] = b*b/(-64.0);
193 
194  if( p[3] < 0. )
195  {
196  CubicRoots(p,r);
197 
198  for( k = 1; k < 4; k++ )
199  {
200  if( r[2][k] == 0. && r[1][k] > 0 )
201  {
202  d = r[1][k]*4;
203  a = a + d;
204 
205  if ( a >= 0. && b >= 0.) { p[1] = std::sqrt(d); }
206  else if( a <= 0. && b <= 0.) { p[1] = std::sqrt(d); }
207  else { p[1] = -std::sqrt(d); }
208 
209  b = 0.5*( a + b/p[1] );
210 
211  p[2] = c/b;
212  QuadRoots(p,r);
213 
214  for( i = 1; i < 3; i++ )
215  {
216  for( j = 1; j < 3; j++ ) { r[j][i+2] = r[j][i]; }
217  }
218  p[1] = -p[1];
219  p[2] = b;
220  QuadRoots(p,r);
221 
222  for( i = 1; i < 5; i++ ) { r[1][i] = r[1][i] - e; }
223 
224  return 4;
225  }
226  }
227  }
228  if( p[2] < 0. )
229  {
230  b = std::sqrt(c);
231  d = b + b - a;
232  p[1] = 0.;
233 
234  if( d > 0. ) { p[1] = std::sqrt(d); }
235  }
236  else
237  {
238  if( p[1] > 0.) { b = std::sqrt(p[2])*2.0 + p[1]; }
239  else { b = -std::sqrt(p[2])*2.0 + p[1]; }
240 
241  if( b != 0.) { p[1] = 0; }
242  else
243  {
244  for(k = 1; k < 5; k++ )
245  {
246  r[1][k] = -e;
247  r[2][k] = 0;
248  }
249  return 0;
250  }
251  }
252 
253  p[2] = c/b;
254  QuadRoots(p,r);
255 
256  for( k = 1; k < 3; k++ )
257  {
258  for( j = 1; j < 3; j++ ) { r[j][k+2] = r[j][k]; }
259  }
260  p[1] = -p[1];
261  p[2] = b;
262  QuadRoots(p,r);
263 
264  for( k = 1; k < 5; k++ ) { r[1][k] = r[1][k] - e; }
265 
266  return 4;
267 }
268 
270 
272 {
273  G4double a0, a1, a2, a3, y1;
274  G4double R2, D2, E2, D, E, R = 0.;
275  G4double a, b, c, d, ds;
276 
277  G4double reRoot[4];
278  G4int k, noReRoots = 0;
279 
280  for( k = 0; k < 4; k++ ) { reRoot[k] = DBL_MAX; }
281 
282  if( p[0] != 1.0 )
283  {
284  for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
285  p[0] = 1.;
286  }
287  a3 = p[1];
288  a2 = p[2];
289  a1 = p[3];
290  a0 = p[4];
291 
292  // resolvent cubic equation cofs:
293 
294  p[1] = -a2;
295  p[2] = a1*a3 - 4*a0;
296  p[3] = 4*a2*a0 - a1*a1 - a3*a3*a0;
297 
298  CubicRoots(p,r);
299 
300  for( k = 1; k < 4; k++ )
301  {
302  if( r[2][k] == 0. ) // find a real root
303  {
304  noReRoots++;
305  reRoot[k] = r[1][k];
306  }
307  else reRoot[k] = DBL_MAX; // kInfinity;
308  }
309  y1 = DBL_MAX; // kInfinity;
310  for( k = 1; k < 4; k++ )
311  {
312  if ( reRoot[k] < y1 ) { y1 = reRoot[k]; }
313  }
314 
315  R2 = 0.25*a3*a3 - a2 + y1;
316  b = 0.25*(4*a3*a2 - 8*a1 - a3*a3*a3);
317  c = 0.75*a3*a3 - 2*a2;
318  a = c - R2;
319  d = 4*y1*y1 - 16*a0;
320 
321  if( R2 > 0.)
322  {
323  R = std::sqrt(R2);
324  D2 = a + b/R;
325  E2 = a - b/R;
326 
327  if( D2 >= 0. )
328  {
329  D = std::sqrt(D2);
330  r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
331  r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
332  r[2][1] = 0.;
333  r[2][2] = 0.;
334  }
335  else
336  {
337  D = std::sqrt(-D2);
338  r[1][1] = -0.25*a3 + 0.5*R;
339  r[1][2] = -0.25*a3 + 0.5*R;
340  r[2][1] = 0.5*D;
341  r[2][2] = -0.5*D;
342  }
343  if( E2 >= 0. )
344  {
345  E = std::sqrt(E2);
346  r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
347  r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
348  r[2][3] = 0.;
349  r[2][4] = 0.;
350  }
351  else
352  {
353  E = std::sqrt(-E2);
354  r[1][3] = -0.25*a3 - 0.5*R;
355  r[1][4] = -0.25*a3 - 0.5*R;
356  r[2][3] = 0.5*E;
357  r[2][4] = -0.5*E;
358  }
359  }
360  else if( R2 < 0.)
361  {
362  R = std::sqrt(-R2);
363  G4complex CD2(a,-b/R);
364  G4complex CD = std::sqrt(CD2);
365 
366  r[1][1] = -0.25*a3 + 0.5*real(CD);
367  r[1][2] = -0.25*a3 - 0.5*real(CD);
368  r[2][1] = 0.5*R + 0.5*imag(CD);
369  r[2][2] = 0.5*R - 0.5*imag(CD);
370  G4complex CE2(a,b/R);
371  G4complex CE = std::sqrt(CE2);
372 
373  r[1][3] = -0.25*a3 + 0.5*real(CE);
374  r[1][4] = -0.25*a3 - 0.5*real(CE);
375  r[2][3] = -0.5*R + 0.5*imag(CE);
376  r[2][4] = -0.5*R - 0.5*imag(CE);
377  }
378  else // R2=0 case
379  {
380  if(d >= 0.)
381  {
382  D2 = c + std::sqrt(d);
383  E2 = c - std::sqrt(d);
384 
385  if( D2 >= 0. )
386  {
387  D = std::sqrt(D2);
388  r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
389  r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
390  r[2][1] = 0.;
391  r[2][2] = 0.;
392  }
393  else
394  {
395  D = std::sqrt(-D2);
396  r[1][1] = -0.25*a3 + 0.5*R;
397  r[1][2] = -0.25*a3 + 0.5*R;
398  r[2][1] = 0.5*D;
399  r[2][2] = -0.5*D;
400  }
401  if( E2 >= 0. )
402  {
403  E = std::sqrt(E2);
404  r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
405  r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
406  r[2][3] = 0.;
407  r[2][4] = 0.;
408  }
409  else
410  {
411  E = std::sqrt(-E2);
412  r[1][3] = -0.25*a3 - 0.5*R;
413  r[1][4] = -0.25*a3 - 0.5*R;
414  r[2][3] = 0.5*E;
415  r[2][4] = -0.5*E;
416  }
417  }
418  else
419  {
420  ds = std::sqrt(-d);
421  G4complex CD2(c,ds);
422  G4complex CD = std::sqrt(CD2);
423 
424  r[1][1] = -0.25*a3 + 0.5*real(CD);
425  r[1][2] = -0.25*a3 - 0.5*real(CD);
426  r[2][1] = 0.5*R + 0.5*imag(CD);
427  r[2][2] = 0.5*R - 0.5*imag(CD);
428 
429  G4complex CE2(c,-ds);
430  G4complex CE = std::sqrt(CE2);
431 
432  r[1][3] = -0.25*a3 + 0.5*real(CE);
433  r[1][4] = -0.25*a3 - 0.5*real(CE);
434  r[2][3] = -0.5*R + 0.5*imag(CE);
435  r[2][4] = -0.5*R - 0.5*imag(CE);
436  }
437  }
438  return 4;
439 }
440 
441 //
442 //