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

Generated for cloudy by doxygen 1.8.1.2