00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "physconst.h"
00008 #include "cddrive.h"
00009 #include "lines_service.h"
00010 #include "iso.h"
00011 #include "continuum.h"
00012 #include "stopcalc.h"
00013 #include "hyperfine.h"
00014 #include "dense.h"
00015 #include "grainvar.h"
00016 #include "version.h"
00017 #include "rt.h"
00018 #include "he.h"
00019 #include "ionbal.h"
00020 #include "taulines.h"
00021 #include "hydrogenic.h"
00022 #include "lines.h"
00023 #include "trace.h"
00024 #include "hcmap.h"
00025 #include "hmi.h"
00026 #include "save.h"
00027 #include "h2.h"
00028 #include "conv.h"
00029 #include "dynamics.h"
00030 #include "opacity.h"
00031 #include "geometry.h"
00032 #include "elementnames.h"
00033 #include "ca.h"
00034 #include "broke.h"
00035 #include "pressure.h"
00036 #include "mole.h"
00037 #include "co.h"
00038 #include "atoms.h"
00039 #include "abund.h"
00040 #include "rfield.h"
00041 #include "colden.h"
00042 #include "phycon.h"
00043 #include "timesc.h"
00044 #include "hextra.h"
00045 #include "radius.h"
00046 #include "iterations.h"
00047 #include "fudgec.h"
00048 #include "called.h"
00049 #include "magnetic.h"
00050 #include "wind.h"
00051 #include "secondaries.h"
00052 #include "struc.h"
00053 #include "oxy.h"
00054 #include "input.h"
00055 #include "thermal.h"
00056 #include "atmdat.h"
00057 #include "warnings.h"
00058
00059
00060 STATIC void chkCaHeps(double *totwid);
00061
00062
00063 STATIC void prt_smooth_predictions(void);
00064
00065 void PrtComment(void)
00066 {
00067 char chLbl[11],
00068 chLine[INPUT_LINE_LENGTH];
00069
00070 bool lgAbort_flag,
00071 lgThick,
00072 lgLots_of_moles;
00073
00074 long int dum1,
00075 dum2,
00076 i,
00077 imas,
00078 ipLo ,
00079 ipHi ,
00080 ipISO,
00081 nelem,
00082 isav,
00083 j,
00084 nc,
00085 nline,
00086 nn,
00087 nneg,
00088 ns,
00089 nw;
00090
00091 double big_ion_jump,
00092 absint ,
00093 aj,
00094 alpha,
00095 big,
00096 bigm,
00097 comfrc,
00098 differ,
00099 error,
00100 flur,
00101 freqn,
00102 rate,
00103 ratio,
00104 rel,
00105 small,
00106 tauneg,
00107 ts ,
00108 HBeta,
00109 relfl ,
00110 relhm,
00111 fedest,
00112 GetHeat,
00113 SumNeg,
00114 c4363,
00115 t4363,
00116 totwid;
00117
00118 double VolComputed , VolExpected , ConComputed , ConExpected;
00119
00120 bool lgLotsSolids;
00121
00122 DEBUG_ENTRY( "PrtComment()" );
00123
00124 if( 0 && lgAbort )
00125 {
00126 return;
00127 }
00128
00129 if( trace.lgTrace )
00130 {
00131 fprintf( ioQQQ, " PrtComment called.\n" );
00132 }
00133
00134
00135
00136
00137
00138
00139 iterations.lgIterAgain = false;
00140
00141
00142 wcnint();
00143
00144 if( t_version::Inst().nBetaVer > 0 )
00145 {
00146 sprintf( chLine,
00147 " !This is beta test version %ld and is intended for testing only.",
00148 t_version::Inst().nBetaVer );
00149 bangin(chLine);
00150 }
00151
00152
00153
00154
00155 if( broke.lgFixit && !t_version::Inst().lgRelease )
00156 {
00157 sprintf( chLine, " !The code needs to be fixed - search for fixit()." );
00158 bangin(chLine);
00159 }
00160
00161
00162
00163
00164 if( broke.lgCheckit && !t_version::Inst().lgRelease )
00165 {
00166 sprintf( chLine, " !New code needs to be reviewed - search for CodeReview()." );
00167 bangin(chLine);
00168 }
00169
00170
00171 if( conv.lgBadStop )
00172 {
00173
00174 sprintf( warnings.chRgcln[0], " W-Calculation stopped because %s Iteration%3ld of %ld",
00175 StopCalc.chReasonStop, iteration, iterations.itermx + 1 );
00176 sprintf( chLine, " W-Calculation stopped because %s",
00177 StopCalc.chReasonStop );
00178 warnin(chLine);
00179 sprintf( chLine, " W-This was not intended." );
00180 warnin(chLine);
00181 }
00182 else
00183 {
00184
00185 if( (conv.lgAutoIt && iteration != iterations.itermx + 1) &&
00186 iteration > 2 )
00187 {
00188 sprintf( warnings.chRgcln[0],
00189 " Calculation stopped because %s Iteration %ld of %ld, not converged due to %s",
00190 StopCalc.chReasonStop,
00191 iteration,
00192 iterations.itermx + 1,
00193 conv.chNotConverged );
00194 }
00195 else
00196 {
00197 sprintf( warnings.chRgcln[0],
00198 " Calculation stopped because %s Iteration %ld of %ld",
00199 StopCalc.chReasonStop, iteration, iterations.itermx + 1 );
00200 }
00201 }
00202
00203
00204
00205 if( (!geometry.lgZoneSet) && geometry.lgZoneTrp )
00206 {
00207 conv.lgBadStop = true;
00208 sprintf( chLine,
00209 " W-Calculation stopped because default number of zones reached. Was this intended???" );
00210 warnin(chLine);
00211 sprintf( chLine,
00212 " W-Default limit can be increased while retaining this check with the SET NEND command." );
00213 warnin(chLine);
00214 }
00215
00216
00217
00218
00219 if( radius.lgDrMinUsed || radius.lgdR2Small )
00220 {
00221 conv.lgBadStop = true;
00222 sprintf( chLine,
00223 " W-Calculation stopped zone thickness became too thin. This was not intended." );
00224 warnin(chLine);
00225 sprintf( chLine,
00226 " W-The most likely reason was an uncontrolled oscillation." );
00227 warnin(chLine);
00228 ShowMe();
00229 }
00230
00231 if( radius.lgdR2Small )
00232 {
00233 sprintf( chLine,
00234 " W-This happened because the globule scale became very small relative to the depth." );
00235 warnin(chLine);
00236 sprintf( chLine,
00237 " W-This problem is described in Hazy." );
00238 warnin(chLine);
00239 }
00240
00241
00242
00243 ASSERT( nzone < struc.nzlim );
00244
00245 if( struc.testr[nzone-1] == 0. && nzone > 1)
00246 struc.testr[nzone-1] = struc.testr[nzone-2];
00247
00248 if( struc.ednstr[nzone-1] == 0. && nzone > 1)
00249 struc.ednstr[nzone-1] = struc.ednstr[nzone-2];
00250
00251
00252 rel = radius.depth/radius.rinner;
00253 if( rel < 0.1 )
00254 {
00255 sprintf( warnings.chRgcln[1], " The geometry is plane-parallel." );
00256 }
00257 else if( rel >= 0.1 && rel < 3. )
00258 {
00259 sprintf( warnings.chRgcln[1], " The geometry is a thick shell." );
00260 }
00261 else
00262 {
00263 sprintf( warnings.chRgcln[1], " The geometry is spherical." );
00264 }
00265
00266
00267
00268
00269
00270
00271
00272
00273 if( broke.lgBroke )
00274 {
00275 sprintf( chLine, " W-The code is broken - search for broken()." );
00276 warnin(chLine);
00277 }
00278
00279
00280 if( dense.lgEdenBad )
00281 {
00282 if( dense.nzEdenBad == nzone )
00283 {
00284 sprintf( chLine, " C-The assumed electron density was incorrect for the last zone." );
00285 caunin(chLine);
00286 sprintf( chLine, " C-Did a temperature discontinuity occur??" );
00287 caunin(chLine);
00288 }
00289 else
00290 {
00291 sprintf( chLine, " W-The assumed electron density was incorrect during the calculation. This is bad." );
00292 warnin(chLine);
00293 ShowMe();
00294 }
00295 }
00296
00297 if( lgAbort )
00298 {
00299 sprintf( chLine, " W-The calculation aborted. Something REALLY went wrong!" );
00300 warnin(chLine);
00301 }
00302
00303
00304 if( hcmap.lgMapDone && !hcmap.lgMapOK )
00305 {
00306 sprintf( chLine, " !The thermal map had changes in slope - check map output." );
00307 bangin(chLine);
00308 }
00309
00310
00311
00312 if( fudgec.nfudge > 0 || fudgec.lgFudgeUsed )
00313 {
00314 sprintf( chLine, " !Fudge factors were used or were checked. Why?" );
00315 bangin(chLine);
00316 }
00317
00318 if( dense.gas_phase[ipHYDROGEN] > 1.1e13 )
00319 {
00320 if( dense.gas_phase[ipHYDROGEN] > 1e15 )
00321 {
00322 sprintf( chLine, " C-Density greater than 10**15, heavy elements are very uncertain." );
00323 caunin(chLine);
00324 }
00325 else
00326 {
00327 sprintf( chLine, " C-Density greater than 10**13" );
00328 caunin(chLine);
00329 }
00330 }
00331
00332
00333 if( cdLine("Pump",4861.36f,&relfl,&absint)<=0 )
00334 {
00335 fprintf( ioQQQ, " PROBLEM Did not find Pump H-beta, set to unity\n" );
00336 relfl = 1.;
00337 absint = 1.;
00338 }
00339
00340
00341
00342
00343 if( cdLine( "H 1",4861.36f,&HBeta,&absint)<=0 )
00344 {
00345 fprintf( ioQQQ, " NOTE Did not find H 1 H-beta - set intensity to unity, "
00346 "will not check on importance of H 1 pumping.\n" );
00347 HBeta = 1.;
00348 absint = 1.;
00349 }
00350 else
00351 {
00352
00353 if( HBeta>SMALLFLOAT )
00354 {
00355 flur = relfl/HBeta;
00356 if( flur > 0.1 )
00357 {
00358 sprintf( chLine, " !Continuum fluorescent production of H-beta was very important." );
00359 bangin(chLine);
00360 }
00361 else if(flur > 0.01 )
00362 {
00363 sprintf( chLine, " Continuum fluorescent production of H-beta was significant." );
00364 notein(chLine);
00365 }
00366 }
00367 }
00368
00369
00370 if( iteration > 1 && conv.lgAutoIt && iteration<iterations.itermx)
00371 {
00372 sprintf( chLine , " Iteration not converged because %s.",
00373 conv.chNotConverged );
00374 notein(chLine);
00375 }
00376
00377
00378 if( wind.lgBallistic() && ((!wind.lgWindOK) || wind.windv < 1e6) )
00379 {
00380 sprintf( chLine, " C-Wind velocity below sonic point; solution is not valid." );
00381 caunin(chLine);
00382 }
00383
00384
00385 if( !wind.lgStatic() )
00386 {
00387 rel = wind.emdot/(wind.windv*dense.gas_phase[ipHYDROGEN])/radius.r1r0sq;
00388 if( fabs(1.-rel)> 0.02 )
00389 {
00390 sprintf( chLine, " C-Wind mass flux error is %g%%",fabs(1.-rel)*100. );
00391 caunin(chLine);
00392 fprintf(ioQQQ,"DEBUG emdot\t%.3e\t%.3e\t%.3e\t%.3e\n",
00393 wind.emdot , wind.windv*dense.gas_phase[ipHYDROGEN],wind.windv,dense.gas_phase[ipHYDROGEN]);
00394 }
00395 }
00396
00397
00398 if( nzone >= struc.nzlim )
00399 {
00400 TotalInsanity();
00401 }
00402
00403
00404 if( !lgConserveEnergy() )
00405 {
00406 ShowMe();
00407 lgAbort = true;
00408 fprintf(ioQQQ,"\n\n Enable per zone energy conservation check by setting "
00409 "CHECK_ENERGY_EVERY_ZONE=true in cloudy.cpp, recompile, then rerun.\n");
00410 }
00411
00412
00413
00414 if( hextra.cryden*magnetic.lgB > 0. )
00415 {
00416 sprintf( chLine,
00417 " !Magnetic field & cosmic rays both present. Their interactions are not treated." );
00418 bangin(chLine);
00419 }
00420
00421
00422 if( hextra.cryden== 0. && StopCalc.TempLoStopZone < phycon.TEMP_STOP_DEFAULT)
00423 {
00424 sprintf( chLine,
00425 " !Background cosmic rays are not included - is this physical? It affects the chemistry." );
00426 bangin(chLine);
00427 }
00428
00429
00430 if( hextra.cryden > 0. && (colden.colden[ipCOL_H0]/10. + colden.colden[ipCOL_Hp]) > 1e23 )
00431 {
00432 sprintf( chLine,
00433 " C-Model is thick to cosmic rays, which are on." );
00434 caunin(chLine);
00435 }
00436
00437
00438 if( hextra.cryden == 0. && iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc < 1e-17 )
00439 {
00440 sprintf( chLine,
00441 " !Ionization rate fell below background cosmic ray ionization rate. Should this be added too?" );
00442 bangin(chLine);
00443 sprintf( chLine,
00444 " ! Use the COSMIC RAY BACKGROUND command." );
00445 bangin(chLine);
00446 }
00447
00448
00449 if( lgTestCodeCalled )
00450 {
00451 sprintf( chLine, " !Test code is in place." );
00452 bangin(chLine);
00453 }
00454
00455
00456 if( rfield.lgComUndr )
00457 {
00458 sprintf( chLine,
00459 " !Compton cooling rate underflows to zero. Is this important?" );
00460 bangin(chLine);
00461 }
00462
00463
00464 if( input.lgUnderscoreFound )
00465 {
00466 sprintf( chLine,
00467 " !Some input lines contained underscores, these were changed to spaces." );
00468 bangin(chLine);
00469 }
00470
00471
00472 if( input.lgBracketFound )
00473 {
00474 sprintf( chLine,
00475 " !Some input lines contained [ or ], these were changed to spaces." );
00476 bangin(chLine);
00477 }
00478
00479
00480 if( rfield.lgHionRad )
00481 {
00482 sprintf( chLine,
00483 " !There is no hydrogen-ionizing radiation. Was this intended?" );
00484 bangin(chLine);
00485 }
00486
00487
00488 if( thermal.nUnstable > 0 )
00489 {
00490 sprintf( chLine,
00491 " Derivative of net cooling negative and so possibly thermally unstable in%4ld zones.",
00492 thermal.nUnstable );
00493 notein(chLine);
00494 }
00495
00496
00497 if( nzone > 1 &&
00498 (realnum)(thermal.nUnstable)/(realnum)(nzone) > 0.25 )
00499 {
00500 sprintf( chLine,
00501 " !A large fraction of the zones were possibly thermally unstable,%4ld out of%4ld",
00502 thermal.nUnstable, nzone );
00503 bangin(chLine);
00504 }
00505
00506
00507 if( thermal.CoolHeatMax > 0.2 )
00508 {
00509 sprintf( chLine,
00510 " !Negative cooling reached %6.1f%% of the local heating, due to %4.4s %.1f",
00511 thermal.CoolHeatMax*100., thermal.chCoolHeatMax, thermal.wlCoolHeatMax );
00512 bangin(chLine);
00513 }
00514 else if( thermal.CoolHeatMax > 0.05 )
00515 {
00516 sprintf( chLine,
00517 " Negative cooling reached %6.1f%% of the local heating, due to %4.4s %.2f",
00518 thermal.CoolHeatMax*100., thermal.chCoolHeatMax, thermal.wlCoolHeatMax );
00519 notein(chLine);
00520 }
00521
00522
00523 if( dynamics.HeatMax > 0.05 )
00524 {
00525 sprintf( chLine,
00526 " !Advection heating reached %.2f%% of the local heating.",
00527 dynamics.HeatMax*100. );
00528 bangin(chLine);
00529 }
00530 else if( dynamics.HeatMax > 0.005 )
00531 {
00532 sprintf( chLine,
00533 " Advection heating reached %.2f%% of the local heating.",
00534 dynamics.HeatMax*100. );
00535 notein(chLine);
00536 }
00537
00538
00539 if( dynamics.CoolMax > 0.05 )
00540 {
00541 sprintf( chLine,
00542 " !Advection cooling reached %.2f%% of the local cooling.",
00543 dynamics.CoolMax*100. );
00544 bangin(chLine);
00545 }
00546 else if( dynamics.CoolMax > 0.005 )
00547 {
00548 sprintf( chLine,
00549 " Advection cooling reached %.2f%% of the local heating.",
00550 dynamics.CoolMax*100. );
00551 notein(chLine);
00552 }
00553
00554
00555
00556 if( dynamics.lgTimeDependentStatic && dynamics.lgRecom )
00557 {
00558 if( rfield.uh > 1. )
00559 {
00560 sprintf( chLine,
00561 " W-Time dependent ionization front cannot now handle strong-R cases - the ionization parameter is too large." );
00562 warnin(chLine);
00563 }
00564 else if( rfield.uh > 0.1 )
00565 {
00566 sprintf( chLine,
00567 " C-Time dependent ionization front cannot now handle strong-R cases - the ionization parameter is too large." );
00568 caunin(chLine);
00569 }
00570 }
00571
00572
00573 if( hydro.HCollIonMax > 0.10 )
00574 {
00575 sprintf( chLine,
00576 " !Thermal collisional ionization of H reached %.2f%% of the local ionization rate.",
00577 hydro.HCollIonMax*100. );
00578 bangin(chLine);
00579 }
00580 else if( hydro.HCollIonMax > 0.005 )
00581 {
00582 sprintf( chLine,
00583 " Thermal collisional ionization of H reached %.2f%% of the local ionization rate.",
00584 hydro.HCollIonMax*100. );
00585 notein(chLine);
00586 }
00587
00588
00589 if( !atmdat.lgHCaseBOK[1][ipHYDROGEN] )
00590 {
00591 sprintf( chLine,
00592 " Te-ne bounds of Case B lookup table exceeded, H I Case B line intensities set to zero." );
00593 notein(chLine);
00594 }
00595 if( !atmdat.lgHCaseBOK[1][ipHELIUM] )
00596 {
00597 sprintf( chLine,
00598 " Te-ne bounds of Case B lookup table exceeded, He II Case B line intensities set to zero." );
00599 notein(chLine);
00600 }
00601
00602 if( dense.EdenMax>1e8 )
00603 {
00604 sprintf( chLine,
00605 " !The high electron density makes the Nussbaumer/Storey CNO recombination predictions unreliable." );
00606 bangin(chLine);
00607 }
00608
00609
00610 if( secondaries.SecHIonMax > 0.10 )
00611 {
00612 sprintf( chLine,
00613 " !Suprathermal collisional ionization of H reached %.2f%% of the local H ionization rate.",
00614 secondaries.SecHIonMax*100. );
00615 bangin(chLine);
00616 }
00617 else if( secondaries.SecHIonMax > 0.005 )
00618 {
00619 sprintf( chLine,
00620 " Suprathermal collisional ionization of H reached %.2f%% of the local H ionization rate.",
00621 secondaries.SecHIonMax*100. );
00622 notein(chLine);
00623 }
00624
00625
00626 if( hmi.HeatH2DexcMax > 0.05 )
00627 {
00628 sprintf( chLine,
00629 " !H2 vib deexec heating reached %.2f%% of the local heating.",
00630 hmi.HeatH2DexcMax*100. );
00631 bangin(chLine);
00632 }
00633 else if( hmi.HeatH2DexcMax > 0.005 )
00634 {
00635 sprintf( chLine,
00636 " H2 vib deexec heating reached %.2f%% of the local heating.",
00637 hmi.HeatH2DexcMax*100. );
00638 notein(chLine);
00639 }
00640
00641
00642 if( hmi.CoolH2DexcMax > 0.05 )
00643 {
00644 sprintf( chLine,
00645 " !H2 deexec cooling reached %.2f%% of the local heating.",
00646 hmi.CoolH2DexcMax*100. );
00647 bangin(chLine);
00648 }
00649 else if( hmi.CoolH2DexcMax > 0.005 )
00650 {
00651 sprintf( chLine,
00652 " H2 deexec cooling reached %.2f%% of the local heating.",
00653 hmi.CoolH2DexcMax*100. );
00654 notein(chLine);
00655 }
00656
00657
00658 if( atmdat.HIonFracMax > 0.10 )
00659 {
00660 sprintf( chLine,
00661 " !Charge transfer ionization of H reached %.1f%% of the local H ionization rate.",
00662 atmdat.HIonFracMax*100. );
00663 bangin(chLine);
00664 }
00665 else if( atmdat.HIonFracMax > 0.005 )
00666 {
00667 sprintf( chLine,
00668 " Charge transfer ionization of H reached %.2f%% of the local H ionization rate.",
00669 atmdat.HIonFracMax*100. );
00670 notein(chLine);
00671 }
00672
00673
00674 if( atmdat.HCharHeatMax > 0.05 )
00675 {
00676 sprintf( chLine,
00677 " !Charge transfer heating reached %.2f%% of the local heating.",
00678 atmdat.HCharHeatMax*100. );
00679 bangin(chLine);
00680 }
00681 else if( atmdat.HCharHeatMax > 0.005 )
00682 {
00683 sprintf( chLine,
00684 " Charge transfer heating reached %.2f%% of the local heating.",
00685 atmdat.HCharHeatMax*100. );
00686 notein(chLine);
00687 }
00688
00689 if( atmdat.HCharCoolMax > 0.05 )
00690 {
00691 sprintf( chLine,
00692 " !Charge transfer cooling reached %.2f%% of the local heating.",
00693 atmdat.HCharCoolMax*100. );
00694 bangin(chLine);
00695 }
00696 else if( atmdat.HCharCoolMax > 0.005 )
00697 {
00698 sprintf( chLine,
00699 " Charge transfer cooling reached %.2f%% of the local heating.",
00700 atmdat.HCharCoolMax*100. );
00701 notein(chLine);
00702 }
00703
00704
00705 if( atoms.xMg2Max > 0.1 )
00706 {
00707 sprintf( chLine,
00708 " !Photoionization of upper level of Mg II 2798 reached %.1f%% of the total Mg+ photo rate.",
00709 atoms.xMg2Max*100. );
00710 bangin(chLine);
00711 }
00712 else if( atoms.xMg2Max > 0.01 )
00713 {
00714 sprintf( chLine,
00715 " Photoionization of upper level of Mg II 2798 reached %.1f%% of the total Mg+ photo rate.",
00716 atoms.xMg2Max*100. );
00717 notein(chLine);
00718 }
00719
00720
00721 if( oxy.poimax > 0.1 )
00722 {
00723 sprintf( chLine,
00724 " !Photoionization of upper levels of [O I] reached %.1f%% of the total O destruction rate.",
00725 oxy.poimax*100. );
00726 bangin(chLine);
00727 }
00728 else if( oxy.poimax > 0.01 )
00729 {
00730 sprintf( chLine,
00731 " Photoionization of upper levels of [O I] reached %.1f%% of the total O destruction rate.",
00732 oxy.poimax*100. );
00733 notein(chLine);
00734 }
00735
00736
00737 if( (oxy.poiii2Max + oxy.poiii3Max) > 0.1 )
00738 {
00739 sprintf( chLine,
00740 " !Photoionization of upper levels of [O III] reached %.1f%% of the total O++ photo rate.",
00741 (oxy.poiii2Max + oxy.poiii3Max)*100. );
00742 bangin(chLine);
00743 }
00744 else if( (oxy.poiii2Max + oxy.poiii3Max) > 0.01 )
00745 {
00746 sprintf( chLine,
00747 " Photoionization of upper levels of [O III] reached %.1f%% of the total O++ photo rate.",
00748 (oxy.poiii2Max + oxy.poiii3Max)*100. );
00749 notein(chLine);
00750 }
00751
00752
00753 if( he.frac_he0dest_23S > 0.1 )
00754 {
00755 sprintf( chLine,
00756 " !Destruction of He 2TriS reached %.1f%% of the total He0 dest rate"
00757 " at zone %li, %.1f%% of that was photoionization.",
00758 he.frac_he0dest_23S*100,
00759 he.nzone,
00760 he.frac_he0dest_23S_photo*100. );
00761 bangin(chLine);
00762 }
00763 else if( he.frac_he0dest_23S > 0.01 )
00764 {
00765 sprintf( chLine,
00766 " Destruction of He 2TriS reached %.1f%% of the total He0 dest rate"
00767 " at zone %li, %.1f%% of that was photoionization.",
00768 he.frac_he0dest_23S*100,
00769 he.nzone,
00770 he.frac_he0dest_23S_photo*100. );
00771 notein(chLine);
00772 }
00773
00774
00775 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00776 {
00777 if( !iso_ctrl.lgCritDensLMix[ipISO] && dense.lgElmtOn[ipISO] )
00778 {
00779 sprintf( chLine,
00780 " The density is too low to l-mix the lowest %s I collapsed level. "
00781 " More resolved levels are needed for accurate line ratios.",
00782 elementnames.chElementSym[ipISO]);
00783 notein(chLine);
00784 }
00785 }
00786
00787 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00788 {
00789
00790 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00791 {
00792 if( iso_sp[ipISO][nelem].lgLevelsLowered && dense.lgElmtOn[nelem] )
00793 {
00794 sprintf( chLine, " !Continuum was lowered into model %s-like %s due to high density. Highest n is %li",
00795 elementnames.chElementSym[ipISO],
00796 elementnames.chElementSym[nelem],
00797 iso_sp[ipISO][nelem].n_HighestResolved_local+iso_sp[ipISO][nelem].nCollapsed_local);
00798 bangin(chLine);
00799 }
00800 else if( iso_sp[ipISO][nelem].lgLevelsEverLowered && dense.lgElmtOn[nelem] )
00801 {
00802 sprintf( chLine, " !Continuum was lowered into model %s-like %s due to high density at SOME point but NOT at the last zone.",
00803 elementnames.chElementSym[ipISO],
00804 elementnames.chElementNameShort[nelem]);
00805 bangin(chLine);
00806 }
00807 }
00808
00809
00810 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00811 {
00812 if( iso_sp[ipISO][nelem].lgPopsRescaled )
00813 {
00814 ASSERT( dense.lgSetIoniz[nelem] );
00815 sprintf( chLine, " C-Populations were rescaled for %s-like %s due to \"element ionization\" command.",
00816 elementnames.chElementSym[ipISO],
00817 elementnames.chElementSym[nelem] );
00818 caunin(chLine);
00819 }
00820 }
00821 }
00822
00823
00824 if( !rfield.lgMMok )
00825 {
00826 sprintf( chLine,
00827 " C-Continuum not defined in extreme infrared - Compton scat, grain heating, not treated properly?" );
00828 caunin(chLine);
00829 }
00830
00831 if( !rfield.lgHPhtOK )
00832 {
00833 sprintf( chLine,
00834 " C-Continuum not defined at photon energies which ionize excited states of H, important for H- and ff heating." );
00835 caunin(chLine);
00836 }
00837
00838 if( !rfield.lgXRayOK )
00839 {
00840 sprintf( chLine,
00841 " C-Continuum not defined at X-Ray energies - Compton scattering and Auger ionization wrong?" );
00842 caunin(chLine);
00843 }
00844
00845 if( !rfield.lgGamrOK )
00846 {
00847 sprintf( chLine,
00848 " C-Continuum not defined at gamma-ray energies - pair production and Compton scattering OK?" );
00849 caunin(chLine);
00850 }
00851
00852 if( continuum.lgCon0 )
00853 {
00854 sprintf( chLine, " C-Continuum zero at some energies." );
00855 caunin(chLine);
00856 }
00857
00858 if( continuum.lgCoStarInterpolationCaution )
00859 {
00860 sprintf( chLine , " C-CoStarInterpolate interpolated between non-adjoining tracks, this may not be valid." );
00861 caunin(chLine);
00862 }
00863
00864 if( rfield.lgOcc1Hi )
00865 {
00866 sprintf( chLine,
00867 " !The continuum occupation number at 1 Ryd is greater than unity." );
00868 bangin(chLine);
00869 }
00870
00871
00872 if( radius.lgDR2Big )
00873 {
00874 sprintf( chLine,
00875 " C-The thickness of the first zone was set larger than optimal by a SET DR command." );
00876 caunin(chLine);
00877
00878
00879 if( nzone<=1 )
00880 sprintf( chLine,
00881 " C-Consider using the STOP THICKNESS command instead." );
00882 caunin(chLine);
00883 }
00884
00885
00886 if( cdLine("TOTL",4363,&t4363,&absint)<=0 )
00887 {
00888 fprintf( ioQQQ, " PrtComment could not find total O III 4363 with cdLine.\n" );
00889 ShowMe();
00890 cdEXIT(EXIT_FAILURE);
00891 }
00892
00893 if( cdLine("Coll",4363,&c4363,&absint)<=0 )
00894 {
00895 fprintf( ioQQQ, " PrtComment could not find collisional O III 4363 with cdLine.\n" );
00896 ShowMe();
00897 cdEXIT(EXIT_FAILURE);
00898 }
00899
00900
00901 if( HBeta > 1e-35 )
00902 {
00903
00904 if( t4363/HBeta > 1e-5 && dense.gas_phase[ipHYDROGEN] < 1e8 )
00905 {
00906 ratio = (t4363 - c4363)/t4363;
00907 if( ratio > 0.01 )
00908 {
00909 sprintf( chLine,
00910 " !Non-collisional excitation of [O III] 4363 reached %.2f%% of the total.",
00911 ratio*100. );
00912 bangin(chLine);
00913 }
00914 else if( ratio > 0.001 )
00915 {
00916 sprintf( chLine,
00917 " Non-collisional excitation of [O III] 4363 reached %.2f%% of the total.",
00918 ratio*100. );
00919 notein(chLine);
00920 }
00921 }
00922 }
00923
00924
00925 if( rfield.lgPlasNu )
00926 {
00927 sprintf( chLine,
00928 " !The largest plasma frequency was %.2e Ryd = %.2e micron The continuum is set to 0 below this.",
00929 rfield.plsfrqmax,
00930
00931 RYDLAM/rfield.plsfrqmax/1e4);
00932 bangin(chLine);
00933 }
00934
00935 if( rfield.occmax > 0.1 )
00936 {
00937 if( rfield.occmnu > 1e-4 )
00938 {
00939 sprintf( chLine,
00940 " !The largest continuum occupation number was %.3e at %.3e Ryd.",
00941 rfield.occmax, rfield.occmnu );
00942 bangin(chLine);
00943 }
00944 else
00945 {
00946
00947
00948 sprintf( chLine,
00949 " The largest continuum occupation number was %.3e at %.3e Ryd.",
00950 rfield.occmax, rfield.occmnu );
00951 notein(chLine);
00952 }
00953 }
00954
00955 if( rfield.occmax > 1e4 && rfield.occ1nu > 0. )
00956 {
00957
00958 if( rfield.occ1nu < 0.0912 )
00959 {
00960 sprintf( chLine,
00961 " The continuum occupation number fell below 1 at %.3e microns.",
00962 0.0912/rfield.occ1nu );
00963 notein(chLine);
00964 }
00965 else if( rfield.occ1nu < 1. )
00966 {
00967 sprintf( chLine,
00968 " The continuum occupation number fell below 1 at %.3e Angstroms.",
00969 912./rfield.occ1nu );
00970 notein(chLine);
00971 }
00972 else
00973 {
00974 sprintf( chLine,
00975 " The continuum occupation number fell below 1 at %.3e Ryd.",
00976 rfield.occ1nu );
00977 notein(chLine);
00978 }
00979 }
00980
00981 if( rfield.tbrmax > 1e3 )
00982 {
00983 sprintf( chLine,
00984 " !The largest continuum brightness temperature was %.3eK at %.3e Ryd.",
00985 rfield.tbrmax, rfield.tbrmnu );
00986 bangin(chLine);
00987 }
00988
00989 if( rfield.tbrmax > 1e4 )
00990 {
00991
00992 if( rfield.tbr4nu < 0.0912 )
00993 {
00994 sprintf( chLine,
00995 " The continuum brightness temperature fell below 10000K at %.3e microns.",
00996 0.0912/rfield.tbr4nu );
00997 notein(chLine);
00998 }
00999 else if( rfield.tbr4nu < 1. )
01000 {
01001 sprintf( chLine,
01002 " The continuum brightness temperature fell below 10000K at %.3e Angstroms.",
01003 912./rfield.tbr4nu );
01004 notein(chLine);
01005 }
01006 else
01007 {
01008 sprintf( chLine,
01009 " The continuum brightness temperature fell below 10000K at %.3e Ryd.",
01010 rfield.tbr4nu );
01011 notein(chLine);
01012 }
01013 }
01014
01015
01016 if( DoppVel.TurbVel > 0. && strcmp(dense.chDenseLaw,"CPRE") == 0 )
01017 {
01018 sprintf( chLine,
01019 " !Both constant pressure and turbulence makes no physical sense?" );
01020 bangin(chLine);
01021 }
01022
01023
01024 if( geometry.FillFac < 1. && strcmp(dense.chDenseLaw,"CPRE") == 0 )
01025 {
01026 sprintf( chLine,
01027 " !Both constant pressure and a filling factor makes no physical sense?" );
01028 bangin(chLine);
01029 }
01030
01031
01032 if( gv.lgDustOn() && abund.lgAbnSolar )
01033 {
01034 sprintf( chLine,
01035 " !Grains are present, but the gas phase abundances were left at the solar default. This is not physical." );
01036 bangin(chLine);
01037 }
01038
01039
01040 if( abund.lgDepln && !gv.lgDustOn() )
01041 {
01042 sprintf( chLine,
01043 " !Grains are not present, but the gas phase abundances were depleted. This is not physical." );
01044 bangin(chLine);
01045 }
01046
01047 if( gv.lgDustOn() )
01048 {
01049 long nBin=0L, nFail=0L;
01050 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01051 {
01052 if( gv.bin[nd]->QHeatFailures > 0L )
01053 {
01054 ++nBin;
01055 nFail += gv.bin[nd]->QHeatFailures;
01056 }
01057 }
01058 if( nFail > 0 )
01059 {
01060 sprintf( chLine,
01061 " !The grain quantum heating treatment failed to converge %ld time(s) in %ld bin(s).", nFail, nBin );
01062 bangin(chLine);
01063 }
01064 }
01065
01066 #if 0
01067
01068
01070 if( gv.lgDustOn() )
01071 {
01072 bool lgPAHsPresent_and_constant = false;
01073 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01074 {
01075 lgPAHsPresent_and_constant = lgPAHsPresent_and_constant ||
01076
01077 (gv.bin[nd]->lgPAHsInIonizedRegion );
01078 }
01079 if( lgPAHsPresent_and_constant )
01080 {
01081 sprintf( chLine,
01082 " C-PAH's were present in the ionized region, this has never been observed in H II regions." );
01083 caunin(chLine);
01084 }
01085 }
01086 #endif
01087
01088
01089 if( thermal.lgTemperatureConstant && thermal.ConstTemp*1.0001 < phycon.TEnerDen )
01090 {
01091 sprintf( chLine,
01092 " C-The continuum energy density temperature (%g K)"
01093 " is greater than the gas kinetic temperature (%g K).",
01094 phycon.TEnerDen , thermal.ConstTemp);
01095 caunin(chLine);
01096 sprintf( chLine, " C-This is unphysical." );
01097 caunin(chLine);
01098 }
01099
01100
01101 if( !gv.lgDustOn() && phycon.TEnerDen < 800. )
01102 {
01103 sprintf( chLine,
01104 " Grains were not present but might survive in this environment (energy density temperature was %.2eK)",
01105 phycon.TEnerDen );
01106 notein(chLine);
01107 }
01108
01109
01110 AgeCheck();
01111
01112
01113
01114 chkCaHeps(&totwid);
01115 if( totwid > 121. )
01116 {
01117 sprintf( chLine, " H-eps and Ca H overlap." );
01118 notein(chLine);
01119 }
01120
01121
01122 if( !phycon.lgPhysOK )
01123 {
01124 sprintf( chLine, " !A physical process has been disabled." );
01125 bangin(chLine);
01126 }
01127
01128
01129 if( dense.gas_phase[ipHYDROGEN] < 1e8 )
01130 {
01131 if( oxy.r5007Max > 0.0263f )
01132 {
01133 sprintf( chLine,
01134 " !Photoionization of upper level of [O III] 5007 reached %.2e%% of the radiative lifetime.",
01135 oxy.r5007Max*100. );
01136 bangin(chLine);
01137 }
01138 else if( oxy.r5007Max > 0.0263f/10.f )
01139 {
01140 sprintf( chLine,
01141 " Photoionization of upper level of [O III] 5007 reached %.2e%% of the radiative lifetime.",
01142 oxy.r5007Max*100. );
01143 notein(chLine);
01144 }
01145 if( oxy.r4363Max > 1.78f )
01146 {
01147 sprintf( chLine,
01148 " !Photoionization of upper level of [O III] 4363 reached %.2e%% of the radiative lifetime.",
01149 oxy.r4363Max*100. );
01150 bangin(chLine);
01151 }
01152 else if( oxy.r4363Max > 1.78f/10.f )
01153 {
01154 sprintf( chLine,
01155 " Photoionization of upper level of [O III] 4363 reached %.2e%% of the radiative lifetime.",
01156 oxy.r4363Max*100. );
01157 notein(chLine);
01158 }
01159 }
01160
01161
01162
01163 error = fabs(thermal.power-thermal.totcol)/SDIV((thermal.power + thermal.totcol)/2.);
01164 if( thermal.lgTemperatureConstant )
01165 {
01166 if( error > 0.05 )
01167 {
01168 sprintf( chLine,
01169 " !Heating - cooling mismatch =%5.1f%%. Caused by constant temperature assumption. ",
01170 error*100. );
01171 bangin(chLine);
01172 }
01173 }
01174
01175 else
01176 {
01177 if( error > 0.05 && error < 0.2 )
01178 {
01179 sprintf( chLine, " C-Heating - cooling mismatch =%.1f%%. What\'s wrong?",
01180 error*100. );
01181 caunin(chLine);
01182 }
01183 else if( error >= 0.2 )
01184 {
01185 sprintf( chLine, " W-Heating - cooling mismatch =%.2e%%. What\'s wrong????",
01186 error*100. );
01187 warnin(chLine);
01188 }
01189 }
01190
01191
01192 if( ca.Ca2RmLya > 0.01 )
01193 {
01194 sprintf( chLine,
01195 " Photoionization of Ca+ 2D level by Ly-alpha reached %6.1f%% of the total rate out.",
01196 ca.Ca2RmLya*100. );
01197 notein(chLine);
01198 }
01199
01200
01201 if( hydro.nLyaHot > 0 )
01202 {
01203 if( hydro.TLyaMax/hydro.TeLyaMax > 1.05 )
01204 {
01205 sprintf( chLine,
01206 " !The excitation temp of Lya exceeded the electron temp, largest value was %.2eK (gas temp there was %.2eK, zone%4ld)",
01207 hydro.TLyaMax, hydro.TeLyaMax, hydro.nZTLaMax );
01208 bangin(chLine);
01209 }
01210 }
01211
01212
01213
01214
01215
01216 if( cdLine("Line",0,&SumNeg,&absint)<=0 )
01217 {
01218 fprintf( ioQQQ, " did not get sumneg cdLine\n" );
01219 ShowMe();
01220 cdEXIT(EXIT_FAILURE);
01221 }
01222
01223
01224 if( cdLine("TotH",0,&GetHeat,&absint)<=0 )
01225 {
01226 fprintf( ioQQQ, " did not get GetHeat cdLine\n" );
01227 ShowMe();
01228 cdEXIT(EXIT_FAILURE);
01229 }
01230
01231 if( GetHeat > 0. )
01232 {
01233 SumNeg /= GetHeat;
01234 if( SumNeg > 0.1 )
01235 {
01236 sprintf( chLine,
01237 " !Line absorption heating reached %.2f%% of the global heating.",
01238 SumNeg*100. );
01239 bangin(chLine);
01240 }
01241 else if( SumNeg > 0.01 )
01242 {
01243 sprintf( chLine,
01244 " Line absorption heating reached %.2f%% of the global heating.",
01245 SumNeg*100. );
01246 notein(chLine);
01247 }
01248 }
01249
01250
01251 if( input.lgSetNoBuffering )
01252 {
01253 sprintf( chLine,
01254 " !NO BUFFERING command was entered - this increases exec time by LARGE amounts.");
01255 bangin(chLine);
01256 }
01257
01258
01259 if( thermal.GBarMax > 0.1 )
01260 {
01261 ASSERT( thermal.ipMaxExtra > 0 );
01262 strcpy( chLbl, chLineLbl(TauLine2[thermal.ipMaxExtra-1]) );
01263
01264 sprintf( chLine,
01265 " !G-bar cooling lines reached %.2f%% of the local cooling. Line=%.10s",
01266 thermal.GBarMax*100., chLbl );
01267 bangin(chLine);
01268 }
01269
01270 else if( thermal.GBarMax > 0.01 )
01271 {
01272 strcpy( chLbl, chLineLbl(TauLine2[thermal.ipMaxExtra-1]) );
01273
01274 sprintf( chLine,
01275 " G-bar cooling lines reached %.2f%% of the local cooling. Line=%.10s",
01276 thermal.GBarMax*100., chLbl );
01277 notein(chLine);
01278 }
01279
01280
01281 if( hyperfine.cooling_max > 0.1 )
01282 {
01283 sprintf( chLine,
01284 " !Hyperfine structure line cooling reached %.2f%% of the local cooling.",
01285 hyperfine.cooling_max*100.);
01286 bangin(chLine);
01287 }
01288
01289 else if( hyperfine.cooling_max > 0.01 )
01290 {
01291 sprintf( chLine,
01292 " Hyperfine structure line cooling reached %.2f%% of the local cooling.",
01293 hyperfine.cooling_max*100. );
01294 notein(chLine);
01295 }
01296
01297
01298
01299 if( thermal.HeatLineMax > 0.1 )
01300 {
01301 long level = -1;
01302 TransitionProxy t = FndLineHt(&level);
01303 strcpy( chLbl, chLineLbl(t) );
01304 sprintf( chLine,
01305 " !Line absorption heating reached %.2f%% of the local heating - largest by level%2ld line %.10s",
01306 thermal.HeatLineMax*100., level, chLbl );
01307 bangin(chLine);
01308 }
01309
01310 else if( thermal.HeatLineMax > 0.01 )
01311 {
01312 sprintf( chLine,
01313 " Line absorption heating reached %.2f%% of the local heating.",
01314 thermal.HeatLineMax*100. );
01315 notein(chLine);
01316 }
01317
01318 if( ionbal.CompHeating_Max > 0.05 )
01319 {
01320 sprintf( chLine,
01321 " !Bound Compton heating reached %.2f%% of the local heating.",
01322 ionbal.CompHeating_Max*100. );
01323 bangin(chLine);
01324 }
01325 else if( ionbal.CompHeating_Max > 0.01 )
01326 {
01327 sprintf( chLine,
01328 " Bound Compton heating reached %.2f%% of the local heating.",
01329 ionbal.CompHeating_Max*100. );
01330 notein(chLine);
01331 }
01332
01333
01334 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01335 {
01336 for( nelem=ipISO; nelem<LIMELM; ++nelem )
01337 {
01338 if( dense.lgElmtOn[nelem] )
01339 {
01340
01341 long int nmax = iso_sp[ipISO][nelem].numLevels_local;
01342
01343
01344 for( ipHi=1; ipHi < nmax - 1; ++ipHi )
01345 {
01346 for( ipLo=0; ipLo < ipHi; ++ipLo )
01347 {
01348 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
01349 continue;
01350
01351
01352 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() < -0.1 )
01353 {
01354 sprintf( chLine,
01355 " !Some iso-structure lines mased: %s-like %s, line %li-%li had optical depth %.2e",
01356 elementnames.chElementSym[ipISO],
01357 elementnames.chElementNameShort[nelem],
01358 ipHi , ipLo ,
01359 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() );
01360 bangin(chLine);
01361 }
01362 }
01363 }
01364 }
01365 }
01366 }
01367
01368 if( dense.gas_phase[ipHYDROGEN] < 1e7 )
01369 {
01370
01371 lgThick = false;
01372 tauneg = 0.;
01373 alpha = 0.;
01374
01375 for( i=3; i <= nLevel1; i++ )
01376 {
01377
01378 if( TauLines[i].EnergyWN() < 10000. )
01379 {
01380 if( TauLines[i].Emis().TauIn() > 1. )
01381 {
01382 lgThick = true;
01383 alpha = MAX2(alpha,(double)TauLines[i].Emis().TauIn());
01384 }
01385 else if( TauLines[i].Emis().TauIn() < (realnum)tauneg )
01386 {
01387 tauneg = TauLines[i].Emis().TauIn();
01388 strcpy( chLbl, chLineLbl(TauLines[i]) );
01389 }
01390 }
01391 }
01392
01393 if( lgThick )
01394 {
01395 sprintf( chLine,
01396 " !Some infrared fine structure lines are optically thick: largest tau was %.2e",
01397 alpha );
01398 bangin(chLine);
01399 }
01400
01401 if( tauneg < -0.01 )
01402 {
01403 sprintf( chLine,
01404 " !Some fine structure lines mased: line %s had optical depth %.2e",
01405 chLbl, tauneg );
01406 bangin(chLine);
01407 }
01408 }
01409
01410
01411 tauneg = 0.;
01412 alpha = 0.;
01413 for( i=1; i <= nLevel1; i++ )
01414 {
01415
01416 if( TauLines[i].EnergyWN() >= 10000. )
01417 {
01418 if( TauLines[i].Emis().TauIn() < (realnum)tauneg )
01419 {
01420 tauneg = TauLines[i].Emis().TauIn();
01421 strcpy( chLbl, chLineLbl(TauLines[i]) );
01422 }
01423 }
01424 }
01425
01426
01427 if( tauneg < -0.01 )
01428 {
01429 sprintf( chLine,
01430 " !Some level1 lines mased: most negative ion and tau were: %s %.2e",
01431 chLbl, tauneg );
01432 bangin(chLine);
01433 }
01434
01435
01436
01437
01438 if( geometry.lgStatic && iterations.lgLastIt && (iteration == 1) &&
01439 !geometry.lgStaticNoIt)
01440 {
01441 sprintf( chLine, " C-I must iterate when SPHERE STATIC is set." );
01442 caunin(chLine);
01443 iterations.lgIterAgain = true;
01444 }
01445
01446
01447 if( save.lgPunContinuum && iteration == 1 && iterations.lgLastIt)
01448 {
01449 sprintf( chLine, " C-I must iterate when save continuum output is done." );
01450 caunin(chLine);
01451 iterations.lgIterAgain = true;
01452 }
01453
01455
01456 {
01457 two_photon& tnu = iso_sp[ipH_LIKE][ipHYDROGEN].TwoNu[0];
01458 if( tnu.induc_dn_max > 1. )
01459 {
01460 sprintf( chLine, " !Rate of induced H 2-photon emission reached %.2e s^-1",
01461 tnu.induc_dn_max );
01462 bangin(chLine);
01463 }
01464 else if( tnu.induc_dn_max > 0.01 )
01465 {
01466 sprintf( chLine, " Rate of induced H 2-photon emission reached %.2e s^-1",
01467 tnu.induc_dn_max );
01468 notein(chLine);
01469 }
01470 }
01471
01472
01473 if( hydro.FracInd > 0.01 )
01474 {
01475 sprintf( chLine,
01476 " Induced recombination was %5.1f%% of the total for H level%3ld",
01477 hydro.FracInd*100., hydro.ndclev );
01478 notein(chLine);
01479 }
01480
01481 if( hydro.fbul > 0.01 )
01482 {
01483 sprintf( chLine,
01484 " Stimulated emission was%6.1f%% of the total for H transition%3ld -%3ld",
01485 hydro.fbul*100., hydro.nbul + 1, hydro.nbul );
01486 notein(chLine);
01487 }
01488
01489
01490
01491 if( cdLine("Fe 2",1216,&fedest,&absint)<=0 )
01492 {
01493 fprintf( ioQQQ, " Did not find Fe II Lya\n" );
01494 ShowMe();
01495 cdEXIT(EXIT_FAILURE);
01496 }
01497
01498
01499 if( cdLine("TOTL",1215.68f,&relhm,&absint)<=0 )
01500 {
01501 fprintf( ioQQQ, " Did not find Lya\n" );
01502 ShowMe();
01503 cdEXIT(EXIT_FAILURE);
01504 }
01505
01506 if( relhm > 0. )
01507 {
01508 ratio = fedest/(fedest + relhm);
01509 if( ratio > 0.1 )
01510 {
01511 sprintf( chLine, " !Fe II destruction of Ly-a removed %.1f%% of the line.",
01512 ratio *100.);
01513 bangin(chLine);
01514 }
01515 else if( ratio > 0.01 )
01516 {
01517 sprintf( chLine, " Fe II destruction of Ly-a removed %.1f%% of the line.",
01518 ratio );
01519 notein(chLine);
01520 }
01521 }
01522
01523 if( cdLine("H-CT",6563,&relhm,&absint)<=0 )
01524 {
01525 fprintf( ioQQQ, " Comment did not find H-CT H-alpha\n" );
01526 ShowMe();
01527 cdEXIT(EXIT_FAILURE);
01528 }
01529
01530 if( HBeta > 0. )
01531 {
01532 if( relhm/HBeta > 0.01 )
01533 {
01534 sprintf( chLine,
01535 " !Mutual neutralization production of H-alpha was significant." );
01536 bangin(chLine);
01537 }
01538 }
01539
01540
01541 if( hydro.lgHiPop2 )
01542 {
01543 sprintf( chLine,
01544 " The population of H n=2 reached %.2e relative to the ground state.",
01545 hydro.pop2mx );
01546 notein(chLine);
01547 }
01548
01549
01550 for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
01551 {
01552 for( nelem=ipISO; nelem < LIMELM; nelem++ )
01553 {
01554 if( iso_sp[ipISO][nelem].CaseBCheck > 1.5 )
01555 {
01556 sprintf( chLine,
01557 " Ratio of computed diffuse emission to case B reached %g for iso %li element %li",
01558 iso_sp[ipISO][nelem].CaseBCheck , ipISO , nelem+1 );
01559 notein(chLine);
01560 }
01561 }
01562 }
01563
01564
01565 if( thermal.thist > 1e9 )
01566 {
01567
01568
01569
01570 if( thermal.thist > 1.0001e10 )
01571 {
01572 sprintf( chLine, " W-Electrons were relativistic; High TE=%.2e",
01573 thermal.thist );
01574 warnin(chLine);
01575 }
01576 else
01577 {
01578 sprintf( chLine, " C-Electrons were mildly relativistic; High TE=%.2e",
01579 thermal.thist );
01580 caunin(chLine);
01581 }
01582 }
01583
01584
01585 rate = timesc.TimeErode*2e-26;
01586 if( rate > 1e-35 )
01587 {
01588
01589
01590
01591 ts = (1./rate)/3e7;
01592 if( ts < 1e3 )
01593 {
01594 sprintf( chLine, " !Timescale-photoerosion of Fe=%.2e yr",
01595 ts );
01596 bangin(chLine);
01597 }
01598 else if( ts < 1e9 )
01599 {
01600 sprintf( chLine, " Timescale-photoerosion of Fe=%.2e yr",
01601 ts );
01602 notein(chLine);
01603 }
01604 }
01605
01606
01607 comfrc = rfield.comtot/SDIV(thermal.power);
01608 if( comfrc > 0.01 )
01609 {
01610 sprintf( chLine, " Compton heating was %5.1f%% of the total.",
01611 comfrc*100. );
01612 notein(chLine);
01613 }
01614
01615
01616 if( comfrc > 0.01 && rfield.cinrat > 0.05 )
01617 {
01618 sprintf( chLine,
01619 " !Induced Compton heating was %.2e of the total Compton heating.",
01620 rfield.cinrat );
01621 bangin(chLine);
01622 }
01623
01624
01625 if( timesc.tcmptn > 5e17 )
01626 {
01627 if( comfrc > 0.05 )
01628 {
01629 sprintf( chLine,
01630 " C-Compton cooling is significant and the equilibrium timescale (%.2e s) is longer than the Hubble time.",
01631 timesc.tcmptn );
01632 caunin(chLine);
01633 }
01634 else
01635 {
01636 sprintf( chLine,
01637 " Compton cooling equilibrium timescale (%.2e s) is longer than Hubble time.",
01638 timesc.tcmptn );
01639 notein(chLine);
01640 }
01641 }
01642
01643 if( timesc.time_therm_long > 5e17 )
01644 {
01645 sprintf( chLine,
01646 " C-Thermal equilibrium timescale, %.2e s, longer than Hubble time; this cloud is not time-steady.",
01647 timesc.time_therm_long );
01648 caunin(chLine);
01649 }
01650
01651
01652
01653
01654 if( log10(radius.depth) > colden.rjnmin )
01655 {
01656
01657 aj = pow(10.,colden.ajmmin - log10(SOLAR_MASS));
01658 if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
01659 {
01660 sprintf( chLine,
01661 " C-Cloud thicker than smallest Jeans length=%8.2ecm; stability problems? (smallest Jeans mass=%8.2eMo)",
01662 pow((realnum)10.f,colden.rjnmin), aj );
01663 caunin(chLine);
01664 }
01665 else
01666 {
01667 sprintf( chLine,
01668 " Cloud thicker than smallest Jeans length=%8.2ecm; stability problems? (smallest Jeans mass=%8.2eMo)",
01669 pow((realnum)10.f,colden.rjnmin), aj );
01670 notein(chLine);
01671 }
01672 }
01673
01674
01675 for( size_t nd=0; nd < gv.bin.size(); nd++ )
01676 {
01677 if( gv.bin[nd]->TeGrainMax > gv.bin[nd]->Tsublimat )
01678 {
01679 sprintf( chLine,
01680 " W-Maximum temperature of grain%-12.12s was %.2eK, above its sublimation temperature, %.2eK.",
01681 gv.bin[nd]->chDstLab, gv.bin[nd]->TeGrainMax,
01682 gv.bin[nd]->Tsublimat );
01683 warnin(chLine);
01684 }
01685 else if( gv.bin[nd]->TeGrainMax > gv.bin[nd]->Tsublimat* 0.9 )
01686 {
01687 sprintf( chLine,
01688 " C-Maximum temperature of grain%-12.12s was %.2eK, near its sublimation temperature, %.2eK.",
01689 gv.bin[nd]->chDstLab, gv.bin[nd]->TeGrainMax,
01690 gv.bin[nd]->Tsublimat );
01691 caunin(chLine);
01692 }
01693 }
01694
01695 if( gv.lgNegGrnDrg )
01696 {
01697 sprintf( chLine, " !Grain drag force <0." );
01698 bangin(chLine);
01699 }
01700
01701
01702 if( gv.GrnElecDonateMax > 0.05 )
01703 {
01704 sprintf( chLine,
01705 " !Grains donated %5.1f%% of the total electrons in some regions.",
01706 gv.GrnElecDonateMax*100. );
01707 bangin(chLine);
01708 }
01709 else if( gv.GrnElecDonateMax > 0.005 )
01710 {
01711 sprintf( chLine,
01712 " Grains donated %5.1f%% of the total electrons in some regions.",
01713 gv.GrnElecDonateMax*100. );
01714 notein(chLine);
01715 }
01716
01717
01718 if( gv.GrnElecHoldMax > 0.05 )
01719 {
01720 sprintf( chLine,
01721 " !Grains contained %5.1f%% of the total electrons in some regions.",
01722 gv.GrnElecHoldMax*100. );
01723 bangin(chLine);
01724 }
01725 else if( gv.GrnElecHoldMax > 0.005 )
01726 {
01727 sprintf( chLine,
01728 " Grains contained %5.1f%% of the total electrons in some regions.",
01729 gv.GrnElecHoldMax*100. );
01730 notein(chLine);
01731 }
01732
01733
01734 if( gv.dphmax > 0.5 )
01735 {
01736 sprintf( chLine,
01737 " !Local grain-gas photoelectric heating rate reached %5.1f%% of the total.",
01738 gv.dphmax*100. );
01739 bangin(chLine);
01740 }
01741 else if( gv.dphmax > 0.05 )
01742 {
01743 sprintf( chLine,
01744 " Local grain-gas photoelectric heating rate reached %5.1f%% of the total.",
01745 gv.dphmax*100. );
01746 notein(chLine);
01747 }
01748
01749 if( gv.TotalDustHeat/SDIV(thermal.power) > 0.01 )
01750 {
01751 sprintf( chLine,
01752 " Global grain photoelectric heating of gas was%5.1f%% of the total.",
01753 gv.TotalDustHeat/thermal.power*100. );
01754 notein(chLine);
01755 if( gv.TotalDustHeat/thermal.power > 0.25 )
01756 {
01757 sprintf( chLine,
01758 " !Grain photoelectric heating is VERY important." );
01759 bangin(chLine);
01760 }
01761 }
01762
01763
01764 if( gv.dclmax > 0.05 )
01765 {
01766 sprintf( chLine,
01767 " Local grain-gas cooling of gas rate reached %5.1f%% of the total.",
01768 gv.dclmax*100. );
01769 notein(chLine);
01770 }
01771
01772
01773 if( h2.renorm_max > 1.05 )
01774 {
01775 if( h2.renorm_max > 1.2 )
01776 {
01777 sprintf( chLine,
01778 " !The large H2 molecule - main chemistry network renormalization factor reached %.2f.",
01779 h2.renorm_max);
01780 bangin(chLine);
01781 }
01782 else
01783 {
01784 sprintf( chLine,
01785 " The large H2 molecule - main chemistry network renormalization factor reached %.2f.",
01786 h2.renorm_max);
01787 notein(chLine);
01788 }
01789 }
01790 if( h2.renorm_min < 0.95 )
01791 {
01792 if( h2.renorm_min < 0.8 )
01793 {
01794 sprintf( chLine,
01795 " !The large H2 molecule - main chemistry network renormalization factor reached %.2f.",
01796 h2.renorm_min);
01797 bangin(chLine);
01798 }
01799 else
01800 {
01801 sprintf( chLine,
01802 " The large H2 molecule - main chemistry network renormalization factor reached %.2f.",
01803 h2.renorm_min);
01804 notein(chLine);
01805 }
01806 }
01807
01808
01809 if( hmi.h2pmax > 0.10 )
01810 {
01811 sprintf( chLine,
01812 " !The local H2+ photodissociation heating rate reached %5.1f%% of the total heating.",
01813 hmi.h2pmax*100. );
01814 bangin(chLine);
01815 }
01816
01817 else if( hmi.h2pmax > 0.01 )
01818 {
01819 sprintf( chLine,
01820 " The local H2+ photodissociation heating rate reached %.1f%% of the total heating.",
01821 hmi.h2pmax*100. );
01822 notein(chLine);
01823 }
01824
01825
01826 if( hmi.h2dfrc > 0.1 )
01827 {
01828 sprintf( chLine,
01829 " !The local H2 photodissociation heating rate reached %.1f%% of the total heating.",
01830 hmi.h2dfrc*100. );
01831 bangin(chLine);
01832 }
01833 else if( hmi.h2dfrc > 0.01 )
01834 {
01835 sprintf( chLine,
01836 " The local H2 photodissociation heating rate reached %.1f%% of the total heating.",
01837 hmi.h2dfrc*100. );
01838 notein(chLine);
01839 }
01840
01841
01842 if( hmi.h2line_cool_frac > 0.1 )
01843 {
01844 sprintf( chLine,
01845 " !The local H2 cooling rate reached %.1f%% of the local cooling.",
01846 hmi.h2line_cool_frac*100. );
01847 bangin(chLine);
01848 }
01849 else if( hmi.h2line_cool_frac > 0.01 )
01850 {
01851 sprintf( chLine,
01852 " The local H2 cooling rate reached %.1f%% of the local cooling.",
01853 hmi.h2line_cool_frac*100. );
01854 notein(chLine);
01855 }
01856
01857 if( hmi.h2dtot/SDIV(thermal.power) > 0.01 )
01858 {
01859 sprintf( chLine,
01860 " Global H2 photodissociation heating of gas was %.1f%% of the total heating.",
01861 hmi.h2dtot/thermal.power*100. );
01862 notein(chLine);
01863 if( hmi.h2dtot/thermal.power > 0.25 )
01864 {
01865 sprintf( chLine, " H2 photodissociation heating is VERY important." );
01866 notein(chLine);
01867 }
01868 }
01869
01870
01871 if( co.codfrc > 0.25 )
01872 {
01873 sprintf( chLine,
01874 " !Local CO photodissociation heating rate reached %.1f%% of the total.",
01875 co.codfrc*100. );
01876 bangin(chLine);
01877 }
01878 else if( co.codfrc > 0.05 )
01879 {
01880 sprintf( chLine,
01881 " Local CO photodissociation heating rate reached %.1f%% of the total.",
01882 co.codfrc*100. );
01883 notein(chLine);
01884 }
01885
01886 if( co.codtot/SDIV(thermal.power) > 0.01 )
01887 {
01888 sprintf( chLine,
01889 " Global CO photodissociation heating of gas was %.1f%% of the total.",
01890 co.codtot/thermal.power*100. );
01891 notein(chLine);
01892 if( co.codtot/thermal.power > 0.25 )
01893 {
01894 sprintf( chLine, " CO photodissociation heating is VERY important." );
01895 notein(chLine);
01896 }
01897 }
01898
01899 if( thermal.lgEdnGTcm )
01900 {
01901 sprintf( chLine,
01902 " Energy density of radiation field was greater than the Compton temperature. Is this physical?" );
01903 notein(chLine);
01904 }
01905
01906
01907 if( hydro.cintot/SDIV(thermal.power) > 0.01 )
01908 {
01909 sprintf( chLine, " Induced recombination cooling was %.1f%% of the total.",
01910 hydro.cintot/thermal.power*100. );
01911 notein(chLine);
01912 }
01913
01914
01915 if( thermal.FreeFreeTotHeat/SDIV(thermal.power) > 0.1 )
01916 {
01917 sprintf( chLine, " !Free-free heating was %.1f%% of the total.",
01918 thermal.FreeFreeTotHeat/thermal.power*100. );
01919 bangin(chLine);
01920 }
01921 else if( thermal.FreeFreeTotHeat/SDIV(thermal.power) > 0.01 )
01922 {
01923 sprintf( chLine, " Free-free heating was %.1f%% of the total.",
01924 thermal.FreeFreeTotHeat/thermal.power*100. );
01925 notein(chLine);
01926 }
01927
01928
01929 if( hmi.hmitot/SDIV(thermal.power) > 0.01 )
01930 {
01931 sprintf( chLine, " H- absorption heating was %.1f%% of the total.",
01932 hmi.hmitot/SDIV(thermal.power)*100. );
01933 notein(chLine);
01934 }
01935
01936
01937 if( mole_global.lgH2Ozer )
01938 {
01939 sprintf( chLine, " Water destruction rate zero." );
01940 notein(chLine);
01941 }
01942
01943
01944 if( atoms.nNegOI > 0 )
01945 {
01946 sprintf( chLine, " C-O I negative level populations %ld times.",
01947 atoms.nNegOI );
01948 caunin(chLine);
01949 }
01950
01951
01952
01953 small = 0.;
01954 imas = 0;
01955 isav = 0;
01956 j = 0;
01957 for( nelem=0; nelem<LIMELM; ++nelem )
01958 {
01959 if( dense.lgElmtOn[nelem] )
01960 {
01961
01962 for( ipLo=ipH2p; ipLo < (iso_sp[ipH_LIKE][nelem].numLevels_local - 1); ipLo++ )
01963 {
01964
01965 for( ipHi=ipLo + 1; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_local; ipHi++ )
01966 {
01967 if( iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
01968 continue;
01969
01970 if( iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauIn() < (realnum)small )
01971 {
01972 small = iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauIn();
01973 imas = ipHi;
01974 j = ipLo;
01975 isav = nelem;
01976 }
01977 }
01978 }
01979 }
01980 }
01981
01982 if( small < -0.05 )
01983 {
01984 sprintf( chLine,
01985 " !Some hydrogenic lines mased, species was %2s%2ld, smallest tau was %.2e, transition %li-%li",
01986 elementnames.chElementSym[isav],
01987 isav+1,small, imas , j );
01988 bangin(chLine);
01989 }
01990
01991
01992 if( opac.lgOpacNeg )
01993 {
01994 sprintf( chLine, " !Some opacities were negative - the SET NEGOPC command will save which ones." );
01995 bangin(chLine);
01996 }
01997
01998
01999 small = 0.;
02000 imas = 0;
02001 isav = 0;
02002 for( nelem=0; nelem<LIMELM; ++nelem )
02003 {
02004 if( dense.lgElmtOn[nelem] )
02005 {
02006
02007 for( i=0; i < iso_sp[ipH_LIKE][nelem].numLevels_local; i++ )
02008 {
02009 if( opac.TauAbsGeo[0][iso_sp[ipH_LIKE][nelem].fb[i].ipIsoLevNIonCon-1] < -0.001 )
02010 {
02011 small = MIN2(small,(double)opac.TauAbsGeo[0][iso_sp[ipH_LIKE][nelem].fb[i].ipIsoLevNIonCon-1]);
02012 imas = i;
02013 isav = nelem;
02014 }
02015 }
02016 }
02017 }
02018
02019 if( small < -0.05 )
02020 {
02021 sprintf( chLine, " !Some hydrogenic (%2s%2ld) continua optical depths were negative; smallest=%.2e level=%3ld",
02022 elementnames.chElementSym[isav],
02023 isav+1,
02024 small, imas );
02025 bangin(chLine);
02026 }
02027
02028
02029 nneg = 0;
02030 tauneg = 0.;
02031 freqn = 0.;
02032 for( i=0; i < rfield.nflux; i++ )
02033 {
02034 if( opac.TauAbsGeo[0][i] < -0.001 )
02035 {
02036 nneg += 1;
02037
02038 if( nneg == 1 )
02039 freqn = rfield.anu[i];
02040 tauneg = MIN2(tauneg,(double)opac.TauAbsGeo[0][i]);
02041 }
02042 }
02043
02044 if( nneg > 0 )
02045 {
02046 sprintf( chLine, " !Some continuous optical depths <0. The lowest freq was %.3e Ryd, and a total of%4ld",
02047 freqn, nneg );
02048 bangin(chLine);
02049 sprintf( chLine, " !The smallest optical depth was %.2e",
02050 tauneg );
02051 bangin(chLine);
02052 }
02053
02054
02055 if( opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1] > 0.05 )
02056 {
02057 sprintf( chLine, " The Balmer continuum optical depth was %.2e.",
02058 opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1] );
02059 notein(chLine);
02060 }
02061
02062
02063 if( opac.stimax[0] > 0.02 && opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] > 0.2 )
02064 {
02065 sprintf( chLine, " The Lyman continuum stimulated emission correction to optical depths reached %.2e.",
02066 opac.stimax[0] );
02067 notein(chLine);
02068 }
02069 else if( opac.stimax[1] > 0.02 && opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1] > 0.1 )
02070 {
02071 sprintf( chLine, " The Balmer continuum stimulated emission correction to optical depths reached %.2e.",
02072 opac.stimax[1] );
02073 notein(chLine);
02074 }
02075
02076
02077 if( opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[3].ipIsoLevNIonCon-1] > 0.2 )
02078 {
02079 sprintf( chLine,
02080 " The Paschen continuum optical depth was %.2e.",
02081 opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[3].ipIsoLevNIonCon-1] );
02082 notein(chLine);
02083 }
02084
02085
02086 if( opac.TauAbsGeo[0][0] > 1. )
02087 {
02088 sprintf( chLine,
02089 " The continuum optical depth at the lowest energy considered (%.3e Ryd) was %.3e.",
02090 rfield.anu[0], opac.TauAbsGeo[0][0] );
02091 notein(chLine);
02092 }
02093
02094
02095
02096 if( colden.colden[ipCOL_H0]*7e-24 > 0.01 )
02097 {
02098 sprintf( chLine,
02099 " The optical depth to Rayleigh scattering at 1300A is %.2e",
02100 colden.colden[ipCOL_H0]*6.71e-24 );
02101 notein(chLine);
02102 }
02103
02104 if( colden.colden[ipCOL_H2p]*7e-18 > 0.1 )
02105 {
02106 sprintf( chLine,
02107 " !The optical depth to the H2+ molecular ion is %.2e",
02108 colden.colden[ipCOL_H2p]*7e-18 );
02109 bangin(chLine);
02110 }
02111 else if( colden.colden[ipCOL_H2p]*7e-18 > 0.01 )
02112 {
02113 sprintf( chLine,
02114 " The optical depth to the H2+ molecular ion is %.2e",
02115 colden.colden[ipCOL_H2p]*7e-18 );
02116 notein(chLine);
02117 }
02118
02119
02120 if( opac.thmin > 0.1 )
02121 {
02122 sprintf( chLine,
02123 " !Optical depth to negative hydrogen ion is %.2e",
02124 opac.thmin );
02125 bangin(chLine);
02126 }
02127 else if( opac.thmin > 0.01 )
02128 {
02129 sprintf( chLine,
02130 " Optical depth to negative hydrogen ion is %.2e",
02131 opac.thmin );
02132 notein(chLine);
02133 }
02134
02135
02136
02137
02138 if( ionbal.ifail > 0 && ionbal.ifail <= 10 )
02139 {
02140 sprintf( chLine,
02141 " 3 body recombination coefficient outside range %ld", ionbal.ifail );
02142 notein(chLine);
02143 }
02144 else if( ionbal.ifail > 10 )
02145 {
02146 sprintf( chLine,
02147 " C-3 body recombination coefficient outside range %ld", ionbal.ifail );
02148 caunin(chLine);
02149 }
02150
02151
02152 if( phycon.TEnerDen < 2.6 )
02153 {
02154 sprintf( chLine,
02155 " !Incident radiation field energy density is less than 2.7K. Add background with CMB command." );
02156 bangin(chLine);
02157 }
02158
02159
02160 if( !rfield.lgCMB_set )
02161 {
02162 sprintf( chLine,
02163 " !The CMB was not included. This is added with the CMB command." );
02164 bangin(chLine);
02165 }
02166
02167
02168 if( rfield.lgHabing )
02169 {
02170 sprintf( chLine,
02171 " !The intensity of the incident radiation field is less than 10 times the Habing diffuse ISM field. Is this OK?" );
02172 bangin(chLine);
02173 sprintf( chLine,
02174 " ! Consider adding diffuse ISM emission with TABLE ISM command." );
02175 bangin(chLine);
02176 }
02177
02178
02179
02180
02181 if( dense.lgElmtOn[ipOXYGEN] && dense.lgElmtOn[ipCARBON] )
02182 {
02183 if( dense.gas_phase[ipCARBON]/dense.gas_phase[ipOXYGEN] > 1. )
02184 {
02185 sprintf( chLine, " !The C/O abundance ratio, %.1f, is greater than unity. The chemistry will be carbon dominated.",
02186 dense.gas_phase[ipCARBON]/dense.gas_phase[ipOXYGEN] );
02187 bangin(chLine);
02188 }
02189 }
02190
02191 lgLots_of_moles = false;
02192 lgLotsSolids = false;
02193
02194 for( i=0; i<mole_global.num_calc; ++i )
02195 {
02196 if( mole.species[i].location == NULL && ( mole_global.list[i]->parentLabel.empty() || mole_global.list[i]->charge<0 ) )
02197 {
02198 if( mole.species[i].xFracLim > 0.1 )
02199 {
02200 sprintf( chLine, " !The fraction of %s in %s reached %.1f%% at some point in the cloud.",
02201 elementnames.chElementSym[mole.species[i].nAtomLim],
02202 mole_global.list[i]->label.c_str(),
02203 mole.species[i].xFracLim*100. );
02204 bangin(chLine);
02205 lgLots_of_moles = true;
02206
02207 if( !mole_global.list[i]->lgGas_Phase )
02208 lgLotsSolids = true;
02209 }
02210 else if( mole.species[i].xFracLim>0.01 )
02211 {
02212 sprintf( chLine, " The fraction of %s in %s reached %.2f%% at some point in the cloud.",
02213 elementnames.chElementSym[mole.species[i].nAtomLim],
02214 mole_global.list[i]->label.c_str(),
02215 mole.species[i].xFracLim*100. );
02216 notein(chLine);
02217 lgLots_of_moles = true;
02218
02219 if( !mole_global.list[i]->lgGas_Phase )
02220 lgLotsSolids = true;
02221 }
02222 else if( mole.species[i].xFracLim > 1e-3 )
02223 {
02224 sprintf( chLine, " The fraction of %s in %s reached %.3f%% at some point in the cloud.",
02225 elementnames.chElementSym[mole.species[i].nAtomLim],
02226 mole_global.list[i]->label.c_str(),
02227 mole.species[i].xFracLim*100. );
02228 notein(chLine);
02229
02230 if( !mole_global.list[i]->lgGas_Phase )
02231 lgLotsSolids = true;
02232 }
02233 }
02234 }
02235
02236
02237 if( lgLots_of_moles )
02238 {
02239
02240 for(nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
02241 {
02242
02243 if( !dense.lgElmtOn[nelem] )
02244 {
02245
02246 sprintf( chLine,
02247 " C-Molecules are important, but %s, part of the chemistry network, is turned off.",
02248 elementnames.chElementName[nelem] );
02249 caunin(chLine);
02250 }
02251 # if 0
02252
02253 for( i=NUM_HEAVY_MOLEC+NUM_ELEMENTS; i<NUM_COMOLE_CALC; ++i )
02254 {
02255 if( nelem==mole_global.list[i].nelem_den )
02256 {
02257
02258 sprintf( chLine,
02259 " C-Molecules are important, but %s, part of the chemistry network, is turned off.",
02260 elementnames.chElementName[nelem] );
02261 caunin(chLine);
02262 }
02263 }
02264 # endif
02265 }
02266 }
02267
02268
02269
02270 if( lgLotsSolids )
02271 {
02272 sprintf( chLine, " !A significant amount of molecules condensed onto grain surfaces." );
02273 bangin(chLine);
02274 sprintf( chLine, " !These are the molecular species with \"grn\" above." );
02275 bangin(chLine);
02276 }
02277
02278
02279 if( rfield.EnergyBremsThin > 0.09 )
02280 {
02281 sprintf( chLine, " !The cloud is optically thick at optical wavelengths, extending to %.3e Ryd =%.3eA",
02282 rfield.EnergyBremsThin, RYDLAM/rfield.EnergyBremsThin );
02283 bangin(chLine);
02284 }
02285 else if( rfield.EnergyBremsThin > 0.009 )
02286 {
02287 sprintf( chLine, " The continuum of the computed structure may be optically thick in the near infrared." );
02288 notein(chLine);
02289 }
02290
02291
02292 if( radius.Radius > 1e23 && radius.Radius/radius.rinner > 10. )
02293 {
02294 sprintf( chLine, " Is an outer radius of %.2e reasonable?",
02295 radius.Radius );
02296 notein(chLine);
02297 }
02298
02299
02300 if( rt.lgMaserCapHit )
02301 {
02302 sprintf( chLine, " Laser maser optical depths capped in RT_line_one_tauinc." );
02303 notein(chLine);
02304 }
02305
02306
02307 if( rt.lgMaserSetDR )
02308 {
02309 sprintf( chLine, " !Line maser set zone thickness in some zones." );
02310 bangin(chLine);
02311 }
02312
02313
02314
02315 if( (pressure.lgPradCap && (strcmp(dense.chDenseLaw,"CPRE") == 0)) &&
02316 pressure.lgPres_radiation_ON )
02317 {
02318 sprintf( chLine, " Radiation pressure kept below gas pressure on this iteration." );
02319 notein(chLine);
02320 }
02321
02322 if( pressure.RadBetaMax > 0.25 )
02323 {
02324 if( pressure.ipPradMax_line == 0 )
02325 {
02326 sprintf( chLine,
02327 " !The ratio of radiation to gas pressure reached %.2e at zone %li. Caused by Lyman alpha.",
02328 pressure.RadBetaMax,
02329 pressure.ipPradMax_nzone);
02330 bangin(chLine);
02331 }
02332 else
02333 {
02334 sprintf( chLine,
02335 " !The ratio of radiation to gas pressure reached %.2e at zone %li. "
02336 "Caused by line number %ld, label %s",
02337 pressure.RadBetaMax,
02338 pressure.ipPradMax_nzone,
02339 pressure.ipPradMax_line,
02340 pressure.chLineRadPres );
02341 bangin(chLine);
02342 }
02343 }
02344
02345 else if( pressure.RadBetaMax > 0.025 )
02346 {
02347 if( pressure.ipPradMax_line == 0 )
02348 {
02349 sprintf( chLine,
02350 " The ratio of radiation to gas pressure reached %.2e at zone %li. Caused by Lyman alpha.",
02351 pressure.RadBetaMax,
02352 pressure.ipPradMax_nzone);
02353 notein(chLine);
02354 }
02355 else
02356 {
02357 sprintf( chLine,
02358 " The ratio of radiation to gas pressure reached %.2e at zone %li. "
02359 "Caused by line number %ld, label %s",
02360 pressure.RadBetaMax,
02361 pressure.ipPradMax_nzone,
02362 pressure.ipPradMax_line,
02363 pressure.chLineRadPres );
02364 notein(chLine);
02365 }
02366 }
02367
02368 if( opac.telec >= 5. )
02369 {
02370 sprintf( chLine, " W-The model is optically thick to electron "
02371 "scattering; tau=%.2e Cloudy is NOT intended for this regime.",
02372 opac.telec );
02373 warnin(chLine);
02374 }
02375 else if( opac.telec > 2.0 )
02376 {
02377 sprintf( chLine, " C-The model is moderately optically thick to electron scattering; tau=%.1f",
02378 opac.telec );
02379 caunin(chLine);
02380 }
02381 else if( opac.telec > 0.1 )
02382 {
02383 sprintf( chLine, " !The model has modest optical depth to electron scattering; tau=%.2f",
02384 opac.telec );
02385 bangin(chLine);
02386 }
02387 else if( opac.telec > 0.01 )
02388 {
02389 sprintf( chLine, " The optical depth to electron scattering is %.3f",
02390 opac.telec );
02391 notein(chLine);
02392 }
02393
02394
02395 if( HFLines[0].Emis().TauIn() > 0.5 )
02396 {
02397 sprintf( chLine, " !The optical depth in the H I 21 cm line is %.2e",HFLines[0].Emis().TauIn() );
02398 bangin(chLine);
02399 }
02400
02401
02402
02403 if( nWindLine==0 )
02404 {
02405
02406 sprintf( chLine, " !The level2 lines are disabled. UV pumping of excited levels within ground terms is not treated." );
02407 bangin(chLine);
02408 }
02409
02410
02411 for( nelem=0; nelem < LIMELM; nelem++ )
02412 {
02413 if( dense.lgElmtOn[nelem] && !dynamics.lgTimeDependentStatic )
02414 {
02415 if( iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauIn() > 0.2 )
02416 {
02417 differ = fabs(1.-iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauIn()*
02418 rt.DoubleTau/iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauTot())*100.;
02419
02420
02421
02422
02423 if( ((iterations.lgLastIt && iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauIn() > 0.8) &&
02424 differ > 20.) && wind.lgStatic() )
02425 {
02426 sprintf( chLine,
02427 " C-This is the last iteration and %2s%2ld Bal(a) optical depth"
02428 " changed by%6.1f%% (was %.2e). Try another iteration.",
02429 elementnames.chElementSym[nelem],
02430 nelem+1, differ,
02431 iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauTot() );
02432 caunin(chLine);
02433 iterations.lgIterAgain = true;
02434 }
02435
02436
02437 if( iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn() > 0. )
02438 {
02439 differ = fabs(1.-iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn()*
02440 rt.DoubleTau/iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauTot())*100.;
02441
02442
02443
02444
02445 if( ((iterations.lgLastIt && iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn() > 0.8) &&
02446 differ > 25.) && wind.lgStatic() )
02447 {
02448 sprintf( chLine,
02449 " C-This is the last iteration and %2s%2ld Ly(a) optical depth"
02450 " changed by%7.0f%% (was %.2e). Try another iteration.",
02451 elementnames.chElementSym[nelem],
02452 nelem+1,differ, iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauTot() );
02453 caunin(chLine);
02454 iterations.lgIterAgain = true;
02455 }
02456 }
02457 }
02458 }
02459 }
02460
02461
02462 if( radius.Radius/radius.rinner > 2. && !geometry.lgSphere )
02463 {
02464 sprintf( chLine, " C-R(out)/R(in)=%.2e and SPHERE was not set.",
02465 radius.Radius/radius.rinner );
02466 caunin(chLine);
02467 }
02468
02469
02470 if( iterations.lgLastIt && !opac.lgCaseB )
02471 {
02472
02473
02474 if( rfield.nflux > iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon )
02475 {
02476 if( opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] < 2. &&
02477 opac.TauAbsGeo[1][iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] > 2. )
02478 {
02479 sprintf( chLine, " C-The H Lyman continuum is thin, and I assumed"
02480 " that it was thick. Try another iteration." );
02481 caunin(chLine);
02482 iterations.lgIterAgain = true;
02483 }
02484 }
02485
02486
02487 if( rfield.nflux > iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon && dense.lgElmtOn[ipHELIUM] )
02488 {
02489 if( (opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1] < 2. &&
02490 opac.TauAbsGeo[1][iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1] > 2.) )
02491 {
02492 sprintf( chLine,
02493 " C-The He II continuum is thin and I assumed that it was thick."
02494 " Try another iteration." );
02495 caunin(chLine);
02496 iterations.lgIterAgain = true;
02497 }
02498 }
02499
02500 if( rfield.nflux > iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon && dense.lgElmtOn[ipHELIUM] )
02501 {
02502 if( (opac.TauAbsGeo[0][iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1] < 2. &&
02503 opac.TauAbsGeo[1][iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1] > 2.) )
02504 {
02505 sprintf( chLine,
02506 " C-The He I continuum is thin and I assumed that it was thick."
02507 " Try another iteration." );
02508 caunin(chLine);
02509 iterations.lgIterAgain = true;
02510 }
02511 }
02512 }
02513
02514
02515 if( iteration > 1 )
02516 {
02517 if( colden.colden_old[ipCOL_HTOT] <= 0. )
02518 {
02519 fprintf( ioQQQ, " colden_old is insane in PrtComment.\n" );
02520 ShowMe();
02521 cdEXIT(EXIT_FAILURE);
02522 }
02523
02524 differ = fabs(1.-colden.colden[ipCOL_HTOT]/
02525 colden.colden_old[ipCOL_HTOT]);
02526
02527 if( differ > 0.1 && differ <= 0.3 )
02528 {
02529 sprintf( chLine,
02530 " The H column density changed by %.2e%% between this and previous iteration.",
02531 differ*100. );
02532 notein(chLine);
02533 }
02534
02535 else if( differ > 0.3 )
02536 {
02537 if( iterations.lgLastIt )
02538 {
02539 sprintf( chLine,
02540 " C-The H column density changed by %.2e%% and this is the last iteration. What happened?",
02541 differ*100. );
02542 caunin(chLine);
02543 }
02544 else
02545 {
02546 sprintf( chLine,
02547 " !The H column density changed by %.2e%% What happened?",
02548 differ*100. );
02549 bangin(chLine);
02550 }
02551 }
02552
02553
02554 if( (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])/SDIV(colden.colden[ipCOL_HTOT]) > 1e-5 )
02555 {
02556 differ = fabs(1.-colden.colden[ipCOL_H2g]/
02557 SDIV(colden.colden_old[ipCOL_H2g]));
02558
02559 if( differ > 0.1 && differ <= 0.3 )
02560 {
02561 sprintf( chLine,
02562 " The H2 column density changed by %.2e%% between this and previous iteration.",
02563 differ*100. );
02564 notein(chLine);
02565 }
02566
02567 else if( differ > 0.3 )
02568 {
02569 if( iterations.lgLastIt )
02570 {
02571 sprintf( chLine,
02572 " C-The H2 column density changed by %.2e%% and this is the last iteration. What happened?",
02573 differ*100. );
02574 caunin(chLine);
02575 }
02576 else
02577 {
02578 sprintf( chLine,
02579 " !The H2 column density changed by %.2e%% What happened?",
02580 differ*100. );
02581 bangin(chLine);
02582 }
02583 }
02584 }
02585 }
02586
02587
02588 differ = fabs(1.-iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn()/
02589 SDIV(iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot()))*100.;
02590
02591 if( iterations.lgLastIt && (pressure.RadBetaMax > 0.1) &&
02592 (differ > 50.) && (pressure.ipPradMax_line == 1) && (pressure.lgPres_radiation_ON) &&
02593 wind.lgStatic() )
02594 {
02595 sprintf( chLine, " C-This is the last iteration, radiation pressure was significant, and the L-a optical depth changed by %7.2f%% (was %.2e)",
02596 differ, iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot() );
02597 caunin(chLine);
02598 }
02599
02600
02601
02602 if( iterations.lgLastIt &&
02603 ( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot() * 1.02 -
02604 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn() ) < 0. )
02605 {
02606 sprintf( chLine, " C-The Lya optical depth scale was overrun and this is the last iteration - Tspin(21 cm) is not valid." );
02607 caunin(chLine);
02608 sprintf( chLine, " C-Another iteration is needed for Tspin(21 cm) to be valid." );
02609 caunin(chLine);
02610 }
02611
02612
02613 if( pressure.lgPradDen )
02614 {
02615 sprintf( chLine, " Line radiation pressure capped by thermalization length." );
02616 notein(chLine);
02617 }
02618
02619
02620 nline = MIN2(conv.nTeFail,10);
02621 if( conv.nTeFail != 0 )
02622 {
02623 long int _o;
02624 if( conv.failmx < 0.1 )
02625 {
02626 _o = sprintf( chLine, " There were %ld minor temperature failures. zones:",
02627 conv.nTeFail );
02628
02629
02630 for( i=0; i < nline; i++ )
02631 {
02632 _o += sprintf( chLine+_o, " %ld", conv.ifailz[i] );
02633 }
02634 notein(chLine);
02635 }
02636 else
02637 {
02638 _o = sprintf( chLine,
02639 " !There were %ld temperature failures, and some were large. The largest was %.1f%%. What happened?",
02640 conv.nTeFail, conv.failmx*100. );
02641 bangin(chLine);
02642
02643
02644
02645 _o = sprintf( chLine , " !The zones were" );
02646 for( i=0; i < nline; i++ )
02647 {
02648 _o += sprintf( chLine+_o, " %ld", conv.ifailz[i] );
02649 }
02650 bangin(chLine);
02651
02652 if( struc.testr[0] > 8e4 && phycon.te < 6e5 )
02653 {
02654 sprintf( chLine, " !I think they may have been caused by the change from hot to nebular gas phase. The physics of this is unclear." );
02655 bangin(chLine);
02656 }
02657 }
02658 }
02659
02660
02661 big_ion_jump = 0.;
02662 j = 0;
02663 for( i=1; i < nzone; i++ )
02664 {
02665 big = fabs(1.-struc.testr[i-1]/struc.testr[i]);
02666 if( big > big_ion_jump )
02667 {
02668 j = i;
02669 big_ion_jump = big;
02670 }
02671 }
02672
02673 if( big_ion_jump > 0.2 )
02674 {
02675
02676 if( j < 1 )
02677 {
02678 fprintf( ioQQQ, " j too small big jump check\n" );
02679 ShowMe();
02680 cdEXIT(EXIT_FAILURE);
02681 }
02682
02683 if( big_ion_jump > 0.4 )
02684 {
02685 sprintf( chLine, " C-A temperature discontinuity occurred at zone %ld from %.2eK to %.2eK.",
02686 j, struc.testr[j-1], struc.testr[j] );
02687 caunin(chLine);
02688
02689
02690
02691
02692 if( struc.testr[j]>100. && struc.testr[j] < 1000. )
02693 {
02694 sprintf( chLine, " C-This was probably due to a thermal front." );
02695 caunin(chLine);
02696 }
02697 }
02698 else if( big_ion_jump > 0.2 )
02699 {
02700 sprintf( chLine, " !A temperature discontinuity occurred at zone %ld from %.2eK to %.2eK.",
02701 j, struc.testr[j-1], struc.testr[j] );
02702 bangin(chLine);
02703
02704
02705
02706
02707 if( struc.testr[j]>100. && struc.testr[j] < 1000. )
02708 {
02709 sprintf( chLine, " !This was probably due to a thermal front." );
02710 bangin(chLine);
02711 }
02712 }
02713 }
02714
02715
02716 if( fabs(conv.BigEdenError) > conv.EdenErrorAllowed )
02717 {
02718
02719 if( fabs(conv.BigEdenError) > conv.EdenErrorAllowed*20. && dense.nzEdenBad !=
02720 nzone )
02721 {
02722 sprintf( chLine, " W-The local error in the electron density reached %.1f%% at zone %ld",
02723 conv.BigEdenError*100, dense.nzEdenBad );
02724 warnin(chLine);
02725 }
02726 else if( fabs(conv.BigEdenError) > conv.EdenErrorAllowed*5. )
02727 {
02728 sprintf( chLine, " C-The local error in the electron density reached %.1f%% at zone %ld",
02729 conv.BigEdenError*100, dense.nzEdenBad );
02730 caunin(chLine);
02731 }
02732 else
02733 {
02734 sprintf( chLine, " The local error in the electron density reached %.1f%% at zone %ld",
02735 conv.BigEdenError*100, dense.nzEdenBad );
02736 notein(chLine);
02737 }
02738 }
02739
02740
02741 big_ion_jump = 0.;
02742 j = 0;
02743 for( i=1; i < (nzone - 1); i++ )
02744 {
02745 big = fabs( (struc.testr[i-1] - struc.testr[i])/struc.testr[i] );
02746 bigm = fabs( (struc.testr[i] - struc.testr[i+1])/struc.testr[i] );
02747
02748
02749 rel = ( (struc.testr[i-1] - struc.testr[i])/struc.testr[i])*
02750 ( (struc.testr[i] - struc.testr[i+1])/struc.testr[i] );
02751
02752 if( rel < 0. && MIN2( bigm , big ) > big_ion_jump )
02753 {
02754 j = i;
02755 big_ion_jump = MIN2( bigm , big );
02756 }
02757 }
02758
02759 if( big_ion_jump > 0.1 )
02760 {
02761
02762 if( j < 1 )
02763 {
02764 fprintf( ioQQQ, " j too small bigjump2 check\n" );
02765 ShowMe();
02766 cdEXIT(EXIT_FAILURE);
02767 }
02768
02769 if( big_ion_jump > 0.3 )
02770 {
02771 sprintf( chLine,
02772 " C-A temperature oscillation occurred at zone%4ld by%5.0f%% from %.2e to %.2e to %.2e",
02773 j, big_ion_jump*100., struc.testr[j-1], struc.testr[j], struc.testr[j+1] );
02774 caunin(chLine);
02775 }
02776 else if( big_ion_jump > 0.1 )
02777 {
02778 sprintf( chLine,
02779 " !A temperature oscillation occurred at zone%4ld by%5.0f%% from %.2e to %.2e to %.2e",
02780 j, big_ion_jump*100., struc.testr[j-1], struc.testr[j], struc.testr[j+1] );
02781 bangin(chLine);
02782 }
02783 }
02784
02785
02786 if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
02787 {
02788 j = 0;
02789 big_ion_jump = 0.;
02790 for( i=1; i < (nzone - 1); i++ )
02791 {
02792 big = (struc.ednstr[i-1] - struc.ednstr[i])/struc.ednstr[i];
02793 if( fabs(big) < conv.EdenErrorAllowed )
02794 big = 0.;
02795 bigm = (struc.ednstr[i] - struc.ednstr[i+1])/struc.ednstr[i];
02796 if( fabs(bigm) < conv.EdenErrorAllowed )
02797 bigm = 0.;
02798 if( big*bigm < 0. &&
02799 fabs(struc.ednstr[i-1]-struc.ednstr[i])/struc.ednstr[i] > big_ion_jump )
02800 {
02801 j = i;
02802 big_ion_jump = fabs(struc.ednstr[i-1]-struc.ednstr[i])/
02803 struc.ednstr[i];
02804 }
02805 }
02806
02807
02808
02809 if( big_ion_jump > conv.EdenErrorAllowed*3. )
02810 {
02811 if( j < 1 )
02812 {
02813 fprintf( ioQQQ, " j too small bigjump3 check\n" );
02814 ShowMe();
02815 cdEXIT(EXIT_FAILURE);
02816 }
02817
02818 if( big_ion_jump > conv.EdenErrorAllowed*10. )
02819 {
02820 sprintf( chLine, " C-An electron density oscillation occurred at zone%4ld by%5.0f%% from %.2e to %.2e to %.2e",
02821 j, big_ion_jump*100., struc.ednstr[j-1], struc.ednstr[j],
02822 struc.ednstr[j+1] );
02823 caunin(chLine);
02824 }
02825 else if( big_ion_jump > conv.EdenErrorAllowed*3. )
02826 {
02827 sprintf( chLine, " !An electron density oscillation occurred at zone%4ld by%5.0f%% from %.2e to %.2e to %.2e",
02828 j, big_ion_jump*100., struc.ednstr[j-1], struc.ednstr[j],
02829 struc.ednstr[j+1] );
02830 bangin(chLine);
02831 }
02832 }
02833 }
02834
02835
02836
02837 prt_smooth_predictions();
02838
02839
02840
02841
02842
02843
02844
02845
02846 i = cdLine( "Unit" , 1 , &rate , &absint );
02847 ASSERT( i> 0 );
02848
02849
02850 VolComputed = LineSv[i].SumLine[0] / 1e-10;
02851
02852
02853 if( radius.Radius/radius.rinner > 1.0001 )
02854 {
02855
02856
02857
02858
02859 VolExpected = geometry.covaper*geometry.FillFac*radius.rinner/(geometry.iEmissPower+1)*
02860 ( powi( radius.Radius/radius.rinner,geometry.iEmissPower+1 ) - 1. );
02861 }
02862 else
02863 {
02864
02865
02866 VolExpected = geometry.covaper*geometry.FillFac*(radius.depth-DEPTH_OFFSET);
02867 }
02868
02869
02870 error = fabs(VolComputed - VolExpected)/SDIV(VolExpected);
02871
02872
02873
02874 if( radius.lgCylnOn || geometry.filpow!=0. )
02875 {
02876 error = 0.;
02877 }
02878
02879
02880 if( error > 0.001 && !lgAbort )
02881 {
02882 sprintf( chLine,
02883 " W-PrtComment insanity - Line unit integration did not verify \n");
02884 warnin(chLine);
02885 fprintf( ioQQQ,
02886 " PROBLEM PrtComment insanity - Line unit integration did not verify \n");
02887 fprintf( ioQQQ,
02888 " expected, derived vols were %g %g \n",
02889 VolExpected , VolComputed );
02890 fprintf( ioQQQ,
02891 " relative difference is %g, ratio is %g.\n",error,VolComputed/VolExpected);
02892 TotalInsanity();
02893 }
02894
02895
02896
02897
02898
02899
02900
02901
02902
02903
02904
02905
02906 ConComputed = rfield.ConInterOut[rfield.nflux]/ 1e-10;
02907
02908
02909 ConComputed /= ( (1. + geometry.covrt)/2. );
02910
02911
02912 if( radius.Radius/radius.rinner < 1.001 )
02913 {
02914
02915 ConExpected = (radius.depth-DEPTH_OFFSET)*geometry.FillFac;
02916 }
02917 else
02918 {
02919
02920 ConExpected = radius.rinner*geometry.FillFac * (1. - radius.rinner/radius.Radius );
02921 }
02922
02923 ASSERT( ConExpected > 0. );
02924
02925
02926 error = fabs(ConComputed - ConExpected)/ConExpected;
02927
02928
02929
02930 if( radius.lgCylnOn || geometry.filpow!=0. )
02931 {
02932 error = 0.;
02933 }
02934
02935
02936
02937
02938
02939
02940
02941 if( error > 0.001 && !lgAbort)
02942 {
02943 sprintf( chLine,
02944 " W-PrtComment insanity - Continuum unit integration did not verify \n");
02945 warnin(chLine);
02946 fprintf( ioQQQ," PROBLEM PrtComment insanity - Continuum unit integration did not verify \n");
02947 fprintf( ioQQQ," exact vol= %g, derived vol= %g relative difference is %g \n",
02948 ConExpected , ConComputed ,error);
02949 fprintf( ioQQQ," ConInterOut= %g, \n",
02950 rfield.ConInterOut[rfield.nflux]);
02951 TotalInsanity();
02952 }
02953
02954
02955 cdNwcns(&lgAbort_flag,&nw,&nc,&nn,&ns,&i,&j,&dum1,&dum2);
02956
02957 warnings.lgWarngs = ( nw > 0 );
02958 warnings.lgCautns = ( nc > 0 );
02959
02960 if( called.lgTalk )
02961 {
02962
02963 fprintf( ioQQQ, " %s\n", input.chTitle );
02964
02965 cdReasonGeo(ioQQQ);
02966
02967 cdWarnings(ioQQQ);
02968
02969 cdCautions(ioQQQ);
02970
02971 cdSurprises(ioQQQ);
02972
02973 cdNotes(ioQQQ);
02974 }
02975
02976
02977 if( lgPrnErr )
02978 {
02979 cdWarnings(ioPrnErr);
02980 }
02981
02982 if( trace.lgTrace )
02983 {
02984 fprintf( ioQQQ, " PrtComment returns.\n" );
02985 }
02986 return;
02987 }
02988
02989
02990
02991 STATIC void chkCaHeps(double *totwid)
02992 {
02993 double conca,
02994 conalog ,
02995 conhe;
02996
02997 DEBUG_ENTRY( "chkCaHeps()" );
02998
02999 *totwid = 0.;
03000
03001 if( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_local +
03002 iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local >= 6 )
03003 {
03004
03005 long ip6p = iso_sp[ipH_LIKE][ipHYDROGEN].QuantumNumbers2Index[6][1][2];
03006
03007 if( TauLines[ipT3969].Emis().TauIn() > 0. &&
03008 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ip6p,ipH2s).Emis().TauIn() > 0. )
03009 {
03010
03011 conca = pow(6.1e-5* TauLines[ipT3969].Emis().TauIn(),0.5);
03012 conalog = log((double)TauLines[ipT3969].Emis().TauIn());
03013 conalog = sqrt(MAX2(1., conalog));
03014 conca = MAX2(conalog,conca);
03015
03016 conalog = log((double)iso_sp[ipH_LIKE][ipHYDROGEN].trans(ip6p,ipH2s).Emis().TauIn());
03017 conalog = sqrt(MAX2(1.,conalog));
03018 conhe = pow(1.7e-6*iso_sp[ipH_LIKE][ipHYDROGEN].trans(ip6p,ipH2s).Emis().TauIn(),0.5);
03019 conhe = MAX2(conalog, conhe);
03020
03021 *totwid = 10.*conhe + 1.6*conca;
03022 }
03023 }
03024 return;
03025 }
03026
03027
03028 STATIC void prt_smooth_predictions(void)
03029 {
03030 long int i,
03031 nzone_oscillation,
03032 nzone_ion_jump,
03033 nzone_den_jump,
03034 nelem,
03035 ion;
03036 double BigOscillation ,
03037 big_ion_jump,
03038 big_jump,
03039 rel,
03040 big,
03041 bigm;
03042
03043 char chLine[INPUT_LINE_LENGTH];
03044
03045 DEBUG_ENTRY( "prt_smooth_predictions()" );
03046
03047
03048 nzone_oscillation = 0;
03049 nzone_ion_jump = 0;
03050
03051 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
03052 {
03053 if( dense.lgElmtOn[nelem] )
03054 {
03055 for( ion=0; ion<=nelem+1; ++ion)
03056 {
03057 BigOscillation = 0.;
03058 big_ion_jump = -15.;
03059
03060 for( i=1; i < (nzone - 1-(int)lgAbort); i++ )
03061 {
03062
03063
03064 if( struc.xIonDense[nelem][ion][i-1]/struc.gas_phase[nelem][i-1]>struc.dr_ionfrac_limit &&
03065 struc.xIonDense[nelem][ion][i ]/struc.gas_phase[nelem][i ]>struc.dr_ionfrac_limit &&
03066 struc.xIonDense[nelem][ion][i+1]/struc.gas_phase[nelem][i+1]>struc.dr_ionfrac_limit )
03067 {
03068
03069
03070 big = fabs( (struc.xIonDense[nelem][ion][i-1] - struc.xIonDense[nelem][ion][i])/struc.xIonDense[nelem][ion][i] );
03071 bigm = fabs( (struc.xIonDense[nelem][ion][i] - struc.xIonDense[nelem][ion][i+1])/struc.xIonDense[nelem][ion][i] );
03072
03073
03074 rel = ( (struc.xIonDense[nelem][ion][i-1] - struc.xIonDense[nelem][ion][i] )/struc.xIonDense[nelem][ion][i])*
03075 ( (struc.xIonDense[nelem][ion][i] - struc.xIonDense[nelem][ion][i+1])/struc.xIonDense[nelem][ion][i] );
03076
03077 if( rel < 0. && MIN2( bigm , big ) > BigOscillation )
03078 {
03079 nzone_oscillation = i;
03080 BigOscillation = MIN2( bigm , big );
03081 }
03082
03083
03084
03085
03086 rel = -log10( (struc.xIonDense[nelem][ion][i]/struc.gas_phase[nelem][i]) /
03087 (struc.xIonDense[nelem][ion][i+1]/struc.gas_phase[nelem][i+1] ) );
03088
03089 if( rel > big_ion_jump )
03090 {
03091 big_ion_jump = rel;
03092 nzone_ion_jump = i;
03093 }
03094 }
03095 }
03096
03097
03098
03099 if( BigOscillation > 0.2 )
03100 {
03101
03102 if( nzone_oscillation < 1 )
03103 {
03104 fprintf( ioQQQ, " nzone_oscillation too small bigjump2 check\n" );
03105 ShowMe();
03106 cdEXIT(EXIT_FAILURE);
03107 }
03108 if( BigOscillation > 3. )
03109 {
03110 sprintf( chLine,
03111 " W-An ionization oscillation occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
03112 nzone_oscillation,
03113 elementnames.chElementSym[nelem], ion+1,
03114 BigOscillation*100.,
03115 struc.xIonDense[nelem][ion][nzone_oscillation-1]/struc.gas_phase[nelem][nzone_oscillation-1],
03116 struc.xIonDense[nelem][ion][nzone_oscillation]/struc.gas_phase[nelem][nzone_oscillation],
03117 struc.xIonDense[nelem][ion][nzone_oscillation+1]/struc.gas_phase[nelem][nzone_oscillation+1] );
03118 warnin(chLine);
03119 }
03120
03121 else if( BigOscillation > 0.7 )
03122 {
03123 sprintf( chLine,
03124 " C-An ionization oscillation occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
03125 nzone_oscillation,
03126 elementnames.chElementSym[nelem], ion+1,
03127 BigOscillation*100.,
03128 struc.xIonDense[nelem][ion][nzone_oscillation-1]/struc.gas_phase[nelem][nzone_oscillation-1],
03129 struc.xIonDense[nelem][ion][nzone_oscillation]/struc.gas_phase[nelem][nzone_oscillation],
03130 struc.xIonDense[nelem][ion][nzone_oscillation+1]/struc.gas_phase[nelem][nzone_oscillation+1] );
03131 caunin(chLine);
03132 }
03133 else if( BigOscillation > 0.2 )
03134 {
03135 sprintf( chLine,
03136 " !An ionization oscillation occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
03137 nzone_oscillation,
03138 elementnames.chElementSym[nelem], ion+1,
03139 BigOscillation*100.,
03140 struc.xIonDense[nelem][ion][nzone_oscillation-1]/struc.gas_phase[nelem][nzone_oscillation-1],
03141 struc.xIonDense[nelem][ion][nzone_oscillation]/struc.gas_phase[nelem][nzone_oscillation],
03142 struc.xIonDense[nelem][ion][nzone_oscillation+1]/struc.gas_phase[nelem][nzone_oscillation+1] );
03143 bangin(chLine);
03144 }
03145 }
03146
03147
03148
03149 big_ion_jump = pow(10., big_ion_jump );
03150 if( big_ion_jump > 1.5 && nzone_ion_jump > 0 )
03151 {
03152 if( big_ion_jump > 10. )
03153 {
03154 sprintf( chLine,
03155 " C-An ionization jump occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
03156 nzone_ion_jump,
03157 elementnames.chElementSym[nelem], ion+1,
03158 big_ion_jump*100.,
03159 struc.xIonDense[nelem][ion][nzone_ion_jump-1]/struc.gas_phase[nelem][nzone_ion_jump-1],
03160 struc.xIonDense[nelem][ion][nzone_ion_jump]/struc.gas_phase[nelem][nzone_ion_jump],
03161 struc.xIonDense[nelem][ion][nzone_ion_jump+1]/struc.gas_phase[nelem][nzone_ion_jump+1] );
03162 caunin(chLine);
03163 }
03164 else
03165 {
03166 sprintf( chLine,
03167 " !An ionization jump occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
03168 nzone_ion_jump,
03169 elementnames.chElementSym[nelem], ion+1,
03170 big_ion_jump*100.,
03171 struc.xIonDense[nelem][ion][nzone_ion_jump-1]/struc.gas_phase[nelem][nzone_ion_jump-1],
03172 struc.xIonDense[nelem][ion][nzone_ion_jump]/struc.gas_phase[nelem][nzone_ion_jump],
03173 struc.xIonDense[nelem][ion][nzone_ion_jump+1]/struc.gas_phase[nelem][nzone_ion_jump+1] );
03174 bangin(chLine);
03175 }
03176 }
03177 }
03178 }
03179 }
03180
03181 big_jump = -15;
03182 nzone_den_jump = 0;
03183
03184
03185 for( i=1; i < (nzone - 1 - (int)lgAbort); i++ )
03186 {
03187
03188 rel = fabs(log10( struc.gas_phase[ipHYDROGEN][i] /
03189 struc.gas_phase[ipHYDROGEN][i+1] ) );
03190
03191 if( rel > big_jump )
03192 {
03193 big_jump = rel;
03194 nzone_den_jump = i;
03195 }
03196 }
03197
03198
03199 big_jump = pow( 10., big_jump );
03200 if( big_jump > 1.2 )
03201 {
03202 if( big_jump > 3. )
03203 {
03204 sprintf( chLine,
03205 " C-The H density jumped at by %.0f%% at zone %ld, from %.2e to %.2e to %.2e",
03206 big_jump*100.,
03207 nzone_den_jump,
03208 struc.gas_phase[ipHYDROGEN][nzone_den_jump-1],
03209 struc.gas_phase[ipHYDROGEN][nzone_den_jump],
03210 struc.gas_phase[ipHYDROGEN][nzone_den_jump+1] );
03211 caunin(chLine);
03212 }
03213 else
03214 {
03215 sprintf( chLine,
03216 " !An H density jump occurred at zone %ld, by %.0f%% from %.2e to %.2e to %.2e",
03217 nzone_den_jump,
03218 big_jump*100.,
03219 struc.gas_phase[ipHYDROGEN][nzone_den_jump-1],
03220 struc.gas_phase[ipHYDROGEN][nzone_den_jump],
03221 struc.gas_phase[ipHYDROGEN][nzone_den_jump+1] );
03222 bangin(chLine);
03223 }
03224 }
03225
03226
03227 big_jump = -15;
03228 nzone_den_jump = 0;
03229
03230
03231
03232
03233 for( i=3; i < (nzone - 2 - (int)lgAbort); i++ )
03234 {
03235
03236 rel = fabs(log10( SDIV(struc.pres_radiation_lines_curr[i]) /
03237 SDIV(0.5*(struc.pres_radiation_lines_curr[i-1]+struc.pres_radiation_lines_curr[i+1])) ) );
03238
03239 if( rel > big_jump )
03240 {
03241 big_jump = rel;
03242 nzone_den_jump = i;
03243 }
03244 }
03245
03246
03247
03248 if( pressure.RadBetaMax > 0.01 )
03249 {
03250 big_jump = pow( 10., big_jump );
03251 if( big_jump > 1.2 )
03252 {
03253
03254
03255 if( big_jump > 3. && strcmp(dense.chDenseLaw,"CPRE") == 0)
03256 {
03257 sprintf( chLine,
03258 " C-The radiation pressure jumped by %.0f%% at zone %ld, from %.2e to %.2e to %.2e",
03259 big_jump*100.,
03260 nzone_den_jump,
03261 struc.pres_radiation_lines_curr[nzone_den_jump-1],
03262 struc.pres_radiation_lines_curr[nzone_den_jump],
03263 struc.pres_radiation_lines_curr[nzone_den_jump+1] );
03264 caunin(chLine);
03265 }
03266 else
03267 {
03268 sprintf( chLine,
03269 " !The radiation pressure jumped by %.0f%% at zone %ld, from %.2e to %.2e to %.2e",
03270 big_jump*100.,
03271 nzone_den_jump,
03272 struc.pres_radiation_lines_curr[nzone_den_jump-1],
03273 struc.pres_radiation_lines_curr[nzone_den_jump],
03274 struc.pres_radiation_lines_curr[nzone_den_jump+1] );
03275 bangin(chLine);
03276 }
03277 }
03278 }
03279
03280
03281 phycon.BigJumpTe = 0.;
03282 phycon.BigJumpne = 0.;
03283 phycon.BigJumpH2 = 0.;
03284 phycon.BigJumpCO = 0.;
03285
03286 for( i=1; i < (nzone - 1 - (int)lgAbort); i++ )
03287 {
03288
03289 rel = fabs(log10( struc.testr[i] / struc.testr[i+1] ) );
03290 if( rel > phycon.BigJumpTe )
03291 {
03292 phycon.BigJumpTe = (realnum)rel;
03293 }
03294
03295
03296 rel = fabs(log10( struc.ednstr[i] / struc.ednstr[i+1] ) );
03297 if( rel > phycon.BigJumpne )
03298 {
03299 phycon.BigJumpne = (realnum)rel;
03300 }
03301
03302
03303 if( (struc.H2_abund[i])>SMALLFLOAT && (struc.H2_abund[i+1]) > SMALLFLOAT
03304
03305 && (struc.H2_abund[i])/struc.gas_phase[ipHYDROGEN][i]>1e-3)
03306 {
03307 rel = fabs(log10( (struc.H2_abund[i]) / SDIV(struc.H2_abund[i+1]) ) );
03308 if( rel > phycon.BigJumpH2 )
03309 {
03310 phycon.BigJumpH2 = (realnum)rel;
03311 }
03312 }
03313
03314 {
03315 int ipCO = findspecies("CO")->index;
03316
03317
03318 if( ipCO != -1 &&
03319 struc.molecules[ipCO][i]>SMALLFLOAT &&
03320 struc.molecules[ipCO][i+1]>SMALLFLOAT &&
03321 struc.molecules[ipCO][i]/SDIV(struc.gas_phase[ipCARBON][i])>1e-3 )
03322 {
03323 rel = fabs(log10( struc.molecules[ipCO][i] / struc.molecules[ipCO][i+1] ) );
03324 if( rel > phycon.BigJumpCO )
03325 {
03326 phycon.BigJumpCO = (realnum)rel;
03327 }
03328 }
03329 }
03330 }
03331
03332
03333 if(phycon.BigJumpTe>0. )
03334 phycon.BigJumpTe = (realnum)pow( 10., (double)phycon.BigJumpTe ) -1.f;
03335
03336 if(phycon.BigJumpne>0. )
03337 phycon.BigJumpne = (realnum)pow( 10., (double)phycon.BigJumpne )-1.f;
03338
03339 if(phycon.BigJumpH2>0. )
03340 phycon.BigJumpH2 = (realnum)pow( 10., (double)phycon.BigJumpH2 )-1.f;
03341
03342 if(phycon.BigJumpCO>0. )
03343 phycon.BigJumpCO = (realnum)pow( 10., (double)phycon.BigJumpCO )-1.f;
03344
03345
03346
03347 if( phycon.BigJumpTe > 0.3 )
03348 {
03349 sprintf( chLine,
03350 " C-The temperature varied by %.1f%% between two zones",
03351 phycon.BigJumpTe*100.);
03352 caunin(chLine);
03353 }
03354 else if( phycon.BigJumpTe > 0.1 )
03355 {
03356 sprintf( chLine,
03357 " !The temperature varied by %.1f%% between two zones",
03358 phycon.BigJumpTe*100.);
03359 bangin(chLine);
03360 }
03361
03362 if( phycon.BigJumpne > 0.3 )
03363 {
03364 sprintf( chLine,
03365 " C-The electron density varied by %.1f%% between two zones",
03366 phycon.BigJumpne*100.);
03367 caunin(chLine);
03368 }
03369 else if( phycon.BigJumpne > 0.1 )
03370 {
03371 sprintf( chLine,
03372 " !The electron density varied by %.1f%% between two zones",
03373 phycon.BigJumpne*100.);
03374 bangin(chLine);
03375 }
03376
03377 if( phycon.BigJumpH2 > 0.8 )
03378 {
03379 sprintf( chLine,
03380 " C-The H2 density varied by %.1f%% between two zones",
03381 phycon.BigJumpH2*100.);
03382 caunin(chLine);
03383 }
03384 else if( phycon.BigJumpH2 > 0.1 )
03385 {
03386 sprintf( chLine,
03387 " !The H2 density varied by %.1f%% between two zones",
03388 phycon.BigJumpH2*100.);
03389 bangin(chLine);
03390 }
03391
03392 if( phycon.BigJumpCO > 0.8 )
03393 {
03394 sprintf( chLine,
03395 " C-The CO density varied by %.1f%% between two zones",
03396 phycon.BigJumpCO*100.);
03397 caunin(chLine);
03398 }
03399 else if( phycon.BigJumpCO > 0.2 )
03400 {
03401 sprintf( chLine,
03402 " !The CO density varied by %.1f%% between two zones",
03403 phycon.BigJumpCO*100.);
03404 bangin(chLine);
03405 }
03406 return;
03407 }