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 ]

Diff markup

Differences between /global/HEPNumerics/src/G4AnalyticalPolSolver.cc (Version 11.3.0) and /global/HEPNumerics/src/G4AnalyticalPolSolver.cc (Version 7.0.p1)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 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(G4doubl    
 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(G4doub    
 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))    
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(    
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    
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(G4dou    
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(G4do    
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     
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 *    
468         r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 *    
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 *    
484         r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 *    
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