42 fprintf(
ioQQQ,
"\n ConvTempEdenIoniz called\n" );
44 if(
trace.nTrConvg >= 2 )
46 fprintf(
ioQQQ,
" ConvTempEdenIoniz called, entering temp loop using solver %s.\n",
50 if( strcmp(
conv.chSolverTemp ,
"vWDB" ) == 0 )
52 conv.lgConvTemp =
false;
84 conv.lgConvTemp =
true;
88 fprintf(
ioQQQ,
" ConvTempEdenIoniz: Te %e C %.4e H %.4e\n",
90 fprintf(
ioQQQ,
" ConvTempEdenIoniz returns ok.\n" );
98 double t1=0, error1=0, t2, error2;
103 for(
int n=0; n < 5 && !
lgAbort; ++n )
105 const int DEF_ITER = 10;
106 const double DEF_FACTOR = 0.2;
107 double step, factor = DEF_FACTOR;
116 for(
int i=0; i < 100 && !
lgAbort; ++i )
121 double maxstep = factor*t1;
125 if( step == 0.0 || step > maxstep )
127 t2 =
max( t1 +
sign( step, -error1 ),
phycon.TEMP_LIMIT_LOW );
129 TeTrack.
add( t2, error2 );
135 if( i >= n && error1*error2 <= 0. )
143 fprintf(
ioQQQ,
" PROBLEM DISASTER - the kinetic temperature appears to be below the lower limit of the code,"
144 " %.3eK. It does not bracket thermal balance.\n",
146 fprintf(
ioQQQ,
" This calculation is aborting.\n Sorry.\n");
152 if(
trace.nTrConvg >= 2 && error1*error2 > 0. && !
lgAbort )
154 fprintf(
ioQQQ,
" ConvTempEdenIoniz: bracket1 fails t1: %e %e t2: %e %e\n",
155 t1, error1, t2, error2 );
164 if( TeTrack.
init_bracket( t1, error1, t2, error2 ) == 0 )
171 TeTrack.
set_tol(2.*DBL_EPSILON*t2);
173 if( error1 != 0.0 || error2 != 0.0 )
174 t2 = (t1*error2-t2*error1)/(error2-error1);
178 for(
int i = 0; i < (1<<(n/2))*DEF_ITER && !
lgAbort; i++ )
185 TeTrack.
add( t2, error2 );
189 if(
conv.lgConvTemp )
194 fprintf(
ioQQQ,
" ConvTempEdenIoniz: brent fails\n" );
208 fprintf(
ioQQQ,
" ConvTempEdenIoniz: Te %e C %.4e H %.4e (C-H)/H %.2f%%"
209 " d(C-H)/dT %.2e +/- %.2e\n",
213 fprintf(
ioQQQ,
" ConvTempEdenIoniz returns converged=%c\n",
TorF(
conv.lgConvTemp) );
218 fprintf(
ioQQQ,
"ConvTempEdenIoniz finds insane solver %s\n",
conv.chSolverTemp );
234 conv.lgConvTemp =
false;
247 ||
thermal.lgTemperatureConstant )
254 conv.lgConvTemp =
true;
261 conv.lgConvTemp =
false;
264 if(
trace.nTrConvg >= 2 )
265 fprintf(
ioQQQ,
" lgConvTemp: C-H rel err %.4e Te rel err %.4e converged=%c\n",
270 return conv.lgConvTemp;
298 conv.hist_temp_temp.clear();
299 conv.hist_temp_heat.clear();
300 conv.hist_temp_cool.clear();
314 if(
trace.nTrConvg >= 2 )
315 fprintf(
ioQQQ,
" CoolHeatError: Te: %.4e C: %.4e H: %.4e (C-H)/H: %.4e\n",
321 if(
thermal.lgTemperatureConstant )
329 multimap<double,string> output;
332 for(
int i=0; i <
thermal.ncltot; ++i )
335 if( abs(
thermal.heatnt[i]) > thres )
338 sprintf( line,
"heat %s %e: %e %e\n",
340 output.insert( pair<const double,string>( fraction,
string(line) ) );
342 if( abs(
thermal.cooling[i]) > thres )
345 sprintf( line,
"cool %s %e: %e %e\n",
347 output.insert( pair<const double,string>( fraction,
string(line) ) );
351 dprintf(
ioQQQ,
" >>>>>>> STARTING COOLING DUMP <<<<<<\n" );
354 for( multimap<double,string>::reverse_iterator i=output.rbegin(); i != output.rend(); ++i )
356 dprintf(
ioQQQ,
" >>>>>>> FINISHED COOLING DUMP <<<<<<\n" );
361 multimap<double,string> output;
364 for(
int nelem=0; nelem <
LIMELM; ++nelem )
366 for(
int i=0; i <
LIMELM; ++i )
369 if( abs(
thermal.heating[nelem][i]) > thres )
371 sprintf( line,
"heating[%i][%i]: %e %e\n",
372 nelem, i,
thermal.heating[nelem][i], fraction );
373 output.insert( pair<const double,string>( fraction,
string(line) ) );
378 dprintf(
ioQQQ,
" >>>>>>> STARTING HEATING DUMP <<<<<<\n" );
381 for( multimap<double,string>::reverse_iterator i=output.rbegin(); i != output.rend(); ++i )
383 dprintf(
ioQQQ,
" >>>>>>> FINISHED HEATING DUMP <<<<<<\n" );
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
int dprintf(FILE *fp, const char *format,...)
bool fp_equal(sys_float x, sys_float y, int n=3)
NORETURN void TotalInsanity(void)
#define DEBUG_ENTRY(funcname)
void add(double x, double fx)
void print_history() const
double bracket_width() const
int init_bracket(double x1, double fx1, double x2, double fx2)
double deriv(int n, double &sigma) const
int ConvTempEdenIoniz(void)
STATIC double CoolHeatError(double temp)
STATIC bool lgConvTemp(const iter_track &TeTrack)
STATIC void DumpCoolStack(double thres)
STATIC void DumpHeatStack(double thres)
UNUSED const double SQRT2
void PresTotCurrent(void)
void TempChange(double TempNew, bool lgForceUpdate)