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