Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4BulirschStoer.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 /geometry/magneticfield/src/G4BulirschStoer.cc (Version 11.3.0) and /geometry/magneticfield/src/G4BulirschStoer.cc (Version 10.5)


  1 // *******************************************      1 // ********************************************************************
  2 // * License and Disclaimer                         2 // * License and Disclaimer                                           *
  3 // *                                                3 // *                                                                  *
  4 // * The  Geant4 software  is  copyright of th      4 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  5 // * the Geant4 Collaboration.  It is provided      5 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  6 // * conditions of the Geant4 Software License      6 // * conditions of the Geant4 Software License,  included in the file *
  7 // * LICENSE and available at  http://cern.ch/      7 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  8 // * include a list of copyright holders.           8 // * include a list of copyright holders.                             *
  9 // *                                                9 // *                                                                  *
 10 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 11 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 12 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 13 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 14 // * use.  Please see the license in the file      14 // * use.  Please see the license in the file  LICENSE  and URL above *
 15 // * for the full disclaimer and the limitatio     15 // * for the full disclaimer and the limitation of liability.         *
 16 // *                                               16 // *                                                                  *
 17 // * This  code  implementation is the result      17 // * This  code  implementation is the result of  the  scientific and *
 18 // * technical work of the GEANT4 collaboratio     18 // * technical work of the GEANT4 collaboration.                      *
 19 // * By using,  copying,  modifying or  distri     19 // * By using,  copying,  modifying or  distributing the software (or *
 20 // * any work based  on the software)  you  ag     20 // * any work based  on the software)  you  agree  to acknowledge its *
 21 // * use  in  resulting  scientific  publicati     21 // * use  in  resulting  scientific  publications,  and indicate your *
 22 // * acceptance of all terms of the Geant4 Sof     22 // * acceptance of all terms of the Geant4 Software license.          *
 23 // *******************************************     23 // ********************************************************************
 24 //                                                 24 //
                                                   >>  25 //
                                                   >>  26 //
 25 // G4BulirschStoer class implementation            27 // G4BulirschStoer class implementation
 26 // Based on bulirsch_stoer.hpp from boost          28 // Based on bulirsch_stoer.hpp from boost
 27 //                                                 29 //
 28 // Author: Dmitry Sorokin, Google Summer of Co <<  30 // Author: Dmitry Sorokin - GSoC 2016
 29 // ------------------------------------------- <<  31 //
                                                   >>  32 ///////////////////////////////////////////////////////////////////////////////
 30                                                    33 
 31 #include "G4BulirschStoer.hh"                      34 #include "G4BulirschStoer.hh"
 32                                                    35 
 33 #include "G4FieldUtils.hh"                         36 #include "G4FieldUtils.hh"
 34                                                    37 
 35 namespace                                          38 namespace
 36 {                                                  39 {
 37   constexpr G4double STEPFAC1 = 0.65;              40   constexpr G4double STEPFAC1 = 0.65;
 38   constexpr G4double STEPFAC2 = 0.94;              41   constexpr G4double STEPFAC2 = 0.94;
 39   constexpr G4double STEPFAC3 = 0.02;              42   constexpr G4double STEPFAC3 = 0.02;
 40   constexpr G4double STEPFAC4 = 4.0;               43   constexpr G4double STEPFAC4 = 4.0;
 41   constexpr G4double KFAC1 = 0.8;                  44   constexpr G4double KFAC1 = 0.8;
 42   constexpr G4double KFAC2 = 0.9;                  45   constexpr G4double KFAC2 = 0.9;
 43   constexpr G4double inv_STEPFAC1 = 1.0 / STEP     46   constexpr G4double inv_STEPFAC1 = 1.0 / STEPFAC1;
 44   constexpr G4double inv_STEPFAC4 = 1.0 / STEP     47   constexpr G4double inv_STEPFAC4 = 1.0 / STEPFAC4;
 45 } // namespace                                     48 } // namespace
 46                                                    49 
 47 G4BulirschStoer::G4BulirschStoer(G4EquationOfM     50 G4BulirschStoer::G4BulirschStoer(G4EquationOfMotion* equation,
 48                                  G4int nvar, G     51                                  G4int nvar, G4double eps_rel, G4double max_dt)
 49   : fnvar(nvar), m_eps_rel(eps_rel), m_midpoin <<  52   : fnvar(nvar), m_eps_rel(eps_rel), m_midpoint(equation,nvar),
                                                   >>  53     m_last_step_rejected(false), m_first(true), m_dt_last(0.0), m_max_dt(max_dt)
 50 {                                                  54 {
 51   /* initialize sequence of stage numbers and      55   /* initialize sequence of stage numbers and work */
 52                                                    56 
 53   for(G4int i = 0; i < m_k_max + 1; ++i)           57   for(G4int i = 0; i < m_k_max + 1; ++i)
 54   {                                                58   {
 55     m_interval_sequence[i] = 2 * (i + 1);          59     m_interval_sequence[i] = 2 * (i + 1);
 56     if (i == 0)                                    60     if (i == 0)
 57     {                                              61     {
 58       m_cost[i] = m_interval_sequence[i];          62       m_cost[i] = m_interval_sequence[i];
 59     }                                              63     }
 60     else                                           64     else
 61     {                                              65     {
 62       m_cost[i] = m_cost[i-1] + m_interval_seq     66       m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
 63     }                                              67     }
 64     for(G4int k = 0; k < i; ++k)                   68     for(G4int k = 0; k < i; ++k)
 65     {                                              69     {
 66       const G4double r = static_cast<G4double>     70       const G4double r = static_cast<G4double>(m_interval_sequence[i])
 67                        / static_cast<G4double>     71                        / static_cast<G4double>(m_interval_sequence[k]);
 68       m_coeff[i][k] = 1.0 / (r * r - 1.0); //      72       m_coeff[i][k] = 1.0 / (r * r - 1.0); // coefficients for extrapolation
 69     }                                              73     }
 70                                                    74 
 71     // crude estimate of optimal order             75     // crude estimate of optimal order
 72     m_current_k_opt = 4;                           76     m_current_k_opt = 4;
 73                                                    77 
 74     // no calculation because log10 might not      78     // no calculation because log10 might not exist for value_type!
 75                                                    79 
 76     //const G4double logfact = -log10(std::max     80     //const G4double logfact = -log10(std::max(eps_rel, 1.0e-12)) * 0.6 + 0.5;
 77     //m_current_k_opt = std::max(1.,               81     //m_current_k_opt = std::max(1.,
 78     //                  std::min(static_cast<G     82     //                  std::min(static_cast<G4double>(m_k_max-1), logfact));
 79   }                                                83   }
 80 }                                                  84 }
 81                                                    85 
 82 G4BulirschStoer::step_result                       86 G4BulirschStoer::step_result
 83 G4BulirschStoer::try_step( const G4double in[]     87 G4BulirschStoer::try_step( const G4double in[], const G4double dxdt[],
 84                            G4double& t, G4doub     88                            G4double& t, G4double out[], G4double& dt)
 85 {                                                  89 {
 86   if(m_max_dt < dt)                                90   if(m_max_dt < dt)
 87   {                                                91   {
 88     // given step size is bigger then max_dt s     92     // given step size is bigger then max_dt set limit and return fail
 89     //                                             93     //
 90     dt = m_max_dt;                                 94     dt = m_max_dt;
 91     return step_result::fail;                      95     return step_result::fail;
 92   }                                                96   }
 93                                                    97 
 94   if (dt != m_dt_last)                             98   if (dt != m_dt_last)
 95   {                                                99   {
 96     reset(); // step size changed from outside    100     reset(); // step size changed from outside -> reset
 97   }                                               101   }
 98                                                   102 
 99   G4bool reject = true;                           103   G4bool reject = true;
100                                                   104 
101   G4double new_h = dt;                            105   G4double new_h = dt;
102                                                   106 
103   /* m_current_k_opt is the estimated current     107   /* m_current_k_opt is the estimated current optimal stage number */
104                                                   108 
105   for(G4int k = 0; k <= m_current_k_opt+1; ++k    109   for(G4int k = 0; k <= m_current_k_opt+1; ++k)
106   {                                               110   {
107     // the stage counts are stored in m_interv    111     // the stage counts are stored in m_interval_sequence
108     //                                            112     //
109     m_midpoint.SetSteps(m_interval_sequence[k]    113     m_midpoint.SetSteps(m_interval_sequence[k]);
110     if(k == 0)                                    114     if(k == 0)
111     {                                             115     {
112       m_midpoint.DoStep(in, dxdt, out, dt);       116       m_midpoint.DoStep(in, dxdt, out, dt);
113       /* the first step, nothing more to do */    117       /* the first step, nothing more to do */
114     }                                             118     }
115     else                                          119     else
116     {                                             120     {
117       m_midpoint.DoStep(in, dxdt, m_table[k-1]    121       m_midpoint.DoStep(in, dxdt, m_table[k-1], dt);
118       extrapolate(k, out);                        122       extrapolate(k, out);
119       // get error estimate                       123       // get error estimate
120       for (G4int i = 0; i < fnvar; ++i)           124       for (G4int i = 0; i < fnvar; ++i)
121       {                                           125       {
122         m_err[i] = out[i] - m_table[0][i];        126         m_err[i] = out[i] - m_table[0][i];
123       }                                           127       }
124       const G4double error =                      128       const G4double error =
125             field_utils::relativeError(out, m_    129             field_utils::relativeError(out, m_err, dt, m_eps_rel);
126       h_opt[k] = calc_h_opt(dt, error, k);        130       h_opt[k] = calc_h_opt(dt, error, k);
127       work[k] = static_cast<G4double>(m_cost[k    131       work[k] = static_cast<G4double>(m_cost[k]) / h_opt[k];
128                                                   132 
129       if( (k == m_current_k_opt-1) || m_first)    133       if( (k == m_current_k_opt-1) || m_first)  // convergence before k_opt ?
130       {                                           134       {
131         if(error < 1.0)                           135         if(error < 1.0)
132         {                                         136         {
133           // convergence                          137           // convergence
134           reject = false;                         138           reject = false;
135           if( (work[k] < KFAC2 * work[k-1]) ||    139           if( (work[k] < KFAC2 * work[k-1]) || (m_current_k_opt <= 2) )
136           {                                       140           {
137             // leave order as is (except we we    141             // leave order as is (except we were in first round)
138             m_current_k_opt = std::min(m_k_max    142             m_current_k_opt = std::min(m_k_max - 1 , std::max(2 , k + 1));
139             new_h = h_opt[k];                     143             new_h = h_opt[k];
140             new_h *= static_cast<G4double>(m_c    144             new_h *= static_cast<G4double>(m_cost[k + 1])
141                    / static_cast<G4double>(m_c    145                    / static_cast<G4double>(m_cost[k]);
142           }                                       146           }
143           else                                    147           else
144           {                                       148           {
145             m_current_k_opt = std::min(m_k_max    149             m_current_k_opt = std::min(m_k_max - 1, std::max(2, k));
146             new_h = h_opt[k];                     150             new_h = h_opt[k];
147           }                                       151           }
148           break;                                  152           break;
149         }                                         153         }
150         if(should_reject(error , k) && !m_firs << 154         else if(should_reject(error , k) && !m_first)
151         {                                         155         {
152           reject = true;                          156           reject = true;
153           new_h = h_opt[k];                       157           new_h = h_opt[k];
154           break;                                  158           break;
155         }                                         159         }
156       }                                           160       }
157       if(k == m_current_k_opt)  // convergence    161       if(k == m_current_k_opt)  // convergence at k_opt ?
158       {                                           162       {
159         if(error < 1.0)                           163         if(error < 1.0)
160         {                                         164         {
161           // convergence                          165           // convergence
162           reject = false;                         166           reject = false;
163           if(work[k-1] < KFAC2 * work[k])         167           if(work[k-1] < KFAC2 * work[k])
164           {                                       168           {
165             m_current_k_opt = std::max( 2 , m_    169             m_current_k_opt = std::max( 2 , m_current_k_opt-1 );
166             new_h = h_opt[m_current_k_opt];       170             new_h = h_opt[m_current_k_opt];
167           }                                       171           }
168           else if( (work[k] < KFAC2 * work[k-1    172           else if( (work[k] < KFAC2 * work[k-1]) && !m_last_step_rejected )
169           {                                       173           {
170             m_current_k_opt = std::min(m_k_max    174             m_current_k_opt = std::min(m_k_max - 1, m_current_k_opt + 1);
171             new_h = h_opt[k];                     175             new_h = h_opt[k];
172             new_h *= static_cast<G4double>(m_c    176             new_h *= static_cast<G4double>(m_cost[m_current_k_opt])
173                    / static_cast<G4double>(m_c    177                    / static_cast<G4double>(m_cost[k]);
174           }                                       178           }
175           else                                    179           else
176           {                                       180           {
177             new_h = h_opt[m_current_k_opt];       181             new_h = h_opt[m_current_k_opt];
178           }                                       182           }
179           break;                                  183           break;
180         }                                         184         }
181         if(should_reject(error, k))            << 185         else if(should_reject(error, k))
182         {                                         186         {
183           reject = true;                          187           reject = true;
184           new_h = h_opt[m_current_k_opt];         188           new_h = h_opt[m_current_k_opt];
185           break;                                  189           break;
186         }                                         190         }
187       }                                           191       }
188       if(k == m_current_k_opt + 1)  // converg    192       if(k == m_current_k_opt + 1)  // convergence at k_opt+1 ?
189       {                                           193       {
190         if(error < 1.0)  // convergence           194         if(error < 1.0)  // convergence
191         {                                         195         {
192           reject = false;                         196           reject = false;
193           if(work[k-2] < KFAC2 * work[k-1])       197           if(work[k-2] < KFAC2 * work[k-1])
194           {                                       198           {
195             m_current_k_opt = std::max(2, m_cu    199             m_current_k_opt = std::max(2, m_current_k_opt - 1);
196           }                                       200           }
197           if((work[k] < KFAC2 * work[m_current    201           if((work[k] < KFAC2 * work[m_current_k_opt]) && !m_last_step_rejected)
198           {                                       202           {
199             m_current_k_opt = std::min(m_k_max    203             m_current_k_opt = std::min(m_k_max - 1 , k);
200           }                                       204           }
201           new_h = h_opt[m_current_k_opt];         205           new_h = h_opt[m_current_k_opt];
202         }                                         206         }
203         else                                      207         else
204         {                                         208         {
205           reject = true;                          209           reject = true;
206           new_h = h_opt[m_current_k_opt];         210           new_h = h_opt[m_current_k_opt];
207         }                                         211         }
208         break;                                    212         break;
209       }                                           213       }
210     }                                             214     }
211   }                                               215   }
212                                                   216 
213   if(!reject)                                     217   if(!reject)
214   {                                               218   {
215     t += dt;                                      219     t += dt;
216   }                                               220   }
217                                                   221 
218   if(!m_last_step_rejected || new_h < dt)         222   if(!m_last_step_rejected || new_h < dt)
219   {                                               223   {
220     // limit step size                            224     // limit step size
221     new_h = std::min(m_max_dt, new_h);            225     new_h = std::min(m_max_dt, new_h);
222     m_dt_last = new_h;                            226     m_dt_last = new_h;
223     dt = new_h;                                   227     dt = new_h;
224   }                                               228   }
225                                                   229 
226   m_last_step_rejected = reject;                  230   m_last_step_rejected = reject;
227   m_first = false;                                231   m_first = false;
228                                                   232 
229   return reject ? step_result::fail : step_res    233   return reject ? step_result::fail : step_result::success;
230 }                                                 234 }
231                                                   235 
232 void G4BulirschStoer::reset()                     236 void G4BulirschStoer::reset()
233 {                                                 237 {
234   m_first = true;                                 238   m_first = true;
235   m_last_step_rejected = false;                   239   m_last_step_rejected = false;
236 }                                                 240 }
237                                                   241 
238 void G4BulirschStoer::extrapolate(std::size_t  << 242 void G4BulirschStoer::extrapolate(size_t k , G4double xest[])
239 {                                                 243 {
240   /* polynomial extrapolation, see http://www.    244   /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
241    * uses the obtained intermediate results to    245    * uses the obtained intermediate results to extrapolate to dt->0 */
242                                                   246 
243   for(std::size_t j = k - 1 ; j > 0; --j)      << 247   for(G4int j = k - 1 ; j > 0; --j)
244   {                                               248   {
245     for (G4int i = 0; i < fnvar; ++i)             249     for (G4int i = 0; i < fnvar; ++i)
246     {                                             250     {
247       m_table[j-1][i] = m_table[j][i] * (1. +     251       m_table[j-1][i] = m_table[j][i] * (1. + m_coeff[k][j])
248                       - m_table[j-1][i] * m_co    252                       - m_table[j-1][i] * m_coeff[k][j];
249     }                                             253     }
250   }                                               254   }
251   for (G4int i = 0; i < fnvar; ++i)               255   for (G4int i = 0; i < fnvar; ++i)
252   {                                               256   {
253     xest[i] = m_table[0][i] * (1. + m_coeff[k]    257     xest[i] = m_table[0][i] * (1. + m_coeff[k][0]) - xest[i] * m_coeff[k][0];
254   }                                               258   }
255 }                                                 259 }
256                                                   260 
257 G4double                                          261 G4double
258 G4BulirschStoer::calc_h_opt(G4double h , G4dou << 262 G4BulirschStoer::calc_h_opt(G4double h , G4double error , size_t k) const
259 {                                                 263 {
260   /* calculates the optimal step size for a gi    264   /* calculates the optimal step size for a given error and stage number */
261                                                   265 
262   const G4double expo =  1.0 / (2 * k + 1);       266   const G4double expo =  1.0 / (2 * k + 1);
263   const G4double facmin = std::pow(STEPFAC3, e    267   const G4double facmin = std::pow(STEPFAC3, expo);
264   G4double fac;                                   268   G4double fac;
265                                                   269 
266   G4double facminInv= 1.0 / facmin;               270   G4double facminInv= 1.0 / facmin;
267   if (error == 0.0)                               271   if (error == 0.0)
268   {                                               272   {
269     fac = facminInv;                              273     fac = facminInv;
270   }                                               274   }
271   else                                            275   else
272   {                                               276   {
273     fac = STEPFAC2 * std::pow(error * inv_STEP    277     fac = STEPFAC2 * std::pow(error * inv_STEPFAC1 , -expo);
274     fac = std::max(facmin * inv_STEPFAC4, std:    278     fac = std::max(facmin * inv_STEPFAC4, std::min( facminInv, fac));
275   }                                               279   }
276                                                   280 
277   return h * fac;                                 281   return h * fac;
278 }                                                 282 }
279                                                   283 
280 //why is not used!!??                             284 //why is not used!!??
281 G4bool G4BulirschStoer::set_k_opt(std::size_t  << 285 G4bool G4BulirschStoer::set_k_opt(size_t k, G4double& dt)
282 {                                                 286 {
283   /* calculates the optimal stage number */       287   /* calculates the optimal stage number */
284                                                   288 
285   if(k == 1)                                      289   if(k == 1)
286   {                                               290   {
287     m_current_k_opt = 2;                          291     m_current_k_opt = 2;
288     return true;                                  292     return true;
289   }                                               293   }
290   if( (work[k-1] < KFAC1 * work[k]) || (k == m    294   if( (work[k-1] < KFAC1 * work[k]) || (k == m_k_max) )   // order decrease
291   {                                               295   {
292     m_current_k_opt = (G4int)k - 1;            << 296     m_current_k_opt = k - 1;
293     dt = h_opt[ m_current_k_opt ];                297     dt = h_opt[ m_current_k_opt ];
294     return true;                                  298     return true;
295   }                                               299   }
296   if( (work[k] < KFAC2 * work[k-1])            << 300   else if( (work[k] < KFAC2 * work[k-1])
297           || m_last_step_rejected || (k == m_k    301           || m_last_step_rejected || (k == m_k_max-1) )
298   {  // same order - also do this if last step    302   {  // same order - also do this if last step got rejected
299     m_current_k_opt = (G4int)k;                << 303     m_current_k_opt = k;
300     dt = h_opt[m_current_k_opt];                  304     dt = h_opt[m_current_k_opt];
301     return true;                                  305     return true;
302   }                                               306   }
303   else {   // order increase - only if last st    307   else {   // order increase - only if last step was not rejected
304     m_current_k_opt = (G4int)k + 1;            << 308     m_current_k_opt = k + 1;
305     dt = h_opt[m_current_k_opt - 1] * m_cost[m    309     dt = h_opt[m_current_k_opt - 1] * m_cost[m_current_k_opt]
306        / m_cost[m_current_k_opt - 1];             310        / m_cost[m_current_k_opt - 1];
307     return true;                                  311     return true;
308   }                                               312   }
309 }                                                 313 }
310                                                   314 
311 G4bool G4BulirschStoer::in_convergence_window(    315 G4bool G4BulirschStoer::in_convergence_window(G4int k) const
312 {                                                 316 {
313   if( (k == m_current_k_opt - 1) && !m_last_st    317   if( (k == m_current_k_opt - 1) && !m_last_step_rejected )
314   {                                               318   {
315     return true; // decrease stepsize only if     319     return true; // decrease stepsize only if last step was not rejected
316   }                                               320   }
317   return (k == m_current_k_opt) || (k == m_cur    321   return (k == m_current_k_opt) || (k == m_current_k_opt + 1);
318 }                                                 322 }
319                                                   323 
320                                                   324 
321 G4bool G4BulirschStoer::should_reject(G4double    325 G4bool G4BulirschStoer::should_reject(G4double error, G4int k) const
322 {                                                 326 {
323   if(k == m_current_k_opt - 1)                    327   if(k == m_current_k_opt - 1)
324   {                                               328   {
325     const auto  d = G4double(m_interval_sequen << 329     const G4double d = G4double(m_interval_sequence[m_current_k_opt]
326                   * m_interval_sequence[m_curr << 330                               * m_interval_sequence[m_current_k_opt+1]);
327     const auto  e = G4double(m_interval_sequen << 331     const G4double e = G4double(m_interval_sequence[0]);
328     const G4double e2 = e*e;                      332     const G4double e2 = e*e; 
329     // step will fail, criterion 17.3.17 in NR    333     // step will fail, criterion 17.3.17 in NR
330     return error * e2 * e2 > d * d;  //  was r    334     return error * e2 * e2 > d * d;  //  was return error > dOld * dOld; (where dOld= d/e; )
331   }                                               335   }
332   if(k == m_current_k_opt)                     << 336   else if(k == m_current_k_opt)
333   {                                               337   {
334     const auto  d = G4double(m_interval_sequen << 338     const G4double d = G4double(m_interval_sequence[m_current_k_opt]);
335     const auto  e = G4double(m_interval_sequen << 339     const G4double e = G4double(m_interval_sequence[0]);
336     return error * e * e > d * d; //  was retu    340     return error * e * e > d * d; //  was return error > dOld * dOld; (where dOld= d/e; )
337   }                                               341   }
338   else                                            342   else
339   {                                               343   {
340     return error > 1.0;                           344     return error > 1.0;
341   }                                               345   }
342 }                                                 346 }
343                                                   347