Ticket #111: patch2

Index: conv_temp_eden_ioniz.cpp
===================================================================
--- conv_temp_eden_ioniz.cpp    (revision 3610)
+++ conv_temp_eden_ioniz.cpp    (working copy)
@@ -96,7 +96,6 @@
                // here starts the standard solver for variable temperature
                iter_track TeTrack;
                double t1=0, error1=0, t2, error2;
-               bool lgHighTeSearch = false;
 
                t2 = phycon.te;
                error2 = CoolHeatError( t2 );
@@ -105,103 +104,43 @@
                {
                        const int DEF_ITER = 10;
                        const double DEF_FACTOR = 0.02;
-                       double step, factor = DEF_FACTOR;
+                       const double BIG_FACTOR = 0.05;
+                       double step, factor;
 
                        TeTrack.clear();
 
-                       // set up an initial guess for the bracket
-                       // t2 was already initialized outside the main loop, or is copied from the
-                       // previous iteration. don't record this evaluation, it may be poorly converged
-                       for( int i=0; i < 2 && !lgAbort; ++i )
-                       {
-                               t1 = t2;
-                               error1 = error2;
+                       t1 = t2;
+                       error1 = error2;
 
-                               // the factor 1.2 creates 20% safety margin
-                               step = safe_div( -1.2*error1, conv.dCmHdT, 0. );
-                               step = sign( min( abs(step), factor*t1 ), step );
-                               t2 = t1 + step;
-                               error2 = CoolHeatError( t2 );
-                               TeTrack.add( t2, error2 );
-                       }
+                       step = abs(safe_div( 2.*error1, conv.dCmHdT, 0. ));
 
-                       int j = 0;
-
-                       // now hunt until we have bracketed the solution
-                       while( error1*error2 > 0. && j++ < 2*DEF_ITER && !lgAbort )
-                       {
-                               t1 = t2;
-                               error1 = error2;
-                               double deriv = TeTrack.deriv(7);
-                               // the factor 1.2 creates 20% safety margin
-                               step = safe_div( -1.2*error1, deriv, 0. );
-                               step = sign( min( abs(step), factor*t1 ), step );
-                               t2 = t1 + step;
-                               error2 = CoolHeatError( t2 );
-                               TeTrack.add( t2, error2 );
-                       }
-
-                       if( trace.nTrConvg >= 2 && error1*error2 > 0. && !lgAbort )
-                       {
-                               fprintf( ioQQQ, "  ConvTempEdenIoniz: bracket1 fails t1: %e %e t2: %e %e\n",
-                                        t1, error1, t2, error2 );
-                               TeTrack.print_history();
-                       }
-
-                       // If we reach this point without having bracketed the solution this means
-                       // that searching using the derivate doesn't work... This typically happens
-                       // when the heating and cooling curve no longer intersect, which implies
-                       // that C-H vs Te now has a local maximum, but doesn't change sign any more.
-                       // Hence we have reached a thermal front and we need to take a discontinuous
-                       // step up or down. In the first while loop we search upwards, in the second
-                       // downwards.
-
-                       int nUP = 2;
-
-                       if( error2 < 0. && !lgHighTeSearch && !lgAbort )
-                       {
-                               lgHighTeSearch = true;
-                               nUP = 25;
-                       }
-
-                       // hunting using the derivative failed, so first start hunting upwards
+                       // set up a bracket for the temperature solution
+                       // t2 was already initialized outside the main loop, or is copied from the
+                       // previous iteration. don't record this evaluation, it may be poorly converged
                        //
-                       // guard against overrunning the maximum temperature limit of the code
-                       // so that we are guaranteed to get a chance to hunt downwards
-                       while( error1*error2 > 0. && j++ < (2+nUP)*DEF_ITER && !lgAbort &&
-                              t2*(1.+factor) <= 0.99999*phycon.TEMP_LIMIT_HIGH )
+                       // if we have reached a thermal front (which happens when the heating and
+                       // cooling curves no longer intersect near the solution from the previous
+                       // call to ConvTempEdenIoniz()) we may need to make a big adjustment to the
+                       // temperature. since we are only allowed to take baby steps, we should allow
+                       // for a large number of iterations here.
+                       for( int i=0; error1*error2 > 0. && i < 250 && !lgAbort; ++i )
                        {
                                t1 = t2;
                                error1 = error2;
-                               if( lgHighTeSearch && t2 > 4000. )
-                                       factor = 0.05;
-                               t2 = t1*(1.+factor);
-                               error2 = CoolHeatError( t2 );
-                               TeTrack.add( t2, error2 );
-                       }
 
-                       // hunting upwards failed, so start hunting downwards
-                       // we may need to take a big plunge downwards, so # of iter should be big
-                       //
-                       // don't guard against dropping below the minimum temperature limit of the
-                       // code so that we get an abort if the solver has completely lost the solution.
-                       // this is a serious condition that needs to be investigated...
-                       while( error1*error2 > 0. && j++ < (2+2*nUP+9)*DEF_ITER && !lgAbort )
-                       {
-                               t1 = t2;
-                               error1 = error2;
-                               if( lgHighTeSearch && t2 <= 4000. )
-                                       factor = DEF_FACTOR;
-                               t2 = t1*(1.-factor);
+                               factor = ( t1 > 4000. ) ? BIG_FACTOR : DEF_FACTOR;
+                               step = min(step,factor*t1);
+                               t2 = t1 + sign( step, -error1 );
                                error2 = CoolHeatError( t2 );
                                TeTrack.add( t2, error2 );
+                               step *= 2.;
                        }
 
                        factor = DEF_FACTOR;
 
                        if( trace.nTrConvg >= 2 && error1*error2 > 0. && !lgAbort )
                        {
-                               fprintf( ioQQQ, "  ConvTempEdenIoniz: bracket2 fails t1: %e %e t2: %e %e\n",
+                               fprintf( ioQQQ, "  ConvTempEdenIoniz: bracket fails t1: %e %e t2: %e %e\n",
                                         t1, error1, t2, error2 );
                                TeTrack.print_history();
                        }