00001
00002
00003
00004 #include "cddefines.h"
00005
00006 #ifdef EPS
00007 # undef EPS
00008 #endif
00009 #define EPS 1.00001
00010 #include "lines.h"
00011 #include "mole.h"
00012 #include "conv.h"
00013 #include "rfield.h"
00014 #include "iterations.h"
00015 #include "trace.h"
00016 #include "dense.h"
00017 #include "colden.h"
00018 #include "taulines.h"
00019 #include "hmi.h"
00020 #include "prt.h"
00021 #include "phycon.h"
00022 #include "geometry.h"
00023 #include "stopcalc.h"
00024 #include "opacity.h"
00025 #include "thermal.h"
00026 #include "cooling.h"
00027 #include "pressure.h"
00028 #include "radius.h"
00029 #include "called.h"
00030 #include "wind.h"
00031 #include "reason.h"
00032 #include "map.h"
00033
00034
00035 static void dmpary(void);
00036
00037 int iter_end_check(void)
00038 {
00039 bool lgDone,
00040 lgEndFun_v,
00041 lgPrinted;
00042 long int i;
00043 double oxy_in_grains;
00044
00045 DEBUG_ENTRY( "iter_end_check()" );
00046
00047
00048
00049
00050 oxy_in_grains = 0.0f;
00051 for(i=0;i<mole.num_comole_calc;++i)
00052 {
00053
00054 oxy_in_grains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];
00055 }
00056
00057
00058
00059
00060 if( trace.lgTrace )
00061 {
00062 fprintf( ioQQQ, " iter_end_check called, zone %li.\n" , nzone);
00063 }
00064 ASSERT( map.MapZone >= 00 || !conv.lgSearch );
00065
00066
00067 if( nzone == 0 )
00068 {
00069 lgEndFun_v = false;
00070
00071 if( trace.lgTrace )
00072 {
00073 fprintf( ioQQQ, " iter_end_check returns, doing nothing since zone 0.\n" );
00074 }
00075 DEBUG_EXIT( "iter_end_check()" );
00076 return( lgEndFun_v );
00077 }
00078
00079
00080 if( (nzone >= trace.nznbug && iteration >= trace.npsbug) && trace.lgTrOvrd )
00081 {
00082 if( trace.nTrConvg==0 )
00083 {
00084 geometry.nprint = 1;
00085 trace.lgTrace = true;
00086 }
00087 else
00088
00089
00090 trace.nTrConvg = abs( trace.nTrConvg );
00091 }
00092
00093
00094 if( prt.lgPrtStart && prt.nstart == nzone )
00095 {
00096 called.lgTalk = true;
00097 }
00098
00099
00100 lgDone = false;
00101
00102
00103 conv.lgBadStop = false;
00104
00105
00106
00107
00108 if( phycon.te <= StopCalc.TeFloor )
00109 {
00110 thermal.lgTSetOn = true;
00111 thermal.ConstTemp = StopCalc.TeFloor;
00112 phycon.te = thermal.ConstTemp;
00113 tfidle(false);
00114 }
00115
00116
00117
00118
00119
00120 if( (pressure.lgPres_radiation_ON && pressure.pbeta > 1.0) &&
00121 (strcmp(dense.chDenseLaw ,"CPRE") == 0) &&
00122
00123 (iterations.lgLastIt||(pressure.pbeta>1000.)) &&
00124
00125
00126
00127 (pressure.lgRadPresAbortOK||(pressure.pbeta>1000.)) )
00128 {
00129
00130 if( called.lgTalk )
00131 {
00132 fprintf( ioQQQ, "\n STOP since P(rad)/P(gas)=%7.3f >1.0\n",
00133 pressure.pbeta );
00134
00135 fprintf( ioQQQ, " Line contributors to radiation pressure follows:\n" );
00136 PrtLinePres();
00137 }
00138 lgDone = true;
00139 conv.lgBadStop = true;
00140 strncpy( reason.chReason, "of radiation pressure.", sizeof(reason.chReason) );
00141 }
00142
00143
00144
00145 else if( !wind.lgVelPos && wind.windv0 > 0.)
00146 {
00147
00148 lgDone = true;
00149 strncpy( reason.chReason, "wind veloc too small.", sizeof(reason.chReason) );
00150 }
00151
00152
00153 else if( StopCalc.lgStop21cm && (HFLines[0].TauCon >= StopCalc.tauend) )
00154 {
00155 lgDone = true;
00156 strncpy( reason.chReason, "21 cm optical depth.", sizeof(reason.chReason) );
00157 }
00158
00159 else if( rfield.extin_mag_V_extended >= StopCalc.AV_extended )
00160 {
00161
00162 lgDone = true;
00163 strncpy( reason.chReason, "A_V reached.", sizeof(reason.chReason) );
00164 }
00165
00166 else if( rfield.extin_mag_V_point >= StopCalc.AV_point )
00167 {
00168
00169 lgDone = true;
00170 strncpy( reason.chReason, "A_V reached.", sizeof(reason.chReason) );
00171 }
00172
00173 else if( StopCalc.xMass!=0. &&
00174 log10(SDIV(dense.xMassTotal))+1.0992099+ 2.*log10(radius.rinner) >= StopCalc.xMass )
00175 {
00176
00177 lgDone = true;
00178 strncpy( reason.chReason, "mass reached.", sizeof(reason.chReason) );
00179 }
00180
00181
00182
00183
00184 else if( pressure.lgSonicPoint && pressure.lgSonicPointAbortOK )
00185 {
00186
00187 lgDone = true;
00188 strncpy( reason.chReason, "sonic point reached.", sizeof(reason.chReason) );
00189 }
00190
00191 else if( dense.EdenTrue==0 )
00192 {
00193
00194 conv.lgBadStop = true;
00195 lgDone = true;
00196 strncpy( reason.chReason, "zero electron density.", sizeof(reason.chReason) );
00197 }
00198
00199 else if( radius.lgdR2Small )
00200 {
00201 lgDone = true;
00202 conv.lgBadStop = true;
00203 strncpy( reason.chReason, "DR small rel to thick.", sizeof(reason.chReason) );
00204 }
00205
00206 else if( dense.eden < StopCalc.StopElecDensity )
00207 {
00208 lgDone = true;
00209 strncpy( reason.chReason, "lowest EDEN reached.", sizeof(reason.chReason) );
00210 }
00211
00212 else if( dense.eden/dense.gas_phase[ipHYDROGEN] < StopCalc.StopElecFrac )
00213 {
00214 lgDone = true;
00215 strncpy( reason.chReason, "low electron fraction.", sizeof(reason.chReason) );
00216 }
00217
00218
00219
00220
00221 else if( dense.lgElmtOn[ipOXYGEN] &&
00222 (oxy_in_grains/MAX2(SMALLFLOAT,dense.gas_phase[ipOXYGEN])) > StopCalc.StopDepleteFrac )
00223 {
00224 lgDone = true;
00225 strncpy( reason.chReason, "freeze out fraction.", sizeof(reason.chReason) );
00226 }
00227
00228
00229 else if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > StopCalc.StopH2MoleFrac )
00230 {
00231 lgDone = true;
00232 strncpy( reason.chReason, "large H_2/H fraction.", sizeof(reason.chReason) );
00233 }
00234
00235 else if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] <
00236 StopCalc.StopHPlusFrac )
00237 {
00238 lgDone = true;
00239 strncpy( reason.chReason, "low H_+/H fraction.", sizeof(reason.chReason) );
00240 }
00241
00242 else if( radius.lgDrMinUsed )
00243 {
00244
00245 conv.lgBadStop = true;
00246 lgDone = true;
00247 strncpy( reason.chReason, "DRAD small- set DRMIN.", sizeof(reason.chReason) );
00248 }
00249
00250 else if( radius.depth >= radius.router[iteration-1]/EPS || radius.lgDrNeg )
00251 {
00252 lgDone = true;
00253 strncpy( reason.chReason, "outer radius reached.", sizeof(reason.chReason) );
00254 }
00255
00256 else if( StopCalc.iptnu >= 0 &&
00257 opac.TauAbsGeo[0][StopCalc.iptnu-1] >= StopCalc.tauend/EPS )
00258 {
00259 lgDone = true;
00260 strncpy( reason.chReason, "optical depth reached.", sizeof(reason.chReason) );
00261 }
00262
00263 else if( colden.colden[ipCOL_HTOT] >= StopCalc.HColStop/EPS )
00264 {
00265 lgDone = true;
00266 strncpy( reason.chReason, "column dens reached.", sizeof(reason.chReason) );
00267 }
00268
00269 else if( colden.colden[ipCOL_Hp] >= StopCalc.colpls/EPS )
00270 {
00271 lgDone = true;
00272 strncpy( reason.chReason, "column dens reached.", sizeof(reason.chReason) );
00273 }
00274
00275 else if( (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) >= StopCalc.col_h2/EPS )
00276 {
00277
00278 lgDone = true;
00279 strncpy( reason.chReason, "column dens reached.", sizeof(reason.chReason) );
00280 }
00281
00282 else if( (2.*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) + colden.colden[ipCOL_H0]) >= StopCalc.col_h2_nut/EPS )
00283 {
00284
00285 lgDone = true;
00286 strncpy( reason.chReason, "column dens reached.", sizeof(reason.chReason) );
00287 }
00288
00289 else if( colden.H0_ov_Tspin >= StopCalc.col_H0_ov_Tspin/EPS )
00290 {
00291
00292 lgDone = true;
00293 strncpy( reason.chReason, "column dens reached.", sizeof(reason.chReason) );
00294 }
00295
00296 else if( findspecies("CO")->hevcol >= StopCalc.col_monoxco/EPS )
00297 {
00298
00299 lgDone = true;
00300 strncpy( reason.chReason, "column dens reached.", sizeof(reason.chReason) );
00301 }
00302
00303 else if( colden.colden[ipCOL_H0] >= StopCalc.colnut/EPS )
00304 {
00305 lgDone = true;
00306 strncpy( reason.chReason, "column dens reached.", sizeof(reason.chReason) );
00307 }
00308
00309
00310
00311 else if( hmi.lgNoH2Mole &&
00312 ( (colden.colden[ipCOL_H0]+2.*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])) >= StopCalc.colnut/EPS) )
00313 {
00314 lgDone = true;
00315 strncpy( reason.chReason, "column dens reached.", sizeof(reason.chReason) );
00316 }
00317
00318 else if( phycon.te > StopCalc.T2High )
00319 {
00320 lgDone = true;
00321 strncpy( reason.chReason, "highest Te reached.", sizeof(reason.chReason) );
00322 }
00323
00324 else if( phycon.te < StopCalc.tend )
00325 {
00326 lgDone = true;
00327 strncpy( reason.chReason, "lowest Te reached.", sizeof(reason.chReason) );
00328 }
00329
00330 else if( nzone >= geometry.nend[iteration-1] )
00331 {
00332 lgDone = true;
00333 geometry.lgZoneTrp = true;
00334 strncpy( reason.chReason, "NZONE reached.", sizeof(reason.chReason) );
00335 }
00336
00337
00338
00339 else if( StopCalc.nstpl > 0 )
00340 {
00341
00342
00343 for( i=0; i < StopCalc.nstpl; i++ )
00344 {
00345
00346 if( LineSv[StopCalc.ipStopLin2[i]].sumlin[LineSave.lgLineEmergent] > 0. )
00347 {
00348 char chString[10];
00349 if( LineSv[StopCalc.ipStopLin1[i]].sumlin[LineSave.lgLineEmergent]/
00350 LineSv[StopCalc.ipStopLin2[i]].sumlin[LineSave.lgLineEmergent] > StopCalc.stpint[i] )
00351 {
00352 lgDone = true;
00353 sprt_wl( chString , StopCalc.StopLineWl[i] );
00354 sprintf( reason.chReason, "line %.9s reached.", chString );
00355 }
00356 }
00357 }
00358 }
00359
00360 else if( radius.drNext <= 0. )
00361 {
00362
00363 if( called.lgTalk )
00364 {
00365 fprintf( ioQQQ, " drNext=%10.2e STOP\n", radius.drNext );
00366 }
00367 lgDone = true;
00368 conv.lgBadStop = true;
00369 strncpy( reason.chReason, "internal error - DRAD.", sizeof(reason.chReason) );
00370 ShowMe();
00371 }
00372
00373 else if( lgAbort )
00374 {
00375
00376 conv.lgBadStop = true;
00377 lgDone = true;
00378 strncpy( reason.chReason, "calculation aborted.", sizeof(reason.chReason) );
00379 }
00380
00381 lgPrinted = false;
00382 if( lgDone )
00383 {
00384
00385 lgEndFun_v = true;
00386 PrtZone();
00387 lgPrinted = true;
00388 }
00389
00390 else
00391 {
00392
00393
00394 if( ((nzone/geometry.nprint)*geometry.nprint == nzone ||
00395 nzone == 1) || trace.nTrConvg )
00396 {
00397 PrtZone();
00398 lgPrinted = true;
00399 }
00400
00401 lgEndFun_v = false;
00402 }
00403
00404
00405 if( prt.nzdump == nzone || prt.nzdump == 0 )
00406 dmpary();
00407
00408
00409
00410
00411
00412 if( !map.lgMapDone && (map.MapZone == 0 || nzone == map.MapZone) )
00413 {
00414
00415 if( !lgPrinted )
00416 {
00417 PrtZone();
00418 }
00419
00420
00421 map.lgMapBeingDone = true;
00422
00423
00424
00425 if( ioMAP != NULL )
00426 map_do(ioMAP, " map");
00427 else
00428 map_do(ioQQQ, " map");
00429
00430
00431 lgEndFun_v = true;
00432 strncpy( reason.chReason, "MAP command used-stop.", sizeof(reason.chReason) );
00433
00434
00435
00436 iterations.itermx = 0;
00437
00438
00439 reason.chReason[sizeof(reason.chReason)-1] = '\0';
00440
00441 if( trace.lgTrace )
00442 {
00443 fprintf( ioQQQ, " iter_end_check returns after map.\n" );
00444 }
00445 DEBUG_EXIT( "iter_end_check()" );
00446 return( lgEndFun_v );
00447 }
00448
00449 if( lgEndFun_v && prt.lgOnlyZone )
00450 {
00451 puts( "[Stop in lgendfun since only zone printout desired.]" );
00452 cdEXIT(EXIT_FAILURE);
00453 }
00454
00455
00456 reason.chReason[sizeof(reason.chReason)-1] = '\0';
00457
00458 if( trace.lgTrace )
00459 {
00460 fprintf( ioQQQ, " iter_end_check bottom return.\n" );
00461 }
00462 DEBUG_EXIT( "iter_end_check()" );
00463 return( lgEndFun_v );
00464 }
00465
00466 #ifdef EPS
00467 # undef EPS
00468 #endif
00469 #define EPS 0.005
00470
00471 static void dmpary(void)
00472 {
00473 long int i;
00474 float ratio;
00475
00476 DEBUG_ENTRY( "dmpary()" );
00477
00478 fprintf( ioQQQ,
00479 " This is a print out of the cooling array for zone number %3ld\n",
00480 nzone );
00481
00482 fprintf( ioQQQ,
00483 " The total heating was HTOT=%10.2e erg/s/cm3, the total cooling was CTOT=%10.2e, and the temperature was%10.3eK.\n",
00484 thermal.htot, thermal.ctot, phycon.te );
00485
00486 fprintf( ioQQQ,
00487 " All coolants greater than%6.2f%% of the total will be printed.\n",
00488 EPS*100. );
00489
00490
00491 coolpr(ioQQQ,"ZERO",1,0.,"ZERO");
00492 for( i=0; i < thermal.ncltot; i++ )
00493 {
00494 ratio = (float)(thermal.cooling[i]/thermal.ctot);
00495 if( fabs(ratio) > EPS )
00496 {
00497 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
00498 ratio,"DOIT");
00499 }
00500
00501 ratio = (float)(thermal.heatnt[i]/thermal.ctot);
00502 if( fabs(ratio) > EPS )
00503 {
00504 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
00505 ratio,"DOIT");
00506 }
00507 }
00508 coolpr(ioQQQ,"DONE",1,0.,"DONE");
00509
00510 DEBUG_EXIT( "dmpary()" );
00511 return;
00512 }
00513