00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "cddefines.h"
00010 #include "hmi.h"
00011 #include "thermal.h"
00012 #include "iso.h"
00013 #include "hydrogenic.h"
00014 #include "colden.h"
00015 #include "h2.h"
00016 #include "stopcalc.h"
00017 #include "pressure.h"
00018 #include "dense.h"
00019 #include "trace.h"
00020 #include "phycon.h"
00021 #include "conv.h"
00022
00023
00024
00025
00026
00027
00028 static const int LIMDEF = 60;
00029
00030
00031 static double CoolHeatError( double temp );
00032
00033
00034 static double TeBrent(
00035 double x1,
00036 double x2);
00037
00038
00039 static void MakeDeriv(
00040
00041 const char *job,
00042
00043 double *DerivNumer);
00044
00045
00046 static void PutHetCol(double te,
00047 double htot,
00048 double ctot);
00049
00050
00051
00052
00053 int ConvTempEdenIoniz( void )
00054 {
00055 char chDEType;
00056 long int limtry,
00057 nTemperature_loop;
00058
00059 double CoolMHeat=0.,
00060 Damper,
00061 dtl,
00062 fneut;
00063 int kase;
00064 double DerivNumer;
00065 static double OldCmHdT = 0.;
00066 long int nTemp_oscil;
00067
00068
00069 double dCoolmHeatdT;
00070
00071 DEBUG_ENTRY( "ConvTempEdenIoniz()" );
00072
00073
00074
00075
00076
00077 if( trace.lgTrace )
00078 {
00079 fprintf( ioQQQ, "\n ConvTempEdenIoniz called\n" );
00080 }
00081 if( trace.nTrConvg>=2 )
00082 {
00083 fprintf( ioQQQ, " ConvTempEdenIoniz1 called. Te:%.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n",
00084 phycon.te, thermal.htot, thermal.ctot, (thermal.htot - thermal.ctot)*
00085 100./MAX2(1e-36,thermal.htot), dense.eden );
00086 }
00087
00088 if( strcmp( conv.chSolverTemp , "new" ) == 0 )
00089 {
00090
00091 double t1 , t2 ,
00092 error1 , error2;
00093 double TeNew;
00094 double factor = 1.02;
00095
00096 t1 = phycon.te;
00097 error1 = CoolHeatError( t1 );
00098 t2 = t1 * factor;
00099 error2 = CoolHeatError( t2 );
00100 TeNew = TeBrent( t1 , t2 );
00101 fprintf(ioQQQ , "DEBUG new temp solver error1 %.2e error2 %.2e TeNew %.2e\n",
00102 error1 , error2 , TeNew );
00103 }
00104
00105 else if( strcmp( conv.chSolverTemp , "simple") == 0 )
00106 {
00107
00108
00109
00110
00111 conv.lgTOscl = false;
00112
00113 conv.lgOscilOTS = false;
00114
00115
00116 Damper = 1.;
00117
00118
00119
00120 if( nzone < 1 )
00121 {
00122 conv.lgCmHOsc = false;
00123 }
00124
00125
00126 MakeDeriv("zero",&DerivNumer);
00127
00128
00129 limtry = LIMDEF;
00130
00131
00132 dtl = 0.;
00133
00134
00135
00136
00137 nTemperature_loop = 0;
00138 nTemp_oscil = 0;
00139
00140
00141 thermal.dTemper = 1e-3f*phycon.te;
00142
00143
00144
00145 while ( true )
00146 {
00147
00148
00149
00150
00151
00152
00153 if( ConvEdenIoniz() )
00154 {
00155
00156 DEBUG_EXIT( "ConvTempEdenIoniz()" );
00157 return 1;
00158 }
00159 PresTotCurrent();
00160
00161
00162 CoolMHeat = thermal.ctot - thermal.htot;
00163
00164
00165
00166
00167
00168 if( thermal.lgTSetOn )
00169 {
00170
00171 CoolMHeat = 0.;
00172 }
00173
00174 if( thermal.lgTLaw )
00175 {
00176
00177 if( thermal.lgTeBD96 )
00178 {
00179
00180 phycon.te = thermal.T0BD96 / (1.f + thermal.SigmaBD96 * colden.colden[ipCOL_HTOT] );
00181 }
00182 else if( thermal.lgTeSN99 )
00183 {
00184
00185
00186
00187 phycon.te = thermal.T0SN99 /
00188
00189 (float)(1.f + 9.*pow(2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN], 4.) );
00190 }
00191 else
00192 TotalInsanity();
00193 tfidle(false);
00194 }
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206 if( (lgConvTemp() )
00207 && conv.lgConvIoniz )
00208 {
00209
00210 if( trace.lgTrace )
00211 fprintf( ioQQQ, " ConvTempEdenIoniz returns ok.\n" );
00212
00213
00214
00215
00216
00217 if( thermal.dCooldT < 0. )
00218 {
00219 thermal.lgUnstable = true;
00220 }
00221 else
00222 {
00223 thermal.lgUnstable = false;
00224 }
00225
00226
00227 thermal.thist = (float)MAX2(phycon.te,thermal.thist);
00228 thermal.tlowst = (float)MIN2(phycon.te,thermal.tlowst);
00229
00230
00231 conv.lgConvTemp = true;
00232
00233 if( trace.nTrConvg>=2 )
00234 {
00235 fprintf( ioQQQ, " ConvTempEdenIoniz2 converged. Te:%11.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n",
00236 phycon.te, thermal.htot, thermal.ctot, (thermal.htot - thermal.ctot)*
00237 100./MAX2(1e-36,thermal.htot), dense.eden );
00238 }
00239
00240 DEBUG_EXIT( "ConvTempEdenIoniz()" );
00241
00242
00243
00244
00245
00246 return 0;
00247 }
00248 else if(nTemperature_loop >= limtry || lgAbort )
00249 {
00250
00251
00252
00253
00254 if( nTemperature_loop == LIMDEF && !conv.lgTOscl &&
00255 conv.lgConvIoniz )
00256 {
00257
00258
00259
00260 limtry = LIMDEF*4;
00261 if( trace.nTrConvg>2 )
00262 {
00263 fprintf( ioQQQ, " ConvTempEdenIoniz9 increases limtry to%3ld\n",
00264 limtry );
00265 }
00266 }
00267 else
00268 {
00269
00270 DEBUG_EXIT( "ConvTempEdenIoniz()" );
00271 break;
00272 }
00273 }
00274
00275
00276 MakeDeriv("incr",&DerivNumer);
00277
00278
00279 PutHetCol(phycon.te,thermal.htot,thermal.ctot);
00280
00281
00282
00283
00284
00285 if( !conv.lgConvEden )
00286 {
00287
00288
00289 CoolMHeat = 0.;
00290 if( trace.nTrConvg>=2 )
00291 {
00292 fprintf( ioQQQ,
00293 " ConvTempEdenIoniz3 breaks since lgConvEden returns failed electron density.\n");
00294 }
00295 break;
00296 }
00297
00298
00299 if( nzone < 1 )
00300 {
00301
00302 OldCmHdT = thermal.ctot/phycon.te;
00303 }
00304
00305 # define USENUMER false
00306 if( USENUMER )
00307 {
00308 dCoolmHeatdT = DerivNumer;
00309 chDEType = 'n';
00310 }
00311 else if( phycon.te < 1e4 && thermal.dCooldT < 0. )
00312 {
00313
00314
00315 dCoolmHeatdT = (double)(thermal.htot)/phycon.te*2.;
00316 chDEType = 'd';
00317 }
00318 else
00319 {
00320
00321
00322 dCoolmHeatdT = thermal.dCooldT - thermal.dHeatdT;
00323 chDEType = 'a';
00324 }
00325
00326
00327
00328 if( OldCmHdT*dCoolmHeatdT < 0. )
00329 {
00330 conv.lgCmHOsc = true;
00331
00332
00333
00334
00335
00336 dCoolmHeatdT = thermal.ctot/phycon.te;
00337 chDEType = 'o';
00338 }
00339 OldCmHdT = dCoolmHeatdT;
00340
00341
00342
00343 thermal.dTemper = (float)(CoolMHeat/dCoolmHeatdT);
00344 kase = 0;
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 if( dtl*thermal.dTemper < 0. && nTemperature_loop > 3)
00357 {
00358 conv.lgTOscl = true;
00359
00360
00361
00362 Damper *= 0.5;
00363 nTemp_oscil = nTemperature_loop;
00364 }
00365 else if( Damper < 1. && nTemperature_loop-nTemp_oscil > 10 )
00366 {
00367
00368
00369
00370 conv.lgTOscl = false;
00371 if( !((nTemperature_loop-nTemp_oscil)%5) )
00372 Damper = MIN2( 1., Damper*1.5 );
00373 }
00374
00375
00376
00377
00378 fneut = (dense.xIonDense[ipHYDROGEN][0] + 2.*hmi.H2_total)/dense.gas_phase[ipHYDROGEN];
00379
00380
00381
00382
00383
00384
00385
00386 if( 0 && h2.lgH2ON && fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > conv.HeatCoolRelErrorAllowed )
00387 {
00388
00389
00390
00391
00392
00393
00394 double absdt = fabs(thermal.dTemper);
00395 thermal.dTemper = (sign(MIN2(absdt,0.0025*phycon.te),
00396 thermal.dTemper));
00397 kase = 1;
00398 }
00399 else if( fneut > 0.08 && hydro.lgHColionImp )
00400 {
00401
00402
00403 double absdt = fabs(thermal.dTemper);
00404 thermal.dTemper = (sign(MIN2(absdt,0.005*phycon.te),
00405 thermal.dTemper));
00406 kase = 2;
00407 }
00408 # if 0
00409
00410
00411 else if( 0 && nTemperature_loop <= 3)
00412 {
00413
00414
00415 double absdt = fabs(thermal.dTemper);
00416 thermal.dTemper = (sign(MIN2(absdt,0.002*phycon.te),
00417 thermal.dTemper));
00418 kase = 3;
00419 }
00420 # endif
00421 else
00422 {
00423
00424
00425 double absdt = fabs(thermal.dTemper);
00426 thermal.dTemper = (sign(MIN2(absdt,0.02*phycon.te),
00427 thermal.dTemper));
00428 kase = 4;
00429 }
00430
00431
00432
00433
00434
00435 if( 0 && h2.lgH2ON && fabs(thermal.dTemper)/ phycon.te > 0.01 )
00436 {
00437 double absdt = fabs(thermal.dTemper);
00438 thermal.dTemper = (sign(MIN2(absdt,0.01*phycon.te),
00439 thermal.dTemper));
00440 kase = 5;
00441 }
00442
00443
00444
00445
00446
00447
00448
00449 if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] < 0.9 &&
00450
00451 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.01 )
00452 {
00453 if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 )
00454
00455
00456
00457
00458 {
00459 double absdt = fabs(thermal.dTemper);
00460
00461
00462
00463
00464
00465
00466
00467
00468 thermal.dTemper = (sign(MIN2(absdt,0.001*
00469 phycon.te),thermal.dTemper));
00470 kase = 6;
00471 }
00472 }
00473
00474
00475 thermal.dTemper *= Damper;
00476
00477
00478
00479 if( fabs(thermal.dTemper)/phycon.te <10.*DBL_EPSILON )
00480 {
00481 kase = 7;
00482 thermal.dTemper = (sign( 10.*DBL_EPSILON*phycon.te , thermal.dTemper));
00483 }
00484
00485
00486 dtl = thermal.dTemper;
00487
00488
00489 if( !thermal.lgTLaw )
00490 {
00491
00492 phycon.te = phycon.te - thermal.dTemper;
00493 tfidle(false);
00494 }
00495
00496 if( trace.nTrConvg>=2 )
00497 {
00498 fprintf( ioQQQ,
00499 " ConvTempEdenIoniz4 a %4li T:%.4e dt:%10.2e%7.2f%% C:%10.2e H:%10.2e"
00500 " CoolMHeat/h:%7.2f%% dC-H/dT:%.2e kase:%i%c damp%.2f\n",
00501 nTemperature_loop,
00502 phycon.te,
00503 thermal.dTemper,
00504 thermal.dTemper/(phycon.te+thermal.dTemper)*100,
00505 thermal.ctot,
00506 thermal.htot,
00507 CoolMHeat/MIN2(thermal.ctot , thermal.htot )*100. ,
00508 dCoolmHeatdT,
00509 kase,
00510
00511 chDEType,
00512 Damper);
00513
00514
00515 }
00516
00517 if( trace.lgTrace )
00518 {
00519 fprintf( ioQQQ,
00520 " ConvTempEdenIoniz TE:%.5e dT%.3e HTOT:%.5e CTOT:%.5e dCoolmHeatdT:%c%.4e C-H%.4e HDEN%.2e kase %i\n",
00521 phycon.te,
00522 thermal.dTemper,
00523 thermal.htot,
00524 thermal.ctot,
00525 chDEType,
00526 dCoolmHeatdT,
00527 CoolMHeat,
00528 dense.gas_phase[ipHYDROGEN],
00529 kase );
00530 }
00531
00532 # if 0
00533
00534
00535
00536
00537
00538
00539
00540 if( ConvEdenIoniz() )
00541 {
00542
00543 DEBUG_EXIT( "ConvTempEdenIoniz()" );
00544 return 1;
00545 }
00546
00547
00548
00549 PresTotCurrent();
00550
00551 CoolMHeat = thermal.ctot - thermal.htot;
00552
00553
00554
00555
00556
00557 if( phycon.te <= 3. )
00558 {
00559
00560
00561 phycon.te = 2.9998f;
00562 tfidle(false);
00563
00564
00565 PresTotCurrent();
00566 CoolMHeat = 0.;
00567 conv.lgConvIoniz = true;
00568 }
00569
00570 if( trace.nTrConvg>=2 )
00571 {
00572 fprintf( ioQQQ,
00573 " ConvTempEdenIoniz5 b %4li T:%.4e dt:%9.2e%7.1f%% C:%10.2e H:%10.2e CoolMHeat/h:%7.1f%% dC-H/dT:%.2e kase:%i damp:%.2f\n",
00574 nTemperature_loop,
00575 phycon.te,
00576 thermal.dTemper,
00577 thermal.dTemper/(phycon.te+thermal.dTemper)*100,
00578 thermal.ctot,
00579 thermal.htot,
00580 CoolMHeat/MIN2(thermal.ctot , thermal.htot )*100. ,
00581 dCoolmHeatdT,
00582 kase ,
00583 Damper);
00584
00585
00586 }
00587
00588 if( trace.lgTrace )
00589 {
00590 fprintf( ioQQQ,
00591 " ConvTempEdenIoniz TE:%11.5e dT%10.3e HTOT:%11.5e CTOT:%11.5e dCoolmHeatdT:%c%.4e C-H%.4e kase %li\n",
00592 phycon.te, thermal.dTemper, thermal.htot, thermal.ctot, chDEType,
00593 dCoolmHeatdT, CoolMHeat,
00594 kase );
00595 }
00596
00597
00598
00599
00600
00601
00602 if( phycon.te <= 3. ||
00603 ((lgConvTemp() || fabs(thermal.dTemper) < 1e-6*phycon.te)
00604 && conv.lgConvIoniz) )
00605 {
00606
00607 if( trace.lgTrace )
00608 fprintf( ioQQQ, " ConvTempEdenIoniz returns ok.\n" );
00609
00610
00611 thermal.thist = (float)MAX2(phycon.te,thermal.thist);
00612 thermal.tlowst = (float)MIN2(phycon.te,thermal.tlowst);
00613
00614
00615
00616 if( thermal.dCooldT < 0. )
00617 {
00618 thermal.lgUnstable = true;
00619 }
00620 else
00621 {
00622 thermal.lgUnstable = false;
00623 }
00624
00625
00626 conv.lgConvTemp = true;
00627
00628 if( trace.nTrConvg>=2 )
00629 {
00630 fprintf( ioQQQ, " ConvTempEdenIoniz6 converged. Te:%11.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n",
00631 phycon.te, thermal.htot, thermal.ctot, (thermal.htot - thermal.ctot)*
00632 100./MAX2(1e-36,thermal.htot), dense.eden );
00633 }
00634
00635 DEBUG_EXIT( "ConvTempEdenIoniz()" );
00636
00637
00638
00639
00640
00641 return 0;
00642 }
00643 # endif
00644
00645
00646 ++nTemperature_loop;
00647 }
00648
00649
00650
00651
00652 if( fabs(CoolMHeat/thermal.htot ) > conv.HeatCoolRelErrorAllowed )
00653 {
00654 conv.lgConvTemp = false;
00655 if( trace.nTrConvg>=2 )
00656 {
00657 fprintf( ioQQQ, " ConvTempEdenIoniz7 fell through loop, heating cooling not converged, setting conv.lgConvTemp = false\n");
00658 }
00659 }
00660 else
00661 {
00662
00663
00664 conv.lgConvTemp = true;
00665 if( trace.nTrConvg>=2 )
00666 {
00667 fprintf( ioQQQ, " ConvTempEdenIoniz8 fell through loop, heating cooling is converged (ion not?), setting conv.lgConvTemp = true\n");
00668 }
00669 }
00670 }
00671 else
00672 {
00673 fprintf( ioQQQ, "ConvTempEdenIoniz finds insane solver%s \n" , conv.chSolverTemp );
00674 ShowMe();
00675 puts( "[Stop in ConvTempEdenIoniz]" );
00676 cdEXIT(EXIT_FAILURE);
00677 }
00678
00679 DEBUG_EXIT( "ConvTempEdenIoniz()" );
00680 return 0;
00681 }
00682
00683
00684 static void MakeDeriv(
00685 const char *job,
00686 double *DerivNumer)
00687 {
00688 static long int nstore;
00689 static double OldTe , OldCool , OldHeat;
00690
00691 DEBUG_ENTRY( "MakeDeriv()" );
00692
00693 if( strcmp(job,"zero") == 0 )
00694 {
00695 nstore = 0;
00696 OldTe = 0.;
00697 OldCool = 0.;
00698 OldHeat = 0.;
00699 }
00700
00701 else if( strcmp(job,"incr") == 0 )
00702 {
00703
00704 if( nstore > 0 && !thermal.lgTSetOn && fabs(phycon.te - OldTe)>SMALLFLOAT )
00705 {
00706
00707 double dCdT = (thermal.ctot - OldCool)/(phycon.te - OldTe);
00708 double dHdT = (thermal.htot - OldHeat)/(phycon.te - OldTe);
00709 *DerivNumer = dCdT - dHdT;
00710 # if 0
00711 fprintf(ioQQQ,"derivderiv\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00712 nzone, dCdT , thermal.dCooldT ,
00713 thermal.ctot , OldCool,phycon.te , OldTe
00714 );
00715 # endif
00716 }
00717
00718 else
00719 {
00720
00721 *DerivNumer = 0.;
00722 }
00723
00724 OldTe = phycon.te;
00725 OldCool = thermal.ctot;
00726 OldHeat = thermal.htot;
00727 ++ nstore;
00728
00729 }
00730 else
00731 {
00732 fprintf( ioQQQ, " MakeDeriv called with insane job=%4.4s\n",
00733 job );
00734 puts( "[Stop in makederiv]" );
00735 cdEXIT(EXIT_FAILURE);
00736 }
00737
00738 DEBUG_EXIT( "MakeDeriv()" );
00739 return;
00740 }
00741
00742
00743 static void PutHetCol(double te,
00744 double htot,
00745 double ctot)
00746 {
00747
00748 DEBUG_ENTRY( "PutHetCol()" );
00749
00750 if( nzone == 0 )
00751 {
00752
00753 thermal.ipGrid = 0;
00754 }
00755 else
00756 {
00757 if( thermal.ipGrid >= NGRID )
00758 thermal.ipGrid = 0;
00759
00760 ASSERT( thermal.ipGrid >= 0 );
00761 ASSERT( thermal.ipGrid < NGRID );
00762 ASSERT( nzone>=0 );
00763
00764 thermal.TeGrid[thermal.ipGrid] = (float)te;
00765 thermal.HtGrid[thermal.ipGrid] = (float)htot;
00766 thermal.ClGrid[thermal.ipGrid] = (float)ctot;
00767 thermal.nZonGrid[thermal.ipGrid] = nzone;
00768
00769 ++thermal.ipGrid;
00770 }
00771
00772 DEBUG_EXIT( "PutHetCol()" );
00773 return;
00774 }
00775
00776
00777 static double CoolHeatError( double temp )
00778 {
00779
00780 phycon.te = temp;
00781 tfidle(false);
00782
00783
00784
00785
00786 if( ConvEdenIoniz() )
00787 lgAbort = true;
00788
00789
00790
00791
00792 PresTotCurrent();
00793
00794
00795 PutHetCol(phycon.te,thermal.htot,thermal.ctot);
00796
00797 return thermal.ctot - thermal.htot;
00798 }
00799
00800
00801 #define ITERMAX 100
00802
00803 static double TeBrent(
00804 double x1,
00805 double x2)
00806 {
00807 long int iter;
00808 double c,
00809 d,
00810 e,
00811 fx1,
00812 fx2,
00813 fc,
00814 p,
00815 q,
00816 r,
00817 s,
00818 xm;
00819
00820 DEBUG_ENTRY( "TeBrent()" );
00821
00822
00823
00824
00825
00826 fx1 = CoolHeatError(x1);
00827 fx2 = CoolHeatError(x2);
00828
00829 if( fx1*fx2 > 0. )
00830 {
00831
00832 fprintf( ioQQQ, " TeBrent called without proper bracket - this is impossible\n" );
00833 fprintf( ioQQQ, " Sorry.\n" );
00834 puts( "[Stop in TeBrent]" );
00835 cdEXIT(EXIT_FAILURE);
00836 }
00837
00838
00839 c = x2;
00840 fc = fx2;
00841 iter = 0;
00842
00843 e = DBL_MAX;
00844 d = DBL_MAX;
00845 while( iter < ITERMAX )
00846 {
00847 if( (fx2 > 0. && fc > 0.) || (fx2 < 0. && fc < 0.) )
00848 {
00849 c = x1;
00850 fc = fx1;
00851 d = x2 - x1;
00852 e = d;
00853 }
00854 if( fabs(fc) < fabs(fx2) )
00855 {
00856 x1 = x2;
00857 x2 = c;
00858 c = x1;
00859 fx1 = fx2;
00860 fx2 = fc;
00861 fc = fx1;
00862 }
00863
00864 xm = .5*(c - x2);
00865
00866
00867 if( fabs(xm) <= thermal.htot*conv.HeatCoolRelErrorAllowed || fx2 == 0. )
00868 break;
00869
00870 if( fabs(e) >= thermal.htot*conv.HeatCoolRelErrorAllowed && fabs(fx1) > fabs(fx2) )
00871 {
00872 double aa , bb;
00873 s = fx2/fx1;
00874
00875 if( x1 == c )
00876
00877 {
00878 p = 2.*xm*s;
00879 q = 1. - s;
00880 }
00881 else
00882 {
00883 q = fx1/fc;
00884 r = fx2/fc;
00885 p = s*(2.*xm*q*(q - r) - (x2 - x1)*(r - 1.));
00886 q = (q - 1.)*(r - 1.)*(s - 1.);
00887 }
00888 if( p > 0. )
00889 q = -q;
00890
00891 p = fabs(p);
00892 aa = fabs(thermal.htot*conv.HeatCoolRelErrorAllowed*q);
00893 bb = fabs(e*q);
00894 if( 2.*p < MIN2(3.*xm*q-aa,bb) )
00895 {
00896 e = d;
00897 d = p/q;
00898 }
00899 else
00900 {
00901 d = xm;
00902 e = d;
00903 }
00904 }
00905 else
00906 {
00907 d = xm;
00908 e = d;
00909 }
00910 x1 = x2;
00911 fx1 = fx2;
00912 if( fabs(d) > thermal.htot*conv.HeatCoolRelErrorAllowed )
00913 {
00914 x2 += d;
00915 }
00916 else
00917 {
00918 x2 += sign(thermal.htot*conv.HeatCoolRelErrorAllowed,xm);
00919 }
00920 fx2 = CoolHeatError(x2);
00921
00922 ++iter;
00923 }
00924
00925
00926
00927 if( iter == ITERMAX )
00928 {
00929
00930 fprintf( ioQQQ, " TeBrent did not converge in %i iterations\n",ITERMAX );
00931 fprintf( ioQQQ, " Sorry.\n" );
00932 puts( "[Stop in TeBrent]" );
00933 cdEXIT(EXIT_FAILURE);
00934 }
00935
00936 DEBUG_EXIT( "TeBrent()" );
00937 return( x2 );
00938
00939 }
00940
00941
00942 bool lgConvTemp(void)
00943 {
00944 bool lgRet;
00945
00946 if( fabs(thermal.htot - thermal.ctot)/thermal.htot <= conv.HeatCoolRelErrorAllowed ||
00947 thermal.lgTSetOn || phycon.te <= StopCalc.TeLowest )
00948 {
00949
00950
00951
00952 lgRet = true;
00953 conv.lgConvTemp = true;
00954 }
00955 else
00956 {
00957
00958 lgRet = false;
00959 conv.lgConvTemp = false;
00960 }
00961 return lgRet;
00962 }