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