cloudy  trunk
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 */
60 STATIC void chkCaHeps(double *totwid);
61 
62 /*prt_smooth_predictions check whether fluctuations in any predicted quantities occurred */
63 STATIC void prt_smooth_predictions(void);
64 
65 void 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  {
400  TotalInsanity();
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();
890  cdEXIT(EXIT_FAILURE);
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();
897  cdEXIT(EXIT_FAILURE);
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();
1220  cdEXIT(EXIT_FAILURE);
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();
1228  cdEXIT(EXIT_FAILURE);
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();
1495  cdEXIT(EXIT_FAILURE);
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();
1503  cdEXIT(EXIT_FAILURE);
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();
1527  cdEXIT(EXIT_FAILURE);
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();
2521  cdEXIT(EXIT_FAILURE);
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();
2680  cdEXIT(EXIT_FAILURE);
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();
2766  cdEXIT(EXIT_FAILURE);
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();
2815  cdEXIT(EXIT_FAILURE);
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 */
2837  prt_smooth_predictions();
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*/
2965  cdReasonGeo(ioQQQ);
2966  /* print all warnings */
2967  cdWarnings(ioQQQ);
2968  /* all cautions */
2969  cdCautions(ioQQQ);
2970  /* surprises, beginning with a ! */
2971  cdSurprises(ioQQQ);
2972  /* notes about the calculations */
2973  cdNotes(ioQQQ);
2974  }
2975 
2976  /* option to print warnings on special io */
2977  if( lgPrnErr )
2978  {
2979  cdWarnings(ioPrnErr);
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 */
2991 STATIC 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 */
3028 STATIC void prt_smooth_predictions(void)
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();
3106  cdEXIT(EXIT_FAILURE);
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_elementnames::chElementName
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
chLineLbl
char * chLineLbl(const TransitionProxy &t)
Definition: transition.cpp:237
t_atoms::xMg2Max
realnum xMg2Max
Definition: atoms.h:260
Wind::lgWindOK
bool lgWindOK
Definition: wind.h:42
colden.h
t_conv::nTeFail
long int nTeFail
Definition: conv.h:208
thermal.h
t_radius::lgdR2Small
bool lgdR2Small
Definition: radius.h:112
t_conv::lgBadStop
bool lgBadStop
Definition: conv.h:253
t_hydro::lgHiPop2
bool lgHiPop2
Definition: hydrogenic.h:55
co
t_co co
Definition: co.cpp:5
t_struc::gas_phase
realnum ** gas_phase
Definition: struc.h:75
caunin
void caunin(char *chLine)
Definition: warnings.cpp:96
ipCOL_H2s
#define ipCOL_H2s
Definition: colden.h:18
t_pressure::lgPradCap
bool lgPradCap
Definition: pressure.h:153
t_rfield::lgHionRad
bool lgHionRad
Definition: rfield.h:469
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:312
StopCalc
t_StopCalc StopCalc
Definition: stopcalc.cpp:5
t_timesc::TimeErode
realnum TimeErode
Definition: timesc.h:48
t_thermal::CoolHeatMax
realnum CoolHeatMax
Definition: thermal.h:105
lines.h
h2
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_geometry::lgZoneSet
bool lgZoneSet
Definition: geometry.h:87
notein
void notein(char *chLine)
Definition: warnings.cpp:50
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
t_abund::lgDepln
bool lgDepln
Definition: abund.h:103
t_version::lgRelease
bool lgRelease
Definition: version.h:28
t_struc::xIonDense
realnum *** xIonDense
Definition: struc.h:64
t_dense::chDenseLaw
char chDenseLaw[5]
Definition: dense.h:158
struc.h
dense
t_dense dense
Definition: dense.cpp:24
t_opac::lgCaseB
bool lgCaseB
Definition: opacity.h:161
t_secondaries::SecHIonMax
realnum SecHIonMax
Definition: secondaries.h:30
elementnames.h
secondaries.h
t_opac::lgOpacNeg
bool lgOpacNeg
Definition: opacity.h:179
atoms.h
t_input::lgSetNoBuffering
bool lgSetNoBuffering
Definition: input.h:70
Singleton< t_version >::Inst
static t_version & Inst()
Definition: cddefines.h:175
t_conv::BigEdenError
realnum BigEdenError
Definition: conv.h:220
rfield
t_rfield rfield
Definition: rfield.cpp:8
t_atmdat::HCharCoolMax
double HCharCoolMax
Definition: atmdat.h:155
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
geometry.h
nLevel1
long int nLevel1
Definition: taulines.cpp:28
t_continuum::lgCon0
bool lgCon0
Definition: continuum.h:93
t_ionbal::CompHeating_Max
double CompHeating_Max
Definition: ionbal.h:190
fudgec.h
t_thermal::wlCoolHeatMax
realnum wlCoolHeatMax
Definition: thermal.h:106
t_co::codfrc
realnum codfrc
Definition: co.h:15
t_dynamics::CoolMax
double CoolMax
Definition: dynamics.h:68
realnum
float realnum
Definition: cddefines.h:103
iterations
t_iterations iterations
Definition: iterations.cpp:5
conv.h
rfield.h
t_phycon::BigJumpTe
realnum BigJumpTe
Definition: phycon.h:106
nWindLine
long nWindLine
Definition: cdinit.cpp:19
t_conv::EdenErrorAllowed
double EdenErrorAllowed
Definition: conv.h:267
STATIC
#define STATIC
Definition: cddefines.h:97
SOLAR_MASS
const UNUSED double SOLAR_MASS
Definition: physconst.h:71
t_atmdat::lgHCaseBOK
bool lgHCaseBOK[2][HS_NZ]
Definition: atmdat.h:193
t_rfield::occ1nu
realnum occ1nu
Definition: rfield.h:477
t_thermal::ipMaxExtra
long int ipMaxExtra
Definition: thermal.h:143
ipCARBON
const int ipCARBON
Definition: cddefines.h:310
mole.h
t_input::chTitle
char chTitle[INPUT_LINE_LENGTH]
Definition: input.h:37
molecule::index
int index
Definition: mole.h:169
t_colden::ajmmin
realnum ajmmin
Definition: colden.h:89
t_struc::testr
realnum * testr
Definition: struc.h:25
abund.h
t_thermal::thist
realnum thist
Definition: thermal.h:56
t_iterations::lgIterAgain
bool lgIterAgain
Definition: iterations.h:39
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
t_thermal::GBarMax
realnum GBarMax
Definition: thermal.h:142
t_struc::dr_ionfrac_limit
realnum dr_ionfrac_limit
Definition: struc.h:84
cdReasonGeo
void cdReasonGeo(FILE *ioOUT)
Definition: cddrive.cpp:175
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
t_dense::nzEdenBad
long int nzEdenBad
Definition: dense.h:200
chkCaHeps
STATIC void chkCaHeps(double *totwid)
Definition: prt_comment.cpp:2991
phycon
t_phycon phycon
Definition: phycon.cpp:6
FndLineHt
const TransitionProxy FndLineHt(long int *level)
Definition: lines_service.cpp:729
trace.h
t_he::frac_he0dest_23S_photo
double frac_he0dest_23S_photo
Definition: he.h:18
t_hydro::ndclev
long int ndclev
Definition: hydrogenic.h:105
t_mole_global::list
MoleculeList list
Definition: mole.h:317
t_phycon::BigJumpH2
realnum BigJumpH2
Definition: phycon.h:106
t_rfield::occmnu
realnum occmnu
Definition: rfield.h:473
two_photon::induc_dn_max
double induc_dn_max
Definition: two_photon.h:55
DoppVel
t_DoppVel DoppVel
Definition: doppvel.cpp:5
lgTestCodeCalled
bool lgTestCodeCalled
Definition: cddefines.cpp:11
t_conv::failmx
realnum failmx
Definition: conv.h:211
t_oxy::r4363Max
realnum r4363Max
Definition: oxy.h:17
ca
t_ca ca
Definition: ca.cpp:5
t_fudgec::lgFudgeUsed
bool lgFudgeUsed
Definition: fudgec.h:19
t_broke::lgFixit
bool lgFixit
Definition: broke.h:15
t_hydro::nbul
long int nbul
Definition: hydrogenic.h:107
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
t_hmi::h2dtot
realnum h2dtot
Definition: hmi.h:50
t_iso_sp::TwoNu
vector< two_photon > TwoNu
Definition: iso.h:586
wcnint
void wcnint(void)
Definition: warnings.cpp:13
t_timesc::time_therm_long
double time_therm_long
Definition: timesc.h:19
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
t_radius::depth
double depth
Definition: radius.h:38
lines_service.h
t_atmdat::HIonFracMax
double HIonFracMax
Definition: atmdat.h:169
t_opac::TauAbsGeo
realnum ** TauAbsGeo
Definition: opacity.h:82
t_broke::lgCheckit
bool lgCheckit
Definition: broke.h:18
t_pressure::ipPradMax_nzone
long int ipPradMax_nzone
Definition: pressure.h:146
dynamics.h
lgPrnErr
bool lgPrnErr
Definition: cddefines.cpp:13
input
t_input input
Definition: input.cpp:12
t_magnetic::lgB
bool lgB
Definition: magnetic.h:27
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
t_hydro::TeLyaMax
realnum TeLyaMax
Definition: hydrogenic.h:74
struc
t_struc struc
Definition: struc.cpp:6
GrainVar::GrnElecDonateMax
realnum GrnElecDonateMax
Definition: grainvar.h:530
iso.h
t_thermal::HeatLineMax
realnum HeatLineMax
Definition: thermal.h:164
t_oxy::poiii3Max
realnum poiii3Max
Definition: oxy.h:14
version.h
cdNotes
void cdNotes(FILE *ioOUT)
Definition: cddrive.cpp:309
t_opac::stimax
realnum stimax[2]
Definition: opacity.h:297
t_dynamics::lgRecom
bool lgRecom
Definition: dynamics.h:102
Wind::lgStatic
bool lgStatic(void) const
Definition: wind.h:24
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
t_struc::H2_abund
realnum * H2_abund
Definition: struc.h:72
t_StopCalc::chReasonStop
char chReasonStop[nCHREASONSTOP]
Definition: stopcalc.h:130
t_geometry::lgStatic
bool lgStatic
Definition: geometry.h:54
opac
t_opac opac
Definition: opacity.cpp:5
wind
Wind wind
Definition: wind.cpp:5
atmdat.h
EmissionProxy::TauIn
realnum & TauIn() const
Definition: emission.h:423
lgConserveEnergy
bool lgConserveEnergy()
Definition: energy.cpp:314
ipCOL_Hp
#define ipCOL_Hp
Definition: colden.h:26
t_colden::colden_old
realnum colden_old[NCOLD]
Definition: colden.h:40
t_thermal::lgTemperatureConstant
bool lgTemperatureConstant
Definition: thermal.h:32
t_colden::rjnmin
realnum rjnmin
Definition: colden.h:88
MIN2
#define MIN2
Definition: cddefines.h:761
TransitionProxy
Definition: transition.h:23
t_thermal::nUnstable
long int nUnstable
Definition: thermal.h:52
hyperfine
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
t_rfield::cinrat
double cinrat
Definition: rfield.h:294
t_conv::lgAutoIt
bool lgAutoIt
Definition: conv.h:256
t_struc::nzlim
long int nzlim
Definition: struc.h:19
cddrive.h
broke.h
nzone
long int nzone
Definition: cddefines.cpp:14
timesc
t_timesc timesc
Definition: timesc.cpp:5
t_input::lgBracketFound
bool lgBracketFound
Definition: input.h:67
t_hextra::cryden
realnum cryden
Definition: hextra.h:14
ipT3969
long ipT3969
Definition: atmdat_readin.cpp:82
t_geometry::covrt
realnum covrt
Definition: geometry.h:51
t_dense::EdenMax
double EdenMax
Definition: dense.h:193
radius
t_radius radius
Definition: radius.cpp:5
t_pressure::lgPradDen
bool lgPradDen
Definition: pressure.h:154
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:488
t_phycon::TEMP_STOP_DEFAULT
const double TEMP_STOP_DEFAULT
Definition: phycon.h:109
t_dynamics::lgTimeDependentStatic
bool lgTimeDependentStatic
Definition: dynamics.h:96
t_hmi::CoolH2DexcMax
realnum CoolH2DexcMax
Definition: hmi.h:48
t_rfield::uh
realnum uh
Definition: rfield.h:364
t_rt::DoubleTau
realnum DoubleTau
Definition: rt.h:247
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
diatomics::renorm_max
double renorm_max
Definition: h2_priv.h:336
broke
t_broke broke
Definition: broke.cpp:5
dense.h
t_rfield::lgOcc1Hi
bool lgOcc1Hi
Definition: rfield.h:480
ioPrnErr
FILE * ioPrnErr
Definition: cddefines.cpp:9
t_iterations::itermx
long int itermx
Definition: iterations.h:26
t_thermal::ConstTemp
realnum ConstTemp
Definition: thermal.h:44
t_hcmap::lgMapOK
bool lgMapOK
Definition: hcmap.h:30
t_thermal::power
double power
Definition: thermal.h:152
t_timesc::tcmptn
double tcmptn
Definition: timesc.h:16
mole
t_mole_local mole
Definition: mole.cpp:7
t_geometry::covaper
realnum covaper
Definition: geometry.h:44
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:371
t_hyperfine::cooling_max
realnum cooling_max
Definition: hyperfine.h:56
t_opac::thmin
realnum thmin
Definition: opacity.h:176
trace
t_trace trace
Definition: trace.cpp:5
t_hydro::pop2mx
realnum pop2mx
Definition: hydrogenic.h:56
t_phycon::lgPhysOK
bool lgPhysOK
Definition: phycon.h:101
t_DoppVel::TurbVel
realnum TurbVel
Definition: doppvel.h:12
t_geometry::lgZoneTrp
bool lgZoneTrp
Definition: geometry.h:90
cddefines.h
t_co::codtot
realnum codtot
Definition: co.h:16
hextra
t_hextra hextra
Definition: hextra.cpp:5
t_pressure::RadBetaMax
realnum RadBetaMax
Definition: pressure.h:136
ipH2s
const int ipH2s
Definition: iso.h:28
t_hcmap::lgMapDone
bool lgMapDone
Definition: hcmap.h:36
t_atoms::nNegOI
long int nNegOI
Definition: atoms.h:252
atoms
t_atoms atoms
Definition: atoms.cpp:5
t_radius::lgDrMinUsed
bool lgDrMinUsed
Definition: radius.h:180
t_radius::Radius
double Radius
Definition: radius.h:25
thermal
t_thermal thermal
Definition: thermal.cpp:5
t_radius::rinner
double rinner
Definition: radius.h:22
t_rfield::lgMMok
bool lgMMok
Definition: rfield.h:459
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
t_ionbal::ifail
long int ifail
Definition: ionbal.h:249
t_hmi::HeatH2DexcMax
realnum HeatH2DexcMax
Definition: hmi.h:46
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
cdNwcns
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
t_abund::lgAbnSolar
bool lgAbnSolar
Definition: abund.h:57
fudgec
t_fudgec fudgec
Definition: fudgec.cpp:5
cdSurprises
void cdSurprises(FILE *ioOUT)
Definition: cddrive.cpp:282
TauLine2
TransitionList TauLine2("TauLine2", &AnonStates)
t_called::lgTalk
bool lgTalk
Definition: called.h:12
t_tag_LineSv::SumLine
double SumLine[4]
Definition: lines.h:125
t_broke::lgBroke
bool lgBroke
Definition: broke.h:12
hyperfine.h
radius.h
t_struc::pres_radiation_lines_curr
realnum * pres_radiation_lines_curr
Definition: struc.h:43
colden
t_colden colden
Definition: colden.cpp:5
t_geometry::filpow
realnum filpow
Definition: geometry.h:20
t_phycon::BigJumpne
realnum BigJumpne
Definition: phycon.h:106
t_isoCTRL::lgCritDensLMix
bool lgCritDensLMix[NISO]
Definition: iso.h:395
t_geometry::lgSphere
bool lgSphere
Definition: geometry.h:24
t_rfield::nflux
long int nflux
Definition: rfield.h:43
t_hydro::nLyaHot
long int nLyaHot
Definition: hydrogenic.h:69
t_hmi::h2pmax
realnum h2pmax
Definition: hmi.h:120
t_hydro::cintot
double cintot
Definition: hydrogenic.h:92
hmi.h
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
t_rfield::lgXRayOK
bool lgXRayOK
Definition: rfield.h:461
t_colden::colden
realnum colden[NCOLD]
Definition: colden.h:38
MAX2
#define MAX2
Definition: cddefines.h:782
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
pressure.h
t_rfield::lgHPhtOK
bool lgHPhtOK
Definition: rfield.h:460
LIMELM
const int LIMELM
Definition: cddefines.h:258
RYDLAM
const UNUSED double RYDLAM
Definition: physconst.h:176
t_warnings::chRgcln
char chRgcln[2][INPUT_LINE_LENGTH]
Definition: warnings.h:46
t_rfield::lgCMB_set
bool lgCMB_set
Definition: rfield.h:124
DEPTH_OFFSET
const double DEPTH_OFFSET
Definition: cddefines.h:272
t_geometry::iEmissPower
long int iEmissPower
Definition: geometry.h:61
t_hmi::h2line_cool_frac
realnum h2line_cool_frac
Definition: hmi.h:52
t_oxy::poimax
realnum poimax
Definition: oxy.h:19
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
GrainVar::TotalDustHeat
realnum TotalDustHeat
Definition: grainvar.h:559
t_rfield::lgHabing
bool lgHabing
Definition: rfield.h:376
t_rfield::occmax
realnum occmax
Definition: rfield.h:472
warnings
t_warnings warnings
Definition: warnings.cpp:11
t_dynamics::HeatMax
double HeatMax
Definition: dynamics.h:68
t_StopCalc::TempLoStopZone
realnum TempLoStopZone
Definition: stopcalc.h:42
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1228
save.h
Wind::lgBallistic
bool lgBallistic(void) const
Definition: wind.h:31
t_elementnames::chElementNameShort
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
co.h
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
he.h
t_phycon::TEnerDen
double TEnerDen
Definition: phycon.h:98
iteration
long int iteration
Definition: cddefines.cpp:16
t_radius::lgDR2Big
bool lgDR2Big
Definition: radius.h:168
TauLines
TransitionList TauLines("TauLines", &AnonStates)
t_radius::r1r0sq
double r1r0sq
Definition: radius.h:49
t_oxy::poiii2Max
realnum poiii2Max
Definition: oxy.h:13
abund
t_abund abund
Definition: abund.cpp:5
t_pressure::lgPres_radiation_ON
bool lgPres_radiation_ON
Definition: pressure.h:130
hydrogenic.h
t_he::frac_he0dest_23S
double frac_he0dest_23S
Definition: he.h:16
t_save::lgPunContinuum
bool lgPunContinuum
Definition: save.h:251
he
t_he he
Definition: he.cpp:5
t_rfield::lgPlasNu
bool lgPlasNu
Definition: rfield.h:445
t_oxy::r5007Max
realnum r5007Max
Definition: oxy.h:16
t_iterations::lgLastIt
bool lgLastIt
Definition: iterations.h:36
t_continuum::lgCoStarInterpolationCaution
bool lgCoStarInterpolationCaution
Definition: continuum.h:94
ipH2p
const int ipH2p
Definition: iso.h:29
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:254
grainvar.h
t_hydro::fbul
realnum fbul
Definition: hydrogenic.h:106
t_conv::ifailz
long int ifailz[12]
Definition: conv.h:241
timesc.h
rt.h
t_thermal::chCoolHeatMax
char chCoolHeatMax[NCOLNT_LAB_LEN+1]
Definition: thermal.h:107
powi
double powi(double, long int)
Definition: service.cpp:604
t_conv::chNotConverged
char chNotConverged[INPUT_LINE_LENGTH]
Definition: conv.h:135
t_pressure::ipPradMax_line
long int ipPradMax_line
Definition: pressure.h:143
t_rfield::lgGamrOK
bool lgGamrOK
Definition: rfield.h:462
t_rfield::anu
double * anu
Definition: rfield.h:58
wind.h
ionbal.h
PrtComment
void PrtComment(void)
Definition: prt_comment.cpp:65
GrainVar::lgNegGrnDrg
bool lgNegGrnDrg
Definition: grainvar.h:482
ipH3p
const int ipH3p
Definition: iso.h:31
magnetic.h
t_fudgec::nfudge
long int nfudge
Definition: fudgec.h:17
cdCautions
void cdCautions(FILE *ioOUT)
Definition: cddrive.cpp:226
t_dense::lgEdenBad
bool lgEdenBad
Definition: dense.h:227
hmi
t_hmi hmi
Definition: hmi.cpp:5
t_hydro::HCollIonMax
realnum HCollIonMax
Definition: hydrogenic.h:86
GrainVar::dclmax
realnum dclmax
Definition: grainvar.h:561
t_rfield::tbrmnu
realnum tbrmnu
Definition: rfield.h:475
physconst.h
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
dynamics
t_dynamics dynamics
Definition: dynamics.cpp:44
t_struc::ednstr
realnum * ednstr
Definition: struc.h:30
ipCOL_H2g
#define ipCOL_H2g
Definition: colden.h:16
called
t_called called
Definition: called.cpp:5
conv
t_conv conv
Definition: conv.cpp:5
gv
GrainVar gv
Definition: grainvar.cpp:5
t_warnings::lgCautns
bool lgCautns
Definition: warnings.h:57
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:814
t_phycon::BigJumpCO
realnum BigJumpCO
Definition: phycon.h:106
magnetic
t_magnetic magnetic
Definition: magnetic.cpp:17
t_opac::telec
realnum telec
Definition: opacity.h:175
two_photon
Definition: two_photon.h:9
t_ca::Ca2RmLya
realnum Ca2RmLya
Definition: ca.h:11
rt
t_rt rt
Definition: rt.cpp:5
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
t_rfield::ConInterOut
realnum * ConInterOut
Definition: rfield.h:164
t_rfield::lgComUndr
bool lgComUndr
Definition: rfield.h:298
ipCOL_HTOT
#define ipCOL_HTOT
Definition: colden.h:12
t_rt::lgMaserCapHit
bool lgMaserCapHit
Definition: rt.h:274
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
t_input::lgUnderscoreFound
bool lgUnderscoreFound
Definition: input.h:63
hextra.h
t_rfield::EnergyBremsThin
realnum EnergyBremsThin
Definition: rfield.h:246
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
cdWarnings
void cdWarnings(FILE *ioPNT)
Definition: cddrive.cpp:198
t_rfield::tbrmax
realnum tbrmax
Definition: rfield.h:474
t_rt::lgMaserSetDR
bool lgMaserSetDR
Definition: rt.h:270
t_warnings::lgWarngs
bool lgWarngs
Definition: warnings.h:56
Wind::emdot
realnum emdot
Definition: wind.h:39
t_dense::lgSetIoniz
bool lgSetIoniz[LIMELM]
Definition: dense.h:149
hcmap.h
AgeCheck
void AgeCheck(void)
Definition: age_check.cpp:14
pressure
t_pressure pressure
Definition: pressure.cpp:5
phycon.h
geometry
t_geometry geometry
Definition: geometry.cpp:5
t_iso_sp::n_HighestResolved_local
long int n_HighestResolved_local
Definition: iso.h:507
t_rfield::tbr4nu
realnum tbr4nu
Definition: rfield.h:476
t_rfield::comtot
double comtot
Definition: rfield.h:291
t_pressure::chLineRadPres
char chLineRadPres[101]
Definition: pressure.h:149
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
ShowMe
void ShowMe(void)
Definition: service.cpp:181
EmissionProxy::TauTot
realnum & TauTot() const
Definition: emission.h:433
ipH1s
const int ipH1s
Definition: iso.h:27
Wind::windv
realnum windv
Definition: wind.h:18
iterations.h
bangin
void bangin(char *chLine)
Definition: warnings.cpp:73
continuum
t_continuum continuum
Definition: continuum.cpp:5
ipCOL_H2p
#define ipCOL_H2p
Definition: colden.h:20
ipCOL_H0
#define ipCOL_H0
Definition: colden.h:22
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
called.h
t_thermal::FreeFreeTotHeat
double FreeFreeTotHeat
Definition: thermal.h:161
t_hmi::hmitot
double hmitot
Definition: hmi.h:25
continuum.h
oxy.h
h2.h
LineSv
LinSv * LineSv
Definition: cdinit.cpp:70
hcmap
t_hcmap hcmap
Definition: hcmap.cpp:21
t_mole_global::lgH2Ozer
bool lgH2Ozer
Definition: mole.h:283
t_struc::molecules
realnum ** molecules
Definition: struc.h:71
t_geometry::FillFac
realnum FillFac
Definition: geometry.h:19
stopcalc.h
t_geometry::lgStaticNoIt
bool lgStaticNoIt
Definition: geometry.h:74
t_phycon::te
double te
Definition: phycon.h:11
t_hydro::nZTLaMax
long int nZTLaMax
Definition: hydrogenic.h:77
t_rfield::plsfrqmax
realnum plsfrqmax
Definition: rfield.h:449
NISO
const int NISO
Definition: cddefines.h:261
warnings.h
t_atmdat::HCharHeatMax
double HCharHeatMax
Definition: atmdat.h:154
t_thermal::lgEdnGTcm
bool lgEdnGTcm
Definition: thermal.h:65
prt_smooth_predictions
STATIC void prt_smooth_predictions(void)
Definition: prt_comment.cpp:3028
t_radius::lgCylnOn
bool lgCylnOn
Definition: radius.h:121
GrainVar::lgDustOn
bool lgDustOn() const
Definition: grainvar.h:471
GrainVar::GrnElecHoldMax
realnum GrnElecHoldMax
Definition: grainvar.h:531
input.h
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
t_he::nzone
long nzone
Definition: he.h:20
t_thermal::totcol
double totcol
Definition: thermal.h:110
t_mole_global::num_calc
int num_calc
Definition: mole.h:314
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
t_hydro::TLyaMax
realnum TLyaMax
Definition: hydrogenic.h:72
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
oxy
t_oxy oxy
Definition: oxy.cpp:5
save
t_save save
Definition: save.cpp:5
warnin
void warnin(char *chLine)
Definition: warnings.cpp:27
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
HFLines
TransitionList HFLines("HFLines", &AnonStates)
t_hydro::FracInd
realnum FracInd
Definition: hydrogenic.h:104
GrainVar::dphmax
realnum dphmax
Definition: grainvar.h:560
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
t_hmi::h2dfrc
realnum h2dfrc
Definition: hmi.h:49
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
diatomics::renorm_min
double renorm_min
Definition: h2_priv.h:337
ca.h
Generated by   doxygen 1.8.17