cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
prt_final.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*PrtFinal create PrtFinal pages of printout, emission line intensities, etc */
4 /*StuffComment routine to stuff comments into the stack of comments, def in lines.h */
5 /*gett2o3 analyze computed [OIII] spectrum to get t^2 */
6 /*gett2 analyze computed structure to get structural t^2 */
7 #include "cddefines.h"
8 #include "iso.h"
9 #include "cddrive.h"
10 #include "dynamics.h"
11 #include "lines.h"
12 #include "taulines.h"
13 #include "warnings.h"
14 #include "phycon.h"
15 #include "dense.h"
16 #include "grainvar.h"
17 #include "h2.h"
18 #include "hmi.h"
19 #include "thermal.h"
20 #include "hydrogenic.h"
21 #include "rt.h"
22 #include "atmdat.h"
23 #include "timesc.h"
24 #include "opacity.h"
25 #include "struc.h"
26 #include "pressure.h"
27 #include "conv.h"
28 #include "geometry.h"
29 #include "called.h"
30 #include "iterations.h"
31 #include "version.h"
32 #include "colden.h"
33 #include "input.h"
34 #include "radius.h"
35 #include "peimbt.h"
36 #include "ipoint.h"
37 #include "mean.h"
38 #include "wind.h"
39 #include "prt.h"
40 #include "rfield.h"
41 #include "freebound.h"
42 #include "lines_service.h"
43 
44 // helper routine to center a line in the output
45 STATIC void PrintCenterLine(FILE* io, // file pointer
46  const char chLine[], // string to print, should not end with '\n'
47  size_t ArrLen, // length of chLine
48  size_t LineLen) // width of the line
49 {
50  unsigned long StrLen = min(strlen(chLine),ArrLen);
51  ASSERT( StrLen < LineLen );
52  unsigned long pad = (LineLen-StrLen)/2;
53  for( unsigned long i=0; i < pad; ++i )
54  fprintf( io, " " );
55  fprintf( io, "%s\n", chLine );
56 }
57 
58 /*gett2o3 analyze computed [OIII] spectrum to get t^2 */
59 STATIC void gett2o3(realnum *tsqr);
60 
61 /*gett2 analyze computed structure to get structural t^2 */
62 STATIC void gett2(realnum *tsqr);
63 
64 /* helper routine for printing averaged quantities */
65 inline void PrintRatio(double q1, double q2)
66 {
67  double ratio = ( q2 > SMALLFLOAT ) ? q1/q2 : 0.;
68  fprintf( ioQQQ, " " );
69  fprintf( ioQQQ, PrintEfmt("%9.2e", ratio) );
70  return;
71 }
72 
74 {
75  vector<realnum> sclsav, scaled;
76  vector<long int> ipSortLines;
77  vector<double> xLog_line_lumin;
78  vector<short int> lgPrt;
79  vector<long> Slines;
80 
81  long int
82  i,
83  ipEmType,
84  iprnt,
85  nline,
86  j,
87  ipNegIntensity[33],
88  nNegIntenLines;
89 
90  double a,
91  /* N.B. 8 is used for following preset in loop */
92  d[8],
93  snorm;
94 
95 
96  DEBUG_ENTRY( "PrintSpectrum()" );
97 
98 
99  /*----------------------------------------------------------
100  *
101  * first set scaled lines */
102 
103  /* get space for scaled */
104  scaled.resize(LineSave.nsum);
105 
106  /* get space for xLog_line_lumin */
107  xLog_line_lumin.resize(LineSave.nsum);
108 
109  /* this is option to not print certain contributions */
110  /* gjf 98 jun 10, get space for array lgPrt */
111  lgPrt.resize(LineSave.nsum);
112 
113  /* get space for sclsav */
114  sclsav.resize(LineSave.nsum );
115 
116  Slines.resize(LineSave.nsum);
117 
118  /* get space for array of indices for lines, for possible sorting */
119  ipSortLines.resize(LineSave.nsum );
120 
121  ASSERT( LineSave.ipNormWavL >= 0 );
122 
123  /* option to also print usual first two sets of line arrays
124  * but for two sets of cumulative arrays for time-dependent sims too */
125  int nEmType = 2;
127  nEmType = 4;
128 
129  for( ipEmType=0; ipEmType<nEmType; ++ipEmType )
130  {
131  if( ! prt.lgPrintBlockIntrinsic && (ipEmType == 0 || ipEmType == 2) )
132  continue;
133  if( ! prt.lgPrintBlockEmergent && (ipEmType == 1 || ipEmType == 3) )
134  continue;
135 
136  if( ipEmType > 1 && strcmp( rfield.chCumuType, "NONE" ) == 0 )
137  continue;
138 
139  /* this is the intensity of the line spectrum will be normalized to */
140  snorm = LineSave.lines[LineSave.ipNormWavL].SumLine(ipEmType);
141 
142  /* check that this line has positive intensity */
143  if( ((snorm <= SMALLDOUBLE ) || (LineSave.ipNormWavL < 0)) || (LineSave.ipNormWavL > LineSave.nsum) )
144  {
145  char wl_str[20] = "";
146  sprt_wl( wl_str, LineSave.lines[LineSave.ipNormWavL].wavelength() );
147  fprintf(ioQQQ,
148  "\n\n"
149  " >>PROBLEM Normalization line (\"%s\" %s) has small or zero intensity, its value was %.2e and its intensity was set to 1.\n"
150  " >>Please consider using another normalization line (this is set with the NORMALIZE command).\n",
151  LineSave.lines[LineSave.ipNormWavL].chALab(), wl_str, snorm);
152  fprintf( ioQQQ, " >>The relative intensities will be meaningless, and many lines may appear too faint.\n" );
153  snorm = 1.;
154  }
155  for( i=0; i < LineSave.nsum; i++ )
156  {
157  /* Do only for emergent lines */
158  if( ipEmType == 1 || ipEmType == 3 )
159  LineSave.lines[i].checkEmergent( ipEmType );
160 
161  /* when normalization line is off-scale small (generally a model
162  * with mis-set parameters) the scaled intensity can be larger than
163  * a realnum - this is not actually a problem since the number will
164  * overflow the format and hence be unreadable */
165  double scale = LineSave.lines[i].SumLine(ipEmType)/snorm*LineSave.ScaleNormLine;
166  /* this will become a realnum, so limit dynamic range */
167  scale = MIN2(BIGFLOAT , scale );
168  scale = MAX2( -BIGFLOAT , scale );
169 
170  /* find logs of ALL line intensities/luminosities */
171  scaled[i] = (realnum)scale;
172 
173  // the keyword volatile works around a bug in g++
174  // see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65425
175  volatile double line_lumin = LineSave.lines[i].SumLine(ipEmType) * radius.Conv2PrtInten;
176  xLog_line_lumin[i] = ( line_lumin > 0. ) ? log10(line_lumin) : 0.;
177  }
178 
179  /* now find which lines to print, which to ignore because they are the wrong type */
180  for( i=0; i < LineSave.nsum; i++ )
181  {
182  /* never print unit normalization check, at least in main list */
183  if( LineSave.lines[i].isUnit() || LineSave.lines[i].isUnitD() )
184  lgPrt[i] = false;
185  else if( !prt.lgPrnColl && LineSave.lines[i].isCollisional() )
186  lgPrt[i] = false;
187  else if( !prt.lgPrnPump && LineSave.lines[i].isPump() )
188  lgPrt[i] = false;
189  else if( !prt.lgPrnInwd && ( LineSave.lines[i].isInward() ||
190  LineSave.lines[i].isInwardContinuum() ||
191  LineSave.lines[i].isInwardTotal() )
192  )
193  lgPrt[i] = false;
194  else if( !prt.lgPrnHeat && LineSave.lines[i].isHeat() )
195  lgPrt[i] = false;
196  else
197  lgPrt[i] = true;
198  }
199 
200  /* do not print relatively faint lines unless requested */
201  nNegIntenLines = 0;
202 
203  /* set ipNegIntensity to bomb to make sure set in following */
204  for(i=0; i< 32; i++ )
205  {
206  ipNegIntensity[i] = LONG_MAX;
207  }
208 
209  for(i=0;i<8;++i)
210  {
211  d[i] = -DBL_MAX;
212  }
213 
214  /* create header for blocks of emission line intensities */
215  const char chIntensityType[4][100]=
216  {" Intrinsic" , " Emergent" , "Cumulative intrinsic" , "Cumulative emergent" };
217  ASSERT( ipEmType==0 || ipEmType==1 || ipEmType==2 || ipEmType==3 );
218  /* if true then printing in 4 columns of lines, this is offset to
219  * center the title */
220  fprintf( ioQQQ, "\n" );
221  if( prt.lgPrtLineArray )
222  fprintf( ioQQQ, " " );
223  fprintf( ioQQQ, "%s" , chIntensityType[ipEmType] );
224  fprintf( ioQQQ, " line intensities\n" );
225  // caution about emergent intensities when outward optical
226  // depths are not yet known
227  if( ipEmType==1 && iteration==1 )
228  fprintf(ioQQQ," Caution: emergent intensities are not reliable on the "
229  "first iteration.\n");
230 
231  /* option to only print brighter lines */
232  if( prt.lgFaintOn )
233  {
234  iprnt = 0;
235  for( i=0; i < LineSave.nsum; i++ )
236  {
237  /* this insanity can happen when arrays overrun */
238  ASSERT( iprnt <= i);
239  if( scaled[i] >= prt.TooFaint && lgPrt[i] )
240  {
241  if( prt.lgPrtLineLog )
242  {
243  xLog_line_lumin[iprnt] = log10(LineSave.lines[i].SumLine(ipEmType) * radius.Conv2PrtInten);
244  }
245  else
246  {
247  xLog_line_lumin[iprnt] = LineSave.lines[i].SumLine(ipEmType) * radius.Conv2PrtInten;
248  }
249  sclsav[iprnt] = scaled[i];
250  /* check that null is correct, string overruns have
251  * been a problem in the past */
252  Slines[iprnt] = i;
253  ++iprnt;
254  }
255  else if( -scaled[i] > prt.TooFaint && lgPrt[i] )
256  {
257  /* negative intensities occur if line absorbs continuum */
258  ipNegIntensity[nNegIntenLines] = i;
259  nNegIntenLines = MIN2(32,nNegIntenLines+1);
260  }
261  /* special labels to give id for blocks of lines
262  * do not add these labels when sorting by wavelength since invalid */
263  else if( LineSave.lines[i].isSeparator() &&!prt.lgSortLines )
264  {
265  xLog_line_lumin[iprnt] = 0.;
266  sclsav[iprnt] = 0.;
267  Slines[iprnt] = i;
268  ++iprnt;
269  }
270  }
271  }
272 
273  else
274  {
275  /* print everything */
276  iprnt = LineSave.nsum;
277  for( i=0; i < LineSave.nsum; i++ )
278  {
279  if( LineSave.lines[i].isSeparator() )
280  {
281  xLog_line_lumin[i] = 0.;
282  sclsav[i] = 0.;
283  }
284  else
285  {
286  sclsav[i] = scaled[i];
287  }
288  Slines[i] = i;
289  if( scaled[i] < 0. )
290  {
291  ipNegIntensity[nNegIntenLines] = i;
292  nNegIntenLines = MIN2(32,nNegIntenLines+1);
293  }
294  }
295  }
296 
297  /* reorder lines according to wavelength for comparison with spectrum
298  * including writing out the results */
299  if( prt.lgSortLines )
300  {
301  int lgFlag;
303  {
304  /* first check if wavelength range specified */
305  if( prt.wlSort1 >-0.1 )
306  {
307  j = 0;
308  /* skip over those lines not desired */
309  for( i=0; i<iprnt; ++i )
310  {
311  realnum wavelength=LineSave.lines[Slines[i]].wavelength();
312  if( wavelength>= prt.wlSort1 && wavelength<= prt.wlSort2 )
313  {
314  if( j!=i )
315  {
316  sclsav[j] = sclsav[i];
317  /* >>chng 05 feb 03, add this, had been left out,
318  * thanks to Marcus Copetti for discovering the bug */
319  xLog_line_lumin[j] = xLog_line_lumin[i];
320  Slines[j] = Slines[i];
321  }
322  ++j;
323  }
324  }
325  iprnt = j;
326  }
327 
328  vector<realnum> wavelength(iprnt);
329  for( i=0; i<iprnt; ++i )
330  {
331  wavelength[i] = LineSave.lines[Slines[i]].wavelength();
332  }
333 
334  /*spsort netlib routine to sort array returning sorted indices */
335  if( iprnt > 0 )
336  {
337  spsort(&wavelength[0],
338  iprnt,
339  &ipSortLines[0],
340  /* flag saying what to do - 1 sorts into increasing order, not changing
341  * the original routine */
342  1,
343  &lgFlag);
344  if( lgFlag )
345  {
346  fprintf(ioQQQ," >>> PROBLEM spsort reports error\n");
347  TotalInsanity();
348  }
349  }
350  }
351  else if( prt.lgSortLineIntensity )
352  {
353  if( iprnt > 0 )
354  {
355  /*spsort netlib routine to sort array returning sorted indices */
356  spsort(&sclsav[0],
357  iprnt,
358  &ipSortLines[0],
359  /* flag saying what to do - -1 sorts into decreasing order, not changing
360  * the original routine */
361  -1,
362  &lgFlag);
363  if( lgFlag )
364  TotalInsanity();
365  }
366  }
367  else
368  {
369  /* only to keep lint happen, or in case vars clobbered */
370  TotalInsanity();
371  }
372  }
373  else
374  {
375  /* do not sort lines (usual case) simply print in incoming order */
376  for( i=0; i<iprnt; ++i )
377  {
378  ipSortLines[i] = i;
379  }
380  }
381 
382  /* the number of lines to print better be positive */
383  if (iprnt == 0)
384  fprintf(ioQQQ," >>>> No strong lines found <<<<\n");
385 
386  /* print out all lines which made it through the faint filter,
387  * there are iprnt lines to print
388  * can print in either 4 column (the default ) or one long
389  * column of lines */
390  if( prt.lgPrtLineArray )
391  {
392  /* four or three columns ? - depends on how many sig figs */
393  if( LineSave.sig_figs >= 5 )
394  {
395  nline = (iprnt + 2)/3;
396  }
397  else
398  {
399  /* nline will be the number of horizontal lines -
400  * the 4 represents the 4 columns */
401  nline = (iprnt + 3)/4;
402  }
403  }
404  else
405  {
406  /* this option print a single column of emission lines */
407  nline = iprnt;
408  }
409 
410  /* now loop over the spectrum, printing lines */
411  for( i=0; i < nline; i++ )
412  {
413 
414  /* this skips over nline per step, to go to the next column in
415  * the output */
416  for( j=i; j<iprnt; j = j + nline)
417  {
418  /* index with possibly reordered set of lines */
419  long ipLin = ipSortLines[j];
420  /* this is the actual print statement for the emission line
421  * spectrum*/
422  if( LineSave.lines[Slines[ipLin]].isSeparator() )
423  {
424  /* special labels */
425  /*fprintf( ioQQQ, "1111222223333333444444444 " ); */
426  /* this string was set in */
428  (int)LineSave.lines[Slines[ipLin]].wavelength()] );
429  }
430  else
431  {
432  /* the label for the line */
433  LineSave.lines[Slines[ipLin]].prt(ioQQQ);
434 
435  /* print the log of the intensity/luminosity of the
436  * lines, the usual case */
437  fprintf( ioQQQ, " " );
438  if( prt.lgPrtLineLog )
439  {
440  fprintf( ioQQQ, "%*.3f", prt_linecol.absint_len, xLog_line_lumin[ipLin] );
441  }
442  else
443  {
444  /* print linear quantity instead */
445  fprintf( ioQQQ, " %*.2e", prt_linecol.absint_len, xLog_line_lumin[ipLin] );
446  }
447 
448  /* print scaled intensity with various formats,
449  * depending on order of magnitude. value is
450  * always the same but the format changes. */
451  fprintf( ioQQQ, " " );
452  if( sclsav[ipLin] < 9999.5 )
453  {
454  fprintf( ioQQQ, "%*.4f", prt_linecol.relint_len, sclsav[ipLin] );
455  }
456  else if( sclsav[ipLin] < 99999.5 )
457  {
458  fprintf( ioQQQ, "%*.3f", prt_linecol.relint_len, sclsav[ipLin] );
459  }
460  else if( sclsav[ipLin] < 999999.5 )
461  {
462  fprintf( ioQQQ, "%*.2f", prt_linecol.relint_len, sclsav[ipLin] );
463  }
464  else if( sclsav[ipLin] < 9999999.5 )
465  {
466  fprintf( ioQQQ, "%*.1f", prt_linecol.relint_len, sclsav[ipLin] );
467  }
468  else
469  {
470  fprintf( ioQQQ, "%s", prt_linecol.relint_outrange.c_str() );
471  }
472  }
473 
474  /* now end the block with some spaces to set next one off */
475  if( j+nline < iprnt )
476  fprintf( ioQQQ, "%s", prt_linecol.col_gap.c_str() );
477  }
478  /* now end the lines */
479  fprintf( ioQQQ, "\n" );
480  }
481 
482  if( nNegIntenLines > 0 )
483  {
484  fprintf( ioQQQ, " Lines with negative intensities - Linear intensities relative to normalize line.\n" );
485  fprintf( ioQQQ, " " );
486 
487  for( i=0; i < nNegIntenLines; ++i )
488  {
489  fprintf( ioQQQ, "%ld %s %.1e, ",
490  ipNegIntensity[i],
491  LineSave.lines[ipNegIntensity[i]].label().c_str(),
492  scaled[ipNegIntensity[i]] );
493  }
494  fprintf( ioQQQ, "\n" );
495  }
496  }
497 
498 
499  /* now find which were the very strongest, more that 5% of cooling */
500  j = 0;
501  for( i=0; i < LineSave.nsum; i++ )
502  {
503  a = (double)LineSave.lines[i].SumLine(0)/(double)thermal.totcol;
504  if( (a >= 0.05) && LineSave.lines[i].chSumTyp() == 'c' )
505  {
506  ipNegIntensity[j] = i;
507  d[j] = a;
508  j = MIN2(j+1,7);
509  }
510  }
511 
512  fprintf( ioQQQ, "\n\n\n %s\n", input.chTitle );
513  if( j != 0 )
514  {
515  fprintf( ioQQQ, " Cooling: " );
516  for( i=0; i < j; i++ )
517  {
518  fprintf( ioQQQ, " ");
519  LineSave.lines[ipNegIntensity[i]].prt(ioQQQ);
520  fprintf( ioQQQ, ":%5.3f",
521  d[i] );
522  }
523  fprintf( ioQQQ, " \n" );
524  }
525 
526  /* now find strongest heating, more that 5% of total */
527  j = 0;
528  for( i=0; i < LineSave.nsum; i++ )
529  {
530  a = safe_div((double)LineSave.lines[i].SumLine(0),(double)thermal.power,0.0);
531  if( (a >= 0.05) && LineSave.lines[i].chSumTyp() == 'h' )
532  {
533  ipNegIntensity[j] = i;
534  d[j] = a;
535  j = MIN2(j+1,7);
536  }
537  }
538 
539  if( j != 0 )
540  {
541  fprintf( ioQQQ, " Heating: " );
542  for( i=0; i < j; i++ )
543  {
544  fprintf( ioQQQ, " ");
545  LineSave.lines[ipNegIntensity[i]].prt(ioQQQ);
546  fprintf( ioQQQ, ":%5.3f",
547  d[i] );
548  }
549  fprintf( ioQQQ, " \n" );
550  }
551 
552  return;
553 }
554 
555 void PrtFinal(void)
556 {
557  char chCKey[5],
558  chGeo[7],
559  chPlaw[21];
560  const char* chUnit;
561 
562  long int
563  i,
564  ip2500,
565  ip2kev,
566  j;
567  double o4363,
568  bacobs,
569  absint,
570  bacthn,
571  hbcab,
572  hbeta,
573  o5007;
574 
575  double a,
576  ajmass,
577  ajmin,
578  alfox,
579  bot,
580  bottom,
581  bremtk,
582  coleff,
583  dmean,
584  ferode,
585  he4686,
586  he5876,
587  heabun,
588  hescal,
589  pion,
590  plow,
591  powerl,
592  qion,
593  qlow,
594  rate,
595  ratio,
596  reclin,
597  rjeans,
598  tmean,
599  top,
600  THI,/* HI 21 cm spin temperature */
601  uhel,
602  uhl,
603  usp,
604  wmean,
605  znit;
606 
607  double bac,
608  tel,
609  x;
610 
611  DEBUG_ENTRY( "PrtFinal()" );
612 
613  /* return if not talking */
614  /* >>chng 05 nov 11, also if nzone is zero, sign of abort before model got started */
615  if( !called.lgTalk || (lgAbort && nzone< 1) )
616  {
617  return;
618  }
619 
620  /* print out header, or parts that were saved */
621 
622  /* this would be a major logical error */
623  ASSERT( LineSave.nsum > 1 );
624 
625  /* print name and version number */
626  fprintf( ioQQQ, "\f\n");
627  fprintf( ioQQQ, "%23c", ' ' );
628  int len = (int)strlen(t_version::Inst().chVersion);
629  int repeat = (72-len)/2;
630  for( i=0; i < repeat; ++i )
631  fprintf( ioQQQ, "*" );
632  fprintf( ioQQQ, "> Cloudy %s <", t_version::Inst().chVersion );
633  for( i=0; i < 72-repeat-len; ++i )
634  fprintf( ioQQQ, "*" );
635  fprintf( ioQQQ, "\n" );
636 
637  for( i=0; i <= input.nSave; i++ )
638  {
639  char chCard[INPUT_LINE_LENGTH];
640  /* print command line, unless it is a continue line */
641 
642  /* copy start of command to key, making it into capitals */
643  cap4(chCKey,input.chCardSav[i]);
644 
645  /* copy input to CAPS to make sure hide not on it */
646  strcpy( chCard , input.chCardSav[i] );
647  caps( chCard );
648 
649  /* print if not continue or hide */
650  /* >>chng 04 jan 21, add hide option */
651  if( (strcmp(chCKey,"CONT")!= 0) && !nMatch( "HIDE" , chCard) )
652  fprintf( ioQQQ, "%23c* %-80s*\n", ' ', input.chCardSav[i] );
653  }
654 
655  /* print log of ionization parameter if greater than zero - U==0 for PDR calcs */
656  if( rfield.uh > 0. )
657  {
658  a = log10(rfield.uh);
659  }
660  else
661  {
662  a = -37.;
663  }
664 
665  fprintf( ioQQQ,
666  " *********************************> Log(U):%6.2f <*********************************\n",
667  a );
668 
669  if( t_version::Inst().nBetaVer > 0 )
670  {
671  fprintf( ioQQQ,
672  "\n This is a beta test version of the code, and is intended for testing only.\n\n" );
673  }
674 
675  if( warnings.lgWarngs )
676  {
677  fprintf( ioQQQ, " \n" );
678  fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
679  fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
680  fprintf( ioQQQ, " >>>>>>>>>> Warning! Warnings exist, this calculation has serious problems.\n" );
681  fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
682  fprintf( ioQQQ, " >>>>>>>>>> Warning!\n" );
683  fprintf( ioQQQ, " \n" );
684  }
685  else if( warnings.lgCautns )
686  {
687  fprintf( ioQQQ,
688  " >>>>>>>>>> Cautions are present.\n" );
689  }
690 
691  if( conv.lgBadStop )
692  {
693  fprintf( ioQQQ,
694  " >>>>>>>>>> The calculation stopped for unintended reasons.\n" );
695  }
696 
697  if( iterations.lgIterAgain )
698  {
699  fprintf( ioQQQ,
700  " >>>>>>>>>> Another iteration is needed. Use the ITERATE command.\n" );
701  }
702 
703  /* open or closed geometry?? */
704  if( geometry.lgSphere )
705  {
706  strcpy( chGeo, "Closed" );
707  }
708  else
709  {
710  strcpy( chGeo, " Open" );
711  }
712 
713  /* now give description of pressure law */
714  if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
715  {
716  strcpy( chPlaw, " Constant Pressure " );
717  }
718 
719  else if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
720  {
721  strcpy( chPlaw, " Constant Density " );
722  }
723 
724  else if( (strcmp(dense.chDenseLaw,"POWD") == 0 || strcmp(dense.chDenseLaw
725  ,"POWR") == 0) || strcmp(dense.chDenseLaw,"POWC") == 0 )
726  {
727  strcpy( chPlaw, " Power Law Density " );
728  }
729 
730  else if( strcmp(dense.chDenseLaw,"SINE") == 0 )
731  {
732  strcpy( chPlaw, " Rapid Fluctuations " );
733  }
734 
735  else if( strncmp(dense.chDenseLaw , "DLW" , 3) == 0 )
736  {
737  strcpy( chPlaw, " Special Density Lw " );
738  }
739 
740  else if( strcmp(dense.chDenseLaw,"HYDR") == 0 )
741  {
742  strcpy( chPlaw, " Hydrostatic Equlib " );
743  }
744 
745  else if( strcmp(dense.chDenseLaw,"WIND") == 0 )
746  {
747  strcpy( chPlaw, " Radia Driven Wind " );
748  }
749 
750  else if( strcmp(dense.chDenseLaw,"DYNA") == 0 )
751  {
752  strcpy( chPlaw, " Dynamical Flow " );
753  }
754 
755  else if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
756  {
757  strcpy( chPlaw, " Globule " );
758  }
759 
760  else
761  {
762  strcpy( chPlaw, " UNRECOGNIZED CPRES " );
763  }
764 
765  fprintf( ioQQQ,
766  "\n Emission Line Spectrum. %20.20sModel. %6.6s geometry. Iteration %ld of %ld.\n",
767  chPlaw, chGeo, iteration, iterations.itermx + 1 );
768 
769  /* emission lines as logs of total luminosity */
771  {
772  char chLine[INPUT_LINE_LENGTH];
773  if( geometry.iEmissPower == 1 && !geometry.lgSizeSet )
774  chUnit = "/arcsec";
775  else if( geometry.iEmissPower == 0 && !geometry.lgSizeSet )
776  chUnit = "/arcsec^2";
777  else
778  chUnit = "";
779 
780  sprintf( chLine, "Flux observed at the Earth (erg/s/cm^2%s).", chUnit );
781  PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
782  }
783  else if( prt.lgSurfaceBrightness )
784  {
785  char chLine[INPUT_LINE_LENGTH];
787  chUnit = "sr";
788  else
789  chUnit = "arcsec^2";
790 
791  sprintf( chLine, "Surface brightness (erg/s/cm^2/%s).", chUnit );
792  PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
793  }
794  else if( radius.lgPredLumin )
795  {
796  char chLine[INPUT_LINE_LENGTH];
797  if( geometry.iEmissPower == 2 )
798  chUnit = "erg/s";
799  else if( geometry.iEmissPower == 1 )
800  chUnit = "erg/s/cm";
801  else if( geometry.iEmissPower == 0 )
802  chUnit = "erg/s/cm^2";
803  else
804  TotalInsanity();
805 
806  char chCoverage[INPUT_LINE_LENGTH];
807  if( fp_equal( geometry.covgeo, realnum(1.) ) )
808  sprintf( chCoverage, "with full coverage" );
809  else
810  sprintf( chCoverage, "with a covering factor of %.1f%%", geometry.covgeo*100. );
811 
812  if( radius.lgCylnOn )
813  sprintf( chLine, "Luminosity (%s) emitted by a partial cylinder %s.", chUnit, chCoverage );
814  else
815  sprintf( chLine, "Luminosity (%s) emitted by a shell %s.", chUnit, chCoverage );
816  PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
817 
818  if( geometry.iEmissPower != 2 )
819  {
820  const char* chAper;
821  if( geometry.iEmissPower == 1 )
822  chAper = "long slit";
823  else if( geometry.iEmissPower == 0 )
824  chAper = "pencil beam";
825  else
826  TotalInsanity();
827 
828  sprintf( chLine, "Observed through a %s with aperture covering factor %.1f%%.",
829  chAper, geometry.covaper*100. );
830  PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
831  }
832  }
833  else
834  {
835  char chLine[INPUT_LINE_LENGTH];
836  sprintf( chLine, "Intensity (erg/s/cm^2)." );
837  PrintCenterLine( ioQQQ, chLine, sizeof(chLine), 130 );
838  }
839 
840  fprintf( ioQQQ, "\n" );
841 
842  /******************************************************************
843  * kill Hummer & Storey case b predictions if outside their table *
844  ******************************************************************/
845 
846  /* >>chng 02 aug 29, from lgHCaseBOK to following - caught by Ryan Porter */
847  if( !atmdat.lgHCaseBOK[1][ipHYDROGEN] )
848  {
849  static const int NWLH = 21;
850  /* list of all H case b lines */
851  realnum wlh[NWLH] = { 6562.85e0f, 4861.36e0f, 4340.49e0f, 4101.76e0f, 1.87511e4f, 1.28181e4f,
852  1.09381e4f, 1.00494e4f, 2.62515e4f, 2.16553e4f, 1.94456e4f, 7.45781e4f,
853  4.6525e4f, 3.73953e4f, 4.05116e4f, 7.45781e4f, 3.29609e4f, 1215.68e0f,
854  1025.73e0f, 972.543e0f, 949.749e0f };
855 
856  /* table exceeded - kill all H case b predictions*/
857  for( i=0; i < LineSave.nsum; i++ )
858  {
859  /* >>chng 04 jun 21, kill both case a and b at same time,
860  * actually lgHCaseBOK has separate A and B flags, but
861  * this is simpler */
862  if( LineSave.lines[i].isCaseB() ||
863  LineSave.lines[i].isCaseA() )
864  {
865  realnum errorwave;
866  /* this logic must be kept parallel with which H lines are
867  * added as case B in lines_hydro - any automatic hosing of
868  * case b would kill both H I and He II */
869  errorwave = WavlenErrorGet( LineSave.lines[i].wavelength(), LineSave.sig_figs );
870  for( j=0; j<NWLH; ++j )
871  {
872  if( fabs(LineSave.lines[i].wavelength()-wlh[j] ) <= errorwave )
873  {
874  LineSave.lines[i].SumLineZero();
875  break;
876  }
877  }
878  }
879  }
880  }
881 
882  if( !atmdat.lgHCaseBOK[1][ipHELIUM] )
883  {
884  /* table exceeded - kill all He case b predictions*/
885  static const int NWLHE = 20;
886  realnum wlhe[NWLHE] = {1640.e0f, 1215.23e0f, 1085.03e0f, 1025.35e0f, 4686.01e0f, 3203.30e0f,
887  2733.46e0f, 2511.36e0f, 1.01242e4f, 6560.44e0f, 5411.80e0f, 4859.57e0f,
888  1.86377e4f, 1.16270e4f, 9345.37e0f, 8237.17e0f, 303.808e0f, 256.338e0f,
889  243.046e0f, 237.350e0f};
890  for( i=0; i < LineSave.nsum; i++ )
891  {
892  if( LineSave.lines[i].isCaseB() ||
893  LineSave.lines[i].isCaseA() )
894  {
895  realnum errorwave;
896  /* this logic must be kept parallel with which H lines are
897  * added as case B in lines_hydro - any automatic hosing of
898  * case b would kill both H I and He II */
899  errorwave = WavlenErrorGet( LineSave.lines[i].wavelength(), LineSave.sig_figs );
900  for( j=0; j<NWLHE; ++j )
901  {
902  if( fabs(LineSave.lines[i].wavelength()-wlhe[j] ) <= errorwave )
903  {
904  LineSave.lines[i].SumLineZero();
905  break;
906  }
907  }
908  }
909  }
910  }
911 
912  /**********************************************************
913  *analyse line arrays for abundances, temp, etc *
914  **********************************************************/
915 
916  /* find apparent helium abundance, will not find these if helium is turned off */
917 
918  if( cdLine("H 1",wlAirVac(4861.36),&hbeta,&absint)<=0 )
919  {
920  if( dense.lgElmtOn[ipHYDROGEN] )
921  {
922  /* this is a major logical error if hydrogen is turned on */
923  fprintf( ioQQQ, " PrtFinal could not find H 1 4861.36 with cdLine.\n" );
925  }
926  else
927  {
928  hbeta = 0;
929  absint = -37.;
930  }
931  }
932 
933  if( dense.lgElmtOn[ipHELIUM] )
934  {
935  /* get HeI 5876 */
936  /* >>chng 06 may 15, changed this so that it works for up to six sig figs. */
937  if( cdLine("He 1",wlAirVac(5875.64),&he5876,&absint)<=0 )
938  {
939  /* 06 aug 28, from numLevels_max to _local. */
940  if( iso_sp[ipHE_LIKE][ipHELIUM].numLevels_local >= 14 )
941  {
942  /* this is a major logical error if helium is turned on */
943  fprintf( ioQQQ, " PrtFinal could not find He 1 5876 with cdLine.\n" );
945  }
946  }
947 
948  /* get HeII 4686 */
949  /* >>chng 06 may 15, changed this so that it works for up to six sig figs. */
950  if( cdLine("He 2",wlAirVac(4686.01),&he4686,&absint)<=0 )
951  {
952  /* 06 aug 28, from numLevels_max to _local. */
953  if( iso_sp[ipH_LIKE][ipHELIUM].numLevels_local > 5 )
954  {
955  /* this is a major logical error if helium is turned on
956  * and the model atom has enough levels */
957  fprintf( ioQQQ, " PrtFinal could not find He 2 4686 with cdLine.\n" );
959  }
960  }
961  }
962  else
963  {
964  he5876 = 0.;
965  absint = 0.;
966  he4686 = 0.;
967  }
968 
969  if( hbeta > 0. )
970  {
971  heabun = (he4686*0.078 + he5876*0.739)/hbeta;
972  }
973  else
974  {
975  heabun = 0.;
976  }
977 
978  if( dense.lgElmtOn[ipHELIUM] )
979  {
980  hescal = heabun/(dense.gas_phase[ipHELIUM]/dense.gas_phase[ipHYDROGEN]);
981  }
982  else
983  {
984  hescal = 0.;
985  }
986 
987  /* get temperature from [OIII] spectrum, o may be turned off,
988  * but lines still dumped into stack */
990  {
991 
992  if( cdLine("blnd",5007,&o5007,&absint)<=0 )
993  {
994  if( dense.lgElmtOn[ipOXYGEN] )
995  {
996  /* this is a major logical error if hydrogen is turned on */
997  fprintf( ioQQQ, " PrtFinal could not find O 3 5007 with cdLine.\n" );
999  }
1000  }
1001 
1002  if( cdLine("blnd",4363,&o4363,&absint)<=0 )
1003  {
1004  if( dense.lgElmtOn[ipOXYGEN] )
1005  {
1006  /* this is a major logical error if hydrogen is turned on */
1007  fprintf( ioQQQ, " PrtFinal could not find H 1 4363 with cdLine.\n" );
1009  }
1010  }
1011 
1012  // Energy ratio of O3 lines -- better to get this from transition structures
1013  realnum o3enro = (realnum)(4.56e-12/3.98e-12);
1014  /* first find low density limit OIII zone temp */
1015  if( (o4363 != 0.) && (o5007 != 0.) )
1016  {
1017  /* following will assume coll excitation only, so correct
1018  * for 4363's that cascade as 5007 */
1019  bot = o5007 - o4363/o3enro;
1020 
1021  if( bot > 0. )
1022  {
1023  ratio = o4363/bot;
1024  }
1025  else
1026  {
1027  ratio = 0.178;
1028  }
1029 
1030  if( ratio > 0.177 )
1031  {
1032  /* ratio was above infinite temperature limit */
1033  peimbt.toiiir = 1e10;
1034  }
1035  else
1036  {
1037  /* o iii 5007+4959, As 96 NIST */
1038  /*The cs of the transitions 3P0,1,2 to 1D2 are added together to give oiii_cs3P1D2 */
1039  /*the cs of the transition 1D2-1S0 is mentioned as oiii_cs1D21S0*/
1040  /*The cs of the transitions 3P0,1,2 to 1S0 are added together to give oiii_cs3P1S0*/
1041  realnum o3cs12 = 2.2347f;
1042  realnum o3cs13 = 0.29811f;
1043  double a31 = 0.2262;
1044  double a32 = 1.685;
1045  realnum o3ex23 = 32947.;
1046  realnum o3br32 = (realnum)(a32/(a31 + a32));
1047 
1048  /* data for following set in OIII cooling routine
1049  * ratio of collision strengths, factor of 4/3 makes 5007+4959
1050  * >>chng 96 sep 07, reset cs to values at 10^4K,
1051  * had been last temp in model */
1052  /*>>chng 06 jul 25, update to recent data */
1053  ratio = ratio/1.3333/(o3cs13/o3cs12);
1054  /* ratio of energies and branching ratio for 4363 */
1055  ratio = ratio/o3enro/o3br32;
1056  peimbt.toiiir = (realnum)(-o3ex23/log(ratio));
1057  }
1058  }
1059  else
1060  {
1061  peimbt.toiiir = 0.;
1062  }
1063  }
1064 
1065  if( geometry.iEmissPower == 2 )
1066  {
1067  /* find temperature from Balmer continuum */
1068  if( cdLine("Bac ",3646.,&bacobs,&absint)<=0 )
1069  {
1070  fprintf( ioQQQ, " PrtFinal could not find Bac 3646 with cdLine.\n" );
1072  }
1073 
1074  /* we pulled hbeta out of stack with cdLine above */
1075  if( hbeta > 0. )
1076  {
1077  bac = bacobs/hbeta;
1078  }
1079  else
1080  {
1081  bac = 0.;
1082  }
1083  }
1084  else
1085  {
1086  bac = 0.;
1087  }
1088 
1089  if( bac > 0. )
1090  {
1091  /*----------------------------------------------------------*
1092  ***** TableCurve c:\tcwin2\balcon.for Sep 6, 1994 11:13:38 AM
1093  ***** log bal/Hbet
1094  ***** X= log temp
1095  ***** Y=
1096  ***** Eqn# 6503 y=a+b/x+c/x^2+d/x^3+e/x^4+f/x^5
1097  ***** r2=0.9999987190883581
1098  ***** r2adj=0.9999910336185065
1099  ***** StdErr=0.001705886369042427
1100  ***** Fval=312277.1895753243
1101  ***** a= -4.82796940090671
1102  ***** b= 33.08493347410885
1103  ***** c= -56.08886262205931
1104  ***** d= 52.44759454803217
1105  ***** e= -25.07958990094209
1106  ***** f= 4.815046760060006
1107  *----------------------------------------------------------*
1108  * bac is double precision!!!! */
1109  x = 1.e0/log10(bac);
1110  tel = -4.827969400906710 + x*(33.08493347410885 + x*(-56.08886262205931 +
1111  x*(52.44759454803217 + x*(-25.07958990094209 + x*(4.815046760060006)))));
1112 
1113  if( tel > 1. && tel < 5. )
1114  {
1115  peimbt.tbac = (realnum)exp10(tel);
1116  }
1117  else
1118  {
1119  peimbt.tbac = 0.;
1120  }
1121  }
1122  else
1123  {
1124  peimbt.tbac = 0.;
1125  }
1126 
1127  if( geometry.iEmissPower == 2 )
1128  {
1129  /* find temperature from optically thin Balmer continuum and case B H-beta */
1130  if( cdLine("thin",3646.,&bacthn,&absint)<=0 )
1131  {
1132  fprintf( ioQQQ, " PrtFinal could not find thin 3646 with cdLine.\n" );
1134  }
1135 
1136  /* >>chng 06 may 15, changed this so that it works for up to six sig figs. */
1137  if( cdLine("Ca B",wlAirVac(4861.36),&hbcab,&absint)<=0 )
1138  {
1139  fprintf( ioQQQ, " PrtFinal could not find Ca B 4861.36 with cdLine.\n" );
1141  }
1142 
1143  if( hbcab > 0. )
1144  {
1145  bacthn /= hbcab;
1146  }
1147  else
1148  {
1149  bacthn = 0.;
1150  }
1151  }
1152  else
1153  {
1154  bacthn = 0.;
1155  }
1156 
1157  if( bacthn > 0. )
1158  {
1159  /*----------------------------------------------------------*
1160  ***** TableCurve c:\tcwin2\balcon.for Sep 6, 1994 11:13:38 AM
1161  ***** log bal/Hbet
1162  ***** X= log temp
1163  ***** Y=
1164  ***** Eqn# 6503 y=a+b/x+c/x^2+d/x^3+e/x^4+f/x^5
1165  ***** r2=0.9999987190883581
1166  ***** r2adj=0.9999910336185065
1167  ***** StdErr=0.001705886369042427
1168  ***** Fval=312277.1895753243
1169  ***** a= -4.82796940090671
1170  ***** b= 33.08493347410885
1171  ***** c= -56.08886262205931
1172  ***** d= 52.44759454803217
1173  ***** e= -25.07958990094209
1174  ***** f= 4.815046760060006
1175  *----------------------------------------------------------*
1176  * bac is double precision!!!! */
1177  x = 1.e0/log10(bacthn);
1178  tel = -4.827969400906710 + x*(33.08493347410885 + x*(-56.08886262205931 +
1179  x*(52.44759454803217 + x*(-25.07958990094209 + x*(4.815046760060006)))));
1180 
1181  if( tel > 1. && tel < 5. )
1182  {
1183  peimbt.tbcthn = (realnum)exp10(tel);
1184  }
1185  else
1186  {
1187  peimbt.tbcthn = 0.;
1188  }
1189  }
1190  else
1191  {
1192  peimbt.tbcthn = 0.;
1193  }
1194 
1195  /* we now have temps from OIII ratio and BAC ratio, now
1196  * do Peimbert analysis, getting To and t^2 */
1197 
1198  peimbt.tohyox = (realnum)((8.5*peimbt.toiiir - 7.5*peimbt.tbcthn - 228200. +
1199  sqrt(POW2(8.5*peimbt.toiiir-7.5*peimbt.tbcthn-228200.)+9.128e5*
1200  peimbt.tbcthn))/2.);
1201 
1202  if( peimbt.tohyox > 0. )
1203  {
1205  }
1206  else
1207  {
1208  peimbt.t2hyox = 0.;
1209  }
1210 
1211  if( prt.lgPrintBlock )
1212  PrintSpectrum();
1213 
1214  // don't print this text twice...
1215  if( !prt.lgPrtCitations )
1216  {
1217  fprintf( ioQQQ, "\n" );
1219  }
1220 
1221  /* print out ionization parameters and radiation making it through */
1222  qlow = 0.;
1223  plow = 0.;
1224  for( i=0; i < (iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon - 1); i++ )
1225  {
1226  /* N.B. in following en1ryd prevents overflow */
1227  plow += (rfield.flux[0][i] + rfield.ConInterOut[i]+ rfield.outlin[0][i] + rfield.outlin_noplot[i])*
1228  EN1RYD*rfield.anu(i);
1229  qlow += rfield.flux[0][i] + rfield.ConInterOut[i]+ rfield.outlin[0][i] + rfield.outlin_noplot[i];
1230  }
1231 
1232  qlow *= radius.r1r0sq;
1233  plow *= radius.r1r0sq;
1234  if( plow > 0. )
1235  {
1236  qlow = log10(qlow * radius.Conv2PrtInten);
1237  plow = log10(plow * radius.Conv2PrtInten);
1238  }
1239  else
1240  {
1241  qlow = 0.;
1242  plow = 0.;
1243  }
1244 
1245  qion = 0.;
1246  pion = 0.;
1247  for( i=iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1; i < rfield.nflux; i++ )
1248  {
1249  /* N.B. in following en1ryd prevents overflow */
1250  pion += (rfield.flux[0][i] + rfield.ConInterOut[i]+ rfield.outlin[0][i] + rfield.outlin_noplot[i])*
1251  EN1RYD*rfield.anu(i);
1252  qion += rfield.flux[0][i] + rfield.ConInterOut[i]+ rfield.outlin[0][i] + rfield.outlin_noplot[i];
1253  }
1254 
1255  qion *= radius.r1r0sq;
1256  pion *= radius.r1r0sq;
1257 
1258  if( pion > 0. )
1259  {
1260  qion = log10(qion * radius.Conv2PrtInten);
1261  pion = log10(pion * radius.Conv2PrtInten);
1262  }
1263  else
1264  {
1265  qion = 0.;
1266  pion = 0.;
1267  }
1268 
1269  /* derive ionization parameter for spherical geometry */
1270  if( rfield.qhtot > 0. )
1271  {
1272  if( rfield.lgUSphON )
1273  {
1274  /* RSTROM is stromgren radius */
1276  2.998e10/dense.gas_phase[ipHYDROGEN];
1277  usp = log10(usp);
1278  }
1279  else
1280  {
1281  /* no stromgren radius found, use outer radius */
1283  usp = log10(usp);
1284  }
1285  }
1286 
1287  else
1288  {
1289  usp = -37.;
1290  }
1291 
1292  if( rfield.uh > 0. )
1293  {
1294  uhl = log10(rfield.uh);
1295  }
1296  else
1297  {
1298  uhl = -37.;
1299  }
1300 
1301  if( rfield.uheii > 0. )
1302  {
1303  uhel = log10(rfield.uheii);
1304  }
1305  else
1306  {
1307  uhel = -37.;
1308  }
1309 
1310  fprintf( ioQQQ,
1311  "\n IONIZE PARMET: U(1-)%8.4f U(4-):%8.4f U(sp):%6.2f "
1312  "Q(ion): %7.3f L(ion)%7.3f Q(low):%7.3f L(low)%7.3f\n",
1313  uhl, uhel, usp, qion, pion, qlow, plow );
1314 
1315  a = safe_div(fabs((thermal.power-thermal.totcol)*100.),thermal.power,0.0);
1316  /* output power and total cooling; can be neg for shocks, collisional ioniz */
1317  if( thermal.power > 0. )
1318  {
1319  powerl = log10(thermal.power * radius.Conv2PrtInten);
1320  }
1321  else
1322  {
1323  powerl = 0.;
1324  }
1325 
1326  if( thermal.totcol > 0. )
1327  {
1329  }
1330  else
1331  {
1332  thermal.totcol = 0.;
1333  }
1334 
1335  if( thermal.FreeFreeTotHeat > 0. )
1336  {
1338  }
1339  else
1340  {
1341  thermal.FreeFreeTotHeat = 0.;
1342  }
1343 
1344  /* following is recombination line intensity */
1345  reclin = totlin('r');
1346  if( reclin > 0. )
1347  {
1348  reclin = log10(reclin * radius.Conv2PrtInten);
1349  }
1350  else
1351  {
1352  reclin = 0.;
1353  }
1354 
1355  fprintf( ioQQQ,
1356  " ENERGY BUDGET: Heat:%8.3f Coolg:%8.3f Error:%5.1f%% Rec Lin:%8.3f F-F H%7.3f P(rad/tot)max ",
1357  powerl,
1358  thermal.totcol,
1359  a,
1360  reclin,
1363  // resolving power in fine continuum for this run
1364  fprintf(ioQQQ," R(F Con):%.3e",
1365  SPEEDLIGHT / rfield.fine_opac_velocity_width );
1366 
1367  fprintf( ioQQQ, "\n");
1368 
1369  /* effective x-ray column density, from 0.5keV attenuation, no scat
1370  * IPXRY is pointer to 73.5 Ryd */
1371  coleff = opac.TauAbsGeo[0][rt.ipxry-1] - rt.tauxry;
1372  coleff /= 2.14e-22;
1373 
1374  /* find t^2 from H II structure */
1375  gett2(&peimbt.t2hstr);
1376 
1377  /* find t^2 from OIII structure */
1379 
1380  fprintf( ioQQQ, "\n Col(Heff): ");
1381  PrintE93(ioQQQ, coleff);
1382  fprintf( ioQQQ, " snd travl time ");
1384  fprintf( ioQQQ, " sec Te-low: ");
1386  fprintf( ioQQQ, " Te-hi: ");
1388 
1389  /* this is the intensity of the UV continuum at the illuminated face, relative to the Habing value, as
1390  * defined by Tielens & Hollenbach 1985 */
1391  fprintf( ioQQQ, " G0TH85: ");
1393  /* this is the intensity of the UV continuum at the illuminated face, relative to the Habing value, as
1394  * defined by Draine & Bertoldi 1985 */
1395  fprintf( ioQQQ, " G0DB96:");
1397 
1398  fprintf( ioQQQ, "\n");
1399 
1400  fprintf( ioQQQ, " Emiss Measure n(e)n(p) dl ");
1402  fprintf( ioQQQ, " n(e)n(He+)dl ");
1404  fprintf( ioQQQ, " En(e)n(He++) dl ");
1406  fprintf( ioQQQ, " ne nC+:");
1408  fprintf( ioQQQ, "\n");
1409 
1410  /* following is wl where gas goes thick to bremsstrahlung, in cm */
1411  if( rfield.EnergyBremsThin > 0. )
1412  {
1413  bremtk = RYDLAM*1e-8/rfield.EnergyBremsThin;
1414  }
1415  else
1416  {
1417  bremtk = 1e30;
1418  }
1419 
1420  /* apparent helium abundance */
1421  fprintf( ioQQQ, " He/Ha:");
1422  PrintE82( ioQQQ, heabun);
1423 
1424  /* he/h relative to correct helium abundance */
1425  fprintf(ioQQQ, " =%7.2f*true Lthin:",hescal);
1426 
1427  /* wavelength were structure is optically thick to bremsstrahlung absorption */
1428  PrintE82( ioQQQ, bremtk);
1429 
1430  /* this is ratio of conv.nTotalIoniz, the number of times ConvBase
1431  * was called during the model, over the number of zones.
1432  * This is a measure of the convergence of the model -
1433  * includes initial search so worse for short calculations.
1434  * It is a measure of how hard the model was to converge */
1435  if( nzone > 0 )
1436  {
1437  /* >>chng 07 feb 23, subtract n calls to do initial solution
1438  * so this is the number of calls needed to converge the zones.
1439  * different is to discount careful approach to molecular solutions
1440  * in one zone models */
1441  znit = (double)(conv.nTotalIoniz-conv.nTotalIoniz_start)/(double)(nzone);
1442  }
1443  else
1444  {
1445  znit = 0.;
1446  }
1447  /* >>chng 02 jan 09, format from 5.3f to 5.2f */
1448  fprintf(ioQQQ, " itr/zn:%5.2f",znit);
1449 
1450  /* sort-of same thing for large H2 molecule - number is number of level pop solutions per zone */
1451  fprintf(ioQQQ, " H2 itr/zn:%6.2f",h2.H2_itrzn());
1452 
1453  /* say whether we used stored opacities (T) or derived them from scratch (F) */
1454  fprintf(ioQQQ, " File Opacity: F" );
1455 
1456  /* log of total mass in grams if spherical, or gm/cm2 if plane parallel */
1457  /* this is mass per unit inner area */
1458  double xmass;
1459  // if spherical change to total mass, if pp leave gm/cm2
1460  if( radius.pirsq > 0. )
1461  {
1462  chUnit = "(gm)";
1463  xmass = dense.xMassTotal * exp10( (double)radius.pirsq ) ;
1464  }
1465  else
1466  {
1467  chUnit = "(gm/cm^2)";
1468  xmass = dense.xMassTotal;
1469  }
1470  fprintf(ioQQQ," MassTot %.2e %s", xmass, chUnit );
1471  fprintf(ioQQQ,"\n");
1472 
1473  /* this block is a series of prints dealing with 21 cm quantities
1474  * first this is the temperature derived from Lya - 21 cm optical depths
1475  * get harmonic mean HI temperature, as needed for 21 cm spin temperature */
1476  if( cdTemp( "opti",0,&THI,"volume" ) )
1477  {
1478  fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting 21cm opt.\n");
1479  TotalInsanity();
1480  }
1481  fprintf( ioQQQ, " Temps(21 cm) T(21cm/Ly a) ");
1482  PrintE82(ioQQQ, THI );
1483 
1484  /* get harmonic mean HI gas kin temperature, as needed for 21 cm spin temperature
1485  * >>chng 06 jul 21, this was over volume but hazy said radius - change to radius,
1486  * bug discovered by Eric Pellegrini & Jack Baldwin */
1487  /*if( cdTemp( "21cm",0,&THI,"volume" ) )*/
1488  if( cdTemp( "21cm",0,&THI,"radius" ) )
1489  {
1490  fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting 21cm temp.\n");
1491  TotalInsanity();
1492  }
1493  fprintf(ioQQQ, " T(<nH/Tkin>) ");
1494  PrintE82(ioQQQ, THI);
1495 
1496  /* get harmonic mean HI 21 21 spin temperature, as needed for 21 cm spin temperature
1497  * for this, always weighted by radius, volume would be ignored were it present */
1498  if( cdTemp( "spin",0,&THI,"radius" ) )
1499  {
1500  fprintf(ioQQQ,"\n>>>>>>>>>>>>>>>>> PrtFinal, impossible problem getting 21cm spin.\n");
1501  TotalInsanity();
1502  }
1503  fprintf(ioQQQ, " T(<nH/Tspin>) ");
1504  PrintE82(ioQQQ, THI);
1505 
1506  /* now convert this HI spin temperature into a brightness temperature */
1507  THI *= (1. - sexp( HFLines[0].Emis().TauIn() ) );
1508  fprintf( ioQQQ, " TB21cm:");
1509  PrintE82(ioQQQ, THI);
1510  fprintf( ioQQQ, "\n");
1511 
1512  fprintf( ioQQQ, " N(H0/Tspin) ");
1514  fprintf( ioQQQ, " N(OH/Tkin) ");
1516 
1517  /* mean B weighted wrt 21 cm absorption */
1518  bot = cdB21cm();
1519  fprintf( ioQQQ, " B(21cm) ");
1520  PrintE82(ioQQQ, bot );
1521 
1522  /* end prints for 21 cm */
1523  fprintf(ioQQQ, "\n");
1524 
1525  /* timescale (sec here) for photoerosion of Fe-Ni */
1526  rate = timesc.TimeErode*2e-26;
1527  if( rate > SMALLFLOAT )
1528  {
1529  ferode = 1./rate;
1530  }
1531  else
1532  {
1533  ferode = 0.;
1534  }
1535 
1536  /* mean acceleration */
1537  if( wind.acldr > 0. )
1538  {
1539  wind.AccelAver /= wind.acldr;
1540  }
1541  else
1542  {
1543  wind.AccelAver = 0.;
1544  }
1545 
1546  /* DMEAN is mean density (gm per cc); mean temp is weighted by mass density */
1547  /* >>chng 02 aug 21, from radius.depth_x_fillfac to integral of radius times fillfac */
1549  tmean = colden.tmas/SDIV(colden.TotMassColl);
1550  /* mean mass per particle */
1551  wmean = colden.wmas/SDIV(colden.TotMassColl);
1552 
1553  fprintf( ioQQQ, " <a>:");
1555  fprintf( ioQQQ, " erdeFe");
1556  PrintE71(ioQQQ , ferode);
1557  fprintf( ioQQQ, " Tcompt");
1559  fprintf( ioQQQ, " Tthr");
1561  fprintf( ioQQQ, " <Tden>: ");
1562  PrintE82(ioQQQ , tmean);
1563  fprintf( ioQQQ, " <dens>:");
1564  PrintE82(ioQQQ , dmean);
1565  fprintf( ioQQQ, " <MolWgt>");
1566  PrintE82(ioQQQ , wmean);
1567  fprintf(ioQQQ,"\n");
1568 
1569  /* log of Jeans length and mass - this is mean over model */
1570  if( tmean > 0. )
1571  {
1572  rjeans = (log10(JEANS)+log10(tmean) - log10(dense.wmole) - log10(dense.xMassDensity*
1573  geometry.FillFac))/2.;
1574  }
1575  else
1576  {
1577  /* tmean undefined, set rjeans to large value so 0 printed below */
1578  rjeans = 40.f;
1579  }
1580 
1581  if( rjeans < 36. )
1582  {
1583  rjeans = exp10(rjeans);
1584  /* AJMASS is Jeans mass in solar units */
1585  ajmass = 3.*log10(rjeans/2.) + log10(4.*PI/3.*dense.xMassDensity*
1586  geometry.FillFac) - log10(SOLAR_MASS);
1587  if( ajmass < 37 )
1588  {
1589  ajmass = exp10(ajmass);
1590  }
1591  else
1592  {
1593  ajmass = 0.;
1594  }
1595  }
1596  else
1597  {
1598  rjeans = 0.;
1599  ajmass = 0.;
1600  }
1601 
1602  /* Jeans length and mass - this is smallest over model */
1603  ajmin = colden.ajmmin - log10(SOLAR_MASS);
1604  if( ajmin < 37 )
1605  {
1606  ajmin = exp10(ajmin);
1607  }
1608  else
1609  {
1610  ajmin = 0.;
1611  }
1612 
1613  /* estimate alpha (o-x) */
1614  if( rfield.anu(rfield.nflux-1) > 150. )
1615  {
1616  /* generate pointers to energies where alpha ox will be evaluated */
1617  ip2kev = ipoint(147.);
1618  ip2500 = ipoint(0.3648);
1619 
1620  /* now get fluxes there */
1621  bottom = rfield.flux[0][ip2500-1]*
1622  rfield.anu(ip2500-1)/rfield.widflx(ip2500-1);
1623 
1624  top = rfield.flux[0][ip2kev-1]*
1625  rfield.anu(ip2kev-1)/rfield.widflx(ip2kev-1);
1626 
1627  /* generate alpha ox if denominator is non-zero */
1628  if( bottom > 1e-30 && top > 1e-30 )
1629  {
1630  ratio = log10(top) - log10(bottom);
1631  if( ratio < 36. && ratio > -36. )
1632  {
1633  ratio = exp10(ratio);
1634  }
1635  else
1636  {
1637  ratio = 0.;
1638  }
1639  }
1640 
1641  else
1642  {
1643  ratio = 0.;
1644  }
1645 
1646  if( ratio > 0. )
1647  {
1648  // the separate variable freq_ratio is needed to work around a bug in icc 10.0
1649  double freq_ratio = rfield.anu(ip2kev-1)/rfield.anu(ip2500-1);
1650  alfox = log(ratio)/log(freq_ratio);
1651  }
1652  else
1653  {
1654  alfox = 0.;
1655  }
1656  }
1657  else
1658  {
1659  alfox = 0.;
1660  }
1661 
1662  if( colden.rjnmin < 37 )
1663  {
1665  }
1666  else
1667  {
1668  colden.rjnmin = FLT_MAX;
1669  }
1670 
1671  fprintf( ioQQQ, " Mean Jeans l(cm)");
1672  PrintE82(ioQQQ, rjeans);
1673  fprintf( ioQQQ, " M(sun)");
1674  PrintE82(ioQQQ, ajmass);
1675  fprintf( ioQQQ, " smallest: len(cm):");
1677  fprintf( ioQQQ, " M(sun):");
1678  PrintE82(ioQQQ, ajmin);
1679  fprintf( ioQQQ, " a_ox tran:%6.2f\n" , alfox);
1680 
1681  fprintf( ioQQQ, " Rion:");
1682  double R_ion;
1683  if( rfield.lgUSphON )
1684  R_ion = rfield.rstrom;
1685  else
1686  R_ion = radius.Radius;
1687  PrintE93(ioQQQ, R_ion);
1688  fprintf( ioQQQ, " Dist:");
1690  fprintf( ioQQQ, " Diam:");
1691  if( radius.distance > 0. )
1692  PrintE93(ioQQQ, 2.*R_ion*AS1RAD/radius.distance);
1693  else
1694  PrintE93(ioQQQ, 0.);
1695  fprintf( ioQQQ, "\n");
1696 
1697  if( prt.lgPrintTime )
1698  {
1699  /* print execution time by default */
1700  alfox = cdExecTime();
1701  }
1702  else
1703  {
1704  /* flag set false with no time command, so that different runs can
1705  * compare exactly */
1706  alfox = 0.;
1707  }
1708 
1709  /* some details about the hydrogen and helium ions */
1710  fprintf( ioQQQ,
1711  " Hatom level%3ld HatomType:%4.4s HInducImp %2c"
1712  " He tot level:%3ld He2 level: %3ld ExecTime%8.1f\n",
1713  /* 06 aug 28, from numLevels_max to _local. */
1714  iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local,
1715  hydro.chHTopType,
1717  /* 06 aug 28, from numLevels_max to _local. */
1719  /* 06 aug 28, from numLevels_max to _local. */
1720  iso_sp[ipH_LIKE][ipHELIUM].numLevels_local ,
1721  alfox );
1722 
1723  /* now give an indication of the convergence error budget */
1724  fprintf( ioQQQ,
1725  " ConvrgError(%%) <eden>%7.3f MaxEden%7.3f <H-C>%7.2f Max(H-C)%8.2f <Press>%8.3f MaxPrs er%7.3f\n",
1726  conv.AverEdenError/nzone*100. ,
1727  conv.BigEdenError*100. ,
1728  conv.AverHeatCoolError/nzone*100. ,
1729  conv.BigHeatCoolError*100. ,
1730  conv.AverPressError/nzone*100. ,
1731  conv.BigPressError*100. );
1732 
1733  fprintf(ioQQQ,
1734  " Continuity(%%) chng Te%6.1f elec den%6.1f n(H2)%7.1f n(CO) %7.1f\n",
1735  phycon.BigJumpTe*100. ,
1736  phycon.BigJumpne*100. ,
1737  phycon.BigJumpH2*100. ,
1738  phycon.BigJumpCO*100. );
1739 
1740  /* print out some average quantities */
1741  fprintf( ioQQQ, "\n Averaged Quantities\n" );
1742  fprintf( ioQQQ, " Te Te(Ne) Te(NeNp) Te(NeHe+)Te(NeHe2+) Te(NeO+) Te(NeO2+)"
1743  " Te(H2) N(H) Ne(O2+) Ne(Np)\n" );
1744  static const char* weight[3] = { "Radius", "Area", "Volume" };
1745  int dmax = geometry.lgGeoPP ? 1 : 3;
1746  for( int dim=0; dim < dmax; ++dim )
1747  {
1748  fprintf( ioQQQ, " %6s:", weight[dim] );
1749  // <Te>/<1>
1750  PrintRatio( mean.TempMean[dim][0], mean.TempMean[dim][1] );
1751  // <Te*ne>/<ne>
1752  PrintRatio( mean.TempEdenMean[dim][0], mean.TempEdenMean[dim][1] );
1753  // <Te*ne*nion>/<ne*nion>
1754  PrintRatio( mean.TempIonEdenMean[dim][ipHYDROGEN][1][0], mean.TempIonEdenMean[dim][ipHYDROGEN][1][1] );
1755  PrintRatio( mean.TempIonEdenMean[dim][ipHELIUM][1][0], mean.TempIonEdenMean[dim][ipHELIUM][1][1] );
1756  PrintRatio( mean.TempIonEdenMean[dim][ipHELIUM][2][0], mean.TempIonEdenMean[dim][ipHELIUM][2][1] );
1757  PrintRatio( mean.TempIonEdenMean[dim][ipOXYGEN][1][0], mean.TempIonEdenMean[dim][ipOXYGEN][1][1] );
1758  PrintRatio( mean.TempIonEdenMean[dim][ipOXYGEN][2][0], mean.TempIonEdenMean[dim][ipOXYGEN][2][1] );
1759  // <Te*nH2>/<nH2>
1760  PrintRatio( mean.TempIonMean[dim][ipHYDROGEN][2][0], mean.TempIonMean[dim][ipHYDROGEN][2][1] );
1761  // <nH>/<1>
1762  PrintRatio( mean.xIonMean[dim][ipHYDROGEN][0][1], mean.TempMean[dim][1] );
1763  // <ne*nion>/<nion>
1764  PrintRatio( mean.TempIonEdenMean[dim][ipOXYGEN][2][1], mean.TempIonMean[dim][ipOXYGEN][2][1] );
1765  PrintRatio( mean.TempIonEdenMean[dim][ipHYDROGEN][1][1], mean.TempIonMean[dim][ipHYDROGEN][1][1] );
1766  fprintf( ioQQQ, "\n" );
1767  }
1768 
1769  /* print out Peimbert analysis, tsqden default 1e7, changed
1770  * with set tsqden command */
1772  {
1773  fprintf( ioQQQ, " \n" );
1774 
1775  /* temperature from the [OIII] 5007/4363 ratio */
1776  fprintf( ioQQQ, " Peimbert T(OIIIr)");
1778 
1779  /* temperature from predicted Balmer jump relative to Hbeta */
1780  fprintf( ioQQQ, " T(Bac)");
1781  PrintE82( ioQQQ , peimbt.tbac);
1782 
1783  /* temperature predicted from optically thin Balmer jump rel to Hbeta */
1784  fprintf( ioQQQ, " T(Hth)");
1786 
1787  /* t^2 predicted from the structure, weighted by H */
1788  fprintf( ioQQQ, " t2(Hstrc)");
1789  fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2hstr));
1790 
1791  /* temperature from both [OIII] and the Balmer jump rel to Hbeta */
1792  fprintf( ioQQQ, " T(O3-BAC)");
1794 
1795  /* t2 from both [OIII] and the Balmer jump rel to Hbeta */
1796  fprintf( ioQQQ, " t2(O3-BC)");
1797  fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2hyox));
1798 
1799  /* structural t2 from the O+2 predicted radial dependence */
1800  fprintf( ioQQQ, " t2(O3str)");
1801  fprintf( ioQQQ,PrintEfmt("%9.2e", peimbt.t2o3str));
1802 
1803  fprintf( ioQQQ, "\n");
1804 
1805  if( gv.lgDustOn() )
1806  {
1807  fprintf( ioQQQ, " Be careful: grains exist. This spectrum was not corrected for reddening before analysis.\n" );
1808  }
1809  }
1810 
1811  /* print average (over radius) grain dust parameters if lgDustOn() is true,
1812  * average quantities are incremented in radius_increment, zeroed in IterStart */
1813  if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
1814  {
1815  char chQHeat;
1816  double AV , AB;
1817  double total_dust2gas = 0.;
1818 
1819  fprintf( ioQQQ, "\n Average Grain Properties (over radius):\n" );
1820 
1821  for( size_t i0=0; i0 < gv.bin.size(); i0 += 10 )
1822  {
1823  if( i0 > 0 )
1824  fprintf( ioQQQ, "\n" );
1825 
1826  /* this is upper limit to how many grain species we will print across line */
1827  size_t i1 = min(i0+10,gv.bin.size());
1828 
1829  fprintf( ioQQQ, " " );
1830  for( size_t nd=i0; nd < i1; nd++ )
1831  {
1832  chQHeat = (char)(gv.bin[nd]->lgEverQHeat ? '*' : ' ');
1833  fprintf( ioQQQ, "%-12.12s%c", gv.bin[nd]->chDstLab, chQHeat );
1834  }
1835  fprintf( ioQQQ, "\n" );
1836 
1837  fprintf( ioQQQ, " nd:" );
1838  for( size_t nd=i0; nd < i1; nd++ )
1839  {
1840  if( nd != i0 ) fprintf( ioQQQ," " );
1841  fprintf( ioQQQ, "%7ld ", (unsigned long)nd );
1842  }
1843  fprintf( ioQQQ, "\n" );
1844 
1845  fprintf( ioQQQ, " <Tgr>:" );
1846  for( size_t nd=i0; nd < i1; nd++ )
1847  {
1848  if( nd != i0 ) fprintf( ioQQQ," " );
1849  fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdust/radius.depth_x_fillfac));
1850  }
1851  fprintf( ioQQQ, "\n" );
1852 
1853  fprintf( ioQQQ, " <Vel>:" );
1854  for( size_t nd=i0; nd < i1; nd++ )
1855  {
1856  if( nd != i0 ) fprintf( ioQQQ," " );
1857  fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdft/radius.depth_x_fillfac));
1858  }
1859  fprintf( ioQQQ, "\n" );
1860 
1861  fprintf( ioQQQ, " <Pot>:" );
1862  for( size_t nd=i0; nd < i1; nd++ )
1863  {
1864  if( nd != i0 ) fprintf( ioQQQ," " );
1865  fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avdpot/radius.depth_x_fillfac));
1866  }
1867  fprintf( ioQQQ, "\n" );
1868 
1869  fprintf( ioQQQ, " <D/G>:" );
1870  for( size_t nd=i0; nd < i1; nd++ )
1871  {
1872  if( nd != i0 ) fprintf( ioQQQ," " );
1873  fprintf( ioQQQ,PrintEfmt("%10.3e", gv.bin[nd]->avDGRatio/radius.depth_x_fillfac));
1874  /* add up total dust to gas mass ratio */
1875  total_dust2gas += gv.bin[nd]->avDGRatio/radius.depth_x_fillfac;
1876  }
1877  fprintf( ioQQQ, "\n" );
1878  }
1879 
1880  fprintf(ioQQQ," Dust to gas ratio (by mass):");
1881  fprintf(ioQQQ,PrintEfmt("%10.3e", total_dust2gas));
1882 
1883  /* total extinction (conv to mags) at V and B per hydrogen, this includes
1884  * forward scattering as an extinction process, so is what would be measured
1885  * for a star, but not for an extended source where forward scattering
1886  * should be discounted */
1889  /* print A_V/N(Htot) for point and extended sources */
1890  fprintf(ioQQQ,", A(V)/N(H)(pnt):%.3e, (ext):%.3e",
1891  AV,
1893 
1894  /* ratio of total to selective extinction */
1895  fprintf(ioQQQ,", R:");
1896 
1897  /* gray grains have AB - AV == 0 */
1898  if( fabs(AB-AV)>SMALLFLOAT )
1899  {
1900  fprintf(ioQQQ,"%.3e", AV/(AB-AV) );
1901  }
1902  else
1903  {
1904  fprintf(ioQQQ,"%.3e", 0. );
1905  }
1906  fprintf(ioQQQ," AV(ext):%.3e (pnt):%.3e\n",
1909  }
1910 
1911  /* option to make short printout */
1912  if( !prt.lgPrtShort && called.lgTalk )
1913  {
1914  /* print log of optical depths,
1915  * calls prtmet if print line optical depths command entered */
1916  PrtAllTau();
1917 
1918  /* only print mean ionization and emergent continuum on last iteration */
1919  if( iterations.lgLastIt )
1920  {
1921  /* option to print column densities, set with print column density command */
1922  if( prt.lgPrintColumns )
1923  PrtColumns(ioQQQ);
1924  /* print mean ionization fractions for all elements */
1925  PrtMeanIon('i', false, ioQQQ);
1926  /* print mean ionization fractions for all elements with density weighting*/
1927  PrtMeanIon('i', true , ioQQQ);
1928  /* print mean temperature fractions for all elements */
1929  PrtMeanIon('t', false, ioQQQ);
1930  /* print mean temperature fractions for all elements with density weighting */
1931  PrtMeanIon('t', true , ioQQQ);
1932  }
1933  }
1934 
1935  /* print input title for model */
1936  fprintf( ioQQQ, "%s\n\n", input.chTitle );
1937  fflush(ioQQQ);
1938  return;
1939 }
1940 
1941 /* routine to stuff comments into the stack of comments,
1942  * return is index to this comment */
1943 long int StuffComment( const char * chComment )
1944 {
1945  int n , i;
1946 
1947  DEBUG_ENTRY( "StuffComment()" );
1948 
1949  /* only do this one time per core load */
1950  if( LineSave.ipass == 0 )
1951  {
1953  {
1954  fprintf( ioQQQ, " Too many comments have been entered into line array; increase the value of NHOLDCOMMENTS.\n" );
1956  }
1957 
1958  strcpy( LineSave.chHoldComments[LineSave.nComment], chComment );
1959 
1960  /* current length of this string */
1961  n = (long)strlen( LineSave.chHoldComments[LineSave.nComment] );
1962  for( i=0; i<prt_linecol.column_len-n; ++i )
1963  {
1964  strcat( LineSave.chHoldComments[LineSave.nComment], ".");
1965  }
1966 
1967 #if 0
1968  for( i=0; i<6; ++i )
1969  {
1970  strcat( LineSave.chHoldComments[LineSave.nComment], " ");
1971  }
1972 #endif
1973  }
1974 
1975  ++LineSave.nComment;
1976  return( LineSave.nComment-1 );
1977 }
1978 
1979 /*gett2 analyze computed structure to get structural t^2 */
1980 STATIC void gett2(realnum *tsqr)
1981 {
1982  long int i;
1983 
1984  double tmean;
1985  double a,
1986  as,
1987  b;
1988 
1989  DEBUG_ENTRY( "gett2()" );
1990 
1991  /* get T, t^2 */
1992  a = 0.;
1993  b = 0.;
1994 
1995  ASSERT( nzone < struc.nzlim);
1996  // struc.volstr[] is radius.dVeffAper saved as a function of nzone
1997  // we need this version of radius.dVeff since we want to compare to
1998  // emission lines that react to the APERTURE command
1999  for( i=0; i < nzone; i++ )
2000  {
2001  as = (double)(struc.volstr[i])*(double)(struc.hiistr[i])*
2002  (double)(struc.ednstr[i]);
2003  a += (double)(struc.testr[i])*as;
2004  /* B is used twice below */
2005  b += as;
2006  }
2007 
2008  if( b <= 0. )
2009  {
2010  *tsqr = 0.;
2011  }
2012  else
2013  {
2014  /* following is H+ weighted mean temp over vol */
2015  tmean = a/b;
2016  a = 0.;
2017 
2018  ASSERT( nzone < struc.nzlim );
2019  for( i=0; i < nzone; i++ )
2020  {
2021  as = (double)(struc.volstr[i])*(double)(struc.hiistr[i])*
2022  struc.ednstr[i];
2023  a += (POW2((double)(struc.testr[i]-tmean)))*as;
2024  }
2025  *tsqr = (realnum)(a/(b*(POW2(tmean))));
2026  }
2027 
2028  return;
2029 }
2030 
2031 /*gett2o3 analyze computed [OIII] spectrum to get t^2 */
2033 {
2034  long int i;
2035  double tmean;
2036  double a,
2037  as,
2038  b;
2039 
2040  DEBUG_ENTRY( "gett2o3()" );
2041 
2042  /* get T, t^2 */ a = 0.;
2043  b = 0.;
2045  // struc.volstr[] is radius.dVeffAper saved as a function of nzone
2046  // we need this version of radius.dVeff since we want to compare to
2047  // emission lines that react to the APERTURE command
2048  for( i=0; i < nzone; i++ )
2049  {
2050  as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
2051  (double)(struc.ednstr[i]);
2052  a += (double)(struc.testr[i])*as;
2053 
2054  /* B is used twice below */
2055  b += as;
2056  }
2057 
2058  if( b <= 0. )
2059  {
2060  *tsqr = 0.;
2061  }
2062 
2063  else
2064  {
2065  /* following is H+ weighted mean temp over vol */
2066  tmean = a/b;
2067  a = 0.;
2068  ASSERT( nzone < struc.nzlim );
2069  for( i=0; i < nzone; i++ )
2070  {
2071  as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
2072  struc.ednstr[i];
2073  a += (POW2((double)(struc.testr[i]-tmean)))*as;
2074  }
2075  *tsqr = (realnum)(a/(b*(POW2(tmean))));
2076  }
2077  return;
2078 }
long int nSave
Definition: input.h:62
realnum toiiir
Definition: peimbt.h:18
bool lgIterAgain
Definition: iterations.h:53
realnum tbac
Definition: peimbt.h:18
realnum wlSort2
Definition: prt.h:150
void PrtFinal(void)
Definition: prt_final.cpp:555
realnum BigJumpTe
Definition: phycon.h:116
realnum t2hyox
Definition: peimbt.h:18
double Radius
Definition: radius.h:31
realnum t2hstr
Definition: peimbt.h:18
t_atmdat atmdat
Definition: atmdat.cpp:6
t_line_col prt_linecol
Definition: prt.cpp:15
t_thermal thermal
Definition: thermal.cpp:6
void PrintE93(FILE *, double)
Definition: service.cpp:1019
double power
Definition: thermal.h:169
t_colden colden
Definition: colden.cpp:5
realnum TimeErode
Definition: timesc.h:61
realnum WavlenErrorGet(realnum wavelength, long sig_figs)
double exp10(double x)
Definition: cddefines.h:1383
const int ipHE_LIKE
Definition: iso.h:65
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
double widflx(size_t i) const
Definition: mesh.h:147
t_input input
Definition: input.cpp:12
t_opac opac
Definition: opacity.cpp:5
realnum ** flux
Definition: rfield.h:70
t_struc struc
Definition: struc.cpp:6
realnum UV_Cont_rel2_Draine_DB96_face
Definition: hmi.h:84
double FreeFreeTotHeat
Definition: thermal.h:178
bool lgPrnHeat
Definition: prt.h:188
const realnum SMALLFLOAT
Definition: cpu.h:246
bool lgPrintLineCumulative
Definition: prt.h:268
realnum EnergyBremsThin
Definition: rfield.h:229
bool lgCylnOn
Definition: radius.h:126
realnum * outlin_noplot
Definition: rfield.h:191
STATIC void PrintCenterLine(FILE *io, const char chLine[], size_t ArrLen, size_t LineLen)
Definition: prt_final.cpp:45
char TorF(bool l)
Definition: cddefines.h:753
const int ipOXYGEN
Definition: cddefines.h:355
#define PrintEfmt(F, V)
Definition: cddefines.h:1364
realnum fine_opac_velocity_width
Definition: rfield.h:369
double tcmptn
Definition: timesc.h:26
long int iEmissPower
Definition: geometry.h:71
const double SMALLDOUBLE
Definition: cpu.h:250
realnum BigJumpH2
Definition: phycon.h:116
t_warnings warnings
Definition: warnings.cpp:11
long nMatch(const char *chKey, const char *chCard)
Definition: service.cpp:537
realnum BigHeatCoolError
Definition: conv.h:176
t_conv conv
Definition: conv.cpp:5
realnum * ednstr
Definition: struc.h:25
realnum tmas
Definition: colden.h:61
STATIC void gett2(realnum *tsqr)
Definition: prt_final.cpp:1980
long int nTotalIoniz_start
Definition: conv.h:164
double distance
Definition: radius.h:70
TransitionList HFLines("HFLines",&AnonStates)
t_phycon phycon
Definition: phycon.cpp:6
t_LineSave LineSave
Definition: lines.cpp:9
string relint_outrange
Definition: prt.h:325
realnum AccelAver
Definition: wind.h:46
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
Definition: cddrive.cpp:1296
sys_float sexp(sys_float x)
Definition: service.cpp:1095
realnum AverPressError
Definition: conv.h:181
char chVersion[INPUT_LINE_LENGTH]
Definition: version.h:19
realnum * volstr
Definition: struc.h:25
realnum BigPressError
Definition: conv.h:180
realnum ajmmin
Definition: colden.h:59
bool lgHInducImp
Definition: hydrogenic.h:131
realnum covgeo
Definition: geometry.h:45
bool lgHCaseBOK[2][HS_NZ]
Definition: atmdat.h:342
realnum ** outlin
Definition: rfield.h:191
FILE * ioQQQ
Definition: cddefines.cpp:7
realnum AverHeatCoolError
Definition: conv.h:177
void PrtMeanIon(char chType, bool lgDensity, FILE *)
Definition: prt_meanion.cpp:11
realnum FillFac
Definition: geometry.h:29
char chTitle[INPUT_LINE_LENGTH]
Definition: input.h:48
long int nzone
Definition: cddefines.cpp:14
bool lgTalk
Definition: called.h:12
realnum rjnmin
Definition: colden.h:59
bool lgPrtCitations
Definition: prt.h:294
t_dynamics dynamics
Definition: dynamics.cpp:42
realnum wlSort1
Definition: prt.h:150
realnum t2o3str
Definition: peimbt.h:18
realnum BigJumpne
Definition: phycon.h:116
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:807
void cap4(char *chCAP, const char *chLab)
Definition: service.cpp:264
double anu(size_t i) const
Definition: mesh.h:111
void PrtAllTau(void)
Definition: prt_alltau.cpp:15
vector< LinSv > lines
Definition: lines.h:132
bool lgSortLineIntensity
Definition: prt.h:146
t_dense dense
Definition: global.cpp:15
realnum tohyox
Definition: peimbt.h:18
static t_version & Inst()
Definition: cddefines.h:209
bool lgBadStop
Definition: conv.h:248
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
double depth_x_fillfac
Definition: radius.h:79
bool lgPrintBlockIntrinsic
Definition: prt.h:134
realnum tauxry
Definition: rt.h:184
Wind wind
Definition: wind.cpp:5
bool lgSphere
Definition: geometry.h:34
long int iteration
Definition: cddefines.cpp:16
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
realnum wmas
Definition: colden.h:61
double rinner
Definition: radius.h:31
t_geometry geometry
Definition: geometry.cpp:5
bool lgdBaseSourceExists[LIMELM][LIMELM+1]
Definition: atmdat.h:452
void PrintE71(FILE *, double)
Definition: service.cpp:969
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:15
realnum TooFaint
Definition: prt.h:244
#define POW2
Definition: cddefines.h:983
double dlnenCp
Definition: colden.h:49
const int ipH1s
Definition: iso.h:29
#define STATIC
Definition: cddefines.h:118
multi_arr< double, 2 > TempEdenMean
Definition: mean.h:38
realnum tlowst
Definition: thermal.h:68
bool lgSurfaceBrightness
Definition: prt.h:179
realnum qhtot
Definition: rfield.h:339
realnum BigJumpCO
Definition: phycon.h:116
bool lgPrintBlockEmergent
Definition: prt.h:138
t_pressure pressure
Definition: pressure.cpp:9
t_rfield rfield
Definition: rfield.cpp:9
double dlnenHep
Definition: colden.h:43
t_mean mean
Definition: mean.cpp:16
realnum * ConInterOut
Definition: rfield.h:156
STATIC void gett2o3(realnum *tsqr)
Definition: prt_final.cpp:2032
float realnum
Definition: cddefines.h:124
realnum uh
Definition: rfield.h:347
#define EXIT_FAILURE
Definition: cddefines.h:168
int relint_len
Definition: prt.h:314
double H0_ov_Tspin
Definition: colden.h:52
const realnum BIGFLOAT
Definition: cpu.h:244
double OH_ov_Tspin
Definition: colden.h:55
const int INPUT_LINE_LENGTH
Definition: cddefines.h:301
multi_arr< double, 4 > TempIonEdenMean
Definition: mean.h:26
realnum uheii
Definition: rfield.h:350
bool lgFaintOn
Definition: prt.h:245
realnum thist
Definition: thermal.h:68
bool lgSurfaceBrightness_SR
Definition: prt.h:179
t_hydro hydro
Definition: hydrogenic.cpp:5
#define cdEXIT(FAIL)
Definition: cddefines.h:484
realnum RadBetaMax
Definition: pressure.h:96
bool lgRadiusKnown
Definition: radius.h:121
double totcol
Definition: thermal.h:130
long min(int a, long b)
Definition: cddefines.h:766
bool lgGrainPhysicsOn
Definition: grainvar.h:481
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
bool lgGeoPP
Definition: geometry.h:21
realnum wlAirVac(double wlAir)
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1015
STATIC void PrintSpectrum()
Definition: prt_final.cpp:73
bool lgPrintTime
Definition: prt.h:161
t_iterations iterations
Definition: iterations.cpp:6
realnum rstrom
Definition: rfield.h:355
void PrtColumns(FILE *ioMEAN)
Definition: prt_columns.cpp:14
long int sig_figs
Definition: lines.h:109
bool lgPrtShort
Definition: prt.h:219
realnum UV_Cont_rel2_Habing_TH85_face
Definition: hmi.h:74
double cdExecTime()
Definition: cddrive.cpp:483
t_radius radius
Definition: radius.cpp:5
t_timesc timesc
Definition: timesc.cpp:7
realnum tsqden
Definition: peimbt.h:18
t_prt prt
Definition: prt.cpp:14
realnum pirsq
Definition: radius.h:148
long int nTotalIoniz
Definition: conv.h:159
bool lgPrtLineLog
Definition: prt.h:264
bool lgElmtOn[LIMELM]
Definition: dense.h:160
realnum TotMassColl
Definition: colden.h:61
double extin_mag_V_point
Definition: rfield.h:260
bool lgPredLumin
Definition: radius.h:144
string col_gap
Definition: prt.h:328
long int nzlim
Definition: struc.h:19
realnum gas_phase[LIMELM]
Definition: dense.h:76
bool lgPrintBlock
Definition: prt.h:130
double totlin(int chInfo)
long int itermx
Definition: iterations.h:37
double dlnenp
Definition: colden.h:40
double Conv2PrtInten
Definition: radius.h:152
realnum wmole
Definition: dense.h:111
bool lgSortLines
Definition: prt.h:142
#define ASSERT(exp)
Definition: cddefines.h:617
void sprt_wl(char *chString, realnum wl)
Definition: prt.cpp:56
realnum covaper
Definition: geometry.h:54
bool lgLastIt
Definition: iterations.h:47
double sound
Definition: timesc.h:40
t_peimbt peimbt
Definition: peimbt.cpp:5
double cdB21cm()
Definition: cddrive.cpp:1255
char chDenseLaw[5]
Definition: dense.h:176
long int nComment
Definition: lines.h:90
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1069
double extin_mag_B_point
Definition: rfield.h:260
double extin_mag_V_extended
Definition: rfield.h:264
long int ipNormWavL
Definition: lines.h:102
const int ipH_LIKE
Definition: iso.h:64
multi_arr< double, 4 > TempIonMean
Definition: mean.h:24
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
const int ipHELIUM
Definition: cddefines.h:349
void PrintRatio(double q1, double q2)
Definition: prt_final.cpp:65
realnum xMassDensity
Definition: dense.h:101
static vector< realnum > wavelength
vector< GrainBin * > bin
Definition: grainvar.h:585
realnum AverEdenError
Definition: conv.h:173
multi_arr< double, 2 > TempMean
Definition: mean.h:36
bool lgPrnPump
Definition: prt.h:188
void DatabasePrintReference()
Definition: service.cpp:1907
char chHTopType[5]
Definition: hydrogenic.h:116
realnum ** TauAbsGeo
Definition: opacity.h:90
static const int NHOLDCOMMENTS
Definition: lines.h:72
bool lgPrintFluxEarth
Definition: prt.h:175
realnum xMassTotal
Definition: dense.h:117
#define MAX2(a, b)
Definition: cddefines.h:828
bool lgPrintColumns
Definition: prt.h:157
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
multi_arr< double, 4 > xIonMean
Definition: mean.h:17
int absint_len
Definition: prt.h:311
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
long int n_initial_relax
Definition: dynamics.h:132
double dlnenHepp
Definition: colden.h:46
long int nsum
Definition: lines.h:87
realnum * testr
Definition: struc.h:25
bool lgSizeSet
Definition: geometry.h:80
void caps(char *chCard)
Definition: service.cpp:304
char chCardSav[NKRD][INPUT_LINE_LENGTH]
Definition: input.h:48
GrainVar gv
Definition: grainvar.cpp:5
t_hmi hmi
Definition: hmi.cpp:5
double r1r0sq
Definition: radius.h:31
bool lgCautns
Definition: warnings.h:44
long int ipxry
Definition: rt.h:181
bool lgUSphON
Definition: rfield.h:353
char chHoldComments[NHOLDCOMMENTS][INPUT_LINE_LENGTH]
Definition: lines.h:99
void PrintE82(FILE *, double)
Definition: service.cpp:920
char chCumuType[5]
Definition: rfield.h:333
double time_therm_long
Definition: timesc.h:29
bool lgPrnColl
Definition: prt.h:188
bool lgDustOn() const
Definition: grainvar.h:477
const int ipHYDROGEN
Definition: cddefines.h:348
realnum BigEdenError
Definition: conv.h:215
long int nflux
Definition: rfield.h:48
realnum * hiistr
Definition: struc.h:25
realnum acldr
Definition: wind.h:46
realnum colden[NCOLD]
Definition: colden.h:32
realnum tbcthn
Definition: peimbt.h:18
double H2_itrzn(void)
Definition: mole_h2.cpp:263
long int numLevels_local
Definition: iso.h:529
bool lgWarngs
Definition: warnings.h:44
long int ipass
Definition: lines.h:96
t_called called
Definition: called.cpp:4
realnum * o3str
Definition: struc.h:25
bool lgPrnInwd
Definition: prt.h:188
bool lgAbort
Definition: cddefines.cpp:10
bool lgSortLineWavelength
Definition: prt.h:146
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition: service.cpp:1318
int column_len
Definition: prt.h:318
double ScaleNormLine
Definition: lines.h:117
bool lgPrtLineArray
Definition: prt.h:260
long int StuffComment(const char *chComment)
Definition: prt_final.cpp:1943
t_rt rt
Definition: rt.cpp:5