cloudy trunk
Loading...
Searching...
No Matches
prt_comment.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*PrtComment analyze model, generating comments on its features */
4/*chkCaHeps check whether CaII K and H epsilon overlap */
5/*prt_smooth_predictions check whether fluctuations in any predicted quantities occurred */
6#include "cddefines.h"
7#include "physconst.h"
8#include "cddrive.h"
9#include "lines_service.h"
10#include "iso.h"
11#include "continuum.h"
12#include "stopcalc.h"
13#include "hyperfine.h"
14#include "dense.h"
15#include "grainvar.h"
16#include "version.h"
17#include "rt.h"
18#include "he.h"
19#include "ionbal.h"
20#include "taulines.h"
21#include "hydrogenic.h"
22#include "lines.h"
23#include "trace.h"
24#include "hcmap.h"
25#include "hmi.h"
26#include "save.h"
27#include "h2.h"
28#include "conv.h"
29#include "dynamics.h"
30#include "opacity.h"
31#include "geometry.h"
32#include "elementnames.h"
33#include "ca.h"
34#include "broke.h"
35#include "pressure.h"
36#include "mole.h"
37#include "co.h"
38#include "atoms.h"
39#include "abund.h"
40#include "rfield.h"
41#include "colden.h"
42#include "phycon.h"
43#include "timesc.h"
44#include "hextra.h"
45#include "radius.h"
46#include "iterations.h"
47#include "fudgec.h"
48#include "called.h"
49#include "magnetic.h"
50#include "wind.h"
51#include "secondaries.h"
52#include "struc.h"
53#include "oxy.h"
54#include "input.h"
55#include "thermal.h"
56#include "atmdat.h"
57#include "warnings.h"
58
59/*chkCaHeps check whether CaII K and H epsilon overlap */
60STATIC void chkCaHeps(double *totwid);
61
62/*prt_smooth_predictions check whether fluctuations in any predicted quantities occurred */
64
65void PrtComment(void)
66{
67 char chLbl[11],
68 chLine[INPUT_LINE_LENGTH];
69
70 bool lgAbort_flag,
71 lgThick,
72 lgLots_of_moles;
73
74 long int dum1,
75 dum2,
76 i,
77 imas,
78 ipLo ,
79 ipHi ,
80 ipISO,
81 nelem,
82 isav,
83 j,
84 nc,
85 nline,
86 nn,
87 nneg,
88 ns,
89 nw;
90
91 double big_ion_jump,
92 absint ,
93 aj,
94 alpha,
95 big,
96 bigm,
97 comfrc,
98 differ,
99 error,
100 flur,
101 freqn,
102 rate,
103 ratio,
104 rel,
105 small,
106 tauneg,
107 ts ,
108 HBeta,
109 relfl ,
110 relhm,
111 fedest,
112 GetHeat,
113 SumNeg,
114 c4363,
115 t4363,
116 totwid;
117
118 double VolComputed , VolExpected , ConComputed , ConExpected;
119
120 bool lgLotsSolids;
121
122 DEBUG_ENTRY( "PrtComment()" );
123
124 if( 0 && lgAbort )
125 {
126 return;
127 }
128
129 if( trace.lgTrace )
130 {
131 fprintf( ioQQQ, " PrtComment called.\n" );
132 }
133
134 /*
135 * enter all comments cautions warnings and surprises into a large
136 * stack of statements
137 * at end of this routine is master call to printing routines
138 */
139 iterations.lgIterAgain = false;
140
141 /* initialize the line saver */
142 wcnint();
143
144 if( t_version::Inst().nBetaVer > 0 )
145 {
146 sprintf( chLine,
147 " !This is beta test version %ld and is intended for testing only.",
148 t_version::Inst().nBetaVer );
149 bangin(chLine);
150 }
151
152 /* this flag set by call to fixit routine,
153 * to show that parts of the code need repair.
154 * lgRelease is true if this is release version */
155 if( broke.lgFixit && !t_version::Inst().lgRelease )
156 {
157 sprintf( chLine, " !The code needs to be fixed - search for fixit()." );
158 bangin(chLine);
159 }
160
161 /* this flag set by call to CodeReview routine,
162 * to show that parts of the code need to be reviewed.
163 * lgRelease is true if this is release version */
164 if( broke.lgCheckit && !t_version::Inst().lgRelease )
165 {
166 sprintf( chLine, " !New code needs to be reviewed - search for CodeReview()." );
167 bangin(chLine);
168 }
169
170 /* say why calculation stopped */
171 if( conv.lgBadStop )
172 {
173 /* this case stop probably was not intended */
174 sprintf( warnings.chRgcln[0], " W-Calculation stopped because %s Iteration%3ld of %ld",
175 StopCalc.chReasonStop, iteration, iterations.itermx + 1 );
176 sprintf( chLine, " W-Calculation stopped because %s",
177 StopCalc.chReasonStop );
178 warnin(chLine);
179 sprintf( chLine, " W-This was not intended." );
180 warnin(chLine);
181 }
182 else
183 {
184 /* for iterate to convergence, print reason why it was not converged on 3rd and higher iterations */
185 if( (conv.lgAutoIt && iteration != iterations.itermx + 1) &&
186 iteration > 2 )
187 {
188 sprintf( warnings.chRgcln[0],
189 " Calculation stopped because %s Iteration %ld of %ld, not converged due to %s",
190 StopCalc.chReasonStop,
191 iteration,
192 iterations.itermx + 1,
193 conv.chNotConverged );
194 }
195 else
196 {
197 sprintf( warnings.chRgcln[0],
198 " Calculation stopped because %s Iteration %ld of %ld",
199 StopCalc.chReasonStop, iteration, iterations.itermx + 1 );
200 }
201 }
202
203 /* check whether stopped because default number of zones hit,
204 * not intended?? */
205 if( (!geometry.lgZoneSet) && geometry.lgZoneTrp )
206 {
207 conv.lgBadStop = true;
208 sprintf( chLine,
209 " W-Calculation stopped because default number of zones reached. Was this intended???" );
210 warnin(chLine);
211 sprintf( chLine,
212 " W-Default limit can be increased while retaining this check with the SET NEND command." );
213 warnin(chLine);
214 }
215
216 /* check whether stopped because zones too thin - only happens when glob set
217 * and depth + dr = depth
218 * not intended */
219 if( radius.lgDrMinUsed || radius.lgdR2Small )
220 {
221 conv.lgBadStop = true;
222 sprintf( chLine,
223 " W-Calculation stopped zone thickness became too thin. This was not intended." );
224 warnin(chLine);
225 sprintf( chLine,
226 " W-The most likely reason was an uncontrolled oscillation." );
227 warnin(chLine);
228 ShowMe();
229 }
230
231 if( radius.lgdR2Small )
232 {
233 sprintf( chLine,
234 " W-This happened because the globule scale became very small relative to the depth." );
235 warnin(chLine);
236 sprintf( chLine,
237 " W-This problem is described in Hazy." );
238 warnin(chLine);
239 }
240
241 /* possible that last zone does not have stored temp - if so
242 * then make it up - this happens for some bad stops */
243 ASSERT( nzone < struc.nzlim );
244
245 if( struc.testr[nzone-1] == 0. && nzone > 1)
246 struc.testr[nzone-1] = struc.testr[nzone-2];
247
248 if( struc.ednstr[nzone-1] == 0. && nzone > 1)
249 struc.ednstr[nzone-1] = struc.ednstr[nzone-2];
250
251 /* give indication of geometry */
252 rel = radius.depth/radius.rinner;
253 if( rel < 0.1 )
254 {
255 sprintf( warnings.chRgcln[1], " The geometry is plane-parallel." );
256 }
257 else if( rel >= 0.1 && rel < 3. )
258 {
259 sprintf( warnings.chRgcln[1], " The geometry is a thick shell." );
260 }
261 else
262 {
263 sprintf( warnings.chRgcln[1], " The geometry is spherical." );
264 }
265
266 /* levels of warnings: Warning (possibly major problems)
267 * Caution (not likely to invalidate the results)
268 * [ !] surprise, but not a problem
269 * [nothing] interesting note
270 */
271 /* this flag set by call to routine broken ( );
272 * and show that the code is broken. */
273 if( broke.lgBroke )
274 {
275 sprintf( chLine, " W-The code is broken - search for broken()." );
276 warnin(chLine);
277 }
278
279 /* incorrect electron density detected */
280 if( dense.lgEdenBad )
281 {
282 if( dense.nzEdenBad == nzone )
283 {
284 sprintf( chLine, " C-The assumed electron density was incorrect for the last zone." );
285 caunin(chLine);
286 sprintf( chLine, " C-Did a temperature discontinuity occur??" );
287 caunin(chLine);
288 }
289 else
290 {
291 sprintf( chLine, " W-The assumed electron density was incorrect during the calculation. This is bad." );
292 warnin(chLine);
293 ShowMe();
294 }
295 }
296
297 if( lgAbort )
298 {
299 sprintf( chLine, " W-The calculation aborted. Something REALLY went wrong!" );
300 warnin(chLine);
301 }
302
303 /* thermal map was done but results were not ok */
304 if( hcmap.lgMapDone && !hcmap.lgMapOK )
305 {
306 sprintf( chLine, " !The thermal map had changes in slope - check map output." );
307 bangin(chLine);
308 }
309
310 /* first is greater than zero if fudge factors were entered, second is
311 * true if fudge ever evaluated, even to see if fudge factors are in place */
312 if( fudgec.nfudge > 0 || fudgec.lgFudgeUsed )
313 {
314 sprintf( chLine, " !Fudge factors were used or were checked. Why?" );
315 bangin(chLine);
316 }
317
318 if( dense.gas_phase[ipHYDROGEN] > 1.1e13 )
319 {
320 if( dense.gas_phase[ipHYDROGEN] > 1e15 )
321 {
322 sprintf( chLine, " C-Density greater than 10**15, heavy elements are very uncertain." );
323 caunin(chLine);
324 }
325 else
326 {
327 sprintf( chLine, " C-Density greater than 10**13" );
328 caunin(chLine);
329 }
330 }
331
332 /* HBeta is used later in the code to check on line intensities */
333 if( cdLine("Pump",4861.36f,&relfl,&absint)<=0 )
334 {
335 fprintf( ioQQQ, " PROBLEM Did not find Pump H-beta, set to unity\n" );
336 relfl = 1.;
337 absint = 1.;
338 }
339
340 /* now find total Hbeta */
341 /* >>chng from "totl" Hbeta which was a special entry, to "H 1" Hbeta, which
342 * is the general case */
343 if( cdLine( "H 1",4861.36f,&HBeta,&absint)<=0 )
344 {
345 fprintf( ioQQQ, " NOTE Did not find H 1 H-beta - set intensity to unity, "
346 "will not check on importance of H 1 pumping.\n" );
347 HBeta = 1.;
348 absint = 1.;
349 }
350 else
351 {
352 /* check on continuum pumping of Balmer lines */
353 if( HBeta>SMALLFLOAT )
354 {
355 flur = relfl/HBeta;
356 if( flur > 0.1 )
357 {
358 sprintf( chLine, " !Continuum fluorescent production of H-beta was very important." );
359 bangin(chLine);
360 }
361 else if(flur > 0.01 )
362 {
363 sprintf( chLine, " Continuum fluorescent production of H-beta was significant." );
364 notein(chLine);
365 }
366 }
367 }
368
369 // iterate to convergence - status of this iteration
370 if( iteration > 1 && conv.lgAutoIt && iteration<iterations.itermx)
371 {
372 sprintf( chLine , " Iteration not converged because %s.",
373 conv.chNotConverged );
374 notein(chLine);
375 }
376
377 /* check if there were problems with initial wind velocity */
378 if( wind.lgBallistic() && ((!wind.lgWindOK) || wind.windv < 1e6) )
379 {
380 sprintf( chLine, " C-Wind velocity below sonic point; solution is not valid." );
381 caunin(chLine);
382 }
383
384 /* now confirm that mass flux here is correct */
385 if( !wind.lgStatic() )
386 {
387 rel = wind.emdot/(wind.windv*dense.gas_phase[ipHYDROGEN])/radius.r1r0sq;
388 if( fabs(1.-rel)> 0.02 )
389 {
390 sprintf( chLine, " C-Wind mass flux error is %g%%",fabs(1.-rel)*100. );
391 caunin(chLine);
392 fprintf(ioQQQ,"DEBUG emdot\t%.3e\t%.3e\t%.3e\t%.3e\n",
393 wind.emdot , wind.windv*dense.gas_phase[ipHYDROGEN],wind.windv,dense.gas_phase[ipHYDROGEN]);
394 }
395 }
396
397 /* check that we didn't overrun zone scale */
398 if( nzone >= struc.nzlim )
399 {
401 }
402
403 // check on energy conservation
404 if( !lgConserveEnergy() )
405 {
406 ShowMe();
407 lgAbort = true;
408 fprintf(ioQQQ,"\n\n Enable per zone energy conservation check by setting "
409 "CHECK_ENERGY_EVERY_ZONE=true in cloudy.cpp, recompile, then rerun.\n");
410 }
411
412 /* comments having to do with cosmic rays */
413 /* comment if cosmic rays and magnetic field both present */
414 if( hextra.cryden*magnetic.lgB > 0. )
415 {
416 sprintf( chLine,
417 " !Magnetic field & cosmic rays both present. Their interactions are not treated." );
418 bangin(chLine);
419 }
420
421 /* comment if cosmic rays are not included and stop temp has been lowered to go into neutral gas */
422 if( hextra.cryden== 0. && StopCalc.TempLoStopZone < phycon.TEMP_STOP_DEFAULT)
423 {
424 sprintf( chLine,
425 " !Background cosmic rays are not included - is this physical? It affects the chemistry." );
426 bangin(chLine);
427 }
428
429 /* check whether cosmic rays on, but model thick to them */
430 if( hextra.cryden > 0. && (colden.colden[ipCOL_H0]/10. + colden.colden[ipCOL_Hp]) > 1e23 )
431 {
432 sprintf( chLine,
433 " C-Model is thick to cosmic rays, which are on." );
434 caunin(chLine);
435 }
436
437 /* was ionization rate less than cosmic ray ionization rate in ISM? */
438 if( hextra.cryden == 0. && iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc < 1e-17 )
439 {
440 sprintf( chLine,
441 " !Ionization rate fell below background cosmic ray ionization rate. Should this be added too?" );
442 bangin(chLine);
443 sprintf( chLine,
444 " ! Use the COSMIC RAY BACKGROUND command." );
445 bangin(chLine);
446 }
447
448 /* PrtComment if test code is in place */
449 if( lgTestCodeCalled && !t_version::Inst().lgReleaseBranch && !t_version::Inst().lgRelease )
450 {
451 sprintf( chLine, " !Test code is in place." );
452 bangin(chLine);
453 }
454
455 /* lgComUndr set to .true. if Compton cooling rate underflows to 0 */
456 if( rfield.lgComUndr )
457 {
458 sprintf( chLine,
459 " !Compton cooling rate underflows to zero. Is this important?" );
460 bangin(chLine);
461 }
462
463 /* make note if input stream contained an underscore, which was converted into a space */
464 if( input.lgUnderscoreFound )
465 {
466 sprintf( chLine,
467 " !Some input lines contained underscores, these were changed to spaces." );
468 bangin(chLine);
469 }
470
471 /* make note if input stream contained a left or right bracket, which was converted into a space */
472 if( input.lgBracketFound )
473 {
474 sprintf( chLine,
475 " !Some input lines contained [ or ], these were changed to spaces." );
476 bangin(chLine);
477 }
478
479 /* lgHionRad set to .true. if no hydrogen ionizing radiation */
480 if( rfield.lgHionRad )
481 {
482 sprintf( chLine,
483 " !There is no hydrogen-ionizing radiation. Was this intended?" );
484 bangin(chLine);
485 }
486
487 /* check whether certain zones were thermally unstable */
488 if( thermal.nUnstable > 0 )
489 {
490 sprintf( chLine,
491 " Derivative of net cooling negative and so possibly thermally unstable in%4ld zones.",
492 thermal.nUnstable );
493 notein(chLine);
494 }
495
496 /* generate a bang if a large fraction of the zones were unstable */
497 if( nzone > 1 &&
498 (realnum)(thermal.nUnstable)/(realnum)(nzone) > 0.25 )
499 {
500 sprintf( chLine,
501 " !A large fraction of the zones were possibly thermally unstable,%4ld out of%4ld",
502 thermal.nUnstable, nzone );
503 bangin(chLine);
504 }
505
506 /* comment if negative coolants were ever significant */
507 if( thermal.CoolHeatMax > 0.2 )
508 {
509 sprintf( chLine,
510 " !Negative cooling reached %6.1f%% of the local heating, due to %4.4s %.1f",
511 thermal.CoolHeatMax*100., thermal.chCoolHeatMax, thermal.wlCoolHeatMax );
512 bangin(chLine);
513 }
514 else if( thermal.CoolHeatMax > 0.05 )
515 {
516 sprintf( chLine,
517 " Negative cooling reached %6.1f%% of the local heating, due to %4.4s %.2f",
518 thermal.CoolHeatMax*100., thermal.chCoolHeatMax, thermal.wlCoolHeatMax );
519 notein(chLine);
520 }
521
522 /* check if advection heating was important */
523 if( dynamics.HeatMax > 0.05 )
524 {
525 sprintf( chLine,
526 " !Advection heating reached %.2f%% of the local heating.",
527 dynamics.HeatMax*100. );
528 bangin(chLine);
529 }
530 else if( dynamics.HeatMax > 0.005 )
531 {
532 sprintf( chLine,
533 " Advection heating reached %.2f%% of the local heating.",
534 dynamics.HeatMax*100. );
535 notein(chLine);
536 }
537
538 /* check if advection cooling was important */
539 if( dynamics.CoolMax > 0.05 )
540 {
541 sprintf( chLine,
542 " !Advection cooling reached %.2f%% of the local cooling.",
543 dynamics.CoolMax*100. );
544 bangin(chLine);
545 }
546 else if( dynamics.CoolMax > 0.005 )
547 {
548 sprintf( chLine,
549 " Advection cooling reached %.2f%% of the local heating.",
550 dynamics.CoolMax*100. );
551 notein(chLine);
552 }
553
554 /* >>chng 06 mar 22, add this comment
555 * check if time dependent ionization front being done with too large a U */
556 if( dynamics.lgTimeDependentStatic && dynamics.lgRecom )
557 {
558 if( rfield.uh > 1. )
559 {
560 sprintf( chLine,
561 " W-Time dependent ionization front cannot now handle strong-R cases - the ionization parameter is too large." );
562 warnin(chLine);
563 }
564 else if( rfield.uh > 0.1 )
565 {
566 sprintf( chLine,
567 " C-Time dependent ionization front cannot now handle strong-R cases - the ionization parameter is too large." );
568 caunin(chLine);
569 }
570 }
571
572 /* check if thermal ionization of ground state of hydrogen was important */
573 if( hydro.HCollIonMax > 0.10 )
574 {
575 sprintf( chLine,
576 " !Thermal collisional ionization of H reached %.2f%% of the local ionization rate.",
577 hydro.HCollIonMax*100. );
578 bangin(chLine);
579 }
580 else if( hydro.HCollIonMax > 0.005 )
581 {
582 sprintf( chLine,
583 " Thermal collisional ionization of H reached %.2f%% of the local ionization rate.",
584 hydro.HCollIonMax*100. );
585 notein(chLine);
586 }
587
588 /* check if lookup table for Hummer & Storey case B was exceeded */
589 if( !atmdat.lgHCaseBOK[1][ipHYDROGEN] )
590 {
591 sprintf( chLine,
592 " Te-ne bounds of Case B lookup table exceeded, H I Case B line intensities set to zero." );
593 notein(chLine);
594 }
595 if( !atmdat.lgHCaseBOK[1][ipHELIUM] )
596 {
597 sprintf( chLine,
598 " Te-ne bounds of Case B lookup table exceeded, He II Case B line intensities set to zero." );
599 notein(chLine);
600 }
601
602 if( dense.EdenMax>1e8 )
603 {
604 sprintf( chLine,
605 " !The high electron density makes the Nussbaumer/Storey CNO recombination predictions unreliable." );
606 bangin(chLine);
607 }
608
609 /* check if secondary ionization of hydrogen was important */
610 if( secondaries.SecHIonMax > 0.10 )
611 {
612 sprintf( chLine,
613 " !Suprathermal collisional ionization of H reached %.2f%% of the local H ionization rate.",
614 secondaries.SecHIonMax*100. );
615 bangin(chLine);
616 }
617 else if( secondaries.SecHIonMax > 0.005 )
618 {
619 sprintf( chLine,
620 " Suprathermal collisional ionization of H reached %.2f%% of the local H ionization rate.",
621 secondaries.SecHIonMax*100. );
622 notein(chLine);
623 }
624
625 /* check if H2 vib-deexcitation heating was important */
626 if( hmi.HeatH2DexcMax > 0.05 )
627 {
628 sprintf( chLine,
629 " !H2 vib deexec heating reached %.2f%% of the local heating.",
630 hmi.HeatH2DexcMax*100. );
631 bangin(chLine);
632 }
633 else if( hmi.HeatH2DexcMax > 0.005 )
634 {
635 sprintf( chLine,
636 " H2 vib deexec heating reached %.2f%% of the local heating.",
637 hmi.HeatH2DexcMax*100. );
638 notein(chLine);
639 }
640
641 /* check if H2 vib-deexcitation heating was important */
642 if( hmi.CoolH2DexcMax > 0.05 )
643 {
644 sprintf( chLine,
645 " !H2 deexec cooling reached %.2f%% of the local heating.",
646 hmi.CoolH2DexcMax*100. );
647 bangin(chLine);
648 }
649 else if( hmi.CoolH2DexcMax > 0.005 )
650 {
651 sprintf( chLine,
652 " H2 deexec cooling reached %.2f%% of the local heating.",
653 hmi.CoolH2DexcMax*100. );
654 notein(chLine);
655 }
656
657 /* check if charge transfer ionization of hydrogen was important */
658 if( atmdat.HIonFracMax > 0.10 )
659 {
660 sprintf( chLine,
661 " !Charge transfer ionization of H reached %.1f%% of the local H ionization rate.",
662 atmdat.HIonFracMax*100. );
663 bangin(chLine);
664 }
665 else if( atmdat.HIonFracMax > 0.005 )
666 {
667 sprintf( chLine,
668 " Charge transfer ionization of H reached %.2f%% of the local H ionization rate.",
669 atmdat.HIonFracMax*100. );
670 notein(chLine);
671 }
672
673 /* check if charge transfer heating cooling was important */
674 if( atmdat.HCharHeatMax > 0.05 )
675 {
676 sprintf( chLine,
677 " !Charge transfer heating reached %.2f%% of the local heating.",
678 atmdat.HCharHeatMax*100. );
679 bangin(chLine);
680 }
681 else if( atmdat.HCharHeatMax > 0.005 )
682 {
683 sprintf( chLine,
684 " Charge transfer heating reached %.2f%% of the local heating.",
685 atmdat.HCharHeatMax*100. );
686 notein(chLine);
687 }
688
689 if( atmdat.HCharCoolMax > 0.05 )
690 {
691 sprintf( chLine,
692 " !Charge transfer cooling reached %.2f%% of the local heating.",
693 atmdat.HCharCoolMax*100. );
694 bangin(chLine);
695 }
696 else if( atmdat.HCharCoolMax > 0.005 )
697 {
698 sprintf( chLine,
699 " Charge transfer cooling reached %.2f%% of the local heating.",
700 atmdat.HCharCoolMax*100. );
701 notein(chLine);
702 }
703
704 /* check whether photo from up level of Mg2 2798 ever important */
705 if( atoms.xMg2Max > 0.1 )
706 {
707 sprintf( chLine,
708 " !Photoionization of upper level of Mg II 2798 reached %.1f%% of the total Mg+ photo rate.",
709 atoms.xMg2Max*100. );
710 bangin(chLine);
711 }
712 else if( atoms.xMg2Max > 0.01 )
713 {
714 sprintf( chLine,
715 " Photoionization of upper level of Mg II 2798 reached %.1f%% of the total Mg+ photo rate.",
716 atoms.xMg2Max*100. );
717 notein(chLine);
718 }
719
720 /* check whether photo from up level of [O I] 6300 ever important */
721 if( oxy.poimax > 0.1 )
722 {
723 sprintf( chLine,
724 " !Photoionization of upper levels of [O I] reached %.1f%% of the total O destruction rate.",
725 oxy.poimax*100. );
726 bangin(chLine);
727 }
728 else if( oxy.poimax > 0.01 )
729 {
730 sprintf( chLine,
731 " Photoionization of upper levels of [O I] reached %.1f%% of the total O destruction rate.",
732 oxy.poimax*100. );
733 notein(chLine);
734 }
735
736 /* check whether photo from up level of [O III] 5007 ever important */
737 if( (oxy.poiii2Max + oxy.poiii3Max) > 0.1 )
738 {
739 sprintf( chLine,
740 " !Photoionization of upper levels of [O III] reached %.1f%% of the total O++ photo rate.",
741 (oxy.poiii2Max + oxy.poiii3Max)*100. );
742 bangin(chLine);
743 }
744 else if( (oxy.poiii2Max + oxy.poiii3Max) > 0.01 )
745 {
746 sprintf( chLine,
747 " Photoionization of upper levels of [O III] reached %.1f%% of the total O++ photo rate.",
748 (oxy.poiii2Max + oxy.poiii3Max)*100. );
749 notein(chLine);
750 }
751
752 /* check whether photoionization of He 2trip S was important */
753 if( he.frac_he0dest_23S > 0.1 )
754 {
755 sprintf( chLine,
756 " !Destruction of He 2TriS reached %.1f%% of the total He0 dest rate"
757 " at zone %li, %.1f%% of that was photoionization.",
758 he.frac_he0dest_23S*100,
759 he.nzone,
760 he.frac_he0dest_23S_photo*100. );
761 bangin(chLine);
762 }
763 else if( he.frac_he0dest_23S > 0.01 )
764 {
765 sprintf( chLine,
766 " Destruction of He 2TriS reached %.1f%% of the total He0 dest rate"
767 " at zone %li, %.1f%% of that was photoionization.",
768 he.frac_he0dest_23S*100,
769 he.nzone,
770 he.frac_he0dest_23S_photo*100. );
771 notein(chLine);
772 }
773
774 /* check for critical density for l-mixing */
775 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
776 {
777 if( !iso_ctrl.lgCritDensLMix[ipISO] && dense.lgElmtOn[ipISO] )
778 {
779 sprintf( chLine,
780 " The density is too low to l-mix the lowest %s I collapsed level. "
781 " More resolved levels are needed for accurate line ratios.",
782 elementnames.chElementSym[ipISO]);
783 notein(chLine);
784 }
785 }
786
787 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
788 {
789 /* report continuum lowering for xx-like iso-sequence. */
790 for( nelem=ipISO; nelem<LIMELM; ++nelem )
791 {
792 if( iso_sp[ipISO][nelem].lgLevelsLowered && dense.lgElmtOn[nelem] )
793 {
794 sprintf( chLine, " !Continuum was lowered into model %s-like %s due to high density. Highest n is %li",
795 elementnames.chElementSym[ipISO],
796 elementnames.chElementSym[nelem],
797 iso_sp[ipISO][nelem].n_HighestResolved_local+iso_sp[ipISO][nelem].nCollapsed_local);
798 bangin(chLine);
799 }
800 else if( iso_sp[ipISO][nelem].lgLevelsEverLowered && dense.lgElmtOn[nelem] )
801 {
802 sprintf( chLine, " !Continuum was lowered into model %s-like %s due to high density at SOME point but NOT at the last zone.",
803 elementnames.chElementSym[ipISO],
804 elementnames.chElementNameShort[nelem]);
805 bangin(chLine);
806 }
807 }
808
809 /* report pop rescaling for xx-like iso-sequence. */
810 for( nelem=ipISO; nelem<LIMELM; ++nelem )
811 {
812 if( iso_sp[ipISO][nelem].lgPopsRescaled )
813 {
814 ASSERT( dense.lgSetIoniz[nelem] );
815 sprintf( chLine, " C-Populations were rescaled for %s-like %s due to \"element ionization\" command.",
816 elementnames.chElementSym[ipISO],
817 elementnames.chElementSym[nelem] );
818 caunin(chLine);
819 }
820 }
821 }
822
823 /* frequency array may not have been defined for all energies */
824 if( !rfield.lgMMok )
825 {
826 sprintf( chLine,
827 " C-Continuum not defined in extreme infrared - Compton scat, grain heating, not treated properly?" );
828 caunin(chLine);
829 }
830
831 if( !rfield.lgHPhtOK )
832 {
833 sprintf( chLine,
834 " C-Continuum not defined at photon energies which ionize excited states of H, important for H- and ff heating." );
835 caunin(chLine);
836 }
837
838 if( !rfield.lgXRayOK )
839 {
840 sprintf( chLine,
841 " C-Continuum not defined at X-Ray energies - Compton scattering and Auger ionization wrong?" );
842 caunin(chLine);
843 }
844
845 if( !rfield.lgGamrOK )
846 {
847 sprintf( chLine,
848 " C-Continuum not defined at gamma-ray energies - pair production and Compton scattering OK?" );
849 caunin(chLine);
850 }
851
852 if( continuum.lgCon0 )
853 {
854 sprintf( chLine, " C-Continuum zero at some energies." );
855 caunin(chLine);
856 }
857
858 if( continuum.lgCoStarInterpolationCaution )
859 {
860 sprintf( chLine , " C-CoStarInterpolate interpolated between non-adjoining tracks, this may not be valid." );
861 caunin(chLine);
862 }
863
864 if( rfield.lgOcc1Hi )
865 {
866 sprintf( chLine,
867 " !The continuum occupation number at 1 Ryd is greater than unity." );
868 bangin(chLine);
869 }
870
871 /* this flag set true it set dr forced first zone to be too big */
872 if( radius.lgDR2Big )
873 {
874 sprintf( chLine,
875 " C-The thickness of the first zone was set larger than optimal by a SET DR command." );
876 caunin(chLine);
877 /* this is case where did one zone of specified thickness - but it
878 * was too large */
879 if( nzone<=1 )
880 sprintf( chLine,
881 " C-Consider using the STOP THICKNESS command instead." );
882 caunin(chLine);
883 }
884
885 /* check whether non-col excitation of 4363 was important */
886 if( cdLine("TOTL",4363,&t4363,&absint)<=0 )
887 {
888 fprintf( ioQQQ, " PrtComment could not find total O III 4363 with cdLine.\n" );
889 ShowMe();
891 }
892
893 if( cdLine("Coll",4363,&c4363,&absint)<=0 )
894 {
895 fprintf( ioQQQ, " PrtComment could not find collisional O III 4363 with cdLine.\n" );
896 ShowMe();
898 }
899
900 /* only print this comment if 4363 is significant and density low */
901 if( HBeta > 1e-35 )
902 {
903 /* >>chng 02 feb 27, lower ratio from -4 to -5, and raise density from 7 to 8 */
904 if( t4363/HBeta > 1e-5 && dense.gas_phase[ipHYDROGEN] < 1e8 )
905 {
906 ratio = (t4363 - c4363)/t4363;
907 if( ratio > 0.01 )
908 {
909 sprintf( chLine,
910 " !Non-collisional excitation of [O III] 4363 reached %.2f%% of the total.",
911 ratio*100. );
912 bangin(chLine);
913 }
914 else if( ratio > 0.001 )
915 {
916 sprintf( chLine,
917 " Non-collisional excitation of [O III] 4363 reached %.2f%% of the total.",
918 ratio*100. );
919 notein(chLine);
920 }
921 }
922 }
923
924 /* check for plasma shielding */
925 if( rfield.lgPlasNu )
926 {
927 sprintf( chLine,
928 " !The largest plasma frequency was %.2e Ryd = %.2e micron The continuum is set to 0 below this.",
929 rfield.plsfrqmax,
930 /* wavelength in microns */
931 RYDLAM/rfield.plsfrqmax/1e4);
932 bangin(chLine);
933 }
934
935 if( rfield.occmax > 0.1 )
936 {
937 if( rfield.occmnu > 1e-4 )
938 {
939 sprintf( chLine,
940 " !The largest continuum occupation number was %.3e at %.3e Ryd.",
941 rfield.occmax, rfield.occmnu );
942 bangin(chLine);
943 }
944 else
945 {
946 /* not surprising if occupation number bigger than 1 around 1e-5 Ryd,
947 * since this is the case for 3K background */
948 sprintf( chLine,
949 " The largest continuum occupation number was %.3e at %.3e Ryd.",
950 rfield.occmax, rfield.occmnu );
951 notein(chLine);
952 }
953 }
954
955 if( rfield.occmax > 1e4 && rfield.occ1nu > 0. )
956 {
957 /* occ1nu is energy (ryd) where continuum occupation number falls below 1 */
958 if( rfield.occ1nu < 0.0912 )
959 {
960 sprintf( chLine,
961 " The continuum occupation number fell below 1 at %.3e microns.",
962 0.0912/rfield.occ1nu );
963 notein(chLine);
964 }
965 else if( rfield.occ1nu < 1. )
966 {
967 sprintf( chLine,
968 " The continuum occupation number fell below 1 at %.3e Angstroms.",
969 912./rfield.occ1nu );
970 notein(chLine);
971 }
972 else
973 {
974 sprintf( chLine,
975 " The continuum occupation number fell below 1 at %.3e Ryd.",
976 rfield.occ1nu );
977 notein(chLine);
978 }
979 }
980
981 if( rfield.tbrmax > 1e3 )
982 {
983 sprintf( chLine,
984 " !The largest continuum brightness temperature was %.3eK at %.3e Ryd.",
985 rfield.tbrmax, rfield.tbrmnu );
986 bangin(chLine);
987 }
988
989 if( rfield.tbrmax > 1e4 )
990 {
991 /* tbr4nu is energy (ryd) where continuum bright temp falls < 1e4 */
992 if( rfield.tbr4nu < 0.0912 )
993 {
994 sprintf( chLine,
995 " The continuum brightness temperature fell below 10000K at %.3e microns.",
996 0.0912/rfield.tbr4nu );
997 notein(chLine);
998 }
999 else if( rfield.tbr4nu < 1. )
1000 {
1001 sprintf( chLine,
1002 " The continuum brightness temperature fell below 10000K at %.3e Angstroms.",
1003 912./rfield.tbr4nu );
1004 notein(chLine);
1005 }
1006 else
1007 {
1008 sprintf( chLine,
1009 " The continuum brightness temperature fell below 10000K at %.3e Ryd.",
1010 rfield.tbr4nu );
1011 notein(chLine);
1012 }
1013 }
1014
1015 /* turbulence AND constant pressure do not make sense */
1016 if( DoppVel.TurbVel > 0. && strcmp(dense.chDenseLaw,"CPRE") == 0 )
1017 {
1018 sprintf( chLine,
1019 " !Both constant pressure and turbulence makes no physical sense?" );
1020 bangin(chLine);
1021 }
1022
1023 /* filling factor AND constant pressure do not make sense */
1024 if( geometry.FillFac < 1. && strcmp(dense.chDenseLaw,"CPRE") == 0 )
1025 {
1026 sprintf( chLine,
1027 " !Both constant pressure and a filling factor makes no physical sense?" );
1028 bangin(chLine);
1029 }
1030
1031 /* grains and solar abundances do not make sense */
1032 if( gv.lgDustOn() && abund.lgAbnSolar )
1033 {
1034 sprintf( chLine,
1035 " !Grains are present, but the gas phase abundances were left at the solar default. This is not physical." );
1036 bangin(chLine);
1037 }
1038
1039 /* check if depletion command set but no grains, another silly thing to do */
1040 if( abund.lgDepln && !gv.lgDustOn() )
1041 {
1042 sprintf( chLine,
1043 " !Grains are not present, but the gas phase abundances were depleted. This is not physical." );
1044 bangin(chLine);
1045 }
1046
1047 if( gv.lgDustOn() )
1048 {
1049 long nBin=0L, nFail=0L;
1050 for( size_t nd=0; nd < gv.bin.size(); nd++ )
1051 {
1052 if( gv.bin[nd]->QHeatFailures > 0L )
1053 {
1054 ++nBin;
1055 nFail += gv.bin[nd]->QHeatFailures;
1056 }
1057 }
1058 if( nFail > 0 )
1059 {
1060 sprintf( chLine,
1061 " !The grain quantum heating treatment failed to converge %ld time(s) in %ld bin(s).", nFail, nBin );
1062 bangin(chLine);
1063 }
1064 }
1065
1066#if 0
1067 /* check if PAHs were present in the ionized region */
1068 /* >>chng 05 jan 01, disabled this code now that PAH's have varying abundances by default, PvH */
1070 if( gv.lgDustOn() )
1071 {
1072 bool lgPAHsPresent_and_constant = false;
1073 for( size_t nd=0; nd < gv.bin.size(); nd++ )
1074 {
1075 lgPAHsPresent_and_constant = lgPAHsPresent_and_constant ||
1076 /* it is ok to have PAHs in the ionized region if the abundances vary */
1077 (gv.bin[nd]->lgPAHsInIonizedRegion /* && !gv.bin[nd]-> lgDustVary */);
1078 }
1079 if( lgPAHsPresent_and_constant )
1080 {
1081 sprintf( chLine,
1082 " C-PAH's were present in the ionized region, this has never been observed in H II regions." );
1083 caunin(chLine);
1084 }
1085 }
1086#endif
1087
1088 /* constant temperature greater than continuum energy density temperature */
1089 if( thermal.lgTemperatureConstant && thermal.ConstTemp*1.0001 < phycon.TEnerDen )
1090 {
1091 sprintf( chLine,
1092 " C-The continuum energy density temperature (%g K)"
1093 " is greater than the gas kinetic temperature (%g K).",
1094 phycon.TEnerDen , thermal.ConstTemp);
1095 caunin(chLine);
1096 sprintf( chLine, " C-This is unphysical." );
1097 caunin(chLine);
1098 }
1099
1100 /* remark that grains not present but energy density was low */
1101 if( !gv.lgDustOn() && phycon.TEnerDen < 800. )
1102 {
1103 sprintf( chLine,
1104 " Grains were not present but might survive in this environment (energy density temperature was %.2eK)",
1105 phycon.TEnerDen );
1106 notein(chLine);
1107 }
1108
1109 /* call routine that will check age of cloud */
1110 AgeCheck();
1111
1112 /* check on Ca H and H-epsilon overlapping
1113 * need to do this since need to refer to lines arrays */
1114 chkCaHeps(&totwid);
1115 if( totwid > 121. )
1116 {
1117 sprintf( chLine, " H-eps and Ca H overlap." );
1118 notein(chLine);
1119 }
1120
1121 /* warning that something was turned off */
1122 if( !phycon.lgPhysOK )
1123 {
1124 sprintf( chLine, " !A physical process has been disabled." );
1125 bangin(chLine);
1126 }
1127
1128 /* check on lifetimes of [O III] against photoionization, only for low den */
1129 if( dense.gas_phase[ipHYDROGEN] < 1e8 )
1130 {
1131 if( oxy.r5007Max > 0.0263f )
1132 {
1133 sprintf( chLine,
1134 " !Photoionization of upper level of [O III] 5007 reached %.2e%% of the radiative lifetime.",
1135 oxy.r5007Max*100. );
1136 bangin(chLine);
1137 }
1138 else if( oxy.r5007Max > 0.0263f/10.f )
1139 {
1140 sprintf( chLine,
1141 " Photoionization of upper level of [O III] 5007 reached %.2e%% of the radiative lifetime.",
1142 oxy.r5007Max*100. );
1143 notein(chLine);
1144 }
1145 if( oxy.r4363Max > 1.78f )
1146 {
1147 sprintf( chLine,
1148 " !Photoionization of upper level of [O III] 4363 reached %.2e%% of the radiative lifetime.",
1149 oxy.r4363Max*100. );
1150 bangin(chLine);
1151 }
1152 else if( oxy.r4363Max > 1.78f/10.f )
1153 {
1154 sprintf( chLine,
1155 " Photoionization of upper level of [O III] 4363 reached %.2e%% of the radiative lifetime.",
1156 oxy.r4363Max*100. );
1157 notein(chLine);
1158 }
1159 }
1160
1161 /* check whether total heating and cooling matched
1162 * >>chng 97 mar 28, added GrossHeat, heat in terms normally heat-cool */
1163 error = fabs(thermal.power-thermal.totcol)/SDIV((thermal.power + thermal.totcol)/2.);
1164 if( thermal.lgTemperatureConstant )
1165 {
1166 if( error > 0.05 )
1167 {
1168 sprintf( chLine,
1169 " !Heating - cooling mismatch =%5.1f%%. Caused by constant temperature assumption. ",
1170 error*100. );
1171 bangin(chLine);
1172 }
1173 }
1174
1175 else
1176 {
1177 if( error > 0.05 && error < 0.2 )
1178 {
1179 sprintf( chLine, " C-Heating - cooling mismatch =%.1f%%. What\'s wrong?",
1180 error*100. );
1181 caunin(chLine);
1182 }
1183 else if( error >= 0.2 )
1184 {
1185 sprintf( chLine, " W-Heating - cooling mismatch =%.2e%%. What\'s wrong????",
1186 error*100. );
1187 warnin(chLine);
1188 }
1189 }
1190
1191 /* say if Ly-alpha photo of Ca+ excited levels was important */
1192 if( ca.Ca2RmLya > 0.01 )
1193 {
1194 sprintf( chLine,
1195 " Photoionization of Ca+ 2D level by Ly-alpha reached %6.1f%% of the total rate out.",
1196 ca.Ca2RmLya*100. );
1197 notein(chLine);
1198 }
1199
1200 /* check if Lya alpha ever hotter than gas */
1201 if( hydro.nLyaHot > 0 )
1202 {
1203 if( hydro.TLyaMax/hydro.TeLyaMax > 1.05 )
1204 {
1205 sprintf( chLine,
1206 " !The excitation temp of Lya exceeded the electron temp, largest value was %.2eK (gas temp there was %.2eK, zone%4ld)",
1207 hydro.TLyaMax, hydro.TeLyaMax, hydro.nZTLaMax );
1208 bangin(chLine);
1209 }
1210 }
1211
1212 /* check if line absorption heating was important */
1213
1214 /* get all negative lines, check if line absorption significant heat source
1215 * this is used in "final" for energy budget print out */
1216 if( cdLine("Line",0,&SumNeg,&absint)<=0 )
1217 {
1218 fprintf( ioQQQ, " did not get sumneg cdLine\n" );
1219 ShowMe();
1221 }
1222
1223 /* this is total heating */
1224 if( cdLine("TotH",0,&GetHeat,&absint)<=0 )
1225 {
1226 fprintf( ioQQQ, " did not get GetHeat cdLine\n" );
1227 ShowMe();
1229 }
1230
1231 if( GetHeat > 0. )
1232 {
1233 SumNeg /= GetHeat;
1234 if( SumNeg > 0.1 )
1235 {
1236 sprintf( chLine,
1237 " !Line absorption heating reached %.2f%% of the global heating.",
1238 SumNeg*100. );
1239 bangin(chLine);
1240 }
1241 else if( SumNeg > 0.01 )
1242 {
1243 sprintf( chLine,
1244 " Line absorption heating reached %.2f%% of the global heating.",
1245 SumNeg*100. );
1246 notein(chLine);
1247 }
1248 }
1249
1250 /* this is check of extra lines added with g-bar */
1251 if( input.lgSetNoBuffering )
1252 {
1253 sprintf( chLine,
1254 " !NO BUFFERING command was entered - this increases exec time by LARGE amounts.");
1255 bangin(chLine);
1256 }
1257
1258 /* this is check of extra lines added with g-bar */
1259 if( thermal.GBarMax > 0.1 )
1260 {
1261 ASSERT( thermal.ipMaxExtra > 0 );
1262 strcpy( chLbl, chLineLbl(TauLine2[thermal.ipMaxExtra-1]) );
1263
1264 sprintf( chLine,
1265 " !G-bar cooling lines reached %.2f%% of the local cooling. Line=%.10s",
1266 thermal.GBarMax*100., chLbl );
1267 bangin(chLine);
1268 }
1269
1270 else if( thermal.GBarMax > 0.01 )
1271 {
1272 strcpy( chLbl, chLineLbl(TauLine2[thermal.ipMaxExtra-1]) );
1273
1274 sprintf( chLine,
1275 " G-bar cooling lines reached %.2f%% of the local cooling. Line=%.10s",
1276 thermal.GBarMax*100., chLbl );
1277 notein(chLine);
1278 }
1279
1280 /* this is check of hyperfine structure lines*/
1281 if( hyperfine.cooling_max > 0.1 )
1282 {
1283 sprintf( chLine,
1284 " !Hyperfine structure line cooling reached %.2f%% of the local cooling.",
1285 hyperfine.cooling_max*100.);
1286 bangin(chLine);
1287 }
1288
1289 else if( hyperfine.cooling_max > 0.01 )
1290 {
1291 sprintf( chLine,
1292 " Hyperfine structure line cooling reached %.2f%% of the local cooling.",
1293 hyperfine.cooling_max*100. );
1294 notein(chLine);
1295 }
1296
1297 /* line absorption heating reached more than 10% of local heating?
1298 * HeatLineMax is largest heating(1,23)/htot */
1299 if( thermal.HeatLineMax > 0.1 )
1300 {
1301 long level = -1;
1302 TransitionProxy t = FndLineHt(&level);
1303 strcpy( chLbl, chLineLbl(t) );
1304 sprintf( chLine,
1305 " !Line absorption heating reached %.2f%% of the local heating - largest by level%2ld line %.10s",
1306 thermal.HeatLineMax*100., level, chLbl );
1307 bangin(chLine);
1308 }
1309
1310 else if( thermal.HeatLineMax > 0.01 )
1311 {
1312 sprintf( chLine,
1313 " Line absorption heating reached %.2f%% of the local heating.",
1314 thermal.HeatLineMax*100. );
1315 notein(chLine);
1316 }
1317
1318 if( ionbal.CompHeating_Max > 0.05 )
1319 {
1320 sprintf( chLine,
1321 " !Bound Compton heating reached %.2f%% of the local heating.",
1322 ionbal.CompHeating_Max*100. );
1323 bangin(chLine);
1324 }
1325 else if( ionbal.CompHeating_Max > 0.01 )
1326 {
1327 sprintf( chLine,
1328 " Bound Compton heating reached %.2f%% of the local heating.",
1329 ionbal.CompHeating_Max*100. );
1330 notein(chLine);
1331 }
1332
1333 /* check whether any lines in the iso sequences mased */
1334 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1335 {
1336 for( nelem=ipISO; nelem<LIMELM; ++nelem )
1337 {
1338 if( dense.lgElmtOn[nelem] )
1339 {
1340 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
1341 long int nmax = iso_sp[ipISO][nelem].numLevels_local;
1342
1343 /* minus one here is to exclude highest level */
1344 for( ipHi=1; ipHi < nmax - 1; ++ipHi )
1345 {
1346 for( ipLo=0; ipLo < ipHi; ++ipLo )
1347 {
1348 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1349 continue;
1350
1351 /* did the line mase */
1352 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() < -0.1 )
1353 {
1354 sprintf( chLine,
1355 " !Some iso-structure lines mased: %s-like %s, line %li-%li had optical depth %.2e",
1356 elementnames.chElementSym[ipISO],
1357 elementnames.chElementNameShort[nelem],
1358 ipHi , ipLo ,
1359 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() );
1360 bangin(chLine);
1361 }
1362 }
1363 }
1364 }
1365 }
1366 }
1367
1368 if( dense.gas_phase[ipHYDROGEN] < 1e7 )
1369 {
1370 /* check on IR fine structure lines - not necessary if dense since will be in LTE */
1371 lgThick = false;
1372 tauneg = 0.;
1373 alpha = 0.;
1374 /* loop from 3, since 0 is dummy, 1 and 2 are spin-flip transitions of H and He */
1375 for( i=3; i <= nLevel1; i++ )
1376 {
1377 /* define IR as anything longward of 1 micron */
1378 if( TauLines[i].EnergyWN() < 10000. )
1379 {
1380 if( TauLines[i].Emis().TauIn() > 1. )
1381 {
1382 lgThick = true;
1383 alpha = MAX2(alpha,(double)TauLines[i].Emis().TauIn());
1384 }
1385 else if( TauLines[i].Emis().TauIn() < (realnum)tauneg )
1386 {
1387 tauneg = TauLines[i].Emis().TauIn();
1388 strcpy( chLbl, chLineLbl(TauLines[i]) );
1389 }
1390 }
1391 }
1392 /* now print results, were any fine structure lines optically thick? */
1393 if( lgThick )
1394 {
1395 sprintf( chLine,
1396 " !Some infrared fine structure lines are optically thick: largest tau was %.2e",
1397 alpha );
1398 bangin(chLine);
1399 }
1400 /* did any fine structure lines mase? */
1401 if( tauneg < -0.01 )
1402 {
1403 sprintf( chLine,
1404 " !Some fine structure lines mased: line %s had optical depth %.2e",
1405 chLbl, tauneg );
1406 bangin(chLine);
1407 }
1408 }
1409
1410 /* were any other lines masing? */
1411 tauneg = 0.;
1412 alpha = 0.;
1413 for( i=1; i <= nLevel1; i++ )
1414 {
1415 /* define UV as anything shortward of 1 micron */
1416 if( TauLines[i].EnergyWN() >= 10000. )
1417 {
1418 if( TauLines[i].Emis().TauIn() < (realnum)tauneg )
1419 {
1420 tauneg = TauLines[i].Emis().TauIn();
1421 strcpy( chLbl, chLineLbl(TauLines[i]) );
1422 }
1423 }
1424 }
1425
1426 /* did any level1 lines mase? */
1427 if( tauneg < -0.01 )
1428 {
1429 sprintf( chLine,
1430 " !Some level1 lines mased: most negative ion and tau were: %s %.2e",
1431 chLbl, tauneg );
1432 bangin(chLine);
1433 }
1434
1435 /* this is check that at least a second iteration was done with sphere static,
1436 * the test is overridden with the (OK) option on the sphere static command,
1437 * which sets geometry.lgStaticNoIt true */
1438 if( geometry.lgStatic && iterations.lgLastIt && (iteration == 1) &&
1439 !geometry.lgStaticNoIt)
1440 {
1441 sprintf( chLine, " C-I must iterate when SPHERE STATIC is set." );
1442 caunin(chLine);
1443 iterations.lgIterAgain = true;
1444 }
1445
1446 /* caution if continuum is punched but only one iteration performed */
1447 if( save.lgPunContinuum && iteration == 1 && iterations.lgLastIt)
1448 {
1449 sprintf( chLine, " C-I must iterate when save continuum output is done." );
1450 caunin(chLine);
1451 iterations.lgIterAgain = true;
1452 }
1453
1455 /* how important was induced two photon?? */
1456 {
1457 two_photon& tnu = iso_sp[ipH_LIKE][ipHYDROGEN].TwoNu[0];
1458 if( tnu.induc_dn_max > 1. )
1459 {
1460 sprintf( chLine, " !Rate of induced H 2-photon emission reached %.2e s^-1",
1461 tnu.induc_dn_max );
1462 bangin(chLine);
1463 }
1464 else if( tnu.induc_dn_max > 0.01 )
1465 {
1466 sprintf( chLine, " Rate of induced H 2-photon emission reached %.2e s^-1",
1467 tnu.induc_dn_max );
1468 notein(chLine);
1469 }
1470 }
1471
1472 /* how important was induced recombination? */
1473 if( hydro.FracInd > 0.01 )
1474 {
1475 sprintf( chLine,
1476 " Induced recombination was %5.1f%% of the total for H level%3ld",
1477 hydro.FracInd*100., hydro.ndclev );
1478 notein(chLine);
1479 }
1480
1481 if( hydro.fbul > 0.01 )
1482 {
1483 sprintf( chLine,
1484 " Stimulated emission was%6.1f%% of the total for H transition%3ld -%3ld",
1485 hydro.fbul*100., hydro.nbul + 1, hydro.nbul );
1486 notein(chLine);
1487 }
1488
1489 /* check whether Fe II destruction of La was important - entry into lines stack
1490 * is in prt_lines_hydro.c */
1491 if( cdLine("Fe 2",1216,&fedest,&absint)<=0 )
1492 {
1493 fprintf( ioQQQ, " Did not find Fe II Lya\n" );
1494 ShowMe();
1496 }
1497
1498 /* find total Lya for comparison */
1499 if( cdLine("TOTL",1215.68f,&relhm,&absint)<=0 )
1500 {
1501 fprintf( ioQQQ, " Did not find Lya\n" );
1502 ShowMe();
1504 }
1505
1506 if( relhm > 0. )
1507 {
1508 ratio = fedest/(fedest + relhm);
1509 if( ratio > 0.1 )
1510 {
1511 sprintf( chLine, " !Fe II destruction of Ly-a removed %.1f%% of the line.",
1512 ratio *100.);
1513 bangin(chLine);
1514 }
1515 else if( ratio > 0.01 )
1516 {
1517 sprintf( chLine, " Fe II destruction of Ly-a removed %.1f%% of the line.",
1518 ratio );
1519 notein(chLine);
1520 }
1521 }
1522
1523 if( cdLine("H-CT",6563,&relhm,&absint)<=0 )
1524 {
1525 fprintf( ioQQQ, " Comment did not find H-CT H-alpha\n" );
1526 ShowMe();
1528 }
1529
1530 if( HBeta > 0. )
1531 {
1532 if( relhm/HBeta > 0.01 )
1533 {
1534 sprintf( chLine,
1535 " !Mutual neutralization production of H-alpha was significant." );
1536 bangin(chLine);
1537 }
1538 }
1539
1540 /* note about very high population in H n=2 rel to ground, set in hydrogenic */
1541 if( hydro.lgHiPop2 )
1542 {
1543 sprintf( chLine,
1544 " The population of H n=2 reached %.2e relative to the ground state.",
1545 hydro.pop2mx );
1546 notein(chLine);
1547 }
1548
1549 /* check where diffuse emission error */
1550 for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
1551 {
1552 for( nelem=ipISO; nelem < LIMELM; nelem++ )
1553 {
1554 if( iso_sp[ipISO][nelem].CaseBCheck > 1.5 )
1555 {
1556 sprintf( chLine,
1557 " Ratio of computed diffuse emission to case B reached %g for iso %li element %li",
1558 iso_sp[ipISO][nelem].CaseBCheck , ipISO , nelem+1 );
1559 notein(chLine);
1560 }
1561 }
1562 }
1563
1564 /* check whether electrons were relativistic */
1565 if( thermal.thist > 1e9 )
1566 {
1567 /* >>chng 06 feb 19, from 5e9 K for warning to 1e10K. add test case at 1e10K
1568 * and don't want warning in test suite. nothing is wrong at this temp - eeff
1569 * is in correctly for relativistic temps and will eventually dominate cooling */
1570 if( thermal.thist > 1.0001e10 )
1571 {
1572 sprintf( chLine, " W-Electrons were relativistic; High TE=%.2e",
1573 thermal.thist );
1574 warnin(chLine);
1575 }
1576 else
1577 {
1578 sprintf( chLine, " C-Electrons were mildly relativistic; High TE=%.2e",
1579 thermal.thist );
1580 caunin(chLine);
1581 }
1582 }
1583
1584 /* check on timescale for photoerosion of elements */
1585 rate = timesc.TimeErode*2e-26;
1586 if( rate > 1e-35 )
1587 {
1588 /* 2E-26 is roughly cross section for photoerosion
1589 * see
1590 * >>refer all photoerode Boyd, R., & Ferland, G.J. ApJ, 318, L21. */
1591 ts = (1./rate)/3e7;
1592 if( ts < 1e3 )
1593 {
1594 sprintf( chLine, " !Timescale-photoerosion of Fe=%.2e yr",
1595 ts );
1596 bangin(chLine);
1597 }
1598 else if( ts < 1e9 )
1599 {
1600 sprintf( chLine, " Timescale-photoerosion of Fe=%.2e yr",
1601 ts );
1602 notein(chLine);
1603 }
1604 }
1605
1606 /* check whether Compton heating was significant */
1607 comfrc = rfield.comtot/SDIV(thermal.power);
1608 if( comfrc > 0.01 )
1609 {
1610 sprintf( chLine, " Compton heating was %5.1f%% of the total.",
1611 comfrc*100. );
1612 notein(chLine);
1613 }
1614
1615 /* check on relative importance of induced Compton heating */
1616 if( comfrc > 0.01 && rfield.cinrat > 0.05 )
1617 {
1618 sprintf( chLine,
1619 " !Induced Compton heating was %.2e of the total Compton heating.",
1620 rfield.cinrat );
1621 bangin(chLine);
1622 }
1623
1624 /* check whether equilibrium timescales are short rel to Hubble time */
1625 if( timesc.tcmptn > 5e17 )
1626 {
1627 if( comfrc > 0.05 )
1628 {
1629 sprintf( chLine,
1630 " C-Compton cooling is significant and the equilibrium timescale (%.2e s) is longer than the Hubble time.",
1631 timesc.tcmptn );
1632 caunin(chLine);
1633 }
1634 else
1635 {
1636 sprintf( chLine,
1637 " Compton cooling equilibrium timescale (%.2e s) is longer than Hubble time.",
1638 timesc.tcmptn );
1639 notein(chLine);
1640 }
1641 }
1642
1643 if( timesc.time_therm_long > 5e17 )
1644 {
1645 sprintf( chLine,
1646 " C-Thermal equilibrium timescale, %.2e s, longer than Hubble time; this cloud is not time-steady.",
1647 timesc.time_therm_long );
1648 caunin(chLine);
1649 }
1650
1651 /* check whether model large relative to Jeans length
1652 * DMEAN is mean density (gm per cc)
1653 * mean temp is weighted by mass density */
1654 if( log10(radius.depth) > colden.rjnmin )
1655 {
1656 /* AJMIN is minimum Jeans mass, log in grams */
1657 aj = pow(10.,colden.ajmmin - log10(SOLAR_MASS));
1658 if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
1659 {
1660 sprintf( chLine,
1661 " C-Cloud thicker than smallest Jeans length=%8.2ecm; stability problems? (smallest Jeans mass=%8.2eMo)",
1662 pow((realnum)10.f,colden.rjnmin), aj );
1663 caunin(chLine);
1664 }
1665 else
1666 {
1667 sprintf( chLine,
1668 " Cloud thicker than smallest Jeans length=%8.2ecm; stability problems? (smallest Jeans mass=%8.2eMo)",
1669 pow((realnum)10.f,colden.rjnmin), aj );
1670 notein(chLine);
1671 }
1672 }
1673
1674 /* check whether grains too hot to survive */
1675 for( size_t nd=0; nd < gv.bin.size(); nd++ )
1676 {
1677 if( gv.bin[nd]->TeGrainMax > gv.bin[nd]->Tsublimat )
1678 {
1679 sprintf( chLine,
1680 " W-Maximum temperature of grain%-12.12s was %.2eK, above its sublimation temperature, %.2eK.",
1681 gv.bin[nd]->chDstLab, gv.bin[nd]->TeGrainMax,
1682 gv.bin[nd]->Tsublimat );
1683 warnin(chLine);
1684 }
1685 else if( gv.bin[nd]->TeGrainMax > gv.bin[nd]->Tsublimat* 0.9 )
1686 {
1687 sprintf( chLine,
1688 " C-Maximum temperature of grain%-12.12s was %.2eK, near its sublimation temperature, %.2eK.",
1689 gv.bin[nd]->chDstLab, gv.bin[nd]->TeGrainMax,
1690 gv.bin[nd]->Tsublimat );
1691 caunin(chLine);
1692 }
1693 }
1694
1695 if( gv.lgNegGrnDrg )
1696 {
1697 sprintf( chLine, " !Grain drag force <0." );
1698 bangin(chLine);
1699 }
1700
1701 /* largest relative number of electrons donated by grains */
1702 if( gv.GrnElecDonateMax > 0.05 )
1703 {
1704 sprintf( chLine,
1705 " !Grains donated %5.1f%% of the total electrons in some regions.",
1706 gv.GrnElecDonateMax*100. );
1707 bangin(chLine);
1708 }
1709 else if( gv.GrnElecDonateMax > 0.005 )
1710 {
1711 sprintf( chLine,
1712 " Grains donated %5.1f%% of the total electrons in some regions.",
1713 gv.GrnElecDonateMax*100. );
1714 notein(chLine);
1715 }
1716
1717 /* largest relative number of electrons on grain surface */
1718 if( gv.GrnElecHoldMax > 0.05 )
1719 {
1720 sprintf( chLine,
1721 " !Grains contained %5.1f%% of the total electrons in some regions.",
1722 gv.GrnElecHoldMax*100. );
1723 bangin(chLine);
1724 }
1725 else if( gv.GrnElecHoldMax > 0.005 )
1726 {
1727 sprintf( chLine,
1728 " Grains contained %5.1f%% of the total electrons in some regions.",
1729 gv.GrnElecHoldMax*100. );
1730 notein(chLine);
1731 }
1732
1733 /* is photoelectric heating of gas by photoionization of grains important */
1734 if( gv.dphmax > 0.5 )
1735 {
1736 sprintf( chLine,
1737 " !Local grain-gas photoelectric heating rate reached %5.1f%% of the total.",
1738 gv.dphmax*100. );
1739 bangin(chLine);
1740 }
1741 else if( gv.dphmax > 0.05 )
1742 {
1743 sprintf( chLine,
1744 " Local grain-gas photoelectric heating rate reached %5.1f%% of the total.",
1745 gv.dphmax*100. );
1746 notein(chLine);
1747 }
1748
1749 if( gv.TotalDustHeat/SDIV(thermal.power) > 0.01 )
1750 {
1751 sprintf( chLine,
1752 " Global grain photoelectric heating of gas was%5.1f%% of the total.",
1753 gv.TotalDustHeat/thermal.power*100. );
1754 notein(chLine);
1755 if( gv.TotalDustHeat/thermal.power > 0.25 )
1756 {
1757 sprintf( chLine,
1758 " !Grain photoelectric heating is VERY important." );
1759 bangin(chLine);
1760 }
1761 }
1762
1763 /* grain-gas collisional cooling of gas */
1764 if( gv.dclmax > 0.05 )
1765 {
1766 sprintf( chLine,
1767 " Local grain-gas cooling of gas rate reached %5.1f%% of the total.",
1768 gv.dclmax*100. );
1769 notein(chLine);
1770 }
1771
1772 /* check how H2 chemistry network performed */
1773 if( h2.renorm_max > 1.05 )
1774 {
1775 if( h2.renorm_max > 1.2 )
1776 {
1777 sprintf( chLine,
1778 " !The large H2 molecule - main chemistry network renormalization factor reached %.2f.",
1779 h2.renorm_max);
1780 bangin(chLine);
1781 }
1782 else
1783 {
1784 sprintf( chLine,
1785 " The large H2 molecule - main chemistry network renormalization factor reached %.2f.",
1786 h2.renorm_max);
1787 notein(chLine);
1788 }
1789 }
1790 if( h2.renorm_min < 0.95 )
1791 {
1792 if( h2.renorm_min < 0.8 )
1793 {
1794 sprintf( chLine,
1795 " !The large H2 molecule - main chemistry network renormalization factor reached %.2f.",
1796 h2.renorm_min);
1797 bangin(chLine);
1798 }
1799 else
1800 {
1801 sprintf( chLine,
1802 " The large H2 molecule - main chemistry network renormalization factor reached %.2f.",
1803 h2.renorm_min);
1804 notein(chLine);
1805 }
1806 }
1807
1808 /* check whether photodissociation of H_2^+ molecular ion was important */
1809 if( hmi.h2pmax > 0.10 )
1810 {
1811 sprintf( chLine,
1812 " !The local H2+ photodissociation heating rate reached %5.1f%% of the total heating.",
1813 hmi.h2pmax*100. );
1814 bangin(chLine);
1815 }
1816
1817 else if( hmi.h2pmax > 0.01 )
1818 {
1819 sprintf( chLine,
1820 " The local H2+ photodissociation heating rate reached %.1f%% of the total heating.",
1821 hmi.h2pmax*100. );
1822 notein(chLine);
1823 }
1824
1825 /* check whether photodissociation of molecular hydrogen (H2)was important */
1826 if( hmi.h2dfrc > 0.1 )
1827 {
1828 sprintf( chLine,
1829 " !The local H2 photodissociation heating rate reached %.1f%% of the total heating.",
1830 hmi.h2dfrc*100. );
1831 bangin(chLine);
1832 }
1833 else if( hmi.h2dfrc > 0.01 )
1834 {
1835 sprintf( chLine,
1836 " The local H2 photodissociation heating rate reached %.1f%% of the total heating.",
1837 hmi.h2dfrc*100. );
1838 notein(chLine);
1839 }
1840
1841 /* check whether cooling by molecular hydrogen (H2) was important */
1842 if( hmi.h2line_cool_frac > 0.1 )
1843 {
1844 sprintf( chLine,
1845 " !The local H2 cooling rate reached %.1f%% of the local cooling.",
1846 hmi.h2line_cool_frac*100. );
1847 bangin(chLine);
1848 }
1849 else if( hmi.h2line_cool_frac > 0.01 )
1850 {
1851 sprintf( chLine,
1852 " The local H2 cooling rate reached %.1f%% of the local cooling.",
1853 hmi.h2line_cool_frac*100. );
1854 notein(chLine);
1855 }
1856
1857 if( hmi.h2dtot/SDIV(thermal.power) > 0.01 )
1858 {
1859 sprintf( chLine,
1860 " Global H2 photodissociation heating of gas was %.1f%% of the total heating.",
1861 hmi.h2dtot/thermal.power*100. );
1862 notein(chLine);
1863 if( hmi.h2dtot/thermal.power > 0.25 )
1864 {
1865 sprintf( chLine, " H2 photodissociation heating is VERY important." );
1866 notein(chLine);
1867 }
1868 }
1869
1870 /* check whether photodissociation of carbon monoxide (co) was important */
1871 if( co.codfrc > 0.25 )
1872 {
1873 sprintf( chLine,
1874 " !Local CO photodissociation heating rate reached %.1f%% of the total.",
1875 co.codfrc*100. );
1876 bangin(chLine);
1877 }
1878 else if( co.codfrc > 0.05 )
1879 {
1880 sprintf( chLine,
1881 " Local CO photodissociation heating rate reached %.1f%% of the total.",
1882 co.codfrc*100. );
1883 notein(chLine);
1884 }
1885
1886 if( co.codtot/SDIV(thermal.power) > 0.01 )
1887 {
1888 sprintf( chLine,
1889 " Global CO photodissociation heating of gas was %.1f%% of the total.",
1890 co.codtot/thermal.power*100. );
1891 notein(chLine);
1892 if( co.codtot/thermal.power > 0.25 )
1893 {
1894 sprintf( chLine, " CO photodissociation heating is VERY important." );
1895 notein(chLine);
1896 }
1897 }
1898
1899 if( thermal.lgEdnGTcm )
1900 {
1901 sprintf( chLine,
1902 " Energy density of radiation field was greater than the Compton temperature. Is this physical?" );
1903 notein(chLine);
1904 }
1905
1906 /* was cooling due to induced recombination important? */
1907 if( hydro.cintot/SDIV(thermal.power) > 0.01 )
1908 {
1909 sprintf( chLine, " Induced recombination cooling was %.1f%% of the total.",
1910 hydro.cintot/thermal.power*100. );
1911 notein(chLine);
1912 }
1913
1914 /* check whether free-free heating was significant */
1915 if( thermal.FreeFreeTotHeat/SDIV(thermal.power) > 0.1 )
1916 {
1917 sprintf( chLine, " !Free-free heating was %.1f%% of the total.",
1918 thermal.FreeFreeTotHeat/thermal.power*100. );
1919 bangin(chLine);
1920 }
1921 else if( thermal.FreeFreeTotHeat/SDIV(thermal.power) > 0.01 )
1922 {
1923 sprintf( chLine, " Free-free heating was %.1f%% of the total.",
1924 thermal.FreeFreeTotHeat/thermal.power*100. );
1925 notein(chLine);
1926 }
1927
1928 /* was heating due to H- absorption important? */
1929 if( hmi.hmitot/SDIV(thermal.power) > 0.01 )
1930 {
1931 sprintf( chLine, " H- absorption heating was %.1f%% of the total.",
1932 hmi.hmitot/SDIV(thermal.power)*100. );
1933 notein(chLine);
1934 }
1935
1936 /* water destruction rate was zero */
1937 if( mole_global.lgH2Ozer )
1938 {
1939 sprintf( chLine, " Water destruction rate zero." );
1940 notein(chLine);
1941 }
1942
1943 /* numerical instability in matrix inversion routine */
1944 if( atoms.nNegOI > 0 )
1945 {
1946 sprintf( chLine, " C-O I negative level populations %ld times.",
1947 atoms.nNegOI );
1948 caunin(chLine);
1949 }
1950
1951 /* check for negative optical depths,
1952 * optical depth in excited state helium lines */
1953 small = 0.;
1954 imas = 0;
1955 isav = 0;
1956 j = 0;
1957 for( nelem=0; nelem<LIMELM; ++nelem )
1958 {
1959 if( dense.lgElmtOn[nelem] )
1960 {
1961 /* >>chng 06 aug 28, from numLevels_max to _local. */
1962 for( ipLo=ipH2p; ipLo < (iso_sp[ipH_LIKE][nelem].numLevels_local - 1); ipLo++ )
1963 {
1964 /* >>chng 06 aug 28, from numLevels_max to _local. */
1965 for( ipHi=ipLo + 1; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_local; ipHi++ )
1966 {
1967 if( iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1968 continue;
1969
1970 if( iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauIn() < (realnum)small )
1971 {
1972 small = iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauIn();
1973 imas = ipHi;
1974 j = ipLo;
1975 isav = nelem;
1976 }
1977 }
1978 }
1979 }
1980 }
1981
1982 if( small < -0.05 )
1983 {
1984 sprintf( chLine,
1985 " !Some hydrogenic lines mased, species was %2s%2ld, smallest tau was %.2e, transition %li-%li",
1986 elementnames.chElementSym[isav],
1987 isav+1,small, imas , j );
1988 bangin(chLine);
1989 }
1990
1991 /* check for negative opacities */
1992 if( opac.lgOpacNeg )
1993 {
1994 sprintf( chLine, " !Some opacities were negative - the SET NEGOPC command will save which ones." );
1995 bangin(chLine);
1996 }
1997
1998 /* now check continua */
1999 small = 0.;
2000 imas = 0;
2001 isav = 0;
2002 for( nelem=0; nelem<LIMELM; ++nelem )
2003 {
2004 if( dense.lgElmtOn[nelem] )
2005 {
2006 /* >>chng 06 aug 28, from numLevels_max to _local. */
2007 for( i=0; i < iso_sp[ipH_LIKE][nelem].numLevels_local; i++ )
2008 {
2009 if( opac.TauAbsGeo[0][iso_sp[ipH_LIKE][nelem].fb[i].ipIsoLevNIonCon-1] < -0.001 )
2010 {
2011 small = MIN2(small,(double)opac.TauAbsGeo[0][iso_sp[ipH_LIKE][nelem].fb[i].ipIsoLevNIonCon-1]);
2012 imas = i;
2013 isav = nelem;
2014 }
2015 }
2016 }
2017 }
2018
2019 if( small < -0.05 )
2020 {
2021 sprintf( chLine, " !Some hydrogenic (%2s%2ld) continua optical depths were negative; smallest=%.2e level=%3ld",
2022 elementnames.chElementSym[isav],
2023 isav+1,
2024 small, imas );
2025 bangin(chLine);
2026 }
2027
2028 /* check whether any continuum optical depths are negative */
2029 nneg = 0;
2030 tauneg = 0.;
2031 freqn = 0.;
2032 for( i=0; i < rfield.nflux; i++ )
2033 {
2034 if( opac.TauAbsGeo[0][i] < -0.001 )
2035 {
2036 nneg += 1;
2037 /* only remember the smallest freq, and most neg optical depth */
2038 if( nneg == 1 )
2039 freqn = rfield.anu[i];
2040 tauneg = MIN2(tauneg,(double)opac.TauAbsGeo[0][i]);
2041 }
2042 }
2043
2044 if( nneg > 0 )
2045 {
2046 sprintf( chLine, " !Some continuous optical depths <0. The lowest freq was %.3e Ryd, and a total of%4ld",
2047 freqn, nneg );
2048 bangin(chLine);
2049 sprintf( chLine, " !The smallest optical depth was %.2e",
2050 tauneg );
2051 bangin(chLine);
2052 }
2053
2054 /* say if Balmer continuum optically thick */
2055 if( opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1] > 0.05 )
2056 {
2057 sprintf( chLine, " The Balmer continuum optical depth was %.2e.",
2058 opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1] );
2059 notein(chLine);
2060 }
2061
2062 /* was correction for stimulated emission significant? */
2063 if( opac.stimax[0] > 0.02 && opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] > 0.2 )
2064 {
2065 sprintf( chLine, " The Lyman continuum stimulated emission correction to optical depths reached %.2e.",
2066 opac.stimax[0] );
2067 notein(chLine);
2068 }
2069 else if( opac.stimax[1] > 0.02 && opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1] > 0.1 )
2070 {
2071 sprintf( chLine, " The Balmer continuum stimulated emission correction to optical depths reached %.2e.",
2072 opac.stimax[1] );
2073 notein(chLine);
2074 }
2075
2076 /* say if Paschen continuum optically thick */
2077 if( opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[3].ipIsoLevNIonCon-1] > 0.2 )
2078 {
2079 sprintf( chLine,
2080 " The Paschen continuum optical depth was %.2e.",
2081 opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[3].ipIsoLevNIonCon-1] );
2082 notein(chLine);
2083 }
2084
2085 /* some comments about near IR total optical depth */
2086 if( opac.TauAbsGeo[0][0] > 1. )
2087 {
2088 sprintf( chLine,
2089 " The continuum optical depth at the lowest energy considered (%.3e Ryd) was %.3e.",
2090 rfield.anu[0], opac.TauAbsGeo[0][0] );
2091 notein(chLine);
2092 }
2093
2094 /* comment if optical depth to Rayleigh scattering is big
2095 * cs from VAL 76 */
2096 if( colden.colden[ipCOL_H0]*7e-24 > 0.01 )
2097 {
2098 sprintf( chLine,
2099 " The optical depth to Rayleigh scattering at 1300A is %.2e",
2100 colden.colden[ipCOL_H0]*6.71e-24 );
2101 notein(chLine);
2102 }
2103
2104 if( colden.colden[ipCOL_H2p]*7e-18 > 0.1 )
2105 {
2106 sprintf( chLine,
2107 " !The optical depth to the H2+ molecular ion is %.2e",
2108 colden.colden[ipCOL_H2p]*7e-18 );
2109 bangin(chLine);
2110 }
2111 else if( colden.colden[ipCOL_H2p]*7e-18 > 0.01 )
2112 {
2113 sprintf( chLine,
2114 " The optical depth to the H2+ molecular ion is %.2e",
2115 colden.colden[ipCOL_H2p]*7e-18 );
2116 notein(chLine);
2117 }
2118
2119 /* warn if optically thick to H- absorption */
2120 if( opac.thmin > 0.1 )
2121 {
2122 sprintf( chLine,
2123 " !Optical depth to negative hydrogen ion is %.2e",
2124 opac.thmin );
2125 bangin(chLine);
2126 }
2127 else if( opac.thmin > 0.01 )
2128 {
2129 sprintf( chLine,
2130 " Optical depth to negative hydrogen ion is %.2e",
2131 opac.thmin );
2132 notein(chLine);
2133 }
2134
2135 /* say if 3-body recombination coefficient function outside range of validity
2136 * tripped if te/z**2 < 100 or approx 10**13: => effect >50% of radiative
2137 * other integers defined in source for da */
2138 if( ionbal.ifail > 0 && ionbal.ifail <= 10 )
2139 {
2140 sprintf( chLine,
2141 " 3 body recombination coefficient outside range %ld", ionbal.ifail );
2142 notein(chLine);
2143 }
2144 else if( ionbal.ifail > 10 )
2145 {
2146 sprintf( chLine,
2147 " C-3 body recombination coefficient outside range %ld", ionbal.ifail );
2148 caunin(chLine);
2149 }
2150
2151 /* check whether energy density less than background */
2152 if( phycon.TEnerDen < 2.6 )
2153 {
2154 sprintf( chLine,
2155 " !Incident radiation field energy density is less than 2.7K. Add background with CMB command." );
2156 bangin(chLine);
2157 }
2158
2159 /* check whether CMB set at all */
2160 if( !rfield.lgCMB_set )
2161 {
2162 sprintf( chLine,
2163 " !The CMB was not included. This is added with the CMB command." );
2164 bangin(chLine);
2165 }
2166
2167 /* incident radiation field is less than background Habing ISM field */
2168 if( rfield.lgHabing )
2169 {
2170 sprintf( chLine,
2171 " !The intensity of the incident radiation field is less than 10 times the Habing diffuse ISM field. Is this OK?" );
2172 bangin(chLine);
2173 sprintf( chLine,
2174 " ! Consider adding diffuse ISM emission with TABLE ISM command." );
2175 bangin(chLine);
2176 }
2177
2178 /* some things dealing with molecules, or molecule formation */
2179
2180 /* if C/O > 1 then chemistry will be carbon dominated rather than oxygen dominated */
2181 if( dense.lgElmtOn[ipOXYGEN] && dense.lgElmtOn[ipCARBON] )
2182 {
2183 if( dense.gas_phase[ipCARBON]/dense.gas_phase[ipOXYGEN] > 1. )
2184 {
2185 sprintf( chLine, " !The C/O abundance ratio, %.1f, is greater than unity. The chemistry will be carbon dominated.",
2186 dense.gas_phase[ipCARBON]/dense.gas_phase[ipOXYGEN] );
2187 bangin(chLine);
2188 }
2189 }
2190
2191 lgLots_of_moles = false;
2192 lgLotsSolids = false;
2193 /* largest fraction in any molecule */
2194 for( i=0; i<mole_global.num_calc; ++i )
2195 {
2196 if( mole.species[i].location == NULL && ( mole_global.list[i]->parentLabel.empty() || mole_global.list[i]->charge<0 ) )
2197 {
2198 if( mole.species[i].xFracLim > 0.1 )
2199 {
2200 sprintf( chLine, " !The fraction of %s in %s reached %.1f%% at some point in the cloud.",
2201 elementnames.chElementSym[mole.species[i].nAtomLim],
2202 mole_global.list[i]->label.c_str(),
2203 mole.species[i].xFracLim*100. );
2204 bangin(chLine);
2205 lgLots_of_moles = true;
2206 /* check whether molecules are on grains */
2207 if( !mole_global.list[i]->lgGas_Phase )
2208 lgLotsSolids = true;
2209 }
2210 else if( mole.species[i].xFracLim>0.01 )
2211 {
2212 sprintf( chLine, " The fraction of %s in %s reached %.2f%% at some point in the cloud.",
2213 elementnames.chElementSym[mole.species[i].nAtomLim],
2214 mole_global.list[i]->label.c_str(),
2215 mole.species[i].xFracLim*100. );
2216 notein(chLine);
2217 lgLots_of_moles = true;
2218 /* check whether molecules are on grains */
2219 if( !mole_global.list[i]->lgGas_Phase )
2220 lgLotsSolids = true;
2221 }
2222 else if( mole.species[i].xFracLim > 1e-3 )
2223 {
2224 sprintf( chLine, " The fraction of %s in %s reached %.3f%% at some point in the cloud.",
2225 elementnames.chElementSym[mole.species[i].nAtomLim],
2226 mole_global.list[i]->label.c_str(),
2227 mole.species[i].xFracLim*100. );
2228 notein(chLine);
2229 /* check whether molecules are on grains */
2230 if( !mole_global.list[i]->lgGas_Phase )
2231 lgLotsSolids = true;
2232 }
2233 }
2234 }
2235
2236 /* generate comment if molecular fraction was significant but some heavy elements are turned off */
2237 if( lgLots_of_moles )
2238 {
2239 /* find all elements that are turned off */
2240 for(nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
2241 {
2242 /* >>chng 05 dec 23, add mole_global.lgElem_in_chemistry */
2243 if( !dense.lgElmtOn[nelem] )
2244 {
2245 /* this triggers if element turned off but it is part of co chem net */
2246 sprintf( chLine,
2247 " C-Molecules are important, but %s, part of the chemistry network, is turned off.",
2248 elementnames.chElementName[nelem] );
2249 caunin(chLine);
2250 }
2251# if 0
2252 /* this element has been turned off - now check if part of chemistry */
2253 for( i=NUM_HEAVY_MOLEC+NUM_ELEMENTS; i<NUM_COMOLE_CALC; ++i )
2254 {
2255 if( nelem==mole_global.list[i].nelem_den )
2256 {
2257 /* this triggers if element turned off but it is part of co chem net */
2258 sprintf( chLine,
2259 " C-Molecules are important, but %s, part of the chemistry network, is turned off.",
2260 elementnames.chElementName[nelem] );
2261 caunin(chLine);
2262 }
2263 }
2264# endif
2265 }
2266 }
2267
2268 /* say if lots of molecules on grains,
2269 * molecules with labels that *GR */
2270 if( lgLotsSolids )
2271 {
2272 sprintf( chLine, " !A significant amount of molecules condensed onto grain surfaces." );
2273 bangin(chLine);
2274 sprintf( chLine, " !These are the molecular species with \"grn\" above." );
2275 bangin(chLine);
2276 }
2277
2278 /* bremsstrahlung optical depth */
2279 if( rfield.EnergyBremsThin > 0.09 )
2280 {
2281 sprintf( chLine, " !The cloud is optically thick at optical wavelengths, extending to %.3e Ryd =%.3eA",
2282 rfield.EnergyBremsThin, RYDLAM/rfield.EnergyBremsThin );
2283 bangin(chLine);
2284 }
2285 else if( rfield.EnergyBremsThin > 0.009 )
2286 {
2287 sprintf( chLine, " The continuum of the computed structure may be optically thick in the near infrared." );
2288 notein(chLine);
2289 }
2290
2291 /* did model run away to very large radius? */
2292 if( radius.Radius > 1e23 && radius.Radius/radius.rinner > 10. )
2293 {
2294 sprintf( chLine, " Is an outer radius of %.2e reasonable?",
2295 radius.Radius );
2296 notein(chLine);
2297 }
2298
2299 /* following set true in RT_line_one_tauinc if maser capped at tau = -1 */
2300 if( rt.lgMaserCapHit )
2301 {
2302 sprintf( chLine, " Laser maser optical depths capped in RT_line_one_tauinc." );
2303 notein(chLine);
2304 }
2305
2306 /* following set true in adius_next if maser cap set dr */
2307 if( rt.lgMaserSetDR )
2308 {
2309 sprintf( chLine, " !Line maser set zone thickness in some zones." );
2310 bangin(chLine);
2311 }
2312
2313 /* lgPradCap is true if radiation pressure was capped on first iteration
2314 * also check that this is a constant total pressure model */
2315 if( (pressure.lgPradCap && (strcmp(dense.chDenseLaw,"CPRE") == 0)) &&
2316 pressure.lgPres_radiation_ON )
2317 {
2318 sprintf( chLine, " Radiation pressure kept below gas pressure on this iteration." );
2319 notein(chLine);
2320 }
2321
2322 if( pressure.RadBetaMax > 0.25 )
2323 {
2324 if( pressure.ipPradMax_line == 0 )
2325 {
2326 sprintf( chLine,
2327 " !The ratio of radiation to gas pressure reached %.2e at zone %li. Caused by Lyman alpha.",
2328 pressure.RadBetaMax,
2329 pressure.ipPradMax_nzone);
2330 bangin(chLine);
2331 }
2332 else
2333 {
2334 sprintf( chLine,
2335 " !The ratio of radiation to gas pressure reached %.2e at zone %li. "
2336 "Caused by line number %ld, label %s",
2337 pressure.RadBetaMax,
2338 pressure.ipPradMax_nzone,
2339 pressure.ipPradMax_line,
2340 pressure.chLineRadPres );
2341 bangin(chLine);
2342 }
2343 }
2344
2345 else if( pressure.RadBetaMax > 0.025 )
2346 {
2347 if( pressure.ipPradMax_line == 0 )
2348 {
2349 sprintf( chLine,
2350 " The ratio of radiation to gas pressure reached %.2e at zone %li. Caused by Lyman alpha.",
2351 pressure.RadBetaMax,
2352 pressure.ipPradMax_nzone);
2353 notein(chLine);
2354 }
2355 else
2356 {
2357 sprintf( chLine,
2358 " The ratio of radiation to gas pressure reached %.2e at zone %li. "
2359 "Caused by line number %ld, label %s",
2360 pressure.RadBetaMax,
2361 pressure.ipPradMax_nzone,
2362 pressure.ipPradMax_line,
2363 pressure.chLineRadPres );
2364 notein(chLine);
2365 }
2366 }
2367
2368 if( opac.telec >= 5. )
2369 {
2370 sprintf( chLine, " W-The model is optically thick to electron "
2371 "scattering; tau=%.2e Cloudy is NOT intended for this regime.",
2372 opac.telec );
2373 warnin(chLine);
2374 }
2375 else if( opac.telec > 2.0 )
2376 {
2377 sprintf( chLine, " C-The model is moderately optically thick to electron scattering; tau=%.1f",
2378 opac.telec );
2379 caunin(chLine);
2380 }
2381 else if( opac.telec > 0.1 )
2382 {
2383 sprintf( chLine, " !The model has modest optical depth to electron scattering; tau=%.2f",
2384 opac.telec );
2385 bangin(chLine);
2386 }
2387 else if( opac.telec > 0.01 )
2388 {
2389 sprintf( chLine, " The optical depth to electron scattering is %.3f",
2390 opac.telec );
2391 notein(chLine);
2392 }
2393
2394 /* optical depth to 21 cm */
2395 if( HFLines[0].Emis().TauIn() > 0.5 )
2396 {
2397 sprintf( chLine, " !The optical depth in the H I 21 cm line is %.2e",HFLines[0].Emis().TauIn() );
2398 bangin(chLine);
2399 }
2400
2401 /* comment if level2 lines are off - they are used to pump excited states
2402 * of ground term by UV light */
2403 if( nWindLine==0 )
2404 {
2405 /* generate comment */
2406 sprintf( chLine, " !The level2 lines are disabled. UV pumping of excited levels within ground terms is not treated." );
2407 bangin(chLine);
2408 }
2409
2410 /* check on optical depth convergence of all hydrogenic lines */
2411 for( nelem=0; nelem < LIMELM; nelem++ )
2412 {
2413 if( dense.lgElmtOn[nelem] && !dynamics.lgTimeDependentStatic )
2414 {
2415 if( iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauIn() > 0.2 )
2416 {
2417 differ = fabs(1.-iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauIn()*
2418 rt.DoubleTau/iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauTot())*100.;
2419
2420 /* check whether H-alpha optical depth changed by much on last iteration
2421 * no tolerance can be finer than autocv, the tolerance on the
2422 * iterate to convergence command. It is 15% */
2423 if( ((iterations.lgLastIt && iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauIn() > 0.8) &&
2424 differ > 20.) && wind.lgStatic() )
2425 {
2426 sprintf( chLine,
2427 " C-This is the last iteration and %2s%2ld Bal(a) optical depth"
2428 " changed by%6.1f%% (was %.2e). Try another iteration.",
2429 elementnames.chElementSym[nelem],
2430 nelem+1, differ,
2431 iso_sp[ipH_LIKE][nelem].trans(ipH3p,ipH2s).Emis().TauTot() );
2432 caunin(chLine);
2433 iterations.lgIterAgain = true;
2434 }
2435
2436 /* only check on Lya convergence if Balmer lines are thick */
2437 if( iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn() > 0. )
2438 {
2439 differ = fabs(1.-iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn()*
2440 rt.DoubleTau/iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauTot())*100.;
2441
2442 /* check whether Lya optical depth changed on last iteration
2443 * no tolerance can be finer than autocv, the tolerance on the
2444 * iterate to convergence command. It is 15% */
2445 if( ((iterations.lgLastIt && iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn() > 0.8) &&
2446 differ > 25.) && wind.lgStatic() )
2447 {
2448 sprintf( chLine,
2449 " C-This is the last iteration and %2s%2ld Ly(a) optical depth"
2450 " changed by%7.0f%% (was %.2e). Try another iteration.",
2451 elementnames.chElementSym[nelem],
2452 nelem+1,differ, iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauTot() );
2453 caunin(chLine);
2454 iterations.lgIterAgain = true;
2455 }
2456 }
2457 }
2458 }
2459 }
2460
2461 /* check whether sphere was set if dr/r large */
2462 if( radius.Radius/radius.rinner > 2. && !geometry.lgSphere )
2463 {
2464 sprintf( chLine, " C-R(out)/R(in)=%.2e and SPHERE was not set.",
2465 radius.Radius/radius.rinner );
2466 caunin(chLine);
2467 }
2468
2469 /* check if thin in hydrogen or helium continua, but assumed to be thick */
2470 if( iterations.lgLastIt && !opac.lgCaseB )
2471 {
2472
2473 /* check if thin in Lyman continuum, and assumed thick */
2474 if( rfield.nflux > iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon )
2475 {
2476 if( opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] < 2. &&
2477 opac.TauAbsGeo[1][iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] > 2. )
2478 {
2479 sprintf( chLine, " C-The H Lyman continuum is thin, and I assumed"
2480 " that it was thick. Try another iteration." );
2481 caunin(chLine);
2482 iterations.lgIterAgain = true;
2483 }
2484 }
2485
2486 /* check on the He+ ionizing continuum */
2487 if( rfield.nflux > iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon && dense.lgElmtOn[ipHELIUM] )
2488 {
2489 if( (opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1] < 2. &&
2490 opac.TauAbsGeo[1][iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1] > 2.) )
2491 {
2492 sprintf( chLine,
2493 " C-The He II continuum is thin and I assumed that it was thick."
2494 " Try another iteration." );
2495 caunin(chLine);
2496 iterations.lgIterAgain = true;
2497 }
2498 }
2499
2500 if( rfield.nflux > iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon && dense.lgElmtOn[ipHELIUM] )
2501 {
2502 if( (opac.TauAbsGeo[0][iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1] < 2. &&
2503 opac.TauAbsGeo[1][iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1] > 2.) )
2504 {
2505 sprintf( chLine,
2506 " C-The He I continuum is thin and I assumed that it was thick."
2507 " Try another iteration." );
2508 caunin(chLine);
2509 iterations.lgIterAgain = true;
2510 }
2511 }
2512 }
2513
2514 /* check whether column density changed by much on this iteration */
2515 if( iteration > 1 )
2516 {
2517 if( colden.colden_old[ipCOL_HTOT] <= 0. )
2518 {
2519 fprintf( ioQQQ, " colden_old is insane in PrtComment.\n" );
2520 ShowMe();
2522 }
2523
2524 differ = fabs(1.-colden.colden[ipCOL_HTOT]/
2525 colden.colden_old[ipCOL_HTOT]);
2526
2527 if( differ > 0.1 && differ <= 0.3 )
2528 {
2529 sprintf( chLine,
2530 " The H column density changed by %.2e%% between this and previous iteration.",
2531 differ*100. );
2532 notein(chLine);
2533 }
2534
2535 else if( differ > 0.3 )
2536 {
2537 if( iterations.lgLastIt )
2538 {
2539 sprintf( chLine,
2540 " C-The H column density changed by %.2e%% and this is the last iteration. What happened?",
2541 differ*100. );
2542 caunin(chLine);
2543 }
2544 else
2545 {
2546 sprintf( chLine,
2547 " !The H column density changed by %.2e%% What happened?",
2548 differ*100. );
2549 bangin(chLine);
2550 }
2551 }
2552
2553 /* check on H2 column density, but only if significant fraction of H is molecular */
2554 if( (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])/SDIV(colden.colden[ipCOL_HTOT]) > 1e-5 )
2555 {
2556 differ = fabs(1.-colden.colden[ipCOL_H2g]/
2557 SDIV(colden.colden_old[ipCOL_H2g]));
2558
2559 if( differ > 0.1 && differ <= 0.3 )
2560 {
2561 sprintf( chLine,
2562 " The H2 column density changed by %.2e%% between this and previous iteration.",
2563 differ*100. );
2564 notein(chLine);
2565 }
2566
2567 else if( differ > 0.3 )
2568 {
2569 if( iterations.lgLastIt )
2570 {
2571 sprintf( chLine,
2572 " C-The H2 column density changed by %.2e%% and this is the last iteration. What happened?",
2573 differ*100. );
2574 caunin(chLine);
2575 }
2576 else
2577 {
2578 sprintf( chLine,
2579 " !The H2 column density changed by %.2e%% What happened?",
2580 differ*100. );
2581 bangin(chLine);
2582 }
2583 }
2584 }
2585 }
2586
2587 /* say if rad pressure caused by la and la optical depth changed too much */
2588 differ = fabs(1.-iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn()/
2589 SDIV(iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot()))*100.;
2590
2591 if( iterations.lgLastIt && (pressure.RadBetaMax > 0.1) &&
2592 (differ > 50.) && (pressure.ipPradMax_line == 1) && (pressure.lgPres_radiation_ON) &&
2593 wind.lgStatic() )
2594 {
2595 sprintf( chLine, " C-This is the last iteration, radiation pressure was significant, and the L-a optical depth changed by %7.2f%% (was %.2e)",
2596 differ, iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot() );
2597 caunin(chLine);
2598 }
2599
2600 /* caution that 21 cm spin temperature is incorrect when Lya optical depth
2601 * scale is overrun */
2602 if( iterations.lgLastIt &&
2603 ( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot() * 1.02 -
2604 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn() ) < 0. )
2605 {
2606 sprintf( chLine, " C-The Lya optical depth scale was overrun and this is the last iteration - Tspin(21 cm) is not valid." );
2607 caunin(chLine);
2608 sprintf( chLine, " C-Another iteration is needed for Tspin(21 cm) to be valid." );
2609 caunin(chLine);
2610 }
2611
2612 /* say if la rad pressure capped by thermalization length */
2613 if( pressure.lgPradDen )
2614 {
2615 sprintf( chLine, " Line radiation pressure capped by thermalization length." );
2616 notein(chLine);
2617 }
2618
2619 /* print te failures */
2620 nline = MIN2(conv.nTeFail,10);
2621 if( conv.nTeFail != 0 )
2622 {
2623 long int _o;
2624 if( conv.failmx < 0.1 )
2625 {
2626 _o = sprintf( chLine, " There were %ld minor temperature failures. zones:",
2627 conv.nTeFail );
2628 /* don't know how many zones we will save, there are nline,
2629 * hence this use of pointer arith */
2630 for( i=0; i < nline; i++ )
2631 {
2632 _o += sprintf( chLine+_o, " %ld", conv.ifailz[i] );
2633 }
2634 notein(chLine);
2635 }
2636 else
2637 {
2638 _o = sprintf( chLine,
2639 " !There were %ld temperature failures, and some were large. The largest was %.1f%%. What happened?",
2640 conv.nTeFail, conv.failmx*100. );
2641 bangin(chLine);
2642
2643 /* don't know how many zones we will save, there are nline,
2644 * hence this use of pointer arith */
2645 _o = sprintf( chLine , " !The zones were" );
2646 for( i=0; i < nline; i++ )
2647 {
2648 _o += sprintf( chLine+_o, " %ld", conv.ifailz[i] );
2649 }
2650 bangin(chLine);
2651
2652 if( struc.testr[0] > 8e4 && phycon.te < 6e5 )
2653 {
2654 sprintf( chLine, " !I think they may have been caused by the change from hot to nebular gas phase. The physics of this is unclear." );
2655 bangin(chLine);
2656 }
2657 }
2658 }
2659
2660 /* check for temperature jumps */
2661 big_ion_jump = 0.;
2662 j = 0;
2663 for( i=1; i < nzone; i++ )
2664 {
2665 big = fabs(1.-struc.testr[i-1]/struc.testr[i]);
2666 if( big > big_ion_jump )
2667 {
2668 j = i;
2669 big_ion_jump = big;
2670 }
2671 }
2672
2673 if( big_ion_jump > 0.2 )
2674 {
2675 /* this is a sanity check, but only do it if jump detected */
2676 if( j < 1 )
2677 {
2678 fprintf( ioQQQ, " j too small big jump check\n" );
2679 ShowMe();
2681 }
2682
2683 if( big_ion_jump > 0.4 )
2684 {
2685 sprintf( chLine, " C-A temperature discontinuity occurred at zone %ld from %.2eK to %.2eK.",
2686 j, struc.testr[j-1], struc.testr[j] );
2687 caunin(chLine);
2688 /* check if the second temperature is between 100 and 1000K */
2689 /* >>chng 05 nov 07, test second not first temperature since second
2690 * will be lower of the two */
2691 /*if( struc.testr[j-1] < 1000. && struc.testr[j-1]>100. )*/
2692 if( struc.testr[j]>100. && struc.testr[j] < 1000. )
2693 {
2694 sprintf( chLine, " C-This was probably due to a thermal front." );
2695 caunin(chLine);
2696 }
2697 }
2698 else if( big_ion_jump > 0.2 )
2699 {
2700 sprintf( chLine, " !A temperature discontinuity occurred at zone %ld from %.2eK to %.2eK.",
2701 j, struc.testr[j-1], struc.testr[j] );
2702 bangin(chLine);
2703 /* check if the second temperature is between 100 and 1000K */
2704 /* >>chng 05 nov 07, test second not first temperature since second
2705 * will be lower of the two */
2706 /*if( struc.testr[j-1] < 1000. && struc.testr[j-1]>100. )*/
2707 if( struc.testr[j]>100. && struc.testr[j] < 1000. )
2708 {
2709 sprintf( chLine, " !This was probably due to a thermal front." );
2710 bangin(chLine);
2711 }
2712 }
2713 }
2714
2715 /* check for largest error in local electron density */
2716 if( fabs(conv.BigEdenError) > conv.EdenErrorAllowed )
2717 {
2718 /* this only produces a warning if not the very last zone */
2719 if( fabs(conv.BigEdenError) > conv.EdenErrorAllowed*20. && dense.nzEdenBad !=
2720 nzone )
2721 {
2722 sprintf( chLine, " W-The local error in the electron density reached %.1f%% at zone %ld",
2723 conv.BigEdenError*100, dense.nzEdenBad );
2724 warnin(chLine);
2725 }
2726 else if( fabs(conv.BigEdenError) > conv.EdenErrorAllowed*5. )
2727 {
2728 sprintf( chLine, " C-The local error in the electron density reached %.1f%% at zone %ld",
2729 conv.BigEdenError*100, dense.nzEdenBad );
2730 caunin(chLine);
2731 }
2732 else
2733 {
2734 sprintf( chLine, " The local error in the electron density reached %.1f%% at zone %ld",
2735 conv.BigEdenError*100, dense.nzEdenBad );
2736 notein(chLine);
2737 }
2738 }
2739
2740 /* check for temperature oscillations or fluctuations*/
2741 big_ion_jump = 0.;
2742 j = 0;
2743 for( i=1; i < (nzone - 1); i++ )
2744 {
2745 big = fabs( (struc.testr[i-1] - struc.testr[i])/struc.testr[i] );
2746 bigm = fabs( (struc.testr[i] - struc.testr[i+1])/struc.testr[i] );
2747
2748 /* this is sign of change in temperature, we are looking for change in sign */
2749 rel = ( (struc.testr[i-1] - struc.testr[i])/struc.testr[i])*
2750 ( (struc.testr[i] - struc.testr[i+1])/struc.testr[i] );
2751
2752 if( rel < 0. && MIN2( bigm , big ) > big_ion_jump )
2753 {
2754 j = i;
2755 big_ion_jump = MIN2( bigm , big );
2756 }
2757 }
2758
2759 if( big_ion_jump > 0.1 )
2760 {
2761 /* only do sanity check if jump detected */
2762 if( j < 1 )
2763 {
2764 fprintf( ioQQQ, " j too small bigjump2 check\n" );
2765 ShowMe();
2767 }
2768
2769 if( big_ion_jump > 0.3 )
2770 {
2771 sprintf( chLine,
2772 " C-A temperature oscillation occurred at zone%4ld by%5.0f%% from %.2e to %.2e to %.2e",
2773 j, big_ion_jump*100., struc.testr[j-1], struc.testr[j], struc.testr[j+1] );
2774 caunin(chLine);
2775 }
2776 else if( big_ion_jump > 0.1 )
2777 {
2778 sprintf( chLine,
2779 " !A temperature oscillation occurred at zone%4ld by%5.0f%% from %.2e to %.2e to %.2e",
2780 j, big_ion_jump*100., struc.testr[j-1], struc.testr[j], struc.testr[j+1] );
2781 bangin(chLine);
2782 }
2783 }
2784
2785 /* check for eden oscillations */
2786 if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
2787 {
2788 j = 0;
2789 big_ion_jump = 0.;
2790 for( i=1; i < (nzone - 1); i++ )
2791 {
2792 big = (struc.ednstr[i-1] - struc.ednstr[i])/struc.ednstr[i];
2793 if( fabs(big) < conv.EdenErrorAllowed )
2794 big = 0.;
2795 bigm = (struc.ednstr[i] - struc.ednstr[i+1])/struc.ednstr[i];
2796 if( fabs(bigm) < conv.EdenErrorAllowed )
2797 bigm = 0.;
2798 if( big*bigm < 0. &&
2799 fabs(struc.ednstr[i-1]-struc.ednstr[i])/struc.ednstr[i] > big_ion_jump )
2800 {
2801 j = i;
2802 big_ion_jump = fabs(struc.ednstr[i-1]-struc.ednstr[i])/
2803 struc.ednstr[i];
2804 }
2805 }
2806
2807 /* only check on j if there was a big jump detected, number must be
2808 * smallest jump */
2809 if( big_ion_jump > conv.EdenErrorAllowed*3. )
2810 {
2811 if( j < 1 )
2812 {
2813 fprintf( ioQQQ, " j too small bigjump3 check\n" );
2814 ShowMe();
2816 }
2817
2818 if( big_ion_jump > conv.EdenErrorAllowed*10. )
2819 {
2820 sprintf( chLine, " C-An electron density oscillation occurred at zone%4ld by%5.0f%% from %.2e to %.2e to %.2e",
2821 j, big_ion_jump*100., struc.ednstr[j-1], struc.ednstr[j],
2822 struc.ednstr[j+1] );
2823 caunin(chLine);
2824 }
2825 else if( big_ion_jump > conv.EdenErrorAllowed*3. )
2826 {
2827 sprintf( chLine, " !An electron density oscillation occurred at zone%4ld by%5.0f%% from %.2e to %.2e to %.2e",
2828 j, big_ion_jump*100., struc.ednstr[j-1], struc.ednstr[j],
2829 struc.ednstr[j+1] );
2830 bangin(chLine);
2831 }
2832 }
2833 }
2834
2835 /*prt_smooth_predictions check whether fluctuations in any predicted quantities occurred */
2836 /* >>chng 03 dec 05, add this test */
2838
2839 /**********************************************************
2840 * check that the volume integrates out ok *
2841 **********************************************************/
2842
2843 /* this was the number 1 fed through the line integrators,
2844 * the number 1e-10 is sent to linadd in lineset1 as follows:*/
2845 /*linadd( 1.e-10 , 1 , "Unit" , 'i' );*/
2846 i = cdLine( "Unit" , 1 , &rate , &absint );
2847 ASSERT( i> 0 );
2848
2849 /* this is now the linear vol, rel to inner radius */
2850 VolComputed = LineSv[i].SumLine[0] / 1e-10;
2851
2852 /* spherical or plane parallel case? */
2853 if( radius.Radius/radius.rinner > 1.0001 )
2854 {
2855 /* spherical case,
2856 * geometry.iEmissPower is usually 2,
2857 * and can be reset to 1 (long slit) or 0 (beam) with
2858 * slit and beam options on aperture */
2859 VolExpected = geometry.covaper*geometry.FillFac*radius.rinner/(geometry.iEmissPower+1)*
2860 ( powi( radius.Radius/radius.rinner,geometry.iEmissPower+1 ) - 1. );
2861 }
2862 else
2863 {
2864 /* plane parallel case */
2865 /* next can be zero for very thin model, depth is always positive */
2866 VolExpected = geometry.covaper*geometry.FillFac*(radius.depth-DEPTH_OFFSET);
2867 }
2868
2869 /* now get the relative difference between computed and expected volumes */
2870 error = fabs(VolComputed - VolExpected)/SDIV(VolExpected);
2871
2872 /* we need to ignore this test if filling factor changes with radius, or
2873 * cylinder geometry in place */
2874 if( radius.lgCylnOn || geometry.filpow!=0. )
2875 {
2876 error = 0.;
2877 }
2878
2879 /* how large is relative error? */
2880 if( error > 0.001 && !lgAbort )
2881 {
2882 sprintf( chLine,
2883 " W-PrtComment insanity - Line unit integration did not verify \n");
2884 warnin(chLine);
2885 fprintf( ioQQQ,
2886 " PROBLEM PrtComment insanity - Line unit integration did not verify \n");
2887 fprintf( ioQQQ,
2888 " expected, derived vols were %g %g \n",
2889 VolExpected , VolComputed );
2890 fprintf( ioQQQ,
2891 " relative difference is %g, ratio is %g.\n",error,VolComputed/VolExpected);
2892 TotalInsanity();
2893 }
2894
2895 /* next do same thing for fake continuum point propagated in highest energy cell, plus 1
2896 * =
2897 * the variable rfield.ConEmitLocal[rfield.nflux]
2898 * are set to
2899 * the number 1.e-10f * Dilution in RT_diffuse. this is the outward
2900 * local emissivity, per unit vol. It is then added to the outward beams
2901 * by the rest of the code, and then checked here.
2902 *
2903 * insanity will be detected if diffuse emission is thrown into the outward beam
2904 * in MadeDiffuse. this happens if the range of ionization encompasses the full
2905 * continuum array, up to nflux. */
2906 ConComputed = rfield.ConInterOut[rfield.nflux]/ 1e-10;
2907 /* correct for fraction that went out, as set in ZoneStart,
2908 * this should now be the volume of the emitting region */
2909 ConComputed /= ( (1. + geometry.covrt)/2. );
2910
2911 /* we expect this to add up to the integral of unity over r^-2 */
2912 if( radius.Radius/radius.rinner < 1.001 )
2913 {
2914 /* plane parallel case, no dilution, use thickness */
2915 ConExpected = (radius.depth-DEPTH_OFFSET)*geometry.FillFac;
2916 }
2917 else
2918 {
2919 /* spherical case */
2920 ConExpected = radius.rinner*geometry.FillFac * (1. - radius.rinner/radius.Radius );
2921 }
2922 /* this is impossible */
2923 ASSERT( ConExpected > 0. );
2924
2925 /* now get the relative difference between computed and expected volumes */
2926 error = fabs(ConComputed - ConExpected)/ConExpected;
2927
2928 /* we need to ignore this test if filling factor changes with radius, or
2929 * cylinder geometry in place */
2930 if( radius.lgCylnOn || geometry.filpow!=0. )
2931 {
2932 error = 0.;
2933 }
2934
2935 /* \todo 2 - These "volumes" seem to be too small by a factor of two.
2936 * rfield.ConInterOut[rfield.nflux] (hence ConComputed) and ConExpected
2937 * should be greater by a factor of 2 if comparison is really of "volume"
2938 * of 1/cc pencil. */
2939
2940 /* how large is relative error? */
2941 if( error > 0.001 && !lgAbort)
2942 {
2943 sprintf( chLine,
2944 " W-PrtComment insanity - Continuum unit integration did not verify \n");
2945 warnin(chLine);
2946 fprintf( ioQQQ," PROBLEM PrtComment insanity - Continuum unit integration did not verify \n");
2947 fprintf( ioQQQ," exact vol= %g, derived vol= %g relative difference is %g \n",
2948 ConExpected , ConComputed ,error);
2949 fprintf( ioQQQ," ConInterOut= %g, \n",
2950 rfield.ConInterOut[rfield.nflux]);
2951 TotalInsanity();
2952 }
2953
2954 /* final printout of warnings, cautions, notes, */
2955 cdNwcns(&lgAbort_flag,&nw,&nc,&nn,&ns,&i,&j,&dum1,&dum2);
2956
2957 warnings.lgWarngs = ( nw > 0 );
2958 warnings.lgCautns = ( nc > 0 );
2959
2960 if( called.lgTalk )
2961 {
2962 /* print the title of the calculation */
2963 fprintf( ioQQQ, " %s\n", input.chTitle );
2964 /* say why the calculation stopped, and indicate the geometry*/
2966 /* print all warnings */
2968 /* all cautions */
2970 /* surprises, beginning with a ! */
2972 /* notes about the calculations */
2973 cdNotes(ioQQQ);
2974 }
2975
2976 /* option to print warnings on special io */
2977 if( lgPrnErr )
2978 {
2980 }
2981
2982 if( trace.lgTrace )
2983 {
2984 fprintf( ioQQQ, " PrtComment returns.\n" );
2985 }
2986 return;
2987}
2988
2989
2990/*chkCaHeps check whether CaII K and H epsilon overlap */
2991STATIC void chkCaHeps(double *totwid)
2992{
2993 double conca,
2994 conalog ,
2995 conhe;
2996
2997 DEBUG_ENTRY( "chkCaHeps()" );
2998
2999 *totwid = 0.;
3000 /* pumping of CaH overlapping with Hy epsilon, 6-2 of H */
3001 if( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_local +
3002 iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local >= 6 )
3003 {
3004 /* this is 6P */
3005 long ip6p = iso_sp[ipH_LIKE][ipHYDROGEN].QuantumNumbers2Index[6][1][2];
3006
3007 if( TauLines[ipT3969].Emis().TauIn() > 0. &&
3008 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ip6p,ipH2s).Emis().TauIn() > 0. )
3009 {
3010 /* casts to double here are to prevent FPE */
3011 conca = pow(6.1e-5* TauLines[ipT3969].Emis().TauIn(),0.5);
3012 conalog = log((double)TauLines[ipT3969].Emis().TauIn());
3013 conalog = sqrt(MAX2(1., conalog));
3014 conca = MAX2(conalog,conca);
3015
3016 conalog = log((double)iso_sp[ipH_LIKE][ipHYDROGEN].trans(ip6p,ipH2s).Emis().TauIn());
3017 conalog = sqrt(MAX2(1.,conalog));
3018 conhe = pow(1.7e-6*iso_sp[ipH_LIKE][ipHYDROGEN].trans(ip6p,ipH2s).Emis().TauIn(),0.5);
3019 conhe = MAX2(conalog, conhe);
3020
3021 *totwid = 10.*conhe + 1.6*conca;
3022 }
3023 }
3024 return;
3025}
3026
3027/*prt_smooth_predictions check whether fluctuations in any predicted quantities occurred */
3029{
3030 long int i,
3031 nzone_oscillation,
3032 nzone_ion_jump,
3033 nzone_den_jump,
3034 nelem,
3035 ion;
3036 double BigOscillation ,
3037 big_ion_jump,
3038 big_jump,
3039 rel,
3040 big,
3041 bigm;
3042
3043 char chLine[INPUT_LINE_LENGTH];
3044
3045 DEBUG_ENTRY( "prt_smooth_predictions()" );
3046
3047 /* check for ionization oscillations or fluctuations and or jumps */
3048 nzone_oscillation = 0;
3049 nzone_ion_jump = 0;
3050
3051 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
3052 {
3053 if( dense.lgElmtOn[nelem] )
3054 {
3055 for( ion=0; ion<=nelem+1; ++ion)
3056 {
3057 BigOscillation = 0.;
3058 big_ion_jump = -15.;
3059 /* >>chng 05 mar 12, add -lgAbort, since in some bad aborts current zone never evaluated */
3060 for( i=1; i < (nzone - 1-(int)lgAbort); i++ )
3061 {
3062
3063 /* only do check if all ions are positive */
3064 if( struc.xIonDense[nelem][ion][i-1]/struc.gas_phase[nelem][i-1]>struc.dr_ionfrac_limit &&
3065 struc.xIonDense[nelem][ion][i ]/struc.gas_phase[nelem][i ]>struc.dr_ionfrac_limit &&
3066 struc.xIonDense[nelem][ion][i+1]/struc.gas_phase[nelem][i+1]>struc.dr_ionfrac_limit )
3067 {
3068
3069 /* this is check for oscillations */
3070 big = fabs( (struc.xIonDense[nelem][ion][i-1] - struc.xIonDense[nelem][ion][i])/struc.xIonDense[nelem][ion][i] );
3071 bigm = fabs( (struc.xIonDense[nelem][ion][i] - struc.xIonDense[nelem][ion][i+1])/struc.xIonDense[nelem][ion][i] );
3072
3073 /* this is sign of change in ionization, we are looking for change in sign */
3074 rel = ( (struc.xIonDense[nelem][ion][i-1] - struc.xIonDense[nelem][ion][i] )/struc.xIonDense[nelem][ion][i])*
3075 ( (struc.xIonDense[nelem][ion][i] - struc.xIonDense[nelem][ion][i+1])/struc.xIonDense[nelem][ion][i] );
3076
3077 if( rel < 0. && MIN2( bigm , big ) > BigOscillation )
3078 {
3079 nzone_oscillation = i;
3080 BigOscillation = MIN2( bigm , big );
3081 }
3082
3083 /* check whether we tripped over an ionization front - a major source
3084 * of instability in a complete linearization code like this one */
3085 /* neg sign picks up only increases in ionization */
3086 rel = -log10( (struc.xIonDense[nelem][ion][i]/struc.gas_phase[nelem][i]) /
3087 (struc.xIonDense[nelem][ion][i+1]/struc.gas_phase[nelem][i+1] ) );
3088 /* only do significant stages of ionization */
3089 if( rel > big_ion_jump )
3090 {
3091 big_ion_jump = rel;
3092 nzone_ion_jump = i;
3093 }
3094 }
3095 }
3096 /* end loop over zones,
3097 * check whether this ion and element underwent fluctuations or jump */
3098
3099 if( BigOscillation > 0.2 )
3100 {
3101 /* only do sanity check if jump detected */
3102 if( nzone_oscillation < 1 )
3103 {
3104 fprintf( ioQQQ, " nzone_oscillation too small bigjump2 check\n" );
3105 ShowMe();
3107 }
3108 if( BigOscillation > 3. )
3109 {
3110 sprintf( chLine,
3111 " W-An ionization oscillation occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
3112 nzone_oscillation,
3113 elementnames.chElementSym[nelem], ion+1,
3114 BigOscillation*100.,
3115 struc.xIonDense[nelem][ion][nzone_oscillation-1]/struc.gas_phase[nelem][nzone_oscillation-1],
3116 struc.xIonDense[nelem][ion][nzone_oscillation]/struc.gas_phase[nelem][nzone_oscillation],
3117 struc.xIonDense[nelem][ion][nzone_oscillation+1]/struc.gas_phase[nelem][nzone_oscillation+1] );
3118 warnin(chLine);
3119 }
3120
3121 else if( BigOscillation > 0.7 )
3122 {
3123 sprintf( chLine,
3124 " C-An ionization oscillation occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
3125 nzone_oscillation,
3126 elementnames.chElementSym[nelem], ion+1,
3127 BigOscillation*100.,
3128 struc.xIonDense[nelem][ion][nzone_oscillation-1]/struc.gas_phase[nelem][nzone_oscillation-1],
3129 struc.xIonDense[nelem][ion][nzone_oscillation]/struc.gas_phase[nelem][nzone_oscillation],
3130 struc.xIonDense[nelem][ion][nzone_oscillation+1]/struc.gas_phase[nelem][nzone_oscillation+1] );
3131 caunin(chLine);
3132 }
3133 else if( BigOscillation > 0.2 )
3134 {
3135 sprintf( chLine,
3136 " !An ionization oscillation occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
3137 nzone_oscillation,
3138 elementnames.chElementSym[nelem], ion+1,
3139 BigOscillation*100.,
3140 struc.xIonDense[nelem][ion][nzone_oscillation-1]/struc.gas_phase[nelem][nzone_oscillation-1],
3141 struc.xIonDense[nelem][ion][nzone_oscillation]/struc.gas_phase[nelem][nzone_oscillation],
3142 struc.xIonDense[nelem][ion][nzone_oscillation+1]/struc.gas_phase[nelem][nzone_oscillation+1] );
3143 bangin(chLine);
3144 }
3145 }
3146
3147 /* big_ion_jump was a log above, convert to linear quantity */
3148 /* if no jump occurred then big_ion_jump is small and nzone_ion_jump is 0 */
3149 big_ion_jump = pow(10., big_ion_jump );
3150 if( big_ion_jump > 1.5 && nzone_ion_jump > 0 )
3151 {
3152 if( big_ion_jump > 10. )
3153 {
3154 sprintf( chLine,
3155 " C-An ionization jump occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
3156 nzone_ion_jump,
3157 elementnames.chElementSym[nelem], ion+1,
3158 big_ion_jump*100.,
3159 struc.xIonDense[nelem][ion][nzone_ion_jump-1]/struc.gas_phase[nelem][nzone_ion_jump-1],
3160 struc.xIonDense[nelem][ion][nzone_ion_jump]/struc.gas_phase[nelem][nzone_ion_jump],
3161 struc.xIonDense[nelem][ion][nzone_ion_jump+1]/struc.gas_phase[nelem][nzone_ion_jump+1] );
3162 caunin(chLine);
3163 }
3164 else
3165 {
3166 sprintf( chLine,
3167 " !An ionization jump occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e",
3168 nzone_ion_jump,
3169 elementnames.chElementSym[nelem], ion+1,
3170 big_ion_jump*100.,
3171 struc.xIonDense[nelem][ion][nzone_ion_jump-1]/struc.gas_phase[nelem][nzone_ion_jump-1],
3172 struc.xIonDense[nelem][ion][nzone_ion_jump]/struc.gas_phase[nelem][nzone_ion_jump],
3173 struc.xIonDense[nelem][ion][nzone_ion_jump+1]/struc.gas_phase[nelem][nzone_ion_jump+1] );
3174 bangin(chLine);
3175 }
3176 }
3177 }
3178 }
3179 }
3180
3181 big_jump = -15;
3182 nzone_den_jump = 0;
3183
3184 /* >>chng 05 mar 12, add -lgAbort, since in some bad aborts current zone never evaluated */
3185 for( i=1; i < (nzone - 1 - (int)lgAbort); i++ )
3186 {
3187 /* this first check is on how the total hydrogen density has changed */
3188 rel = fabs(log10( struc.gas_phase[ipHYDROGEN][i] /
3189 struc.gas_phase[ipHYDROGEN][i+1] ) );
3190 /* only do significant stages of ionization */
3191 if( rel > big_jump )
3192 {
3193 big_jump = rel;
3194 nzone_den_jump = i;
3195 }
3196 }
3197
3198 /* check how stable density was */
3199 big_jump = pow( 10., big_jump );
3200 if( big_jump > 1.2 )
3201 {
3202 if( big_jump > 3. )
3203 {
3204 sprintf( chLine,
3205 " C-The H density jumped at by %.0f%% at zone %ld, from %.2e to %.2e to %.2e",
3206 big_jump*100.,
3207 nzone_den_jump,
3208 struc.gas_phase[ipHYDROGEN][nzone_den_jump-1],
3209 struc.gas_phase[ipHYDROGEN][nzone_den_jump],
3210 struc.gas_phase[ipHYDROGEN][nzone_den_jump+1] );
3211 caunin(chLine);
3212 }
3213 else
3214 {
3215 sprintf( chLine,
3216 " !An H density jump occurred at zone %ld, by %.0f%% from %.2e to %.2e to %.2e",
3217 nzone_den_jump,
3218 big_jump*100.,
3219 struc.gas_phase[ipHYDROGEN][nzone_den_jump-1],
3220 struc.gas_phase[ipHYDROGEN][nzone_den_jump],
3221 struc.gas_phase[ipHYDROGEN][nzone_den_jump+1] );
3222 bangin(chLine);
3223 }
3224 }
3225
3226 /* now do check on smoothness of radiation pressure */
3227 big_jump = -15;
3228 nzone_den_jump = 0;
3229
3230 /* loop starts on zone 3 since dramatic fall in radiation pressure across first
3231 * few zones is normal behavior */
3232 /* >>chng 05 mar 12, add -lgAbort, since in some bad aborts current zone never evaluated */
3233 for( i=3; i < (nzone - 2 - (int)lgAbort); i++ )
3234 {
3235 /* this first check is on how the total hydrogen density has changed */
3236 rel = fabs(log10( SDIV(struc.pres_radiation_lines_curr[i]) /
3237 SDIV(0.5*(struc.pres_radiation_lines_curr[i-1]+struc.pres_radiation_lines_curr[i+1])) ) );
3238 /* only do significant stages of ionization */
3239 if( rel > big_jump )
3240 {
3241 big_jump = rel;
3242 nzone_den_jump = i;
3243 }
3244 }
3245 /* note that changing log big_jump to linear takes place in next branch */
3246
3247 /* check how stable radiation pressure was, but only if significant */
3248 if( pressure.RadBetaMax > 0.01 )
3249 {
3250 big_jump = pow( 10., big_jump );
3251 if( big_jump > 1.2 )
3252 {
3253 /* only make it a caution is pressure jumped, and we were trying
3254 * to do a constant pressure model */
3255 if( big_jump > 3. && strcmp(dense.chDenseLaw,"CPRE") == 0)
3256 {
3257 sprintf( chLine,
3258 " C-The radiation pressure jumped by %.0f%% at zone %ld, from %.2e to %.2e to %.2e",
3259 big_jump*100.,
3260 nzone_den_jump,
3261 struc.pres_radiation_lines_curr[nzone_den_jump-1],
3262 struc.pres_radiation_lines_curr[nzone_den_jump],
3263 struc.pres_radiation_lines_curr[nzone_den_jump+1] );
3264 caunin(chLine);
3265 }
3266 else
3267 {
3268 sprintf( chLine,
3269 " !The radiation pressure jumped by %.0f%% at zone %ld, from %.2e to %.2e to %.2e",
3270 big_jump*100.,
3271 nzone_den_jump,
3272 struc.pres_radiation_lines_curr[nzone_den_jump-1],
3273 struc.pres_radiation_lines_curr[nzone_den_jump],
3274 struc.pres_radiation_lines_curr[nzone_den_jump+1] );
3275 bangin(chLine);
3276 }
3277 }
3278 }
3279
3280 /* these will be used to check on continuity */
3281 phycon.BigJumpTe = 0.;
3282 phycon.BigJumpne = 0.;
3283 phycon.BigJumpH2 = 0.;
3284 phycon.BigJumpCO = 0.;
3285
3286 for( i=1; i < (nzone - 1 - (int)lgAbort); i++ )
3287 {
3288 /* check on how much temperature has changed */
3289 rel = fabs(log10( struc.testr[i] / struc.testr[i+1] ) );
3290 if( rel > phycon.BigJumpTe )
3291 {
3292 phycon.BigJumpTe = (realnum)rel;
3293 }
3294
3295 /* check on how much electron density has changed */
3296 rel = fabs(log10( struc.ednstr[i] / struc.ednstr[i+1] ) );
3297 if( rel > phycon.BigJumpne )
3298 {
3299 phycon.BigJumpne = (realnum)rel;
3300 }
3301
3302 /* check on how much H2 density has changed */
3303 if( (struc.H2_abund[i])>SMALLFLOAT && (struc.H2_abund[i+1]) > SMALLFLOAT
3304 /* only do this test if H2 abund is significant */
3305 && (struc.H2_abund[i])/struc.gas_phase[ipHYDROGEN][i]>1e-3)
3306 {
3307 rel = fabs(log10( (struc.H2_abund[i]) / SDIV(struc.H2_abund[i+1]) ) );
3308 if( rel > phycon.BigJumpH2 )
3309 {
3310 phycon.BigJumpH2 = (realnum)rel;
3311 }
3312 }
3313
3314 {
3315 int ipCO = findspecies("CO")->index;
3316 //fprintf(ioQQQ,"PRTCO %d %ld %d\n",ipCO,i,mole_global.num_calc);
3317 /* check on how much CO density has changed */
3318 if( ipCO != -1 &&
3319 struc.molecules[ipCO][i]>SMALLFLOAT &&
3320 struc.molecules[ipCO][i+1]>SMALLFLOAT &&
3321 struc.molecules[ipCO][i]/SDIV(struc.gas_phase[ipCARBON][i])>1e-3 )
3322 {
3323 rel = fabs(log10( struc.molecules[ipCO][i] / struc.molecules[ipCO][i+1] ) );
3324 if( rel > phycon.BigJumpCO )
3325 {
3326 phycon.BigJumpCO = (realnum)rel;
3327 }
3328 }
3329 }
3330 }
3331
3332 /* convert to linear change - subtract 1 to make it the residual difference */
3333 if(phycon.BigJumpTe>0. )
3334 phycon.BigJumpTe = (realnum)pow( 10., (double)phycon.BigJumpTe ) -1.f;
3335
3336 if(phycon.BigJumpne>0. )
3337 phycon.BigJumpne = (realnum)pow( 10., (double)phycon.BigJumpne )-1.f;
3338
3339 if(phycon.BigJumpH2>0. )
3340 phycon.BigJumpH2 = (realnum)pow( 10., (double)phycon.BigJumpH2 )-1.f;
3341
3342 if(phycon.BigJumpCO>0. )
3343 phycon.BigJumpCO = (realnum)pow( 10., (double)phycon.BigJumpCO )-1.f;
3344 /*fprintf(ioQQQ,"DEBUG continuity large change %.2e %.2e %.2e %.2e \n",
3345 phycon.BigJumpTe , phycon.BigJumpne , phycon.BigJumpH2 , phycon.BigJumpCO );*/
3346
3347 if( phycon.BigJumpTe > 0.3 )
3348 {
3349 sprintf( chLine,
3350 " C-The temperature varied by %.1f%% between two zones",
3351 phycon.BigJumpTe*100.);
3352 caunin(chLine);
3353 }
3354 else if( phycon.BigJumpTe > 0.1 )
3355 {
3356 sprintf( chLine,
3357 " !The temperature varied by %.1f%% between two zones",
3358 phycon.BigJumpTe*100.);
3359 bangin(chLine);
3360 }
3361
3362 if( phycon.BigJumpne > 0.3 )
3363 {
3364 sprintf( chLine,
3365 " C-The electron density varied by %.1f%% between two zones",
3366 phycon.BigJumpne*100.);
3367 caunin(chLine);
3368 }
3369 else if( phycon.BigJumpne > 0.1 )
3370 {
3371 sprintf( chLine,
3372 " !The electron density varied by %.1f%% between two zones",
3373 phycon.BigJumpne*100.);
3374 bangin(chLine);
3375 }
3376
3377 if( phycon.BigJumpH2 > 0.8 )
3378 {
3379 sprintf( chLine,
3380 " C-The H2 density varied by %.1f%% between two zones",
3381 phycon.BigJumpH2*100.);
3382 caunin(chLine);
3383 }
3384 else if( phycon.BigJumpH2 > 0.1 )
3385 {
3386 sprintf( chLine,
3387 " !The H2 density varied by %.1f%% between two zones",
3388 phycon.BigJumpH2*100.);
3389 bangin(chLine);
3390 }
3391
3392 if( phycon.BigJumpCO > 0.8 )
3393 {
3394 sprintf( chLine,
3395 " C-The CO density varied by %.1f%% between two zones",
3396 phycon.BigJumpCO*100.);
3397 caunin(chLine);
3398 }
3399 else if( phycon.BigJumpCO > 0.2 )
3400 {
3401 sprintf( chLine,
3402 " !The CO density varied by %.1f%% between two zones",
3403 phycon.BigJumpCO*100.);
3404 bangin(chLine);
3405 }
3406 return;
3407}
t_abund abund
Definition abund.cpp:5
void AgeCheck(void)
Definition age_check.cpp:14
t_atmdat atmdat
Definition atmdat.cpp:6
long ipT3969
t_atoms atoms
Definition atoms.cpp:5
t_broke broke
Definition broke.cpp:5
t_ca ca
Definition ca.cpp:5
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioPrnErr
Definition cddefines.cpp:9
bool lgTestCodeCalled
Definition cddefines.cpp:11
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgPrnErr
Definition cddefines.cpp:13
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
const int ipCARBON
Definition cddefines.h:310
#define MIN2
Definition cddefines.h:761
const int ipOXYGEN
Definition cddefines.h:312
const double DEPTH_OFFSET
Definition cddefines.h:272
const int LIMELM
Definition cddefines.h:258
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
double powi(double, long int)
Definition service.cpp:604
void cdCautions(FILE *ioOUT)
Definition cddrive.cpp:226
void cdSurprises(FILE *ioOUT)
Definition cddrive.cpp:282
void cdNotes(FILE *ioOUT)
Definition cddrive.cpp:309
void cdNwcns(bool *lgAbort_ret, long int *NumberWarnings, long int *NumberCautions, long int *NumberNotes, long int *NumberSurprises, long int *NumberTempFailures, long int *NumberPresFailures, long int *NumberIonFailures, long int *NumberNeFailures)
Definition cddrive.cpp:1481
void cdWarnings(FILE *ioPNT)
Definition cddrive.cpp:198
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition cddrive.cpp:1228
void cdReasonGeo(FILE *ioOUT)
Definition cddrive.cpp:175
long nWindLine
Definition cdinit.cpp:19
LinSv * LineSv
Definition cdinit.cpp:70
static t_version & Inst()
Definition cddefines.h:175
int index
Definition mole.h:169
double induc_dn_max
Definition two_photon.h:55
t_co co
Definition co.cpp:5
t_colden colden
Definition colden.cpp:5
#define ipCOL_HTOT
Definition colden.h:12
#define ipCOL_H2s
Definition colden.h:18
#define ipCOL_H2g
Definition colden.h:16
#define ipCOL_Hp
Definition colden.h:26
#define ipCOL_H0
Definition colden.h:22
#define ipCOL_H2p
Definition colden.h:20
t_continuum continuum
Definition continuum.cpp:5
t_conv conv
Definition conv.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_DoppVel DoppVel
Definition doppvel.cpp:5
t_dynamics dynamics
Definition dynamics.cpp:44
t_elementnames elementnames
bool lgConserveEnergy()
Definition energy.cpp:314
t_fudgec fudgec
Definition fudgec.cpp:5
t_geometry geometry
Definition geometry.cpp:5
GrainVar gv
Definition grainvar.cpp:5
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_hcmap hcmap
Definition hcmap.cpp:21
t_he he
Definition he.cpp:5
t_hextra hextra
Definition hextra.cpp:5
t_hmi hmi
Definition hmi.cpp:5
t_hydro hydro
Definition hydrogenic.cpp:5
t_hyperfine hyperfine
Definition hyperfine.cpp:5
t_input input
Definition input.cpp:12
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipH3p
Definition iso.h:31
const int ipH1s
Definition iso.h:27
const int ipHE_LIKE
Definition iso.h:63
const int ipH2p
Definition iso.h:29
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
t_iterations iterations
Definition iterations.cpp:5
const TransitionProxy FndLineHt(long int *level)
t_magnetic magnetic
Definition magnetic.cpp:17
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molecule * findspecies(const char buf[])
t_opac opac
Definition opacity.cpp:5
t_oxy oxy
Definition oxy.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double SOLAR_MASS
Definition physconst.h:71
UNUSED const double RYDLAM
Definition physconst.h:176
t_pressure pressure
Definition pressure.cpp:5
STATIC void prt_smooth_predictions(void)
STATIC void chkCaHeps(double *totwid)
void PrtComment(void)
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
t_save save
Definition save.cpp:5
t_secondaries secondaries
t_StopCalc StopCalc
Definition stopcalc.cpp:5
t_struc struc
Definition struc.cpp:6
TransitionList TauLine2("TauLine2", &AnonStates)
TransitionList HFLines("HFLines", &AnonStates)
long int nLevel1
Definition taulines.cpp:28
TransitionList TauLines("TauLines", &AnonStates)
t_thermal thermal
Definition thermal.cpp:5
t_timesc timesc
Definition timesc.cpp:5
t_trace trace
Definition trace.cpp:5
char * chLineLbl(const TransitionProxy &t)
t_warnings warnings
Definition warnings.cpp:11
void wcnint(void)
Definition warnings.cpp:13
void caunin(char *chLine)
Definition warnings.cpp:96
void warnin(char *chLine)
Definition warnings.cpp:27
void bangin(char *chLine)
Definition warnings.cpp:73
void notein(char *chLine)
Definition warnings.cpp:50
Wind wind
Definition wind.cpp:5