Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/src/G4AnalyticalPolSolver.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4AnalyticalPolSolver class implementation
 27 //
 28 // Author: V.Grichine, 24.04.97
 29 // --------------------------------------------------------------------
 30 
 31 #include "globals.hh"
 32 #include <complex>
 33 
 34 #include "G4AnalyticalPolSolver.hh"
 35 
 36 //////////////////////////////////////////////////////////////////////////////
 37 
 38 G4AnalyticalPolSolver::G4AnalyticalPolSolver() { ; }
 39 
 40 //////////////////////////////////////////////////////////////////////////////
 41 
 42 G4AnalyticalPolSolver::~G4AnalyticalPolSolver() { ; }
 43 
 44 //////////////////////////////////////////////////////////////////////////////
 45 //
 46 // Array r[3][5]  p[5]
 47 // Roots of poly p[0] x^2 + p[1] x+p[2]=0
 48 //
 49 // x = r[1][k] + i r[2][k];  k = 1, 2
 50 
 51 G4int G4AnalyticalPolSolver::QuadRoots(G4double p[5], G4double r[3][5])
 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 }
 78 
 79 //////////////////////////////////////////////////////////////////////////////
 80 //
 81 // Array r[3][5]  p[5]
 82 // Roots of poly p[0] x^3 + p[1] x^2...+p[3]=0
 83 // x=r[1][k] + i r[2][k]  k=1,...,3
 84 // Assumes 0<arctan(x)<pi/2 for x>0
 85 
 86 G4int G4AnalyticalPolSolver::CubicRoots(G4double p[5], G4double r[3][5])
 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 }
196 
197 //////////////////////////////////////////////////////////////////////////////
198 //
199 // Array r[3][5]  p[5]
200 // Roots of poly p[0] x^4 + p[1] x^3...+p[4]=0
201 // x=r[1][k] + i r[2][k]  k=1,...,4
202 
203 G4int G4AnalyticalPolSolver::BiquadRoots(G4double p[5], G4double r[3][5])
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 }
336 
337 //////////////////////////////////////////////////////////////////////////////
338 
339 G4int G4AnalyticalPolSolver::QuarticRoots(G4double p[5], G4double r[3][5])
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;
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       reRoot[k] = r[1][k];
379     }
380     else
381     {
382       reRoot[k] = DBL_MAX;  // kInfinity;
383     }
384   }
385   y1 = DBL_MAX;  // kInfinity;
386   for(k = 1; k < 4; k++)
387   {
388     if(reRoot[k] < y1)
389     {
390       y1 = reRoot[k];
391     }
392   }
393 
394   R2 = 0.25 * a3 * a3 - a2 + y1;
395   b  = 0.25 * (4 * a3 * a2 - 8 * a1 - a3 * a3 * a3);
396   c  = 0.75 * a3 * a3 - 2 * a2;
397   a  = c - R2;
398   d  = 4 * y1 * y1 - 16 * a0;
399 
400   if(R2 > 0.)
401   {
402     R  = std::sqrt(R2);
403     D2 = a + b / R;
404     E2 = a - b / R;
405 
406     if(D2 >= 0.)
407     {
408       D       = std::sqrt(D2);
409       r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * D;
410       r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * D;
411       r[2][1] = 0.;
412       r[2][2] = 0.;
413     }
414     else
415     {
416       D       = std::sqrt(-D2);
417       r[1][1] = -0.25 * a3 + 0.5 * R;
418       r[1][2] = -0.25 * a3 + 0.5 * R;
419       r[2][1] = 0.5 * D;
420       r[2][2] = -0.5 * D;
421     }
422     if(E2 >= 0.)
423     {
424       E       = std::sqrt(E2);
425       r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * E;
426       r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * E;
427       r[2][3] = 0.;
428       r[2][4] = 0.;
429     }
430     else
431     {
432       E       = std::sqrt(-E2);
433       r[1][3] = -0.25 * a3 - 0.5 * R;
434       r[1][4] = -0.25 * a3 - 0.5 * R;
435       r[2][3] = 0.5 * E;
436       r[2][4] = -0.5 * E;
437     }
438   }
439   else if(R2 < 0.)
440   {
441     R = std::sqrt(-R2);
442     G4complex CD2(a, -b / R);
443     G4complex CD = std::sqrt(CD2);
444 
445     r[1][1] = -0.25 * a3 + 0.5 * real(CD);
446     r[1][2] = -0.25 * a3 - 0.5 * real(CD);
447     r[2][1] = 0.5 * R + 0.5 * imag(CD);
448     r[2][2] = 0.5 * R - 0.5 * imag(CD);
449     G4complex CE2(a, b / R);
450     G4complex CE = std::sqrt(CE2);
451 
452     r[1][3] = -0.25 * a3 + 0.5 * real(CE);
453     r[1][4] = -0.25 * a3 - 0.5 * real(CE);
454     r[2][3] = -0.5 * R + 0.5 * imag(CE);
455     r[2][4] = -0.5 * R - 0.5 * imag(CE);
456   }
457   else  // R2=0 case
458   {
459     if(d >= 0.)
460     {
461       D2 = c + std::sqrt(d);
462       E2 = c - std::sqrt(d);
463 
464       if(D2 >= 0.)
465       {
466         D       = std::sqrt(D2);
467         r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * D;
468         r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * D;
469         r[2][1] = 0.;
470         r[2][2] = 0.;
471       }
472       else
473       {
474         D       = std::sqrt(-D2);
475         r[1][1] = -0.25 * a3 + 0.5 * R;
476         r[1][2] = -0.25 * a3 + 0.5 * R;
477         r[2][1] = 0.5 * D;
478         r[2][2] = -0.5 * D;
479       }
480       if(E2 >= 0.)
481       {
482         E       = std::sqrt(E2);
483         r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * E;
484         r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * E;
485         r[2][3] = 0.;
486         r[2][4] = 0.;
487       }
488       else
489       {
490         E       = std::sqrt(-E2);
491         r[1][3] = -0.25 * a3 - 0.5 * R;
492         r[1][4] = -0.25 * a3 - 0.5 * R;
493         r[2][3] = 0.5 * E;
494         r[2][4] = -0.5 * E;
495       }
496     }
497     else
498     {
499       ds = std::sqrt(-d);
500       G4complex CD2(c, ds);
501       G4complex CD = std::sqrt(CD2);
502 
503       r[1][1] = -0.25 * a3 + 0.5 * real(CD);
504       r[1][2] = -0.25 * a3 - 0.5 * real(CD);
505       r[2][1] = 0.5 * R + 0.5 * imag(CD);
506       r[2][2] = 0.5 * R - 0.5 * imag(CD);
507 
508       G4complex CE2(c, -ds);
509       G4complex CE = std::sqrt(CE2);
510 
511       r[1][3] = -0.25 * a3 + 0.5 * real(CE);
512       r[1][4] = -0.25 * a3 - 0.5 * real(CE);
513       r[2][3] = -0.5 * R + 0.5 * imag(CE);
514       r[2][4] = -0.5 * R - 0.5 * imag(CE);
515     }
516   }
517   return 4;
518 }
519 
520 //
521 //
522 //////////////////////////////////////////////////////////////////////////////
523