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