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

Generated for cloudy by doxygen 1.8.4